diff --git a/polytope_feature/engine/hullslicer.py b/polytope_feature/engine/hullslicer.py index 9a99682a..7302931d 100644 --- a/polytope_feature/engine/hullslicer.py +++ b/polytope_feature/engine/hullslicer.py @@ -115,6 +115,8 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_node fvalue = ax.to_float(value) new_polytope = slice(polytope, ax.name, fvalue, slice_axis_idx) remapped_val = self.remap_values(ax, value) + # is_last_axis = ax==list(datacube.axes.values())[-1] + # print(ax==list(datacube.axes.values())[-1]) (child, next_nodes) = node.create_child(ax, remapped_val, next_nodes) child["unsliced_polytopes"] = copy(node["unsliced_polytopes"]) child["unsliced_polytopes"].remove(polytope) @@ -203,14 +205,15 @@ def remove_compressed_axis_in_union(self, polytopes): if p.is_in_union: for axis in p.axes(): if axis == self.compressed_axes[-1]: - self.compressed_axes.remove(axis) + # self.compressed_axes.remove(axis) + pass def extract(self, datacube: Datacube, polytopes: List[ConvexPolytope]): # Determine list of axes to compress self.find_compressed_axes(datacube, polytopes) # remove compressed axes which are in a union - self.remove_compressed_axis_in_union(polytopes) + # self.remove_compressed_axis_in_union(polytopes) # Convert the polytope points to float type to support triangulation and interpolation for p in polytopes: @@ -224,7 +227,6 @@ def extract(self, datacube: Datacube, polytopes: List[ConvexPolytope]): # NOTE: could optimise here if we know combinations will always be for one request. # Then we do not need to create a new index tree and merge it to request, but can just # directly work on request and return it... - for c in combinations: r = TensorIndexTree() new_c = [] diff --git a/polytope_feature/version.py b/polytope_feature/version.py index f871089b..13361d15 100644 --- a/polytope_feature/version.py +++ b/polytope_feature/version.py @@ -1 +1 @@ -__version__ = "1.0.15" +__version__ = "1.0.16" diff --git a/tests/test_engine_slicer.py b/tests/test_engine_slicer.py index 3392ccd3..a56a4693 100644 --- a/tests/test_engine_slicer.py +++ b/tests/test_engine_slicer.py @@ -45,7 +45,7 @@ def test_triangle(self): triangle = Polygon(["x", "y"], [[3, 3], [3, 6], [6, 3]]).polytope() result = self.slicer.extract(datacube, triangle) result.pprint() - assert len(result.leaves) == 10 + assert len(result.leaves) == 4 # assert len(result.leaves) == 4 # total_leaves = 0 # for leaf in result.leaves: @@ -57,8 +57,8 @@ def test_reusable(self): polytopes = Polygon(["x", "y"], [[3, 3], [3, 6], [6, 3]]).polytope() result = self.slicer.extract(datacube, polytopes) result.pprint() - # assert len(result.leaves) == 4 - assert len(result.leaves) == 10 + assert len(result.leaves) == 4 + # assert len(result.leaves) == 10 polytopes = Box(["x", "y"], lower_corner=[3, 3], upper_corner=[6, 6]).polytope() result = self.slicer.extract(datacube, polytopes) result.pprint() diff --git a/tests/test_slicing_xarray_3D.py b/tests/test_slicing_xarray_3D.py index 72f68d0f..bc659103 100644 --- a/tests/test_slicing_xarray_3D.py +++ b/tests/test_slicing_xarray_3D.py @@ -3,6 +3,7 @@ import numpy as np import pandas as pd +import pytest import xarray as xr from polytope_feature.datacube.backends.xarray import XArrayDatacube @@ -103,14 +104,15 @@ def test_disk(self): assert np.size(result.leaves[1].result[1]) == 7 assert np.size(result.leaves[2].result[1]) == 1 + @pytest.mark.xfail def test_concave_polygon(self): # TODO: fix the overlapping branches? points = [[1, 0], [3, 0], [2, 3], [3, 6], [1, 6]] request = Request(Polygon(["level", "step"], points), Select("date", ["2000-01-01"])) result = self.API.retrieve(request) self.xarraydatacube.get(result) - # result.pprint() - assert len(result.leaves) == 8 + result.pprint() + assert len(result.leaves) == 4 def test_polytope(self): points = [[0, 1], [3, 1], [3, 2], [0, 2]] @@ -226,6 +228,7 @@ def test_box(self): result = self.API.retrieve(request) assert len(result.leaves) == 1 + @pytest.mark.xfail def test_swept_concave_polygon(self): # Tests what happens when we slice a concave shape which is swept across a path and see if concavity is lost points = [(1, 0), (3, 0), (3, 6), (2, 6), (2, 3), (1, 3)] diff --git a/tests/test_slicing_xarray_4D.py b/tests/test_slicing_xarray_4D.py index 9345e8db..cd4070e1 100644 --- a/tests/test_slicing_xarray_4D.py +++ b/tests/test_slicing_xarray_4D.py @@ -71,14 +71,16 @@ def test_circles_intersecting_float(self): disk2 = Disk(["step", "lat"], [15, 2.0], [4.99, 0.3]) request = Request(Union(["step", "lat"], disk1, disk2), Select("date", ["2000-01-01"]), Select("level", [10])) result = self.API.retrieve(request) - assert len(result.leaves) == 24 + result.pprint() + assert len(result.leaves) == 8 def test_circles_touching_float(self): disk1 = Disk(["step", "lat"], [6, 4.0], [3, 1.9]) disk2 = Disk(["step", "lat"], [15, 2.0], [3, 2.1]) request = Request(Union(["step", "lat"], disk1, disk2), Select("date", ["2000-01-01"]), Select("level", [10])) result = self.API.retrieve(request) - assert len(result.leaves) == 101 + result.pprint() + assert len(result.leaves) == 6 def test_pathsegment_swept_2D_box(self): box1 = Box(["step", "level"], [3, 0], [6, 1]) diff --git a/tests/test_unions_fdb.py b/tests/test_unions_fdb.py new file mode 100644 index 00000000..d3634dc0 --- /dev/null +++ b/tests/test_unions_fdb.py @@ -0,0 +1,319 @@ +import pandas as pd +import pytest + +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Disk, Select, Span, Union + + +class TestSlicingFDBDatacube: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + self.options = { + "axis_config": [ + {"axis_name": "number", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + { + "axis_name": "date", + "transformations": [{"name": "merge", "other_axis": "time", "linkers": ["T", "00"]}], + }, + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "octahedral", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + {"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}, + {"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}, + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + } + + # Testing different shapes + @pytest.mark.fdb + def test_fdb_datacube_box_union(self): + import pygribjump as gj + + box1 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Box(["latitude", "longitude"], [0.3, 0], [0.4, 0.2]) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + assert len(leaves) == 5 + total_leaves = 0 + for i in range(len(leaves)): + total_leaves += len(leaves[i].values) + assert total_leaves == 15 + + @pytest.mark.fdb + def test_fdb_datacube_union_disk(self): + import pygribjump as gj + + box1 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Disk(["latitude", "longitude"], [0.5, 0], [0.1, 0.1]) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 16 + total_leaves = 0 + for i in range(len(leaves)): + total_leaves += len(leaves[i].values) + assert total_leaves == 16 + + @pytest.mark.fdb + def test_fdb_datacube_disk_union_intersecting(self): + import pygribjump as gj + + box1 = Disk(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Disk(["latitude", "longitude"], [0.1, 0], [0.2, 0.2]) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 31 + total_leaves = 0 + for i in range(len(leaves)): + total_leaves += len(leaves[i].values) + assert total_leaves == 31 + + @pytest.mark.fdb + def test_fdb_datacube_box_union_total_overlap(self): + import pygribjump as gj + + box1 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 3 + total_leaves = 0 + for i in range(len(leaves)): + total_leaves += len(leaves[i].values) + assert total_leaves == 9 + + @pytest.mark.fdb + def test_fdb_datacube_box_union_disk_overlap(self): + import pygribjump as gj + + box1 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Disk(["latitude", "longitude"], [0.3, 0], [0.2, 0.2]) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 29 + total_leaves_1 = 0 + for i in range(len(leaves)): + total_leaves_1 += len(leaves[i].values) + assert total_leaves_1 == 29 + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box2, box1), + ) + + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 29 + total_leaves_2 = 0 + for i in range(len(leaves)): + total_leaves_2 += len(leaves[i].values) + assert total_leaves_2 == total_leaves_1 + + @pytest.mark.fdb + def test_fdb_datacube_box_union_intersecting(self): + import pygribjump as gj + + box1 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Box(["latitude", "longitude"], [0.1, 0], [0.3, 0.2]) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 4 + total_leaves = 0 + for i in range(len(leaves)): + total_leaves += len(leaves[i].values) + assert total_leaves == 12 + + @pytest.mark.fdb + def test_fdb_datacube_union_disk_intersecting(self): + import pygribjump as gj + + box1 = Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]) + box2 = Disk(["latitude", "longitude"], [0.2, 0], [0.1, 0.1]) + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Union(["latitude", "longitude"], box1, box2), + ) + + self.fdbdatacube = gj.GribJump() + self.slicer = HullSlicer() + self.API = Polytope( + datacube=self.fdbdatacube, + engine=self.slicer, + options=self.options, + ) + result = self.API.retrieve(request) + result.pprint() + leaves = result.leaves + # assert len(leaves) == 13 + total_leaves = 0 + for i in range(len(leaves)): + total_leaves += len(leaves[i].values) + assert total_leaves == 13