From 5f39743dca8351b94c9e075dc521da4a8b08f6ad Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Wed, 30 Aug 2023 16:44:28 +0200 Subject: [PATCH 01/11] start doing surrounding points option --- polytope/datacube/backends/datacube.py | 8 +-- polytope/datacube/datacube_axis.py | 44 ++++++++++---- polytope/engine/hullslicer.py | 3 +- polytope/polytope.py | 2 +- polytope/shapes.py | 9 ++- tests/test_snapping.py | 79 ++++++++++++++++++++++++++ 6 files changed, 125 insertions(+), 20 deletions(-) create mode 100644 tests/test_snapping.py diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index d26f3f952..6b81e0b4d 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -79,7 +79,7 @@ def fit_path(self, path): path.pop(key) return path - def get_indices(self, path: DatacubePath, axis, lower, upper): + def get_indices(self, path: DatacubePath, axis, lower, upper, method): """ Given a path to a subset of the datacube, return the discrete indexes which exist between two non-discrete values (lower, upper) for a particular axis (given by label) @@ -95,21 +95,21 @@ def get_indices(self, path: DatacubePath, axis, lower, upper): for r in original_search_ranges: offset = axis.offset(r) search_ranges_offset.append(offset) - idx_between = self._look_up_datacube(search_ranges, search_ranges_offset, indexes, axis) + idx_between = self._look_up_datacube(search_ranges, search_ranges_offset, indexes, axis, method) # Remove duplicates even if difference of the order of the axis tolerance if offset is not None: # Note that we can only do unique if not dealing with time values idx_between = unique(idx_between) return idx_between - def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis): + def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, method): idx_between = [] for i in range(len(search_ranges)): r = search_ranges[i] offset = search_ranges_offset[i] low = r[0] up = r[1] - indexes_between = axis.find_indices_between([indexes], low, up, self) + indexes_between = axis.find_indices_between([indexes], low, up, self, method) # Now the indexes_between are values on the cyclic range so need to remap them to their original # values before returning them for j in range(len(indexes_between)): diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 80fec89ae..508495633 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -232,7 +232,8 @@ def unmap_total_path_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube): + def find_indices_between(index_ranges, low, up, datacube, method): + # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: if isinstance(transform, DatacubeMapper): @@ -303,7 +304,8 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube): + def find_indices_between(index_ranges, low, up, datacube, method): + # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: if isinstance(transform, DatacubeAxisMerger): @@ -347,7 +349,8 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube): + def find_indices_between(index_ranges, low, up, datacube, method): + # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: if isinstance(transform, DatacubeAxisReverse): @@ -415,7 +418,8 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube): + def find_indices_between(index_ranges, low, up, datacube, method): + # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: if isinstance(transform, DatacubeAxisTypeChange): @@ -508,20 +512,38 @@ def unmap_total_path_to_datacube(self, path, unmapped_path): def remap_to_requeest(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(self, index_ranges, low, up, datacube): + def find_indices_between(self, index_ranges, low, up, datacube, method): + # TODO: add method for snappping indexes_between_ranges = [] for indexes in index_ranges: if self.name in datacube.complete_axes: # Find the range of indexes between lower and upper # https://pandas.pydata.org/docs/reference/api/pandas.Index.searchsorted.html # Assumes the indexes are already sorted (could sort to be sure) and monotonically increasing - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - indexes_between = indexes[start:end].to_list() - indexes_between_ranges.append(indexes_between) + if method == "snap": + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + start = max(start-1, 0) + end = min(end+1, len(indexes)) + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) + else: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) else: - indexes_between = [i for i in indexes if low <= i <= up] - indexes_between_ranges.append(indexes_between) + if method == "snap": + start = indexes.index(low) + end = indexes.index(up) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + indexes_between = [i for i in indexes if low <= i <= up] + indexes_between_ranges.append(indexes_between) return indexes_between_ranges @staticmethod diff --git a/polytope/engine/hullslicer.py b/polytope/engine/hullslicer.py index 809fa641c..1015f193e 100644 --- a/polytope/engine/hullslicer.py +++ b/polytope/engine/hullslicer.py @@ -46,7 +46,8 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, lower, upper, nex lower = ax.from_float(lower - tol) upper = ax.from_float(upper + tol) flattened = node.flatten() - for value in datacube.get_indices(flattened, ax, lower, upper): + method = polytope.method + for value in datacube.get_indices(flattened, ax, lower, upper, method): # convert to float for slicing fvalue = ax.to_float(value) new_polytope = slice(polytope, ax.name, fvalue) diff --git a/polytope/polytope.py b/polytope/polytope.py index ed1c708aa..5e8fc7133 100644 --- a/polytope/polytope.py +++ b/polytope/polytope.py @@ -42,7 +42,7 @@ def slice(self, polytopes: List[ConvexPolytope]): """Low-level API which takes a polytope geometry object and uses it to slice the datacube""" return self.engine.extract(self.datacube, polytopes) - def retrieve(self, request: Request): + def retrieve(self, request: Request, method="standard"): """Higher-level API which takes a request and uses it to slice the datacube""" request_tree = self.engine.extract(self.datacube, request.polytopes()) self.datacube.get(request_tree) diff --git a/polytope/shapes.py b/polytope/shapes.py index acc456731..e0d0b0bd4 100644 --- a/polytope/shapes.py +++ b/polytope/shapes.py @@ -23,9 +23,11 @@ def axes(self) -> List[str]: class ConvexPolytope(Shape): - def __init__(self, axes, points): + + def __init__(self, axes, points, method=None): self._axes = list(axes) self.points = points + self.method = method def extents(self, axis): slice_axis_idx = self.axes().index(axis) @@ -48,15 +50,16 @@ def polytope(self): class Select(Shape): """Matches several discrete value""" - def __init__(self, axis, values): + def __init__(self, axis, values, method=None): self.axis = axis self.values = values + self.method = method def axes(self): return [self.axis] def polytope(self): - return [ConvexPolytope([self.axis], [[v]]) for v in self.values] + return [ConvexPolytope([self.axis], [[v]], self.method) for v in self.values] class Span(Shape): diff --git a/tests/test_snapping.py b/tests/test_snapping.py new file mode 100644 index 000000000..672e678ed --- /dev/null +++ b/tests/test_snapping.py @@ -0,0 +1,79 @@ +import numpy as np +import xarray as xr + +from polytope.datacube.backends.xarray import XArrayDatacube +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Select + + +class TestSlicing3DXarrayDatacube: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + array = xr.DataArray( + np.random.randn(3, 3), + dims=("level", "step"), + coords={ + "level": [1, 3, 5], + "step": [1, 3, 5], + }, + ) + self.xarraydatacube = XArrayDatacube(array) + self.slicer = HullSlicer() + self.API = Polytope(datacube=array, engine=self.slicer) + + # Testing different shapes + + def test_2D_point(self): + request = Request(Select("level", [2], method="snap"), Select("step", [4], method="snap")) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 4 + for leaf in result.leaves: + path = leaf.flatten() + assert path["level"] in [1, 3] + assert path["step"] in [3, 5] + + def test_2D_point_outside_datacube_left(self): + request = Request(Select("level", [2], method="snap"), Select("step", [0], method="snap")) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 2 + for leaf in result.leaves: + path = leaf.flatten() + assert path["level"] in [1, 3] + assert path["step"] == 1 + + def test_2D_point_outside_datacube_right(self): + request = Request(Select("level", [2], method="snap"), Select("step", [6], method="snap")) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 2 + for leaf in result.leaves: + path = leaf.flatten() + assert path["level"] in [1, 3] + assert path["step"] == 5 + + def test_1D_point_outside_datacube_right(self): + request = Request(Select("level", [1]), Select("step", [6], method="snap")) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 1 + for leaf in result.leaves: + path = leaf.flatten() + assert path["level"] == 1 + assert path["step"] == 5 + + def test_1D_nonexisting_point(self): + request = Request(Select("level", [2]), Select("step", [6], method="snap")) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 1 + assert result.is_root() + + def test_1D_nonexisting_point_v2(self): + request = Request(Select("level", [2], method="snap"), Select("step", [6])) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 1 + assert result.is_root() From 63857e3e14790075c3f25ce82ad79aa9d8bccac3 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 31 Aug 2023 10:26:06 +0200 Subject: [PATCH 02/11] make surrounding method work for all transformations --- polytope/datacube/backends/datacube.py | 26 ++++-- polytope/datacube/backends/mock.py | 2 +- polytope/datacube/datacube_axis.py | 110 +++++++++++++++++++++---- tests/test_cyclic_snapping.py | 36 ++++++++ tests/test_snapping.py | 12 +-- 5 files changed, 157 insertions(+), 29 deletions(-) create mode 100644 tests/test_cyclic_snapping.py diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index 6b81e0b4d..1c02ed423 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -79,7 +79,7 @@ def fit_path(self, path): path.pop(key) return path - def get_indices(self, path: DatacubePath, axis, lower, upper, method): + def get_indices(self, path: DatacubePath, axis, lower, upper, method=None): """ Given a path to a subset of the datacube, return the discrete indexes which exist between two non-discrete values (lower, upper) for a particular axis (given by label) @@ -113,12 +113,24 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, # Now the indexes_between are values on the cyclic range so need to remap them to their original # values before returning them for j in range(len(indexes_between)): - for k in range(len(indexes_between[j])): - if offset is None: - indexes_between[j][k] = indexes_between[j][k] - else: - indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) - idx_between.append(indexes_between[j][k]) + # if we have a special indexes between range that needs additional offset, treat it here + if indexes_between[j][0] == "need_offset": + new_offset = indexes_between[j][1] + for k in range(2, len(indexes_between[j])): + if offset is None: + indexes_between[j][k] = indexes_between[j][k] + else: + offset = offset + new_offset + indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + idx_between.append(indexes_between[j][k]) + else: + # do normal offset if no new offset + for k in range(len(indexes_between[j])): + if offset is None: + indexes_between[j][k] = indexes_between[j][k] + else: + indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + idx_between.append(indexes_between[j][k]) return idx_between def get_mapper(self, axis): diff --git a/polytope/datacube/backends/mock.py b/polytope/datacube/backends/mock.py index 351c0208e..b6bc32d3e 100644 --- a/polytope/datacube/backends/mock.py +++ b/polytope/datacube/backends/mock.py @@ -41,7 +41,7 @@ def get(self, requests: IndexTree): def get_mapper(self, axis): return self.mappers[axis] - def get_indices(self, path: DatacubePath, axis, lower, upper): + def get_indices(self, path: DatacubePath, axis, lower, upper, method=None): if lower == upper == math.ceil(lower): if lower >= 0: return [int(lower)] diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 508495633..01e8eda85 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -135,6 +135,41 @@ def unmap_to_datacube(path, unmapped_path): (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) return (path, unmapped_path) + old_find_indices_between = cls.find_indices_between + + def find_indices_between(index_ranges, low, up, datacube, method=None): + update_range() + range_length = cls.range[1] - cls.range[0] + indexes_between_ranges = [] + if method == "surrounding": + for indexes in index_ranges: + if cls.name in datacube.complete_axes: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + else: + start = indexes.index(low) + end = indexes.index(up) + if start-1 < 0: + start_offset_indicator = "need_offset" + new_offset = - range_length # TODO: not sure about this offset + index_val_found = indexes[-1:][0] + indexes_between_before = [start_offset_indicator, new_offset, index_val_found].to_list() + indexes_between_ranges.append(indexes_between_before) + if end+1 > len(indexes): + start_offset_indicator = "need_offset" + new_offset = range_length + index_val_found = indexes[:1][0] + indexes_between_after = [start_offset_indicator, new_offset, index_val_found].to_list() + indexes_between_ranges.append(indexes_between_after) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + else: + return old_find_indices_between(index_ranges, low, up, datacube, method) + def offset(range): # We first unpad the range by the axis tolerance to make sure that # we find the wanted range of the cyclic axis since we padded by the axis tolerance before. @@ -150,6 +185,7 @@ def offset(range): cls.find_indexes = find_indexes cls.unmap_to_datacube = unmap_to_datacube cls.unmap_total_path_to_datacube = unmap_total_path_to_datacube + cls.find_indices_between = find_indices_between return cls @@ -232,7 +268,7 @@ def unmap_total_path_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -240,8 +276,19 @@ def find_indices_between(index_ranges, low, up, datacube, method): transformation = transform if cls.name in transformation._mapped_axes(): for idxs in index_ranges: - indexes_between = [i for i in idxs if low <= i <= up] - indexes_between_ranges.append(indexes_between) + if method == "surrounding": + start = idxs.index(low) + end = idxs.index(up) + start = max(start-1, 0) + end = min(end+1, len(idxs)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = idxs[start:end] + indexes_between_ranges.append(indexes_between) + else: + indexes_between = [i for i in idxs if low <= i <= up] + indexes_between_ranges.append(indexes_between) + # indexes_between = [i for i in idxs if low <= i <= up] + # indexes_between_ranges.append(indexes_between) return indexes_between_ranges old_remap = cls.remap @@ -304,7 +351,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -312,8 +359,19 @@ def find_indices_between(index_ranges, low, up, datacube, method): transformation = transform if cls.name in transformation._mapped_axes(): for indexes in index_ranges: - indexes_between = [i for i in indexes if low <= i <= up] - indexes_between_ranges.append(indexes_between) + if method == "surrounding": + start = indexes.index(low) + end = indexes.index(up) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + indexes_between = [i for i in indexes if low <= i <= up] + indexes_between_ranges.append(indexes_between) + # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between_ranges.append(indexes_between) return indexes_between_ranges def remap(range): @@ -349,7 +407,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -357,8 +415,19 @@ def find_indices_between(index_ranges, low, up, datacube, method): transformation = transform if cls.name == transformation.name: for indexes in index_ranges: - indexes_between = [i for i in indexes if low <= i <= up] - indexes_between_ranges.append(indexes_between) + if method == "surrounding": + start = indexes.index(low) + end = indexes.index(up) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + indexes_between = [i for i in indexes if low <= i <= up] + indexes_between_ranges.append(indexes_between) + # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between_ranges.append(indexes_between) return indexes_between_ranges def remap(range): @@ -418,7 +487,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -426,8 +495,19 @@ def find_indices_between(index_ranges, low, up, datacube, method): transformation = transform if cls.name == transformation.name: for indexes in index_ranges: - indexes_between = [i for i in indexes if low <= i <= up] - indexes_between_ranges.append(indexes_between) + if method == "surrounding": + start = indexes.index(low) + end = indexes.index(up) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + indexes_between = [i for i in indexes if low <= i <= up] + indexes_between_ranges.append(indexes_between) + # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between_ranges.append(indexes_between) return indexes_between_ranges def remap(range): @@ -512,7 +592,7 @@ def unmap_total_path_to_datacube(self, path, unmapped_path): def remap_to_requeest(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(self, index_ranges, low, up, datacube, method): + def find_indices_between(self, index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for indexes in index_ranges: @@ -520,7 +600,7 @@ def find_indices_between(self, index_ranges, low, up, datacube, method): # Find the range of indexes between lower and upper # https://pandas.pydata.org/docs/reference/api/pandas.Index.searchsorted.html # Assumes the indexes are already sorted (could sort to be sure) and monotonically increasing - if method == "snap": + if method == "surrounding": start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") start = max(start-1, 0) @@ -533,7 +613,7 @@ def find_indices_between(self, index_ranges, low, up, datacube, method): indexes_between = indexes[start:end].to_list() indexes_between_ranges.append(indexes_between) else: - if method == "snap": + if method == "surrounding": start = indexes.index(low) end = indexes.index(up) start = max(start-1, 0) diff --git a/tests/test_cyclic_snapping.py b/tests/test_cyclic_snapping.py new file mode 100644 index 000000000..9f7b045f2 --- /dev/null +++ b/tests/test_cyclic_snapping.py @@ -0,0 +1,36 @@ +import xarray as xr + +from polytope.datacube.backends.xarray import XArrayDatacube +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Select + + +class TestSlicing3DXarrayDatacube: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + array = xr.DataArray( + [0, 1, 2], + dims=("long"), + coords={ + "long": [0, 0.5, 1.0], + }, + ) + options = {"long": {"transformation": {"cyclic": [0, 1.0]}}} + self.xarraydatacube = XArrayDatacube(array) + self.slicer = HullSlicer() + self.API = Polytope(datacube=array, engine=self.slicer, axis_options=options) + + # Testing different shapes + + def test_cyclic_float_axis_across_seam(self): + request = Request( + Select("long", [-0.2], method="surrounding") + ) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 2 + assert result.leaves[0].flatten()["long"] == -0.5 + assert result.leaves[1].flatten()["long"] == 0. + assert result.leaves[0].result == (None, 1) + assert result.leaves[1].result == (None, 0) diff --git a/tests/test_snapping.py b/tests/test_snapping.py index 672e678ed..7d127a877 100644 --- a/tests/test_snapping.py +++ b/tests/test_snapping.py @@ -25,7 +25,7 @@ def setup_method(self, method): # Testing different shapes def test_2D_point(self): - request = Request(Select("level", [2], method="snap"), Select("step", [4], method="snap")) + request = Request(Select("level", [2], method="surrounding"), Select("step", [4], method="surrounding")) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 4 @@ -35,7 +35,7 @@ def test_2D_point(self): assert path["step"] in [3, 5] def test_2D_point_outside_datacube_left(self): - request = Request(Select("level", [2], method="snap"), Select("step", [0], method="snap")) + request = Request(Select("level", [2], method="surrounding"), Select("step", [0], method="surrounding")) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 2 @@ -45,7 +45,7 @@ def test_2D_point_outside_datacube_left(self): assert path["step"] == 1 def test_2D_point_outside_datacube_right(self): - request = Request(Select("level", [2], method="snap"), Select("step", [6], method="snap")) + request = Request(Select("level", [2], method="surrounding"), Select("step", [6], method="surrounding")) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 2 @@ -55,7 +55,7 @@ def test_2D_point_outside_datacube_right(self): assert path["step"] == 5 def test_1D_point_outside_datacube_right(self): - request = Request(Select("level", [1]), Select("step", [6], method="snap")) + request = Request(Select("level", [1]), Select("step", [6], method="surrounding")) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 1 @@ -65,14 +65,14 @@ def test_1D_point_outside_datacube_right(self): assert path["step"] == 5 def test_1D_nonexisting_point(self): - request = Request(Select("level", [2]), Select("step", [6], method="snap")) + request = Request(Select("level", [2]), Select("step", [6], method="surrounding")) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 1 assert result.is_root() def test_1D_nonexisting_point_v2(self): - request = Request(Select("level", [2], method="snap"), Select("step", [6])) + request = Request(Select("level", [2], method="surrounding"), Select("step", [6])) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 1 From 7917b68a0e69135bcc9bb348f2a25b58e7a7675b Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 31 Aug 2023 16:56:12 +0200 Subject: [PATCH 03/11] try to make cyclic surrounding method work --- polytope/datacube/backends/datacube.py | 41 ++-- polytope/datacube/datacube_axis.py | 328 +++++++++++++++++++++---- tests/test_snapping_real_data.py | 59 +++++ 3 files changed, 367 insertions(+), 61 deletions(-) create mode 100644 tests/test_snapping_real_data.py diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index 1c02ed423..ae1ba527b 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -109,28 +109,35 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, offset = search_ranges_offset[i] low = r[0] up = r[1] - indexes_between = axis.find_indices_between([indexes], low, up, self, method) + indexes_between = axis.find_indices_between([indexes], low, up, self, offset, method) # Now the indexes_between are values on the cyclic range so need to remap them to their original # values before returning them for j in range(len(indexes_between)): # if we have a special indexes between range that needs additional offset, treat it here - if indexes_between[j][0] == "need_offset": - new_offset = indexes_between[j][1] - for k in range(2, len(indexes_between[j])): - if offset is None: - indexes_between[j][k] = indexes_between[j][k] - else: - offset = offset + new_offset - indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) - idx_between.append(indexes_between[j][k]) + if len(indexes_between[j]) == 0: + idx_between = idx_between else: - # do normal offset if no new offset - for k in range(len(indexes_between[j])): - if offset is None: - indexes_between[j][k] = indexes_between[j][k] - else: - indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) - idx_between.append(indexes_between[j][k]) + if indexes_between[j][0] == "need_offset": + print(indexes_between[j]) + new_offset = indexes_between[j][1] + for k in range(2, len(indexes_between[j])): + if offset is None: + indexes_between[j][k] = indexes_between[j][k] + else: + # offset = offset + new_offset + offset = offset + new_offset + indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + idx_between.append(indexes_between[j][k]) + else: + # do normal offset if no new offset + for k in range(len(indexes_between[j])): + print("HERE") + print(offset) + if offset is None: + indexes_between[j][k] = indexes_between[j][k] + else: + indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + idx_between.append(indexes_between[j][k]) return idx_between def get_mapper(self, axis): diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 01e8eda85..ed716d8f7 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -137,38 +137,250 @@ def unmap_to_datacube(path, unmapped_path): old_find_indices_between = cls.find_indices_between - def find_indices_between(index_ranges, low, up, datacube, method=None): + def find_indices_between(index_ranges, low, up, datacube, offset, method=None): update_range() range_length = cls.range[1] - cls.range[0] indexes_between_ranges = [] - if method == "surrounding": - for indexes in index_ranges: - if cls.name in datacube.complete_axes: - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - else: - start = indexes.index(low) - end = indexes.index(up) - if start-1 < 0: - start_offset_indicator = "need_offset" - new_offset = - range_length # TODO: not sure about this offset - index_val_found = indexes[-1:][0] - indexes_between_before = [start_offset_indicator, new_offset, index_val_found].to_list() - indexes_between_ranges.append(indexes_between_before) - if end+1 > len(indexes): - start_offset_indicator = "need_offset" - new_offset = range_length - index_val_found = indexes[:1][0] - indexes_between_after = [start_offset_indicator, new_offset, index_val_found].to_list() - indexes_between_ranges.append(indexes_between_after) - start = max(start-1, 0) - end = min(end+1, len(indexes)) - # indexes_between = [i for i in indexes if low <= i <= up] - indexes_between = indexes[start:end].to_list() - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - else: - return old_find_indices_between(index_ranges, low, up, datacube, method) + + if offset < 0: + new_offset = 0 + while low > cls.range[0]: + low = low - range_length + new_offset -= range_length + up = up - range_length + if method == "surrounding": + for indexes in index_ranges: + if cls.name in datacube.complete_axes: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + # TODO: need to include the offset here to understand when start-1< 0 with the offset accounted for + if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) + start_offset_indicator = "need_offset" + # new_offset = new_offset - range_length + # new_offset = offset - range_length + index_val_found = indexes[-1:][0] + indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_before) + if end+1 > len(indexes): + start_offset_indicator = "need_offset" + index_val_found = indexes[:1][0] + indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_after) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end].to_list() + indexes_between = [i - new_offset for i in indexes_between] + indexes_between_ranges.append(indexes_between) + print("IDXS") + print(indexes_between_ranges) + else: + start = indexes.index(low) + end = indexes.index(up) + if start-1 < 0: + start_offset_indicator = "need_offset" + index_val_found = indexes[-1:][0] + indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_before) + if end+1 > len(indexes): + start_offset_indicator = "need_offset" + index_val_found = indexes[:1][0] + indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_after) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + else: + return old_find_indices_between(index_ranges, low, up, datacube, method, offset) + + + + # if offset > 0: + # low = low - offset + # up = up - offset + # else: + # low = low + range_length + offset + # up = up + range_length + offset + # if method == "surrounding": + # for indexes in index_ranges: + # if cls.name in datacube.complete_axes: + # start = indexes.searchsorted(low, "left") + # end = indexes.searchsorted(up, "right") + # # TODO: need to include the offset here to understand when start-1< 0 with the offset accounted for + # if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) + # start_offset_indicator = "need_offset" + # # new_offset = - range_length # TODO: not sure about this offset + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # print(range_length) + # new_offset = range_length + # else: + # new_offset = - range_length + # # new_offset = - range_length + # index_val_found = indexes[-1:][0] + # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_before) + # if end+1 > len(indexes): + # start_offset_indicator = "need_offset" + # # new_offset = range_length + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # new_offset = - range_length + # else: + # new_offset = range_length + # index_val_found = indexes[:1][0] + # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_after) + # start = max(start-1, 0) + # end = min(end+1, len(indexes)) + # # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between = indexes[start:end].to_list() + # indexes_between_ranges.append(indexes_between) + # else: + # start = indexes.index(low) + # end = indexes.index(up) + # if start-1 < 0: + # start_offset_indicator = "need_offset" + # # new_offset = - range_length # TODO: not sure about this offset + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # new_offset = range_length + # else: + # new_offset = - range_length + # index_val_found = indexes[-1:][0] + # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_before) + # if end+1 > len(indexes): + # start_offset_indicator = "need_offset" + # # new_offset = range_length + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # new_offset = - range_length + # else: + # new_offset = range_length + # index_val_found = indexes[:1][0] + # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_after) + # start = max(start-1, 0) + # end = min(end+1, len(indexes)) + # # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between = indexes[start:end] + # indexes_between_ranges.append(indexes_between) + # return indexes_between_ranges + # else: + # return old_find_indices_between(index_ranges, low, up, datacube, method, offset) + + # def find_indices_between(index_ranges, low, up, datacube, offset, method=None): + # update_range() + # range_length = cls.range[1] - cls.range[0] + # indexes_between_ranges = [] + # if offset > 0: + # low = low - offset + # up = up - offset + # else: + # print("HERE") + # print(low) + # # low = low + offset + # low = low + range_length + offset + # print(low) + # # up = up + offset + # up = up + range_length + offset + # if method == "surrounding": + # for indexes in index_ranges: + # # if offset > 0: + # # low = low - offset + # # up = up - offset + # # else: + # # print("HERE") + # # print(low) + # # low = low + offset + # # print(low) + # # up = up + offset + # if cls.name in datacube.complete_axes: + # start = indexes.searchsorted(low, "left") + # end = indexes.searchsorted(up, "right") + # # TODO: need to include the offset here to understand when start-1< 0 with the offset accounted for + # if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) + # start_offset_indicator = "need_offset" + # # new_offset = - range_length # TODO: not sure about this offset + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # print(range_length) + # new_offset = range_length + # else: + # new_offset = - range_length + # # new_offset = - range_length + # index_val_found = indexes[-1:][0] + # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_before) + # if end+1 > len(indexes): + # start_offset_indicator = "need_offset" + # # new_offset = range_length + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # new_offset = - range_length + # else: + # new_offset = range_length + # index_val_found = indexes[:1][0] + # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_after) + # start = max(start-1, 0) + # end = min(end+1, len(indexes)) + # # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between = indexes[start:end].to_list() + # indexes_between_ranges.append(indexes_between) + # else: + # start = indexes.index(low) + # end = indexes.index(up) + # if start-1 < 0: + # start_offset_indicator = "need_offset" + # # new_offset = - range_length # TODO: not sure about this offset + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # new_offset = range_length + # else: + # new_offset = - range_length + # index_val_found = indexes[-1:][0] + # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_before) + # if end+1 > len(indexes): + # start_offset_indicator = "need_offset" + # # new_offset = range_length + # # new_offset = 0 + # if offset == 0: + # new_offset = 0 + # elif offset < 0: + # new_offset = - range_length + # else: + # new_offset = range_length + # index_val_found = indexes[:1][0] + # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + # indexes_between_ranges.append(indexes_between_after) + # start = max(start-1, 0) + # end = min(end+1, len(indexes)) + # # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between = indexes[start:end] + # indexes_between_ranges.append(indexes_between) + # return indexes_between_ranges + # else: + # return old_find_indices_between(index_ranges, low, up, datacube, method, offset) def offset(range): # We first unpad the range by the axis tolerance to make sure that @@ -268,7 +480,7 @@ def unmap_total_path_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method=None): + def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -351,7 +563,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method=None): + def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -407,7 +619,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method=None): + def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -415,17 +627,45 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): transformation = transform if cls.name == transformation.name: for indexes in index_ranges: - if method == "surrounding": - start = indexes.index(low) - end = indexes.index(up) - start = max(start-1, 0) - end = min(end+1, len(indexes)) - # indexes_between = [i for i in indexes if low <= i <= up] - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) + if cls.name in datacube.complete_axes: + # Find the range of indexes between lower and upper + # https://pandas.pydata.org/docs/reference/api/pandas.Index.searchsorted.html + # Assumes the indexes are already sorted (could sort to be sure) and monotonically increasing + if method == "surrounding": + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + start = max(start-1, 0) + end = min(end+1, len(indexes)) + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) + else: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) else: - indexes_between = [i for i in indexes if low <= i <= up] - indexes_between_ranges.append(indexes_between) + if method == "surrounding": + start = indexes.index(low) + end = indexes.index(up) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + # indexes_between = [i for i in indexes if low <= i <= up] + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + indexes_between = [i for i in indexes if low <= i <= up] + indexes_between_ranges.append(indexes_between) + # if method == "surrounding": + # start = indexes.index(low) + # end = indexes.index(up) + # start = max(start-1, 0) + # end = min(end+1, len(indexes)) + # # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between = indexes[start:end] + # indexes_between_ranges.append(indexes_between) + # else: + # indexes_between = [i for i in indexes if low <= i <= up] + # indexes_between_ranges.append(indexes_between) # indexes_between = [i for i in indexes if low <= i <= up] # indexes_between_ranges.append(indexes_between) return indexes_between_ranges @@ -487,7 +727,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, method=None): + def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -592,7 +832,7 @@ def unmap_total_path_to_datacube(self, path, unmapped_path): def remap_to_requeest(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(self, index_ranges, low, up, datacube, method=None): + def find_indices_between(self, index_ranges, low, up, datacube, offset, method=None): # TODO: add method for snappping indexes_between_ranges = [] for indexes in index_ranges: diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py new file mode 100644 index 000000000..9513daa5f --- /dev/null +++ b/tests/test_snapping_real_data.py @@ -0,0 +1,59 @@ +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +from earthkit import data + +from polytope.datacube.backends.xarray import XArrayDatacube +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Box, Select + + +class TestSlicingEra5Data: + def setup_method(self, method): + ds = data.from_source("file", "./tests/data/era5-levels-members.grib") + array = ds.to_xarray().isel(step=0).t + self.xarraydatacube = XArrayDatacube(array) + self.slicer = HullSlicer() + options = {"latitude": {"transformation": {"reverse": {True}}}, + "longitude": {"transformation": {"cyclic": [0, 360.0]}}} + self.API = Polytope(datacube=array, engine=self.slicer, axis_options=options) + + def test_surrounding_on_grid_point(self): + requested_lat = 0 + requested_lon = -720 + request = Request( + Box(["number", "isobaricInhPa"], [6, 500.0], [6, 850.0]), + Select("time", ["2017-01-02T12:00:00"]), + # Box(["latitude", "longitude"], lower_corner=[0.0, 0.0], upper_corner=[10.0, 30.0]), + Select("latitude", [requested_lat], method="surrounding"), + Select("longitude", [requested_lon], method="surrounding"), + Select("step", [np.timedelta64(0, "s")]), + ) + result = self.API.retrieve(request) + result.pprint() + country_points_plotting = [] + lats = [] + longs = [] + temps = [] + for i in range(len(result.leaves)): + cubepath = result.leaves[i].flatten() + lat = cubepath["latitude"] + long = cubepath["longitude"] + latlong_point = [lat, long] + lats.append(lat) + longs.append(long) + t_idx = result.leaves[i].result[1] + temps.append(t_idx) + country_points_plotting.append(latlong_point) + temps = np.array(temps) + + # Plot all the points on a world map + worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) + fig, ax = plt.subplots(figsize=(12, 6)) + worldmap.plot(color="darkgrey", ax=ax) + + plt.scatter(longs, lats, s=16, c=temps, cmap="YlOrRd") + plt.scatter([requested_lon], [requested_lat], s=16, c="blue") + plt.colorbar(label="Temperature") + plt.show() From 3bf4f2db5be6c3ebff3cb36cba74d34cd03221dd Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Fri, 1 Sep 2023 11:14:08 +0200 Subject: [PATCH 04/11] make surrounding method work with cyclic axes --- polytope/datacube/backends/datacube.py | 10 +- polytope/datacube/datacube_axis.py | 245 +++++-------------------- tests/test_snapping_real_data.py | 2 +- 3 files changed, 55 insertions(+), 202 deletions(-) diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index ae1ba527b..45a320d9d 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -118,25 +118,23 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, idx_between = idx_between else: if indexes_between[j][0] == "need_offset": - print(indexes_between[j]) new_offset = indexes_between[j][1] for k in range(2, len(indexes_between[j])): if offset is None: indexes_between[j][k] = indexes_between[j][k] else: - # offset = offset + new_offset offset = offset + new_offset - indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + indexes_between[j][k] = round(indexes_between[j][k] + offset, + int(-math.log10(axis.tol))) idx_between.append(indexes_between[j][k]) else: # do normal offset if no new offset for k in range(len(indexes_between[j])): - print("HERE") - print(offset) if offset is None: indexes_between[j][k] = indexes_between[j][k] else: - indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + indexes_between[j][k] = round(indexes_between[j][k] + offset, + int(-math.log10(axis.tol))) idx_between.append(indexes_between[j][k]) return idx_between diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index ed716d8f7..1c77b1912 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -142,9 +142,11 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): range_length = cls.range[1] - cls.range[0] indexes_between_ranges = [] - if offset < 0: + if offset != 0: + # NOTE that if the offset is not 0, then we need to recenter the low and up + # values to be within the datacube range new_offset = 0 - while low > cls.range[0]: + while low >= cls.range[0] - cls.tol: low = low - range_length new_offset -= range_length up = up - range_length @@ -153,11 +155,8 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): if cls.name in datacube.complete_axes: start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") - # TODO: need to include the offset here to understand when start-1< 0 with the offset accounted for if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) start_offset_indicator = "need_offset" - # new_offset = new_offset - range_length - # new_offset = offset - range_length index_val_found = indexes[-1:][0] indexes_between_before = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_before) @@ -168,12 +167,9 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): indexes_between_ranges.append(indexes_between_after) start = max(start-1, 0) end = min(end+1, len(indexes)) - # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = indexes[start:end].to_list() indexes_between = [i - new_offset for i in indexes_between] indexes_between_ranges.append(indexes_between) - print("IDXS") - print(indexes_between_ranges) else: start = indexes.index(low) end = indexes.index(up) @@ -189,198 +185,56 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): indexes_between_ranges.append(indexes_between_after) start = max(start-1, 0) end = min(end+1, len(indexes)) - # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = indexes[start:end] indexes_between_ranges.append(indexes_between) return indexes_between_ranges else: return old_find_indices_between(index_ranges, low, up, datacube, method, offset) - - - # if offset > 0: - # low = low - offset - # up = up - offset - # else: - # low = low + range_length + offset - # up = up + range_length + offset - # if method == "surrounding": - # for indexes in index_ranges: - # if cls.name in datacube.complete_axes: - # start = indexes.searchsorted(low, "left") - # end = indexes.searchsorted(up, "right") - # # TODO: need to include the offset here to understand when start-1< 0 with the offset accounted for - # if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) - # start_offset_indicator = "need_offset" - # # new_offset = - range_length # TODO: not sure about this offset - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # print(range_length) - # new_offset = range_length - # else: - # new_offset = - range_length - # # new_offset = - range_length - # index_val_found = indexes[-1:][0] - # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_before) - # if end+1 > len(indexes): - # start_offset_indicator = "need_offset" - # # new_offset = range_length - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # new_offset = - range_length - # else: - # new_offset = range_length - # index_val_found = indexes[:1][0] - # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_after) - # start = max(start-1, 0) - # end = min(end+1, len(indexes)) - # # indexes_between = [i for i in indexes if low <= i <= up] - # indexes_between = indexes[start:end].to_list() - # indexes_between_ranges.append(indexes_between) - # else: - # start = indexes.index(low) - # end = indexes.index(up) - # if start-1 < 0: - # start_offset_indicator = "need_offset" - # # new_offset = - range_length # TODO: not sure about this offset - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # new_offset = range_length - # else: - # new_offset = - range_length - # index_val_found = indexes[-1:][0] - # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_before) - # if end+1 > len(indexes): - # start_offset_indicator = "need_offset" - # # new_offset = range_length - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # new_offset = - range_length - # else: - # new_offset = range_length - # index_val_found = indexes[:1][0] - # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_after) - # start = max(start-1, 0) - # end = min(end+1, len(indexes)) - # # indexes_between = [i for i in indexes if low <= i <= up] - # indexes_between = indexes[start:end] - # indexes_between_ranges.append(indexes_between) - # return indexes_between_ranges - # else: - # return old_find_indices_between(index_ranges, low, up, datacube, method, offset) - - # def find_indices_between(index_ranges, low, up, datacube, offset, method=None): - # update_range() - # range_length = cls.range[1] - cls.range[0] - # indexes_between_ranges = [] - # if offset > 0: - # low = low - offset - # up = up - offset - # else: - # print("HERE") - # print(low) - # # low = low + offset - # low = low + range_length + offset - # print(low) - # # up = up + offset - # up = up + range_length + offset - # if method == "surrounding": - # for indexes in index_ranges: - # # if offset > 0: - # # low = low - offset - # # up = up - offset - # # else: - # # print("HERE") - # # print(low) - # # low = low + offset - # # print(low) - # # up = up + offset - # if cls.name in datacube.complete_axes: - # start = indexes.searchsorted(low, "left") - # end = indexes.searchsorted(up, "right") - # # TODO: need to include the offset here to understand when start-1< 0 with the offset accounted for - # if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) - # start_offset_indicator = "need_offset" - # # new_offset = - range_length # TODO: not sure about this offset - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # print(range_length) - # new_offset = range_length - # else: - # new_offset = - range_length - # # new_offset = - range_length - # index_val_found = indexes[-1:][0] - # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_before) - # if end+1 > len(indexes): - # start_offset_indicator = "need_offset" - # # new_offset = range_length - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # new_offset = - range_length - # else: - # new_offset = range_length - # index_val_found = indexes[:1][0] - # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_after) - # start = max(start-1, 0) - # end = min(end+1, len(indexes)) - # # indexes_between = [i for i in indexes if low <= i <= up] - # indexes_between = indexes[start:end].to_list() - # indexes_between_ranges.append(indexes_between) - # else: - # start = indexes.index(low) - # end = indexes.index(up) - # if start-1 < 0: - # start_offset_indicator = "need_offset" - # # new_offset = - range_length # TODO: not sure about this offset - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # new_offset = range_length - # else: - # new_offset = - range_length - # index_val_found = indexes[-1:][0] - # indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_before) - # if end+1 > len(indexes): - # start_offset_indicator = "need_offset" - # # new_offset = range_length - # # new_offset = 0 - # if offset == 0: - # new_offset = 0 - # elif offset < 0: - # new_offset = - range_length - # else: - # new_offset = range_length - # index_val_found = indexes[:1][0] - # indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - # indexes_between_ranges.append(indexes_between_after) - # start = max(start-1, 0) - # end = min(end+1, len(indexes)) - # # indexes_between = [i for i in indexes if low <= i <= up] - # indexes_between = indexes[start:end] - # indexes_between_ranges.append(indexes_between) - # return indexes_between_ranges - # else: - # return old_find_indices_between(index_ranges, low, up, datacube, method, offset) + else: + # If the offset is 0, then the first value found on the left has an offset of range_length + new_offset = 0 + if method == "surrounding": + for indexes in index_ranges: + if cls.name in datacube.complete_axes: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) + start_offset_indicator = "need_offset" + index_val_found = indexes[-1:][0] + new_offset = -range_length + indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_before) + if end+1 > len(indexes): + start_offset_indicator = "need_offset" + index_val_found = indexes[:1][0] + indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_after) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + indexes_between = indexes[start:end].to_list() + indexes_between = [i - new_offset for i in indexes_between] + indexes_between_ranges.append(indexes_between) + else: + start = indexes.index(low) + end = indexes.index(up) + if start-1 < 0: + start_offset_indicator = "need_offset" + index_val_found = indexes[-1:][0] + indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_before) + if end+1 > len(indexes): + start_offset_indicator = "need_offset" + index_val_found = indexes[:1][0] + indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_after) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + else: + return old_find_indices_between(index_ranges, low, up, datacube, method, offset) def offset(range): # We first unpad the range by the axis tolerance to make sure that @@ -630,7 +484,8 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): if cls.name in datacube.complete_axes: # Find the range of indexes between lower and upper # https://pandas.pydata.org/docs/reference/api/pandas.Index.searchsorted.html - # Assumes the indexes are already sorted (could sort to be sure) and monotonically increasing + # Assumes the indexes are already sorted (could sort to be sure) and monotonically + # increasing if method == "surrounding": start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py index 9513daa5f..feb9e01f9 100644 --- a/tests/test_snapping_real_data.py +++ b/tests/test_snapping_real_data.py @@ -21,7 +21,7 @@ def setup_method(self, method): def test_surrounding_on_grid_point(self): requested_lat = 0 - requested_lon = -720 + requested_lon = 720 request = Request( Box(["number", "isobaricInhPa"], [6, 500.0], [6, 850.0]), Select("time", ["2017-01-02T12:00:00"]), From 2b7cd43fc54784d106a2323a25b797531929beba Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Tue, 5 Sep 2023 09:49:31 +0200 Subject: [PATCH 05/11] make all tests work --- polytope/datacube/datacube_axis.py | 12 ++++++---- tests/test_cyclic_simple.py | 37 ++++++++++++++++++++++++++++++ tests/test_cyclic_snapping.py | 2 +- tests/test_datacube_axes_init.py | 2 +- tests/test_snapping_real_data.py | 24 ++++++++++--------- 5 files changed, 60 insertions(+), 17 deletions(-) create mode 100644 tests/test_cyclic_simple.py diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 1c77b1912..db364b4b1 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -142,14 +142,18 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): range_length = cls.range[1] - cls.range[0] indexes_between_ranges = [] + if method != "surrounding": + return old_find_indices_between(index_ranges, low, up, datacube, method, offset) + if offset != 0: # NOTE that if the offset is not 0, then we need to recenter the low and up # values to be within the datacube range new_offset = 0 - while low >= cls.range[0] - cls.tol: - low = low - range_length - new_offset -= range_length - up = up - range_length + if low <= cls.range[0] + cls.tol: + while low >= cls.range[0] - cls.tol: + low = low - range_length + new_offset -= range_length + up = up - range_length if method == "surrounding": for indexes in index_ranges: if cls.name in datacube.complete_axes: diff --git a/tests/test_cyclic_simple.py b/tests/test_cyclic_simple.py new file mode 100644 index 000000000..be48cb8b9 --- /dev/null +++ b/tests/test_cyclic_simple.py @@ -0,0 +1,37 @@ +import numpy as np +import pandas as pd +import xarray as xr + +from polytope.datacube.backends.xarray import XArrayDatacube +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Box, Select + + +class TestSlicing3DXarrayDatacube: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + array = xr.DataArray( + np.random.randn(3, 6, 129, 11), + dims=("date", "step", "level", "long"), + coords={ + "date": pd.date_range("2000-01-01", "2000-01-03", 3), + "step": [0, 3, 6, 9, 12, 15], + "level": range(1, 130), + "long": [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], + }, + ) + options = {"long": {"transformation": {"cyclic": [0, 1.0]}}, "level": {"transformation": {"cyclic": [1, 129]}}} + self.xarraydatacube = XArrayDatacube(array) + self.slicer = HullSlicer() + self.API = Polytope(datacube=array, engine=self.slicer, axis_options=options) + + # Testing different shapes + + def test_cyclic_float_axis_across_seam(self): + request = Request( + Box(["step", "long"], [0, 0.9], [0, 1.2]), Select("date", ["2000-01-01"]), Select("level", [128]) + ) + result = self.API.retrieve(request) + assert len(result.leaves) == 4 + assert [leaf.value for leaf in result.leaves] == [0.9, 1.0, 1.1, 1.2] diff --git a/tests/test_cyclic_snapping.py b/tests/test_cyclic_snapping.py index 9f7b045f2..ac1b11c7f 100644 --- a/tests/test_cyclic_snapping.py +++ b/tests/test_cyclic_snapping.py @@ -31,6 +31,6 @@ def test_cyclic_float_axis_across_seam(self): result.pprint() assert len(result.leaves) == 2 assert result.leaves[0].flatten()["long"] == -0.5 - assert result.leaves[1].flatten()["long"] == 0. + assert result.leaves[1].flatten()["long"] == 0.0 assert result.leaves[0].result == (None, 1) assert result.leaves[1].result == (None, 0) diff --git a/tests/test_datacube_axes_init.py b/tests/test_datacube_axes_init.py index 83461e926..bea08086b 100644 --- a/tests/test_datacube_axes_init.py +++ b/tests/test_datacube_axes_init.py @@ -49,7 +49,7 @@ def test_created_axes(self): assert path == {} assert unmapped_path == {"values": 0} assert lat_ax.find_indices_between([[89.94618771566562, 89.87647835333229]], - 89.87, 90, self.datacube) == [[89.94618771566562, 89.87647835333229]] + 89.87, 90, self.datacube, 0) == [[89.94618771566562, 89.87647835333229]] def test_mapper_transformation_request(self): request = Request( diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py index feb9e01f9..247a276b8 100644 --- a/tests/test_snapping_real_data.py +++ b/tests/test_snapping_real_data.py @@ -1,5 +1,5 @@ -import geopandas as gpd -import matplotlib.pyplot as plt +# import geopandas as gpd +# import matplotlib.pyplot as plt import numpy as np from earthkit import data @@ -21,7 +21,7 @@ def setup_method(self, method): def test_surrounding_on_grid_point(self): requested_lat = 0 - requested_lon = 720 + requested_lon = 0 request = Request( Box(["number", "isobaricInhPa"], [6, 500.0], [6, 850.0]), Select("time", ["2017-01-02T12:00:00"]), @@ -48,12 +48,14 @@ def test_surrounding_on_grid_point(self): country_points_plotting.append(latlong_point) temps = np.array(temps) - # Plot all the points on a world map - worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) - fig, ax = plt.subplots(figsize=(12, 6)) - worldmap.plot(color="darkgrey", ax=ax) + # # Plot all the points on a world map + # worldmap = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) + # fig, ax = plt.subplots(figsize=(12, 6)) + # worldmap.plot(color="darkgrey", ax=ax) - plt.scatter(longs, lats, s=16, c=temps, cmap="YlOrRd") - plt.scatter([requested_lon], [requested_lat], s=16, c="blue") - plt.colorbar(label="Temperature") - plt.show() + # plt.scatter(longs, lats, s=16, c=temps, cmap="YlOrRd") + # plt.scatter([requested_lon], [requested_lat], s=16, c="blue") + # plt.colorbar(label="Temperature") + # plt.show() + for lon in longs: + assert lon in [-3, 0, 3] From ce0661f9cb39e2d3c72c9299fabc56bc05995b74 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Tue, 5 Sep 2023 13:08:29 +0200 Subject: [PATCH 06/11] try to make cyclic without offsets --- polytope/datacube/backends/datacube.py | 15 +- polytope/datacube/datacube_axis.py | 28 ++- polytope/datacube/index_tree.py | 17 +- polytope/engine/hullslicer.py | 10 +- tests/test_cyclic_axis_over_negative_vals.py | 167 ++++++-------- tests/test_cyclic_axis_slicer_not_0.py | 173 +++++++-------- tests/test_cyclic_axis_slicing.py | 215 +++++++------------ tests/test_cyclic_simple.py | 31 ++- tests/test_cyclic_snapping.py | 8 +- 9 files changed, 315 insertions(+), 349 deletions(-) diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index 45a320d9d..0bc710c93 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -1,4 +1,5 @@ import importlib +import itertools import math from abc import ABC, abstractmethod from typing import Any @@ -100,10 +101,15 @@ def get_indices(self, path: DatacubePath, axis, lower, upper, method=None): if offset is not None: # Note that we can only do unique if not dealing with time values idx_between = unique(idx_between) + print("inside get indices") + print(idx_between) return idx_between def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, method): + print("inside look up datacube") idx_between = [] + print(axis.name) + print(search_ranges) for i in range(len(search_ranges)): r = search_ranges[i] offset = search_ranges_offset[i] @@ -118,12 +124,11 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, idx_between = idx_between else: if indexes_between[j][0] == "need_offset": - new_offset = indexes_between[j][1] for k in range(2, len(indexes_between[j])): if offset is None: indexes_between[j][k] = indexes_between[j][k] else: - offset = offset + new_offset + offset = offset indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) idx_between.append(indexes_between[j][k]) @@ -133,9 +138,12 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, if offset is None: indexes_between[j][k] = indexes_between[j][k] else: + offset = offset indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) idx_between.append(indexes_between[j][k]) + print("indexes between") + print(idx_between) return idx_between def get_mapper(self, axis): @@ -149,6 +157,9 @@ def remap_path(self, path: DatacubePath): value = path[key] path[key] = self._axes[key].remap([value, value])[0][0] return path + + # def remap_tree(self, tree: IndexTree): + # for @staticmethod def create(datacube, axis_options: dict): diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index db364b4b1..4c7a66d50 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -149,24 +149,28 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # NOTE that if the offset is not 0, then we need to recenter the low and up # values to be within the datacube range new_offset = 0 - if low <= cls.range[0] + cls.tol: - while low >= cls.range[0] - cls.tol: - low = low - range_length - new_offset -= range_length - up = up - range_length + # if low <= cls.range[0] + cls.tol: + # while low >= cls.range[0] - cls.tol: + # low = low - range_length + # new_offset -= range_length + # up = up - range_length + print("INSIDE THE DATACUBE AXIS DECORATOR") + print((low, up)) if method == "surrounding": for indexes in index_ranges: if cls.name in datacube.complete_axes: start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") + print(start, end) if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) start_offset_indicator = "need_offset" - index_val_found = indexes[-1:][0] + index_val_found = indexes[-2:][0] indexes_between_before = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_before) if end+1 > len(indexes): + print("end exceeded indices length") start_offset_indicator = "need_offset" - index_val_found = indexes[:1][0] + index_val_found = indexes[:2][1] indexes_between_after = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_after) start = max(start-1, 0) @@ -199,25 +203,29 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # If the offset is 0, then the first value found on the left has an offset of range_length new_offset = 0 if method == "surrounding": + print("no offset") for indexes in index_ranges: if cls.name in datacube.complete_axes: start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") + print((start, end)) if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) + print("low start") start_offset_indicator = "need_offset" - index_val_found = indexes[-1:][0] + index_val_found = indexes[-2:][0] + print(index_val_found) new_offset = -range_length indexes_between_before = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_before) if end+1 > len(indexes): start_offset_indicator = "need_offset" - index_val_found = indexes[:1][0] + index_val_found = indexes[:2][1] indexes_between_after = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_after) start = max(start-1, 0) end = min(end+1, len(indexes)) indexes_between = indexes[start:end].to_list() - indexes_between = [i - new_offset for i in indexes_between] + indexes_between = [i for i in indexes_between] indexes_between_ranges.append(indexes_between) else: start = indexes.index(low) diff --git a/polytope/datacube/index_tree.py b/polytope/datacube/index_tree.py index 9833fd240..2d285327e 100644 --- a/polytope/datacube/index_tree.py +++ b/polytope/datacube/index_tree.py @@ -58,7 +58,18 @@ def __hash__(self): def __eq__(self, other): if not isinstance(other, IndexTree): return False - return (self.axis.name, self.value) == (other.axis.name, other.value) + if self.axis.name != other.axis.name: + return False + else: + if other.value == self.value: + return True + if other.value - 2*other.axis.tol <= self.value <= other.value + 2*other.axis.tol: + return True + elif self.value - 2*self.axis.tol <= other.value <= self.value + 2*self.axis.tol: + return True + else: + return False + # return (self.axis.name, self.value) == (other.axis.name, other.value) def __lt__(self, other): return (self.axis.name, self.value) < (other.axis.name, other.value) @@ -73,12 +84,12 @@ def add_child(self, node): self.children.add(node) node._parent = self - def create_child(self, axis, value): + def create_child_not_safe(self, axis, value): node = IndexTree(axis, value) self.add_child(node) return node - def create_child_safe(self, axis, value): + def create_child(self, axis, value): node = IndexTree(axis, value) existing = self.find_child(node) if not existing: diff --git a/polytope/engine/hullslicer.py b/polytope/engine/hullslicer.py index 1015f193e..5553f4e5a 100644 --- a/polytope/engine/hullslicer.py +++ b/polytope/engine/hullslicer.py @@ -1,6 +1,7 @@ from copy import copy from itertools import chain from typing import List +import math import scipy.spatial @@ -52,7 +53,14 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, lower, upper, nex fvalue = ax.to_float(value) new_polytope = slice(polytope, ax.name, fvalue) # store the native type - child = node.create_child(ax, value) + # remapped_val = (ax.remap([value, value])[0][0] + ax.remap([value, value])[0][1])/2 + remapped_val = value + if ax.is_cyclic: + remapped_val = (ax.remap([value, value])[0][0] + ax.remap([value, value])[0][1])/2 + remapped_val = round(remapped_val, int(-math.log10(ax.tol))) + print(remapped_val) + # child = node.create_child(ax, value) + child = node.create_child(ax, remapped_val) child["unsliced_polytopes"] = copy(node["unsliced_polytopes"]) child["unsliced_polytopes"].remove(polytope) if new_polytope is not None: diff --git a/tests/test_cyclic_axis_over_negative_vals.py b/tests/test_cyclic_axis_over_negative_vals.py index 151065ad0..a67333cce 100644 --- a/tests/test_cyclic_axis_over_negative_vals.py +++ b/tests/test_cyclic_axis_over_negative_vals.py @@ -36,29 +36,29 @@ def test_cyclic_float_axis_across_seam(self): Box(["step", "long"], [0, 0.8], [3, 1.7]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) - # result.pprint() + result.pprint() assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2 ] def test_cyclic_float_axis_inside_cyclic_range(self): @@ -69,22 +69,22 @@ def test_cyclic_float_axis_inside_cyclic_range(self): # result.pprint() assert len(result.leaves) == 16 assert [leaf.value for leaf in result.leaves] == [ - 0.0, - 0.1, - 0.2, - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, - 0.0, - 0.1, - 0.2, - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, ] def test_cyclic_float_axis_above_axis_range(self): @@ -94,66 +94,37 @@ def test_cyclic_float_axis_above_axis_range(self): result = self.API.retrieve(request) # result.pprint() assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [1.3, 1.4, 1.5, 1.6, 1.7, 1.3, 1.4, 1.5, 1.6, 1.7] + assert [leaf.value for leaf in result.leaves] == [-0.7, -0.6, -0.5, -0.4, -0.3, + -0.7, -0.6, -0.5, -0.4, -0.3] def test_cyclic_float_axis_two_range_loops(self): request = Request( Box(["step", "long"], [0, 0.3], [3, 2.7]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) - # result.pprint() - assert len(result.leaves) == 50 + result.pprint() + assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, - 2.1, - 2.2, - 2.3, - 2.4, - 2.5, - 2.6, - 2.7, - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, - 2.1, - 2.2, - 2.3, - 2.4, - 2.5, - 2.6, - 2.7, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2 ] def test_cyclic_float_axis_below_axis_range(self): @@ -171,8 +142,11 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): ) result = self.API.retrieve(request) # result.pprint() - assert len(result.leaves) == 22 + assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ + -1.0, + -0.9, + -0.8, -0.7, -0.6, -0.5, @@ -180,10 +154,9 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): -0.3, -0.2, -0.1, - 0.0, - 0.1, - 0.2, - 0.3, + -1.0, + -0.9, + -0.8, -0.7, -0.6, -0.5, @@ -191,8 +164,4 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): -0.3, -0.2, -0.1, - 0.0, - 0.1, - 0.2, - 0.3, ] diff --git a/tests/test_cyclic_axis_slicer_not_0.py b/tests/test_cyclic_axis_slicer_not_0.py index 3173ee77a..33ef062e6 100644 --- a/tests/test_cyclic_axis_slicer_not_0.py +++ b/tests/test_cyclic_axis_slicer_not_0.py @@ -36,28 +36,29 @@ def test_cyclic_float_axis_across_seam(self): Box(["step", "long"], [0, 0.8], [3, 1.7]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) + result.pprint() assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2, ] def test_cyclic_float_axis_inside_cyclic_range(self): @@ -67,22 +68,22 @@ def test_cyclic_float_axis_inside_cyclic_range(self): result = self.API.retrieve(request) assert len(result.leaves) == 16 assert [leaf.value for leaf in result.leaves] == [ - 0.0, - 0.1, - 0.2, - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, - 0.0, - 0.1, - 0.2, - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, ] def test_cyclic_float_axis_above_axis_range(self): @@ -91,65 +92,45 @@ def test_cyclic_float_axis_above_axis_range(self): ) result = self.API.retrieve(request) assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [1.3, 1.4, 1.5, 1.6, 1.7, 1.3, 1.4, 1.5, 1.6, 1.7] + assert [leaf.value for leaf in result.leaves] == [ + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3] def test_cyclic_float_axis_two_range_loops(self): request = Request( Box(["step", "long"], [0, 0.3], [3, 2.7]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) - assert len(result.leaves) == 50 + assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, - 2.1, - 2.2, - 2.3, - 2.4, - 2.5, - 2.6, - 2.7, - 0.3, - 0.4, - 0.5, - 0.6, - 0.7, - 0.8, - 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, - 2.1, - 2.2, - 2.3, - 2.4, - 2.5, - 2.6, - 2.7, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2, + -1.1, + -1.0, + -0.9, + -0.8, + -0.7, + -0.6, + -0.5, + -0.4, + -0.3, + -0.2, ] def test_cyclic_float_axis_below_axis_range(self): @@ -165,8 +146,11 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): Box(["step", "long"], [0, -0.7], [3, 0.3]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) - assert len(result.leaves) == 22 + assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ + -1.0, + -0.9, + -0.8, -0.7, -0.6, -0.5, @@ -174,10 +158,9 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): -0.3, -0.2, -0.1, - 0.0, - 0.1, - 0.2, - 0.3, + -1.0, + -0.9, + -0.8, -0.7, -0.6, -0.5, @@ -185,8 +168,4 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): -0.3, -0.2, -0.1, - 0.0, - 0.1, - 0.2, - 0.3, ] diff --git a/tests/test_cyclic_axis_slicing.py b/tests/test_cyclic_axis_slicing.py index fcf213480..3d3c68465 100644 --- a/tests/test_cyclic_axis_slicing.py +++ b/tests/test_cyclic_axis_slicing.py @@ -34,27 +34,28 @@ def test_cyclic_float_axis_across_seam(self): ) result = self.API.retrieve(request) assert len(result.leaves) == 20 + result.pprint() assert [leaf.value for leaf in result.leaves] == [ + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, 0.8, 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, + 1., + 0.1, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, 0.8, 0.9, - 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, + 1., ] def test_cyclic_float_axis_across_seam_repeated(self): @@ -94,8 +95,8 @@ def test_cyclic_float_axis_across_seam_repeated_twice(self): Box(["step", "long"], [0, 0.0], [3, 2.0]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) - # result.pprint() - assert len(result.leaves) == 22 * 2 - 2 + result.pprint() + assert len(result.leaves) == 22 assert [leaf.value for leaf in result.leaves] == [ 0.0, 0.1, @@ -108,16 +109,6 @@ def test_cyclic_float_axis_across_seam_repeated_twice(self): 0.8, 0.9, 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, 0.0, 0.1, 0.2, @@ -129,16 +120,6 @@ def test_cyclic_float_axis_across_seam_repeated_twice(self): 0.8, 0.9, 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, ] def test_cyclic_float_axis_inside_cyclic_range(self): @@ -176,7 +157,7 @@ def test_cyclic_float_axis_above_axis_range(self): result = self.API.retrieve(request) # result.pprint() assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [1.3, 1.4, 1.5, 1.6, 1.7, 1.3, 1.4, 1.5, 1.6, 1.7] + assert [leaf.value for leaf in result.leaves] == [0.3, 0.4, 0.5, 0.6, 0.7, 0.3, 0.4, 0.5, 0.6, 0.7] def test_cyclic_float_axis_two_range_loops(self): request = Request( @@ -184,8 +165,11 @@ def test_cyclic_float_axis_two_range_loops(self): ) result = self.API.retrieve(request) # result.pprint() - assert len(result.leaves) == 50 + assert len(result.leaves) == 22 assert [leaf.value for leaf in result.leaves] == [ + 0.0, + 0.1, + 0.2, 0.3, 0.4, 0.5, @@ -194,23 +178,9 @@ def test_cyclic_float_axis_two_range_loops(self): 0.8, 0.9, 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, - 2.1, - 2.2, - 2.3, - 2.4, - 2.5, - 2.6, - 2.7, + 0.0, + 0.1, + 0.2, 0.3, 0.4, 0.5, @@ -219,23 +189,6 @@ def test_cyclic_float_axis_two_range_loops(self): 0.8, 0.9, 1.0, - 1.1, - 1.2, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.8, - 1.9, - 2.0, - 2.1, - 2.2, - 2.3, - 2.4, - 2.5, - 2.6, - 2.7, ] def test_cyclic_float_axis_below_axis_range(self): @@ -245,7 +198,7 @@ def test_cyclic_float_axis_below_axis_range(self): result = self.API.retrieve(request) # result.pprint() assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [-0.7, -0.6, -0.5, -0.4, -0.3, -0.7, -0.6, -0.5, -0.4, -0.3] + assert [leaf.value for leaf in result.leaves] == [0.3, 0.4, 0.5, 0.6, 0.7, 0.3, 0.4, 0.5, 0.6, 0.7] def test_cyclic_float_axis_below_axis_range_crossing_seam(self): request = Request( @@ -253,30 +206,28 @@ def test_cyclic_float_axis_below_axis_range_crossing_seam(self): ) result = self.API.retrieve(request) # result.pprint() - assert len(result.leaves) == 22 + assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - -0.7, - -0.6, - -0.5, - -0.4, - -0.3, - -0.2, - -0.1, 0.0, 0.1, 0.2, 0.3, - -0.7, - -0.6, - -0.5, - -0.4, - -0.3, - -0.2, - -0.1, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, 0.0, 0.1, 0.2, 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.8, + 0.9, ] def test_cyclic_float_axis_reversed(self): @@ -286,34 +237,34 @@ def test_cyclic_float_axis_reversed(self): result = self.API.retrieve(request) # result.pprint() assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [1.3, 1.4, 1.5, 1.6, 1.7, 1.3, 1.4, 1.5, 1.6, 1.7] + assert [leaf.value for leaf in result.leaves] == [0.3, 0.4, 0.5, 0.6, 0.7, 0.3, 0.4, 0.5, 0.6, 0.7] def test_two_cyclic_axis_wrong_axis_order(self): request = Request(Box(["step", "long", "level"], [0, 1.3, 131], [3, 1.7, 132]), Select("date", ["2000-01-01"])) result = self.API.retrieve(request) - # result.pprint() + result.pprint() assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, ] def test_two_cyclic_axis(self): @@ -322,26 +273,26 @@ def test_two_cyclic_axis(self): # result.pprint() assert len(result.leaves) == 20 assert [leaf.value for leaf in result.leaves] == [ - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, - 1.3, - 1.4, - 1.5, - 1.6, - 1.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, + 0.3, + 0.4, + 0.5, + 0.6, + 0.7, ] def test_select_cyclic_float_axis_edge(self): diff --git a/tests/test_cyclic_simple.py b/tests/test_cyclic_simple.py index be48cb8b9..ec213e2c8 100644 --- a/tests/test_cyclic_simple.py +++ b/tests/test_cyclic_simple.py @@ -33,5 +33,34 @@ def test_cyclic_float_axis_across_seam(self): Box(["step", "long"], [0, 0.9], [0, 1.2]), Select("date", ["2000-01-01"]), Select("level", [128]) ) result = self.API.retrieve(request) + result.pprint() assert len(result.leaves) == 4 - assert [leaf.value for leaf in result.leaves] == [0.9, 1.0, 1.1, 1.2] + assert [leaf.value for leaf in result.leaves] == [0.1, 0.2, 0.9, 1.0] + + def test_cyclic_float_surrounding(self): + request = Request( + Select("step", [0]), + Select("long", [1.], method="surrounding"), + Select("date", ["2000-01-01"]), + Select("level", [128]) + ) + result = self.API.retrieve(request) + result.pprint() + for leaf in result.leaves: + path = leaf.flatten() + lon_val = path["long"] + assert lon_val in [0., 0.1, 0.9, 1.] + + def test_cyclic_float_surrounding_below_seam(self): + request = Request( + Select("step", [0]), + Select("long", [0.], method="surrounding"), + Select("date", ["2000-01-01"]), + Select("level", [128]) + ) + result = self.API.retrieve(request) + result.pprint() + for leaf in result.leaves: + path = leaf.flatten() + lon_val = path["long"] + assert lon_val in [0., 0.1, 0.9] diff --git a/tests/test_cyclic_snapping.py b/tests/test_cyclic_snapping.py index ac1b11c7f..b0a6580ba 100644 --- a/tests/test_cyclic_snapping.py +++ b/tests/test_cyclic_snapping.py @@ -30,7 +30,7 @@ def test_cyclic_float_axis_across_seam(self): result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 2 - assert result.leaves[0].flatten()["long"] == -0.5 - assert result.leaves[1].flatten()["long"] == 0.0 - assert result.leaves[0].result == (None, 1) - assert result.leaves[1].result == (None, 0) + assert result.leaves[0].flatten()["long"] == 0.0 + assert result.leaves[1].flatten()["long"] == 0.5 + assert result.leaves[0].result == (None, 0) + assert result.leaves[1].result == (None, 1) From abc896b4e4925ec63306d7f587048ce749da177c Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 7 Sep 2023 10:31:08 +0200 Subject: [PATCH 07/11] make surrounding method work with the cyclic axis --- polytope/datacube/backends/datacube.py | 11 ----------- polytope/datacube/datacube_axis.py | 14 ++++---------- polytope/engine/hullslicer.py | 3 +-- tests/test_cyclic_simple.py | 2 +- tests/test_merge_octahedral_one_axis.py | 3 ++- tests/test_snapping.py | 10 ++++++++++ tests/test_snapping_real_data.py | 4 ++-- 7 files changed, 20 insertions(+), 27 deletions(-) diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index 0bc710c93..8510a68f9 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -1,5 +1,4 @@ import importlib -import itertools import math from abc import ABC, abstractmethod from typing import Any @@ -101,15 +100,10 @@ def get_indices(self, path: DatacubePath, axis, lower, upper, method=None): if offset is not None: # Note that we can only do unique if not dealing with time values idx_between = unique(idx_between) - print("inside get indices") - print(idx_between) return idx_between def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, method): - print("inside look up datacube") idx_between = [] - print(axis.name) - print(search_ranges) for i in range(len(search_ranges)): r = search_ranges[i] offset = search_ranges_offset[i] @@ -142,8 +136,6 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) idx_between.append(indexes_between[j][k]) - print("indexes between") - print(idx_between) return idx_between def get_mapper(self, axis): @@ -157,9 +149,6 @@ def remap_path(self, path: DatacubePath): value = path[key] path[key] = self._axes[key].remap([value, value])[0][0] return path - - # def remap_tree(self, tree: IndexTree): - # for @staticmethod def create(datacube, axis_options: dict): diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 4c7a66d50..3efa48573 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -146,6 +146,7 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): return old_find_indices_between(index_ranges, low, up, datacube, method, offset) if offset != 0: + # print("HERE OFFSET NOT 0") # NOTE that if the offset is not 0, then we need to recenter the low and up # values to be within the datacube range new_offset = 0 @@ -154,21 +155,17 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): # low = low - range_length # new_offset -= range_length # up = up - range_length - print("INSIDE THE DATACUBE AXIS DECORATOR") - print((low, up)) if method == "surrounding": for indexes in index_ranges: if cls.name in datacube.complete_axes: start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") - print(start, end) if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) start_offset_indicator = "need_offset" - index_val_found = indexes[-2:][0] + index_val_found = indexes[-1:][0] indexes_between_before = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_before) if end+1 > len(indexes): - print("end exceeded indices length") start_offset_indicator = "need_offset" index_val_found = indexes[:2][1] indexes_between_after = [start_offset_indicator, new_offset, index_val_found] @@ -200,20 +197,17 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): return old_find_indices_between(index_ranges, low, up, datacube, method, offset) else: + # print("HERE OFFSET IS 0") # If the offset is 0, then the first value found on the left has an offset of range_length new_offset = 0 if method == "surrounding": - print("no offset") for indexes in index_ranges: if cls.name in datacube.complete_axes: start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") - print((start, end)) if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) - print("low start") start_offset_indicator = "need_offset" - index_val_found = indexes[-2:][0] - print(index_val_found) + index_val_found = indexes[-1:][0] new_offset = -range_length indexes_between_before = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_before) diff --git a/polytope/engine/hullslicer.py b/polytope/engine/hullslicer.py index 5553f4e5a..39a88d350 100644 --- a/polytope/engine/hullslicer.py +++ b/polytope/engine/hullslicer.py @@ -1,7 +1,7 @@ +import math from copy import copy from itertools import chain from typing import List -import math import scipy.spatial @@ -58,7 +58,6 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, lower, upper, nex if ax.is_cyclic: remapped_val = (ax.remap([value, value])[0][0] + ax.remap([value, value])[0][1])/2 remapped_val = round(remapped_val, int(-math.log10(ax.tol))) - print(remapped_val) # child = node.create_child(ax, value) child = node.create_child(ax, remapped_val) child["unsliced_polytopes"] = copy(node["unsliced_polytopes"]) diff --git a/tests/test_cyclic_simple.py b/tests/test_cyclic_simple.py index ec213e2c8..f7dda6dd8 100644 --- a/tests/test_cyclic_simple.py +++ b/tests/test_cyclic_simple.py @@ -63,4 +63,4 @@ def test_cyclic_float_surrounding_below_seam(self): for leaf in result.leaves: path = leaf.flatten() lon_val = path["long"] - assert lon_val in [0., 0.1, 0.9] + assert lon_val in [0., 0.1, 0.9, 1.] diff --git a/tests/test_merge_octahedral_one_axis.py b/tests/test_merge_octahedral_one_axis.py index 9739a6c79..aec9ee7c6 100644 --- a/tests/test_merge_octahedral_one_axis.py +++ b/tests/test_merge_octahedral_one_axis.py @@ -34,4 +34,5 @@ def test_merge_axis(self): ) result = self.API.retrieve(request) # result.pprint() - assert result.leaves[2].flatten()["longitude"] == 360.0 + assert result.leaves[-1].flatten()["longitude"] == 360.0 + assert result.leaves[0].flatten()["longitude"] == 0.070093457944 diff --git a/tests/test_snapping.py b/tests/test_snapping.py index 7d127a877..83d472203 100644 --- a/tests/test_snapping.py +++ b/tests/test_snapping.py @@ -77,3 +77,13 @@ def test_1D_nonexisting_point_v2(self): result.pprint() assert len(result.leaves) == 1 assert result.is_root() + + def test_1D_nonexisting_point_surrounding(self): + request = Request(Select("level", [0], method="surrounding"), Select("step", [6], method="surrounding")) + result = self.API.retrieve(request) + result.pprint() + assert len(result.leaves) == 1 + for leaf in result.leaves: + path = leaf.flatten() + assert path["level"] == 1 + assert path["step"] == 5 diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py index 247a276b8..bfe32b21d 100644 --- a/tests/test_snapping_real_data.py +++ b/tests/test_snapping_real_data.py @@ -21,7 +21,7 @@ def setup_method(self, method): def test_surrounding_on_grid_point(self): requested_lat = 0 - requested_lon = 0 + requested_lon = -720 request = Request( Box(["number", "isobaricInhPa"], [6, 500.0], [6, 850.0]), Select("time", ["2017-01-02T12:00:00"]), @@ -58,4 +58,4 @@ def test_surrounding_on_grid_point(self): # plt.colorbar(label="Temperature") # plt.show() for lon in longs: - assert lon in [-3, 0, 3] + assert lon in [357, 0, 3] From c03b7a5686dee321d3bf8c7932bd6fad8a900b17 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 7 Sep 2023 10:37:23 +0200 Subject: [PATCH 08/11] clean up --- polytope/datacube/datacube_axis.py | 59 +----------------------------- 1 file changed, 2 insertions(+), 57 deletions(-) diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 3efa48573..5467ec60a 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -139,22 +139,15 @@ def unmap_to_datacube(path, unmapped_path): def find_indices_between(index_ranges, low, up, datacube, offset, method=None): update_range() - range_length = cls.range[1] - cls.range[0] indexes_between_ranges = [] if method != "surrounding": return old_find_indices_between(index_ranges, low, up, datacube, method, offset) - if offset != 0: - # print("HERE OFFSET NOT 0") + if True: # NOTE that if the offset is not 0, then we need to recenter the low and up # values to be within the datacube range new_offset = 0 - # if low <= cls.range[0] + cls.tol: - # while low >= cls.range[0] - cls.tol: - # low = low - range_length - # new_offset -= range_length - # up = up - range_length if method == "surrounding": for indexes in index_ranges: if cls.name in datacube.complete_axes: @@ -185,7 +178,7 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): indexes_between_ranges.append(indexes_between_before) if end+1 > len(indexes): start_offset_indicator = "need_offset" - index_val_found = indexes[:1][0] + index_val_found = indexes[:2][0] indexes_between_after = [start_offset_indicator, new_offset, index_val_found] indexes_between_ranges.append(indexes_between_after) start = max(start-1, 0) @@ -193,54 +186,6 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): indexes_between = indexes[start:end] indexes_between_ranges.append(indexes_between) return indexes_between_ranges - else: - return old_find_indices_between(index_ranges, low, up, datacube, method, offset) - - else: - # print("HERE OFFSET IS 0") - # If the offset is 0, then the first value found on the left has an offset of range_length - new_offset = 0 - if method == "surrounding": - for indexes in index_ranges: - if cls.name in datacube.complete_axes: - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) - start_offset_indicator = "need_offset" - index_val_found = indexes[-1:][0] - new_offset = -range_length - indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_before) - if end+1 > len(indexes): - start_offset_indicator = "need_offset" - index_val_found = indexes[:2][1] - indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_after) - start = max(start-1, 0) - end = min(end+1, len(indexes)) - indexes_between = indexes[start:end].to_list() - indexes_between = [i for i in indexes_between] - indexes_between_ranges.append(indexes_between) - else: - start = indexes.index(low) - end = indexes.index(up) - if start-1 < 0: - start_offset_indicator = "need_offset" - index_val_found = indexes[-1:][0] - indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_before) - if end+1 > len(indexes): - start_offset_indicator = "need_offset" - index_val_found = indexes[:1][0] - indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_after) - start = max(start-1, 0) - end = min(end+1, len(indexes)) - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - else: - return old_find_indices_between(index_ranges, low, up, datacube, method, offset) def offset(range): # We first unpad the range by the axis tolerance to make sure that From 78306d00020bd260a4a2ada7532a270c946d50c5 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 7 Sep 2023 10:50:26 +0200 Subject: [PATCH 09/11] refactor find_indices_between for cyclic axis --- polytope/datacube/datacube_axis.py | 68 ++++++++++++------------------ 1 file changed, 26 insertions(+), 42 deletions(-) diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 5467ec60a..80d6920e1 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -143,49 +143,33 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): if method != "surrounding": return old_find_indices_between(index_ranges, low, up, datacube, method, offset) - - if True: - # NOTE that if the offset is not 0, then we need to recenter the low and up - # values to be within the datacube range + else: new_offset = 0 - if method == "surrounding": - for indexes in index_ranges: - if cls.name in datacube.complete_axes: - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - if start-1 < 0: # NOTE TODO: here the boundaries will not necessarily be 0 or len(indexes) - start_offset_indicator = "need_offset" - index_val_found = indexes[-1:][0] - indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_before) - if end+1 > len(indexes): - start_offset_indicator = "need_offset" - index_val_found = indexes[:2][1] - indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_after) - start = max(start-1, 0) - end = min(end+1, len(indexes)) - indexes_between = indexes[start:end].to_list() - indexes_between = [i - new_offset for i in indexes_between] - indexes_between_ranges.append(indexes_between) - else: - start = indexes.index(low) - end = indexes.index(up) - if start-1 < 0: - start_offset_indicator = "need_offset" - index_val_found = indexes[-1:][0] - indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_before) - if end+1 > len(indexes): - start_offset_indicator = "need_offset" - index_val_found = indexes[:2][0] - indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_after) - start = max(start-1, 0) - end = min(end+1, len(indexes)) - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges + for indexes in index_ranges: + if cls.name in datacube.complete_axes: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + else: + start = indexes.index(low) + end = indexes.index(up) + if start-1 < 0: + start_offset_indicator = "need_offset" + index_val_found = indexes[-1:][0] + indexes_between_before = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_before) + if end+1 > len(indexes): + start_offset_indicator = "need_offset" + index_val_found = indexes[:2][0] + indexes_between_after = [start_offset_indicator, new_offset, index_val_found] + indexes_between_ranges.append(indexes_between_after) + start = max(start-1, 0) + end = min(end+1, len(indexes)) + if cls.name in datacube.complete_axes: + indexes_between = indexes[start:end].to_list() + else: + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges def offset(range): # We first unpad the range by the axis tolerance to make sure that From f04ded57f93c06e01c963ca17a3221b663af3030 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 7 Sep 2023 10:58:22 +0200 Subject: [PATCH 10/11] refactor the surrounding method code in datacube.py _look_up_datacube and find_indices_between in cyclic datacube axis decorator --- polytope/datacube/backends/datacube.py | 27 +++++++------------------- polytope/datacube/datacube_axis.py | 23 +++++++++------------- 2 files changed, 16 insertions(+), 34 deletions(-) diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index 8510a68f9..96554550f 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -109,7 +109,7 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, offset = search_ranges_offset[i] low = r[0] up = r[1] - indexes_between = axis.find_indices_between([indexes], low, up, self, offset, method) + indexes_between = axis.find_indices_between([indexes], low, up, self, method) # Now the indexes_between are values on the cyclic range so need to remap them to their original # values before returning them for j in range(len(indexes_between)): @@ -117,25 +117,12 @@ def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, if len(indexes_between[j]) == 0: idx_between = idx_between else: - if indexes_between[j][0] == "need_offset": - for k in range(2, len(indexes_between[j])): - if offset is None: - indexes_between[j][k] = indexes_between[j][k] - else: - offset = offset - indexes_between[j][k] = round(indexes_between[j][k] + offset, - int(-math.log10(axis.tol))) - idx_between.append(indexes_between[j][k]) - else: - # do normal offset if no new offset - for k in range(len(indexes_between[j])): - if offset is None: - indexes_between[j][k] = indexes_between[j][k] - else: - offset = offset - indexes_between[j][k] = round(indexes_between[j][k] + offset, - int(-math.log10(axis.tol))) - idx_between.append(indexes_between[j][k]) + for k in range(len(indexes_between[j])): + if offset is None: + indexes_between[j][k] = indexes_between[j][k] + else: + indexes_between[j][k] = round(indexes_between[j][k] + offset, int(-math.log10(axis.tol))) + idx_between.append(indexes_between[j][k]) return idx_between def get_mapper(self, axis): diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 80d6920e1..736349d22 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -137,14 +137,13 @@ def unmap_to_datacube(path, unmapped_path): old_find_indices_between = cls.find_indices_between - def find_indices_between(index_ranges, low, up, datacube, offset, method=None): + def find_indices_between(index_ranges, low, up, datacube, method=None): update_range() indexes_between_ranges = [] if method != "surrounding": - return old_find_indices_between(index_ranges, low, up, datacube, method, offset) + return old_find_indices_between(index_ranges, low, up, datacube, method) else: - new_offset = 0 for indexes in index_ranges: if cls.name in datacube.complete_axes: start = indexes.searchsorted(low, "left") @@ -153,15 +152,11 @@ def find_indices_between(index_ranges, low, up, datacube, offset, method=None): start = indexes.index(low) end = indexes.index(up) if start-1 < 0: - start_offset_indicator = "need_offset" index_val_found = indexes[-1:][0] - indexes_between_before = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_before) + indexes_between_ranges.append([index_val_found]) if end+1 > len(indexes): - start_offset_indicator = "need_offset" index_val_found = indexes[:2][0] - indexes_between_after = [start_offset_indicator, new_offset, index_val_found] - indexes_between_ranges.append(indexes_between_after) + indexes_between_ranges.append([index_val_found]) start = max(start-1, 0) end = min(end+1, len(indexes)) if cls.name in datacube.complete_axes: @@ -269,7 +264,7 @@ def unmap_total_path_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, offset, method=None): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -352,7 +347,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, offset, method=None): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -408,7 +403,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, offset, method=None): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -517,7 +512,7 @@ def unmap_to_datacube(path, unmapped_path): def remap_to_requested(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(index_ranges, low, up, datacube, offset, method=None): + def find_indices_between(index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for transform in cls.transformations: @@ -622,7 +617,7 @@ def unmap_total_path_to_datacube(self, path, unmapped_path): def remap_to_requeest(path, unmapped_path): return (path, unmapped_path) - def find_indices_between(self, index_ranges, low, up, datacube, offset, method=None): + def find_indices_between(self, index_ranges, low, up, datacube, method=None): # TODO: add method for snappping indexes_between_ranges = [] for indexes in index_ranges: From f82c8e84c03f5a24337b8fbdaebaa4e068560985 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Thu, 7 Sep 2023 11:01:59 +0200 Subject: [PATCH 11/11] black --- polytope/datacube/backends/FDB_datacube.py | 1 - polytope/datacube/backends/datacube.py | 10 +++-- polytope/datacube/datacube_axis.py | 42 ++++++++++--------- polytope/datacube/index_tree.py | 4 +- .../datacube_transformations.py | 12 +++--- polytope/engine/hullslicer.py | 2 +- polytope/shapes.py | 1 - tests/test_cyclic_axis_over_negative_vals.py | 7 ++-- tests/test_cyclic_axis_slicer_not_0.py | 12 +----- tests/test_cyclic_axis_slicing.py | 4 +- tests/test_cyclic_simple.py | 12 +++--- tests/test_cyclic_snapping.py | 4 +- tests/test_datacube_axes_init.py | 36 +++++++++------- tests/test_snapping_real_data.py | 6 ++- 14 files changed, 76 insertions(+), 77 deletions(-) diff --git a/polytope/datacube/backends/FDB_datacube.py b/polytope/datacube/backends/FDB_datacube.py index 46e2c673c..764ae56c8 100644 --- a/polytope/datacube/backends/FDB_datacube.py +++ b/polytope/datacube/backends/FDB_datacube.py @@ -19,7 +19,6 @@ def update_fdb_dataarray(fdb_dataarray): class FDBDatacube(Datacube): - def __init__(self, config={}, axis_options={}): self.axis_options = axis_options self.grid_mapper = None diff --git a/polytope/datacube/backends/datacube.py b/polytope/datacube/backends/datacube.py index 96554550f..37e7faa8f 100644 --- a/polytope/datacube/backends/datacube.py +++ b/polytope/datacube/backends/datacube.py @@ -29,10 +29,12 @@ def validate(self, axes): def _create_axes(self, name, values, transformation_type_key, transformation_options): # first check what the final axes are for this axis name given transformations - final_axis_names = DatacubeAxisTransformation.get_final_axes(name, transformation_type_key, - transformation_options) - transformation = DatacubeAxisTransformation.create_transform(name, transformation_type_key, - transformation_options) + final_axis_names = DatacubeAxisTransformation.get_final_axes( + name, transformation_type_key, transformation_options + ) + transformation = DatacubeAxisTransformation.create_transform( + name, transformation_type_key, transformation_options + ) for blocked_axis in transformation.blocked_axes(): self.blocked_axes.append(blocked_axis) for axis_name in final_axis_names: diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 736349d22..9582bd0b5 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -151,14 +151,14 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): else: start = indexes.index(low) end = indexes.index(up) - if start-1 < 0: + if start - 1 < 0: index_val_found = indexes[-1:][0] indexes_between_ranges.append([index_val_found]) - if end+1 > len(indexes): + if end + 1 > len(indexes): index_val_found = indexes[:2][0] indexes_between_ranges.append([index_val_found]) - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) if cls.name in datacube.complete_axes: indexes_between = indexes[start:end].to_list() else: @@ -190,6 +190,7 @@ def mapper(cls): from .transformations.datacube_mappers import DatacubeMapper if cls.has_mapper: + def find_indexes(path, datacube): # first, find the relevant transformation object that is a mapping in the cls.transformation dico for transform in cls.transformations: @@ -275,8 +276,8 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): if method == "surrounding": start = idxs.index(low) end = idxs.index(up) - start = max(start-1, 0) - end = min(end+1, len(idxs)) + start = max(start - 1, 0) + end = min(end + 1, len(idxs)) # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = idxs[start:end] indexes_between_ranges.append(indexes_between) @@ -306,6 +307,7 @@ def merge(cls): from .transformations.datacube_merger import DatacubeAxisMerger if cls.has_merger: + def find_indexes(path, datacube): # first, find the relevant transformation object that is a mapping in the cls.transformation dico for transform in cls.transformations: @@ -358,8 +360,8 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): if method == "surrounding": start = indexes.index(low) end = indexes.index(up) - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = indexes[start:end] indexes_between_ranges.append(indexes_between) @@ -387,6 +389,7 @@ def reverse(cls): from .transformations.datacube_reverse import DatacubeAxisReverse if cls.reorder: + def find_indexes(path, datacube): # first, find the relevant transformation object that is a mapping in the cls.transformation dico subarray = datacube.dataarray.sel(path, method="nearest") @@ -419,8 +422,8 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): if method == "surrounding": start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) indexes_between = indexes[start:end].to_list() indexes_between_ranges.append(indexes_between) else: @@ -432,8 +435,8 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): if method == "surrounding": start = indexes.index(low) end = indexes.index(up) - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = indexes[start:end] indexes_between_ranges.append(indexes_between) @@ -471,7 +474,6 @@ def type_change(cls): from .transformations.datacube_type_change import DatacubeAxisTypeChange if cls.type_change: - old_find_indexes = cls.find_indexes def find_indexes(path, datacube): @@ -523,8 +525,8 @@ def find_indices_between(index_ranges, low, up, datacube, method=None): if method == "surrounding": start = indexes.index(low) end = indexes.index(up) - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = indexes[start:end] indexes_between_ranges.append(indexes_between) @@ -628,8 +630,8 @@ def find_indices_between(self, index_ranges, low, up, datacube, method=None): if method == "surrounding": start = indexes.searchsorted(low, "left") end = indexes.searchsorted(up, "right") - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) indexes_between = indexes[start:end].to_list() indexes_between_ranges.append(indexes_between) else: @@ -641,8 +643,8 @@ def find_indices_between(self, index_ranges, low, up, datacube, method=None): if method == "surrounding": start = indexes.index(low) end = indexes.index(up) - start = max(start-1, 0) - end = min(end+1, len(indexes)) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) # indexes_between = [i for i in indexes if low <= i <= up] indexes_between = indexes[start:end] indexes_between_ranges.append(indexes_between) @@ -699,7 +701,7 @@ class FloatDatacubeAxis(DatacubeAxis): tol = 1e-12 range = None transformations = [] - type = 0. + type = 0.0 def parse(self, value: Any) -> Any: return float(value) diff --git a/polytope/datacube/index_tree.py b/polytope/datacube/index_tree.py index 2d285327e..b7a37aff5 100644 --- a/polytope/datacube/index_tree.py +++ b/polytope/datacube/index_tree.py @@ -63,9 +63,9 @@ def __eq__(self, other): else: if other.value == self.value: return True - if other.value - 2*other.axis.tol <= self.value <= other.value + 2*other.axis.tol: + if other.value - 2 * other.axis.tol <= self.value <= other.value + 2 * other.axis.tol: return True - elif self.value - 2*self.axis.tol <= other.value <= self.value + 2*self.axis.tol: + elif self.value - 2 * self.axis.tol <= other.value <= self.value + 2 * self.axis.tol: return True else: return False diff --git a/polytope/datacube/transformations/datacube_transformations.py b/polytope/datacube/transformations/datacube_transformations.py index d0af6dcec..4afac6ce8 100644 --- a/polytope/datacube/transformations/datacube_transformations.py +++ b/polytope/datacube/transformations/datacube_transformations.py @@ -4,7 +4,6 @@ class DatacubeAxisTransformation(ABC): - @staticmethod def create_transform(name, transformation_type_key, transformation_options): transformation_type = _type_to_datacube_transformation_lookup[transformation_type_key] @@ -20,8 +19,9 @@ def create_transform(name, transformation_type_key, transformation_options): @staticmethod def get_final_axes(name, transformation_type_key, transformation_options): - new_transformation = DatacubeAxisTransformation.create_transform(name, transformation_type_key, - transformation_options) + new_transformation = DatacubeAxisTransformation.create_transform( + name, transformation_type_key, transformation_options + ) transformation_axis_names = new_transformation.transformation_axes_final() return transformation_axis_names @@ -102,8 +102,8 @@ def change_val_type(self, axis_name, values): has_transform = { "mapper": "has_mapper", - "cyclic" : "is_cyclic", - "merge" : "has_merger", + "cyclic": "is_cyclic", + "merge": "has_merger", "reverse": "reorder", - "type_change": "type_change" + "type_change": "type_change", } diff --git a/polytope/engine/hullslicer.py b/polytope/engine/hullslicer.py index 39a88d350..5ea0c951e 100644 --- a/polytope/engine/hullslicer.py +++ b/polytope/engine/hullslicer.py @@ -56,7 +56,7 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, lower, upper, nex # remapped_val = (ax.remap([value, value])[0][0] + ax.remap([value, value])[0][1])/2 remapped_val = value if ax.is_cyclic: - remapped_val = (ax.remap([value, value])[0][0] + ax.remap([value, value])[0][1])/2 + remapped_val = (ax.remap([value, value])[0][0] + ax.remap([value, value])[0][1]) / 2 remapped_val = round(remapped_val, int(-math.log10(ax.tol))) # child = node.create_child(ax, value) child = node.create_child(ax, remapped_val) diff --git a/polytope/shapes.py b/polytope/shapes.py index e0d0b0bd4..94562abc5 100644 --- a/polytope/shapes.py +++ b/polytope/shapes.py @@ -23,7 +23,6 @@ def axes(self) -> List[str]: class ConvexPolytope(Shape): - def __init__(self, axes, points, method=None): self._axes = list(axes) self.points = points diff --git a/tests/test_cyclic_axis_over_negative_vals.py b/tests/test_cyclic_axis_over_negative_vals.py index a67333cce..ccb0bbaa0 100644 --- a/tests/test_cyclic_axis_over_negative_vals.py +++ b/tests/test_cyclic_axis_over_negative_vals.py @@ -58,7 +58,7 @@ def test_cyclic_float_axis_across_seam(self): -0.5, -0.4, -0.3, - -0.2 + -0.2, ] def test_cyclic_float_axis_inside_cyclic_range(self): @@ -94,8 +94,7 @@ def test_cyclic_float_axis_above_axis_range(self): result = self.API.retrieve(request) # result.pprint() assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [-0.7, -0.6, -0.5, -0.4, -0.3, - -0.7, -0.6, -0.5, -0.4, -0.3] + assert [leaf.value for leaf in result.leaves] == [-0.7, -0.6, -0.5, -0.4, -0.3, -0.7, -0.6, -0.5, -0.4, -0.3] def test_cyclic_float_axis_two_range_loops(self): request = Request( @@ -124,7 +123,7 @@ def test_cyclic_float_axis_two_range_loops(self): -0.5, -0.4, -0.3, - -0.2 + -0.2, ] def test_cyclic_float_axis_below_axis_range(self): diff --git a/tests/test_cyclic_axis_slicer_not_0.py b/tests/test_cyclic_axis_slicer_not_0.py index 33ef062e6..daba54a86 100644 --- a/tests/test_cyclic_axis_slicer_not_0.py +++ b/tests/test_cyclic_axis_slicer_not_0.py @@ -92,17 +92,7 @@ def test_cyclic_float_axis_above_axis_range(self): ) result = self.API.retrieve(request) assert len(result.leaves) == 10 - assert [leaf.value for leaf in result.leaves] == [ - -0.7, - -0.6, - -0.5, - -0.4, - -0.3, - -0.7, - -0.6, - -0.5, - -0.4, - -0.3] + assert [leaf.value for leaf in result.leaves] == [-0.7, -0.6, -0.5, -0.4, -0.3, -0.7, -0.6, -0.5, -0.4, -0.3] def test_cyclic_float_axis_two_range_loops(self): request = Request( diff --git a/tests/test_cyclic_axis_slicing.py b/tests/test_cyclic_axis_slicing.py index 3d3c68465..554ae7a85 100644 --- a/tests/test_cyclic_axis_slicing.py +++ b/tests/test_cyclic_axis_slicing.py @@ -45,7 +45,7 @@ def test_cyclic_float_axis_across_seam(self): 0.7, 0.8, 0.9, - 1., + 1.0, 0.1, 0.2, 0.3, @@ -55,7 +55,7 @@ def test_cyclic_float_axis_across_seam(self): 0.7, 0.8, 0.9, - 1., + 1.0, ] def test_cyclic_float_axis_across_seam_repeated(self): diff --git a/tests/test_cyclic_simple.py b/tests/test_cyclic_simple.py index f7dda6dd8..bf6660ee7 100644 --- a/tests/test_cyclic_simple.py +++ b/tests/test_cyclic_simple.py @@ -40,27 +40,27 @@ def test_cyclic_float_axis_across_seam(self): def test_cyclic_float_surrounding(self): request = Request( Select("step", [0]), - Select("long", [1.], method="surrounding"), + Select("long", [1.0], method="surrounding"), Select("date", ["2000-01-01"]), - Select("level", [128]) + Select("level", [128]), ) result = self.API.retrieve(request) result.pprint() for leaf in result.leaves: path = leaf.flatten() lon_val = path["long"] - assert lon_val in [0., 0.1, 0.9, 1.] + assert lon_val in [0.0, 0.1, 0.9, 1.0] def test_cyclic_float_surrounding_below_seam(self): request = Request( Select("step", [0]), - Select("long", [0.], method="surrounding"), + Select("long", [0.0], method="surrounding"), Select("date", ["2000-01-01"]), - Select("level", [128]) + Select("level", [128]), ) result = self.API.retrieve(request) result.pprint() for leaf in result.leaves: path = leaf.flatten() lon_val = path["long"] - assert lon_val in [0., 0.1, 0.9, 1.] + assert lon_val in [0.0, 0.1, 0.9, 1.0] diff --git a/tests/test_cyclic_snapping.py b/tests/test_cyclic_snapping.py index b0a6580ba..ab6189ba3 100644 --- a/tests/test_cyclic_snapping.py +++ b/tests/test_cyclic_snapping.py @@ -24,9 +24,7 @@ def setup_method(self, method): # Testing different shapes def test_cyclic_float_axis_across_seam(self): - request = Request( - Select("long", [-0.2], method="surrounding") - ) + request = Request(Select("long", [-0.2], method="surrounding")) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 2 diff --git a/tests/test_datacube_axes_init.py b/tests/test_datacube_axes_init.py index bea08086b..decee5bb5 100644 --- a/tests/test_datacube_axes_init.py +++ b/tests/test_datacube_axes_init.py @@ -29,27 +29,35 @@ def test_created_axes(self): assert self.datacube._axes["longitude"].has_mapper assert type(self.datacube._axes["longitude"]) == FloatDatacubeAxis assert not ("values" in self.datacube._axes.keys()) - assert self.datacube._axes["latitude"].find_indexes({}, - self.datacube)[:5] == [89.94618771566562, - 89.87647835333229, - 89.80635731954224, - 89.73614327160958, - 89.6658939412157] - assert self.datacube._axes["longitude"].find_indexes({"latitude": 89.94618771566562}, - self.datacube)[:8] == [0.0, 18.0, 36.0, 54.0, - 72.0, 90.0, 108.0, 126.0] - assert len(self.datacube._axes["longitude"].find_indexes({"latitude": 89.94618771566562}, - self.datacube)) == 20 + assert self.datacube._axes["latitude"].find_indexes({}, self.datacube)[:5] == [ + 89.94618771566562, + 89.87647835333229, + 89.80635731954224, + 89.73614327160958, + 89.6658939412157, + ] + assert self.datacube._axes["longitude"].find_indexes({"latitude": 89.94618771566562}, self.datacube)[:8] == [ + 0.0, + 18.0, + 36.0, + 54.0, + 72.0, + 90.0, + 108.0, + 126.0, + ] + assert len(self.datacube._axes["longitude"].find_indexes({"latitude": 89.94618771566562}, self.datacube)) == 20 lon_ax = self.datacube._axes["longitude"] lat_ax = self.datacube._axes["latitude"] (path, unmapped_path) = lat_ax.unmap_to_datacube({"latitude": 89.94618771566562}, {}) assert path == {} - assert unmapped_path == {'latitude': 89.94618771566562} + assert unmapped_path == {"latitude": 89.94618771566562} (path, unmapped_path) = lon_ax.unmap_to_datacube({"longitude": 0.0}, {"latitude": 89.94618771566562}) assert path == {} assert unmapped_path == {"values": 0} - assert lat_ax.find_indices_between([[89.94618771566562, 89.87647835333229]], - 89.87, 90, self.datacube, 0) == [[89.94618771566562, 89.87647835333229]] + assert lat_ax.find_indices_between([[89.94618771566562, 89.87647835333229]], 89.87, 90, self.datacube, 0) == [ + [89.94618771566562, 89.87647835333229] + ] def test_mapper_transformation_request(self): request = Request( diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py index bfe32b21d..4cf4e01f7 100644 --- a/tests/test_snapping_real_data.py +++ b/tests/test_snapping_real_data.py @@ -15,8 +15,10 @@ def setup_method(self, method): array = ds.to_xarray().isel(step=0).t self.xarraydatacube = XArrayDatacube(array) self.slicer = HullSlicer() - options = {"latitude": {"transformation": {"reverse": {True}}}, - "longitude": {"transformation": {"cyclic": [0, 360.0]}}} + options = { + "latitude": {"transformation": {"reverse": {True}}}, + "longitude": {"transformation": {"cyclic": [0, 360.0]}}, + } self.API = Polytope(datacube=array, engine=self.slicer, axis_options=options) def test_surrounding_on_grid_point(self):