diff --git a/.gitignore b/.gitignore index 7fcdd7723..0a5171b87 100644 --- a/.gitignore +++ b/.gitignore @@ -13,10 +13,10 @@ site *.grib *.gif *.html +*.req +.mypy_cache example_eo example_mri -.mypy_cache -*.req polytope_python.egg-info new_test_venv_polytope new_venv_polytope @@ -26,3 +26,9 @@ new_updated_numpy_venv newest-polytope-venv serializedTree new_polytope_venv + +newest-polytope-venv +*.pstats +*.profile +new_polytope_venv +*.json \ No newline at end of file diff --git a/examples/3D_shipping_route.py b/examples/3D_shipping_route.py index 50db10aa2..571f3a046 100644 --- a/examples/3D_shipping_route.py +++ b/examples/3D_shipping_route.py @@ -6,9 +6,9 @@ import numpy as np from earthkit import data -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Ellipsoid, Path, Select +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Ellipsoid, Path, Select class Test: diff --git a/examples/3D_shipping_route_wave_model.py b/examples/3D_shipping_route_wave_model.py index e4edabe9f..81463ba51 100644 --- a/examples/3D_shipping_route_wave_model.py +++ b/examples/3D_shipping_route_wave_model.py @@ -4,10 +4,10 @@ import pytest from eccodes import codes_grib_find_nearest, codes_grib_new_from_file -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Ellipsoid, Path, Select +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Ellipsoid, Path, Select from tests.helper_functions import download_test_data diff --git a/examples/4D_flight_path.py b/examples/4D_flight_path.py index d2d1f1eaa..0dc0f8552 100644 --- a/examples/4D_flight_path.py +++ b/examples/4D_flight_path.py @@ -4,10 +4,10 @@ from earthkit import data from PIL import Image -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Path, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Path, Select class Test: diff --git a/examples/country_slicing.py b/examples/country_slicing.py index 4a5ed30db..fa88b3a5f 100644 --- a/examples/country_slicing.py +++ b/examples/country_slicing.py @@ -4,10 +4,10 @@ from earthkit import data from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/examples/cyclic_route_around_earth.py b/examples/cyclic_route_around_earth.py index 4e971d677..26cf5ef81 100644 --- a/examples/cyclic_route_around_earth.py +++ b/examples/cyclic_route_around_earth.py @@ -3,10 +3,10 @@ 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, PathSegment, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, PathSegment, Select class Test: diff --git a/examples/healpix_grid_box_example.py b/examples/healpix_grid_box_example.py index 6cfbb5e78..4e759e28f 100644 --- a/examples/healpix_grid_box_example.py +++ b/examples/healpix_grid_box_example.py @@ -3,10 +3,10 @@ from earthkit import data from eccodes import codes_grib_find_nearest, codes_grib_new_from_file -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 +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select class TestOctahedralGrid: diff --git a/examples/octahedral_grid_box_example.py b/examples/octahedral_grid_box_example.py index a4cddcacd..8293c9eb7 100644 --- a/examples/octahedral_grid_box_example.py +++ b/examples/octahedral_grid_box_example.py @@ -6,10 +6,10 @@ from eccodes import codes_grib_find_nearest, codes_grib_new_from_file from matplotlib import markers -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 +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select def find_nearest_latlon(grib_file, target_lat, target_lon): diff --git a/examples/octahedral_grid_country_example.py b/examples/octahedral_grid_country_example.py index 280701d2b..9d777872b 100644 --- a/examples/octahedral_grid_country_example.py +++ b/examples/octahedral_grid_country_example.py @@ -7,10 +7,10 @@ from matplotlib import markers from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union def find_nearest_latlon(grib_file, target_lat, target_lon): diff --git a/examples/read_me_example.py b/examples/read_me_example.py index 0c13e1274..41e3247e6 100644 --- a/examples/read_me_example.py +++ b/examples/read_me_example.py @@ -1,8 +1,8 @@ import numpy as np from earthkit import data -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select ds = data.from_source("file", "./examples/data/winds.grib") array = ds.to_xarray() diff --git a/examples/slicing_all_ecmwf_countries.py b/examples/slicing_all_ecmwf_countries.py index 7178305a9..dd5e029af 100644 --- a/examples/slicing_all_ecmwf_countries.py +++ b/examples/slicing_all_ecmwf_countries.py @@ -4,10 +4,10 @@ from earthkit import data from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/examples/timeseries_example.py b/examples/timeseries_example.py index bd350cc13..edbed2efc 100644 --- a/examples/timeseries_example.py +++ b/examples/timeseries_example.py @@ -3,10 +3,10 @@ from earthkit import data from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/examples/wind_farms.py b/examples/wind_farms.py index 23699eb36..7c8284c8a 100644 --- a/examples/wind_farms.py +++ b/examples/wind_farms.py @@ -6,10 +6,10 @@ from osgeo import gdal from shapely.geometry import shape -from polytope.datacube.backends.xarray import XArrayDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Polygon, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Polygon, Select, Union class Test: diff --git a/performance/fdb_performance.py b/performance/fdb_performance.py index 78819d462..a8d9e0133 100644 --- a/performance/fdb_performance.py +++ b/performance/fdb_performance.py @@ -2,10 +2,10 @@ import pandas as pd -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select class TestSlicingFDBDatacube: diff --git a/performance/fdb_performance_3D.py b/performance/fdb_performance_3D.py index 88fb85c92..75c87ba62 100644 --- a/performance/fdb_performance_3D.py +++ b/performance/fdb_performance_3D.py @@ -2,10 +2,10 @@ import pandas as pd -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.engine.hullslicer import HullSlicer -from polytope.polytope import Polytope, Request -from polytope.shapes import Box, Select, Span +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select, Span class TestSlicingFDBDatacube: diff --git a/performance/fdb_slice_many_numbers_timeseries.py b/performance/fdb_slice_many_numbers_timeseries.py index 6d62df54c..662b5cc38 100644 --- a/performance/fdb_slice_many_numbers_timeseries.py +++ b/performance/fdb_slice_many_numbers_timeseries.py @@ -3,13 +3,12 @@ import pandas as pd import pygribjump as gj -from polytope.polytope import Polytope, Request -from polytope.shapes import All, Point, Select +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import All, Point, Select time1 = time.time() # Create a dataarray with 3 labelled axes using different index types -# config = {"class": "od", "expver": "0001", "levtype": "sfc", "type": "pf"} options = { "axis_config": [ {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, @@ -62,7 +61,6 @@ Select("class", ["od"]), Select("stream", ["enfo"]), Select("type", ["pf"]), - # Select("latitude", [0.035149384216], method="surrounding"), Point(["latitude", "longitude"], [[0.04, 0]], method="surrounding"), All("number"), ) diff --git a/performance/performance_many_num_steps.py b/performance/performance_many_num_steps.py index eb47b0a06..0cdcefcb2 100644 --- a/performance/performance_many_num_steps.py +++ b/performance/performance_many_num_steps.py @@ -4,9 +4,9 @@ import pandas as pd from earthkit import data -from polytope.datacube.backends.fdb import FDBDatacube -from polytope.polytope import Polytope, Request -from polytope.shapes import All, Point, Select, Span +from polytope_feature.datacube.backends.fdb import FDBDatacube +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import All, Point, Select, Span time1 = time.time() # Create a dataarray with 3 labelled axes using different index types diff --git a/performance/scalability_test.py b/performance/scalability_test.py index ab3c5b299..c0aba173f 100644 --- a/performance/scalability_test.py +++ b/performance/scalability_test.py @@ -3,10 +3,10 @@ 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 Box, Disk, Ellipsoid, Select +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Disk, Ellipsoid, Select class Test: diff --git a/performance/scalability_test_2.py b/performance/scalability_test_2.py index b10e86f4c..7ba0bff4a 100644 --- a/performance/scalability_test_2.py +++ b/performance/scalability_test_2.py @@ -3,10 +3,10 @@ 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 Box, Select, Union +from polytope_feature.datacube.backends.xarray import XArrayDatacube +from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select, Union class Test: diff --git a/polytope_feature/datacube/backends/datacube.py b/polytope_feature/datacube/backends/datacube.py index 97a8c7c77..9f961b174 100644 --- a/polytope_feature/datacube/backends/datacube.py +++ b/polytope_feature/datacube/backends/datacube.py @@ -152,20 +152,40 @@ def remap_path(self, path: DatacubePath): return path @staticmethod - def create(datacube, config={}, axis_options={}, compressed_axes_options=[], alternative_axes=[], context=None): + def create( + datacube, + config={}, + axis_options={}, + compressed_axes_options=[], + point_cloud_options=None, + alternative_axes=[], + context=None, + ): # TODO: get the configs as None for pre-determined value and change them to empty dictionary inside the function if type(datacube).__name__ == "DataArray": from .xarray import XArrayDatacube - xadatacube = XArrayDatacube(datacube, axis_options, compressed_axes_options, context) + xadatacube = XArrayDatacube(datacube, axis_options, compressed_axes_options, point_cloud_options, context) return xadatacube if type(datacube).__name__ == "GribJump": from .fdb import FDBDatacube fdbdatacube = FDBDatacube( - datacube, config, axis_options, compressed_axes_options, alternative_axes, context + datacube, + config, + axis_options, + compressed_axes_options, + point_cloud_options, + alternative_axes, + context, ) return fdbdatacube + if type(datacube).__name__ == "MockDatacube": + return datacube def check_branching_axes(self, request): pass + + @abstractmethod + def find_point_cloud(self): + pass diff --git a/polytope_feature/datacube/backends/fdb.py b/polytope_feature/datacube/backends/fdb.py index bd60665af..858cdafc7 100644 --- a/polytope_feature/datacube/backends/fdb.py +++ b/polytope_feature/datacube/backends/fdb.py @@ -10,7 +10,14 @@ class FDBDatacube(Datacube): def __init__( - self, gj, config=None, axis_options=None, compressed_axes_options=[], alternative_axes=[], context=None + self, + gj, + config=None, + axis_options=None, + compressed_axes_options=[], + point_cloud_options=None, + alternative_axes=[], + context=None, ): if config is None: config = {} @@ -23,6 +30,7 @@ def __init__( self.unwanted_path = {} self.axis_options = axis_options + self.has_point_cloud = point_cloud_options # NOTE: here, will be True/False partial_request = config # Find values in the level 3 FDB datacube @@ -70,7 +78,12 @@ def __init__( val = self._axes[name].type self._check_and_add_axes(options, name, val) - logging.info("Polytope created axes for %s", self._axes.keys()) + logging.info("Polytope created axes for: " + str(self._axes.keys())) + + def find_point_cloud(self): + # TODO: somehow, find the point cloud of irregular grid if it exists + if self.has_point_cloud: + return self.has_point_cloud def check_branching_axes(self, request): polytopes = request.polytopes() @@ -121,6 +134,7 @@ def get(self, requests: TensorIndexTree, context=None): for i, key in enumerate(compressed_request[0].keys()): uncompressed_request[key] = combi[i] complete_uncompressed_request = (uncompressed_request, compressed_request[1], self.grid_md5_hash) + print(self.grid_md5_hash) complete_list_complete_uncompressed_requests.append(complete_uncompressed_request) complete_fdb_decoding_info.append(fdb_requests_decoding_info[j]) @@ -282,6 +296,7 @@ def get_2nd_last_values(self, requests, leaf_path=None): leaf_path_copy = deepcopy(leaf_path) leaf_path_copy.pop("values", None) + leaf_path_copy.pop("index") return (leaf_path_copy, current_start_idxs, fdb_node_ranges, lat_length) def get_last_layer_before_leaf(self, requests, leaf_path, current_idx, fdb_range_n): @@ -290,6 +305,7 @@ def get_last_layer_before_leaf(self, requests, leaf_path, current_idx, fdb_range for i, c in enumerate(requests.children): # now c are the leaves of the initial tree key_value_path = {c.axis.name: c.values} + leaf_path["index"] = c.indexes ax = c.axis (key_value_path, leaf_path, self.unwanted_path) = ax.unmap_path_key( key_value_path, leaf_path, self.unwanted_path diff --git a/polytope_feature/datacube/backends/mock.py b/polytope_feature/datacube/backends/mock.py index c23b80802..f652bfe99 100644 --- a/polytope_feature/datacube/backends/mock.py +++ b/polytope_feature/datacube/backends/mock.py @@ -13,10 +13,10 @@ def __init__(self, dimensions, compressed_axes_options=[]): self.dimensions = dimensions - self.mappers = {} + self._axes = {} for name in self.dimensions: - self.mappers[name] = deepcopy(IntDatacubeAxis()) - self.mappers[name].name = name + self._axes[name] = deepcopy(IntDatacubeAxis()) + self._axes[name].name = name self.stride = {} stride_cumulative = 1 @@ -42,7 +42,7 @@ def get(self, requests: TensorIndexTree, context=None): r.remove_branch() def get_mapper(self, axis): - return self.mappers[axis] + return self._axes[axis] def get_indices(self, path: DatacubePath, axis, lower, upper, method=None): if lower == upper == math.ceil(lower): @@ -59,7 +59,7 @@ def has_index(self, path: DatacubePath, axis, index): @property def axes(self): - return self.mappers + return self._axes def validate(self, axes): return validate_axes(self.axes, axes) @@ -69,3 +69,6 @@ def ax_vals(self, name): def _find_indexes_between(self, axis, indexes, low, up): pass + + def find_point_cloud(self): + pass diff --git a/polytope_feature/datacube/backends/xarray.py b/polytope_feature/datacube/backends/xarray.py index 735c87e25..67f3e1cb5 100644 --- a/polytope_feature/datacube/backends/xarray.py +++ b/polytope_feature/datacube/backends/xarray.py @@ -9,7 +9,14 @@ class XArrayDatacube(Datacube): """Xarray arrays are labelled, axes can be defined as strings or integers (e.g. "time" or 0).""" - def __init__(self, dataarray: xr.DataArray, axis_options=None, compressed_axes_options=[], context=None): + def __init__( + self, + dataarray: xr.DataArray, + axis_options=None, + compressed_axes_options=[], + point_cloud_options=None, + context=None, + ): super().__init__(axis_options, compressed_axes_options) if axis_options is None: axis_options = {} @@ -17,6 +24,7 @@ def __init__(self, dataarray: xr.DataArray, axis_options=None, compressed_axes_o self.axis_counter = 0 self._axes = None self.dataarray = dataarray + self.has_point_cloud = point_cloud_options for name, values in dataarray.coords.variables.items(): options = None @@ -50,9 +58,12 @@ def __init__(self, dataarray: xr.DataArray, axis_options=None, compressed_axes_o val = self._axes[name].type self._check_and_add_axes(options, name, val) + def find_point_cloud(self): + # TODO: somehow, find the point cloud of irregular grid if it exists + if self.has_point_cloud: + return self.has_point_cloud + def get(self, requests, context=None, leaf_path=None, axis_counter=0): - if context is None: - context = {} if leaf_path is None: leaf_path = {} if requests.axis.name == "root": @@ -68,6 +79,8 @@ def get(self, requests, context=None, leaf_path=None, axis_counter=0): if len(requests.children) != 0: # We are not a leaf and we loop over for c in requests.children: + if axis_counter == self.axis_counter - 1: + leaf_path["index"] = c.indexes self.get(c, context, leaf_path, axis_counter + 1) else: if self.axis_counter != axis_counter: @@ -78,7 +91,8 @@ def get(self, requests, context=None, leaf_path=None, axis_counter=0): unmapped_path = {} self.refit_path(leaf_path_copy, unmapped_path, leaf_path) for key in leaf_path_copy: - leaf_path_copy[key] = list(leaf_path_copy[key]) + if isinstance(leaf_path_copy[key], tuple): + leaf_path_copy[key] = list(leaf_path_copy[key]) for key in unmapped_path: if isinstance(unmapped_path[key], tuple): unmapped_path[key] = list(unmapped_path[key]) @@ -108,9 +122,9 @@ def refit_path(self, path_copy, unmapped_path, path): for key in path.keys(): if key not in self.dataarray.dims: path_copy.pop(key) - if key not in self.dataarray.coords.dtypes: + elif key not in self.dataarray.coords.dtypes: unmapped_path.update({key: path[key]}) - path_copy.pop(key) + path_copy.pop(key, None) for key in self.dataarray.coords.dtypes: key_dtype = self.dataarray.coords.dtypes[key] if key_dtype.type is np.str_ and key in path.keys(): diff --git a/polytope_feature/datacube/datacube_axis.py b/polytope_feature/datacube/datacube_axis.py index bf4dd7b8d..d736cdfb6 100644 --- a/polytope_feature/datacube/datacube_axis.py +++ b/polytope_feature/datacube/datacube_axis.py @@ -100,6 +100,11 @@ def _remap_val_to_axis_range(self, value): value = transformation._remap_val_to_axis_range(value, self) return value + def remap_polytopes(self, polytopes): + for transformation in self.transformations[::-1]: + polytopes = transformation.remap_polytopes(self, polytopes) + return polytopes + def find_standard_indices_between(self, indexes, low, up, datacube, method=None): indexes_between_ranges = [] diff --git a/polytope_feature/datacube/quad_tree.py b/polytope_feature/datacube/quad_tree.py new file mode 100644 index 000000000..4dd5d3020 --- /dev/null +++ b/polytope_feature/datacube/quad_tree.py @@ -0,0 +1,175 @@ +from ..utility.slicing_tools import slice, slice_in_two + +""" + + QuadTree data structure + + comes from: https://github.com/mdrasmus/compbio/blob/master/rasmus/quadtree.py + + is specific to 2D + +""" + + +class QuadNode: + def __init__(self, item, index): + self.item = item + self.index = index + + def is_contained_in(self, polygon): + # implement method to check if the node point is inside the polygon + node_x, node_y = self.item + + sliced_vertical_polygon = slice(polygon, polygon._axes[0], node_x, 0) + if sliced_vertical_polygon: + lower_y, upper_y, idx = sliced_vertical_polygon.extents(polygon._axes[1]) + if lower_y <= node_y <= upper_y: + return True + return False + + +class QuadTree: + # TODO: do we need the max_depth? + MAX = 3 + MAX_DEPTH = 20 + + def __init__(self, x=0, y=0, size=[180, 90], depth=0): + self.nodes = [] + self.children = [] + self.center = (x, y) + self.size = tuple(size) + self.depth = depth + self.node_items = set() + + def quadrant_rectangle_points(self): + return set( + [ + (self.center[0] + (self.size[0]), self.center[1] + (self.size[1])), + (self.center[0] + (self.size[0]), self.center[1] - (self.size[1])), + (self.center[0] - (self.size[0]), self.center[1] + (self.size[1])), + (self.center[0] - (self.size[0]), self.center[1] - (self.size[1])), + ] + ) + + def build_point_tree(self, points): + # TODO: SLOW, scales linearly with number of points + for index, p in enumerate(points): + self.insert(tuple(p), index) + + def pprint(self): + if self.depth == 0: + print("\n") + if len(self.children) == 0: + for n in self.nodes: + print("\t" * (self.depth + 1) + "\u21b3" + str(n.item[0]) + " , " + str(n.item[1])) + for child in self.children: + print( + "\t" * (self.depth + 1) + + "\u21b3" + + str(child.center[0] - child.size[0]) + + " , " + + str(child.center[1] - child.size[1]) + + " , " + + str(child.center[0] + child.size[0]) + + " , " + + str(child.center[1] + child.size[1]) + ) + child.pprint() + + def insert(self, item, index): + if not self.children: + node = QuadNode(item, index) + if item not in self.node_items: + self.nodes.append(node) + self.node_items.add(node.item) + + if len(self.nodes) > self.MAX and self.depth < self.MAX_DEPTH: + self.split() + return node + else: + return self.insert_into_children(item, index) + + def insert_into_children(self, item, index): + x, y = item + cx, cy = self.center + # try to insert into children + if x <= cx: + if y <= cy: + self.children[0].insert(item, index) + if y >= cy: + self.children[1].insert(item, index) + if x >= cx: + if y <= cy: + self.children[2].insert(item, index) + if y >= cy: + self.children[3].insert(item, index) + + def split(self): + half_size = [s / 2 for s in self.size] + x_center, y_center = self.center[0], self.center[1] + hx, hy = half_size + + new_centers = [ + (x_center - hx, y_center - hy), + (x_center - hx, y_center + hy), + (x_center + hx, y_center - hy), + (x_center + hx, y_center + hy), + ] + + self.children = [ + QuadTree(new_center[0], new_center[1], half_size, self.depth + 1) for new_center in new_centers + ] + + nodes = self.nodes + self.nodes = [] + for node in nodes: + self.insert_into_children(node.item, node.index) + + def query_polygon(self, polygon, results=None): + # intersect quad tree with a 2D polygon + if results is None: + results = set() + + # intersect the children with the polygon + # TODO: here, we create None polygons... think about how to handle them + if polygon is None: + pass + else: + polygon_points = set([tuple(point) for point in polygon.points]) + # TODO: are these the right points which we are comparing, ie the points on the polygon + # and the points on the rectangle quadrant? + if polygon_points == self.quadrant_rectangle_points(): + for node in self.find_nodes_in(): + results.add(node) + else: + if len(self.children) > 0: + # first slice vertically + left_polygon, right_polygon = slice_in_two(polygon, self.center[0], 0) + + # then slice horizontally + # ie need to slice the left and right polygons each in two to have the 4 quadrant polygons + + q1_polygon, q2_polygon = slice_in_two(left_polygon, self.center[1], 1) + q3_polygon, q4_polygon = slice_in_two(right_polygon, self.center[1], 1) + + # now query these 4 polygons further down the quadtree + self.children[0].query_polygon(q1_polygon, results) + self.children[1].query_polygon(q2_polygon, results) + self.children[2].query_polygon(q3_polygon, results) + self.children[3].query_polygon(q4_polygon, results) + + results.update(node for node in self.nodes if node.is_contained_in(polygon)) + + return results + + def find_nodes_in(self, results=None): + # TODO: find the nodes that are in this subtree + if results is None: + results = set() + if len(self.children) > 0: + # there are children which we need to iterate through + for child in self.children: + child.find_nodes_in(results) + for node in self.nodes: + results.add(node) + return results diff --git a/polytope_feature/datacube/tensor_index_tree.py b/polytope_feature/datacube/tensor_index_tree.py index 604aa7441..b22276af2 100644 --- a/polytope_feature/datacube/tensor_index_tree.py +++ b/polytope_feature/datacube/tensor_index_tree.py @@ -109,6 +109,7 @@ def add_value(self, value): self.values = tuple(new_values) def create_child(self, axis, value, next_nodes): + # TODO: what if we remove the next nodes here? node = TensorIndexTree(axis, (value,)) existing_child = self.find_child(node) if not existing_child: @@ -138,12 +139,9 @@ def is_root(self): def find_child(self, node): index = self.children.bisect_left(node) - if index >= len(self.children): - return None - child = self.children[index] - if not child == node: - return None - return child + if index < len(self.children) and self.children[index] == node: + return self.children[index] + return None def add_node_layer_after(self, ax_name, vals): ax = IntDatacubeAxis() diff --git a/polytope_feature/datacube/transformations/datacube_cyclic/datacube_cyclic.py b/polytope_feature/datacube/transformations/datacube_cyclic/datacube_cyclic.py index 3373dd082..600e3718d 100644 --- a/polytope_feature/datacube/transformations/datacube_cyclic/datacube_cyclic.py +++ b/polytope_feature/datacube/transformations/datacube_cyclic/datacube_cyclic.py @@ -2,6 +2,7 @@ from copy import deepcopy from ....utility.list_tools import unique +from ....utility.slicing_tools import slice_in_two from ..datacube_transformations import DatacubeAxisTransformation @@ -28,6 +29,40 @@ def blocked_axes(self): def unwanted_axes(self): return [] + def remap_polytopes(self, axis, polytopes): + # Find extents of the polytope on that axis + all_sliced_polys = [] + for polytope in polytopes: + lower, upper, slice_axis_idx = polytope.extents(axis.name) + polytope_range = [lower, upper] + + # Find intervals within this range + intervals = self.to_intervals(polytope_range, [], axis) + + slice_vals = [] + + for interval in intervals[:-1]: + slice_vals.append(interval[1]) + + # Successively slice the polytope on each of the slice vals + sliced_polys = self.slice_polytope([polytope], slice_vals, slice_axis_idx) + all_sliced_polys.extend(sliced_polys) + all_sliced_polys = [poly for poly in all_sliced_polys if poly is not None] + + for poly in all_sliced_polys: + for point in poly.points: + point[slice_axis_idx] = self._remap_val_to_axis_range(point[slice_axis_idx], axis) + return all_sliced_polys + + def slice_polytope(self, polytope_list, slice_vals, slice_axis_idx): + for slice_val in slice_vals: + polytope_to_slice = polytope_list[-1] + left_sliced_poly, right_sliced_poly = slice_in_two(polytope_to_slice, slice_val, slice_axis_idx) + polytope_list = polytope_list[:-1] + polytope_list.append(left_sliced_poly) + polytope_list.append(right_sliced_poly) + return polytope_list + def update_range(self, axis): axis.range = self.range diff --git a/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py b/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py index e28566e69..23e891398 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/datacube_mappers.py @@ -86,8 +86,8 @@ def find_second_idx(self, first_val, second_val): def unmap_first_val_to_start_line_idx(self, first_val): return self._final_transformation.unmap_first_val_to_start_line_idx(first_val) - def unmap(self, first_val, second_val): - return self._final_transformation.unmap(first_val, second_val) + def unmap(self, first_val, second_val, unmapped_idx=None): + return self._final_transformation.unmap(first_val, second_val, unmapped_idx) def find_modified_indexes(self, indexes, path, datacube, axis): if axis.name == self._mapped_axes()[0]: @@ -105,10 +105,16 @@ def unmap_path_key(self, key_value_path, leaf_path, unwanted_path, axis): unwanted_path[axis.name] = unwanted_val if axis.name == self._mapped_axes()[1]: first_val = unwanted_path[self._mapped_axes()[0]] - unmapped_idx = [] - for val in value: - unmapped_idx.append(self.unmap(first_val, (val,))) - # unmapped_idx = self.unmap(first_val, value) + unmapped_idx = leaf_path.get("index", None) + if unmapped_idx is not None and len(unmapped_idx) > 0: + # for val in value: + # unmapped_idx.append(self.unmap(first_val, (val,))) + unmapped_idx = list(unmapped_idx) + else: + unmapped_idx = [] + for val in value: + unmapped_idx.append(self.unmap(first_val, (val,))) + # unmapped_idx = self.unmap(first_val, value, unmapped_idx) leaf_path.pop(self._mapped_axes()[0], None) key_value_path.pop(axis.name) key_value_path[self.old_axis] = unmapped_idx @@ -136,5 +142,6 @@ def unmap_tree_node(self, node, unwanted_path): "regular": "RegularGridMapper", "reduced_ll": "ReducedLatLonMapper", "local_regular": "LocalRegularGridMapper", + "irregular": "IrregularGridMapper", "healpix_nested": "NestedHealpixGridMapper", } diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py index 6fdf9d54a..556701312 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix.py @@ -133,7 +133,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): else: return idx - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_value = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_value) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py index c71f59f5b..36a69492a 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/healpix_nested.py @@ -137,7 +137,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): else: return idx - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_value = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_value) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/irregular.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/irregular.py new file mode 100644 index 000000000..0dad07552 --- /dev/null +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/irregular.py @@ -0,0 +1,23 @@ +from ..datacube_mappers import DatacubeMapper + + +class IrregularGridMapper(DatacubeMapper): + # def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + def __init__(self, base_axis, mapped_axes, resolution, md5_hash=None, local_area=[], axis_reversed=None): + self._mapped_axes = mapped_axes + self._base_axis = base_axis + self._resolution = resolution + self._axis_reversed = False + self.compressed_grid_axes = [self._mapped_axes[1]] + if md5_hash is not None: + self.md5_hash = md5_hash + else: + self.md5_hash = _md5_hash.get(resolution, None) + + def unmap(self, first_val, second_val, unmapped_idx=None): + # TODO: But to unmap for the irregular grid, need the request tree + # Suppose we get the idx value somehow from the tree, as an idx input + return unmapped_idx[0] + + +_md5_hash = {} diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py index df611250e..fef0cb1f2 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/local_regular.py @@ -81,7 +81,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): first_idx = self._first_axis_vals.index(first_val) return first_idx * self.second_resolution - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_val = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_val) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py index e03a45b67..415a5ba70 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/octahedral.py @@ -7883,7 +7883,7 @@ def find_second_axis_idx(self, first_val, second_val): second_idx = int(second_val[0] / second_axis_spacing) return (first_idx, second_idx) - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): (first_idx, second_idx) = self.find_second_axis_idx(first_val, second_val) octahedral_index = self.axes_idx_to_octahedral_idx(first_idx, second_idx) return octahedral_index diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py index 391e1bbe0..036060c09 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py @@ -1504,7 +1504,7 @@ def find_second_idx(self, first_val, second_val): second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) return second_idx - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_value = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_value) diff --git a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py index f4dfb6486..7bac1ea92 100644 --- a/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py +++ b/polytope_feature/datacube/transformations/datacube_mappers/mapper_types/regular.py @@ -61,7 +61,7 @@ def unmap_first_val_to_start_line_idx(self, first_val): first_idx = self._first_axis_vals.index(first_val) return first_idx * 4 * self._resolution - def unmap(self, first_val, second_val): + def unmap(self, first_val, second_val, unmapped_idx=None): tol = 1e-8 first_val = [i for i in self._first_axis_vals if first_val[0] - tol <= i <= first_val[0] + tol][0] first_idx = self._first_axis_vals.index(first_val) diff --git a/polytope_feature/datacube/transformations/datacube_transformations.py b/polytope_feature/datacube/transformations/datacube_transformations.py index f174bf1aa..4bb9ca215 100644 --- a/polytope_feature/datacube/transformations/datacube_transformations.py +++ b/polytope_feature/datacube/transformations/datacube_transformations.py @@ -70,6 +70,9 @@ def remap(self, range, ranges, axis): def to_intervals(self, range, intervals, axis): return intervals + def remap_polytopes(self, axis, polytopes): + return polytopes + _type_to_datacube_transformation_lookup = { "mapper": "DatacubeMapper", diff --git a/polytope_feature/engine/engine.py b/polytope_feature/engine/engine.py index c714db0b2..23d7e6238 100644 --- a/polytope_feature/engine/engine.py +++ b/polytope_feature/engine/engine.py @@ -1,19 +1,56 @@ +import math +from abc import abstractmethod from typing import List from ..datacube.backends.datacube import Datacube +from ..datacube.datacube_axis import UnsliceableDatacubeAxis from ..datacube.tensor_index_tree import TensorIndexTree from ..shapes import ConvexPolytope class Engine: - def __init__(self): - pass + def __init__(self, engine_options=None): + if engine_options is None: + engine_options = {} + self.engine_options = engine_options + self.ax_is_unsliceable = {} + + self.axis_values_between = {} + self.sliced_polytopes = {} + self.remapped_vals = {} + self.compressed_axes = [] def extract(self, datacube: Datacube, polytopes: List[ConvexPolytope]) -> TensorIndexTree: + # Delegate to the right slicer that the axes within the polytopes need to use pass + def check_slicer(self, ax): + # Return the slicer instance if ax is sliceable. + # If the ax is unsliceable, return None. + if isinstance(ax, UnsliceableDatacubeAxis): + return None + slicer_type = self.engine_options[ax.name] + slicer = self.generate_slicer(slicer_type) + return slicer + @staticmethod def default(): from .hullslicer import HullSlicer return HullSlicer() + + def remap_values(self, ax, value): + remapped_val = self.remapped_vals.get((value, ax.name), None) + if remapped_val is None: + remapped_val = value + if ax.is_cyclic: + remapped_val_interm = ax.remap([value, value])[0] + remapped_val = (remapped_val_interm[0] + remapped_val_interm[1]) / 2 + if ax.can_round: + remapped_val = round(remapped_val, int(-math.log10(ax.tol))) + self.remapped_vals[(value, ax.name)] = remapped_val + return remapped_val + + @abstractmethod + def _build_branch(self, ax, node, datacube, next_nodes, api): + pass diff --git a/polytope_feature/engine/hullslicer.py b/polytope_feature/engine/hullslicer.py index 9a99682a6..a505a888f 100644 --- a/polytope_feature/engine/hullslicer.py +++ b/polytope_feature/engine/hullslicer.py @@ -1,41 +1,13 @@ -import math from copy import copy -from itertools import chain -from typing import List -import scipy.spatial - -from ..datacube.backends.datacube import Datacube -from ..datacube.datacube_axis import UnsliceableDatacubeAxis -from ..datacube.tensor_index_tree import TensorIndexTree -from ..shapes import ConvexPolytope -from ..utility.combinatorics import group, tensor_product from ..utility.exceptions import UnsliceableShapeError -from ..utility.geometry import lerp -from ..utility.list_tools import argmax, argmin, unique +from ..utility.slicing_tools import slice from .engine import Engine class HullSlicer(Engine): def __init__(self): - self.ax_is_unsliceable = {} - self.axis_values_between = {} - self.has_value = {} - self.sliced_polytopes = {} - self.remapped_vals = {} - self.compressed_axes = [] - - def _unique_continuous_points(self, p: ConvexPolytope, datacube: Datacube): - for i, ax in enumerate(p._axes): - mapper = datacube.get_mapper(ax) - if self.ax_is_unsliceable.get(ax, None) is None: - self.ax_is_unsliceable[ax] = isinstance(mapper, UnsliceableDatacubeAxis) - if self.ax_is_unsliceable[ax]: - break - for j, val in enumerate(p.points): - p.points[j][i] = mapper.to_float(mapper.parse(p.points[j][i])) - # Remove duplicate points - unique(p.points) + super().__init__() def _build_unsliceable_child(self, polytope, ax, node, datacube, lowers, next_nodes, slice_axis_idx): if not polytope.is_flat: @@ -97,21 +69,9 @@ def find_values_between(self, polytope, ax, node, datacube, lower, upper): self.axis_values_between[(flattened_tuple, ax.name, lower, upper, method)] = values return values - def remap_values(self, ax, value): - remapped_val = self.remapped_vals.get((value, ax.name), None) - if remapped_val is None: - remapped_val = value - if ax.is_cyclic: - remapped_val_interm = ax.remap([value, value])[0] - remapped_val = (remapped_val_interm[0] + remapped_val_interm[1]) / 2 - if ax.can_round: - remapped_val = round(remapped_val, int(-math.log10(ax.tol))) - self.remapped_vals[(value, ax.name)] = remapped_val - return remapped_val - - def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_nodes, slice_axis_idx): + def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_nodes, slice_axis_idx, api): for i, value in enumerate(values): - if i == 0 or ax.name not in self.compressed_axes: + if i == 0 or ax.name not in api.compressed_axes: fvalue = ax.to_float(value) new_polytope = slice(polytope, ax.name, fvalue, slice_axis_idx) remapped_val = self.remap_values(ax, value) @@ -125,8 +85,8 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, values, next_node remapped_val = self.remap_values(ax, value) child.add_value(remapped_val) - def _build_branch(self, ax, node, datacube, next_nodes): - if ax.name not in self.compressed_axes: + def _build_branch(self, ax, node, datacube, next_nodes, api): + if ax.name not in api.compressed_axes: parent_node = node.parent right_unsliced_polytopes = [] for polytope in node["unsliced_polytopes"]: @@ -137,7 +97,7 @@ def _build_branch(self, ax, node, datacube, next_nodes): lower, upper, slice_axis_idx = polytope.extents(ax.name) # here, first check if the axis is an unsliceable axis and directly build node if it is # NOTE: we should have already created the ax_is_unsliceable cache before - if self.ax_is_unsliceable[ax.name]: + if api.ax_is_unsliceable[ax.name]: self._build_unsliceable_child(polytope, ax, node, datacube, [lower], next_nodes, slice_axis_idx) else: values = self.find_values_between(polytope, ax, node, datacube, lower, upper) @@ -148,7 +108,7 @@ def _build_branch(self, ax, node, datacube, next_nodes): # we have iterated all polytopes and we can now remove the node if we need to if len(values) == 0 and len(node.children) == 0: node.remove_branch() - self._build_sliceable_child(polytope, ax, node, datacube, values, next_nodes, slice_axis_idx) + self._build_sliceable_child(polytope, ax, node, datacube, values, next_nodes, slice_axis_idx, api) else: all_values = [] all_lowers = [] @@ -164,12 +124,12 @@ def _build_branch(self, ax, node, datacube, next_nodes): lower, upper, slice_axis_idx = polytope.extents(ax.name) if not first_slice_axis_idx: first_slice_axis_idx = slice_axis_idx - if self.ax_is_unsliceable[ax.name]: + if api.ax_is_unsliceable[ax.name]: all_lowers.append(lower) else: values = self.find_values_between(polytope, ax, node, datacube, lower, upper) all_values.extend(values) - if self.ax_is_unsliceable[ax.name]: + if api.ax_is_unsliceable[ax.name]: self._build_unsliceable_child( first_polytope, ax, node, datacube, all_lowers, next_nodes, first_slice_axis_idx ) @@ -177,140 +137,7 @@ def _build_branch(self, ax, node, datacube, next_nodes): if len(all_values) == 0: node.remove_branch() self._build_sliceable_child( - first_polytope, ax, node, datacube, all_values, next_nodes, first_slice_axis_idx + first_polytope, ax, node, datacube, all_values, next_nodes, first_slice_axis_idx, api ) del node["unsliced_polytopes"] - - def find_compressed_axes(self, datacube, polytopes): - # First determine compressable axes from input polytopes - compressable_axes = [] - for polytope in polytopes: - if polytope.is_orthogonal: - for ax in polytope.axes(): - compressable_axes.append(ax) - # Cross check this list with list of compressable axis from datacube - # (should not include any merged or coupled axes) - for compressed_axis in compressable_axes: - if compressed_axis in datacube.compressed_axes: - self.compressed_axes.append(compressed_axis) - # add the last axis of the grid always (longitude) as a compressed axis - k, last_value = _, datacube.axes[k] = datacube.axes.popitem() - self.compressed_axes.append(k) - - def remove_compressed_axis_in_union(self, polytopes): - for p in polytopes: - if p.is_in_union: - for axis in p.axes(): - if axis == self.compressed_axes[-1]: - self.compressed_axes.remove(axis) - - 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) - - # Convert the polytope points to float type to support triangulation and interpolation - for p in polytopes: - self._unique_continuous_points(p, datacube) - - groups, input_axes = group(polytopes) - datacube.validate(input_axes) - request = TensorIndexTree() - combinations = tensor_product(groups) - - # 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 = [] - for combi in c: - if isinstance(combi, list): - new_c.extend(combi) - else: - new_c.append(combi) - r["unsliced_polytopes"] = set(new_c) - current_nodes = [r] - for ax in datacube.axes.values(): - next_nodes = [] - interm_next_nodes = [] - for node in current_nodes: - self._build_branch(ax, node, datacube, interm_next_nodes) - next_nodes.extend(interm_next_nodes) - interm_next_nodes = [] - current_nodes = next_nodes - - request.merge(r) - return request - - -def _find_intersects(polytope, slice_axis_idx, value): - intersects = [] - # Find all points above and below slice axis - above_slice = [p for p in polytope.points if p[slice_axis_idx] >= value] - below_slice = [p for p in polytope.points if p[slice_axis_idx] <= value] - - # Get the intersection of every pair above and below, this will create excess interior points - for a in above_slice: - for b in below_slice: - # edge is incident with slice plane, don't need these points - if a[slice_axis_idx] == b[slice_axis_idx]: - intersects.append(b) - continue - - # Linearly interpolate all coordinates of two points (a,b) of the polytope - interp_coeff = (value - b[slice_axis_idx]) / (a[slice_axis_idx] - b[slice_axis_idx]) - intersect = lerp(a, b, interp_coeff) - intersects.append(intersect) - return intersects - - -def _reduce_dimension(intersects, slice_axis_idx): - temp_intersects = [] - for point in intersects: - point = [p for i, p in enumerate(point) if i != slice_axis_idx] - temp_intersects.append(point) - return temp_intersects - - -def slice(polytope: ConvexPolytope, axis, value, slice_axis_idx): - if polytope.is_flat: - if value in chain(*polytope.points): - intersects = [[value]] - else: - return None - else: - intersects = _find_intersects(polytope, slice_axis_idx, value) - - if len(intersects) == 0: - return None - - # Reduce dimension of intersection points, removing slice axis - intersects = _reduce_dimension(intersects, slice_axis_idx) - - axes = copy(polytope._axes) - axes.remove(axis) - - if len(intersects) < len(intersects[0]) + 1: - return ConvexPolytope(axes, intersects) - # Compute convex hull (removing interior points) - if len(intersects[0]) == 0: - return None - elif len(intersects[0]) == 1: # qhull doesn't like 1D, do it ourselves - amin = argmin(intersects) - amax = argmax(intersects) - vertices = [amin, amax] - else: - try: - hull = scipy.spatial.ConvexHull(intersects) - vertices = hull.vertices - - except scipy.spatial.qhull.QhullError as e: - if "less than" or "flat" in str(e): - return ConvexPolytope(axes, intersects) - # Sliced result is simply the convex hull - return ConvexPolytope(axes, [intersects[i] for i in vertices]) diff --git a/polytope_feature/engine/quadtree_slicer.py b/polytope_feature/engine/quadtree_slicer.py new file mode 100644 index 000000000..6fd0236b2 --- /dev/null +++ b/polytope_feature/engine/quadtree_slicer.py @@ -0,0 +1,74 @@ +from copy import copy + +from ..datacube.quad_tree import QuadTree +from ..datacube.transformations.datacube_mappers import DatacubeMapper +from .engine import Engine + + +class QuadTreeSlicer(Engine): + def __init__(self, points): + # here need to construct quadtree, which is specific to datacube + # NOTE: should this be inside of the datacube instead that we create the quadtree? + super().__init__() + quad_tree = QuadTree() + quad_tree.build_point_tree(points) + self.quad_tree = quad_tree + self.second_axis = None + + def extract_single(self, datacube, polytope): + # extract a single polygon + + # need to find points of the datacube contained within the polytope + # We do this by intersecting the datacube point cloud quad tree with the polytope here + polygon_points = self.quad_tree.query_polygon(polytope) + return polygon_points + + def _build_branch(self, ax, node, datacube, next_nodes, api): + for polytope in node["unsliced_polytopes"]: + if ax.name in polytope._axes: + # here, first check if the axis is an unsliceable axis and directly build node if it is + + # NOTE: here, we only have sliceable children, since the unsliceable children are handled by the + # hullslicer engine? + self._build_sliceable_child(polytope, ax, node, datacube, next_nodes, api) + del node["unsliced_polytopes"] + + def break_up_polytope(self, datacube, polytope): + new_polytopes = [polytope] + for ax in polytope.axes(): + axis = datacube._axes[ax] + if axis.is_cyclic: + new_polytopes = axis.remap_polytopes(new_polytopes) + return new_polytopes + + def _build_sliceable_child(self, polytope, ax, node, datacube, next_nodes, api): + pre_sliced_polytopes = self.break_up_polytope(datacube, polytope) + + if not self.second_axis: + for transform in ax.transformations: + if isinstance(transform, DatacubeMapper): + self.second_axis = datacube._axes[transform.grid_axes[1]] + + node["unsliced_polytopes"].remove(polytope) + + for poly in pre_sliced_polytopes: + node["unsliced_polytopes"].add(poly) + + extracted_points = self.extract_single(datacube, poly) + # add the sliced points as node to the tree and update the next_nodes + if len(extracted_points) == 0: + node.remove_branch() + + for point in extracted_points: + # convert to float for slicing + value = point.index + lat_val = point.item[0] + lon_val = point.item[1] + + # store the native type + (child, _) = node.create_child(ax, lat_val, []) + (grand_child, _) = child.create_child(self.second_axis, lon_val, []) + # NOTE: the index of the point is stashed in the branches' result + grand_child.indexes = [value] + grand_child["unsliced_polytopes"] = copy(node["unsliced_polytopes"]) + grand_child["unsliced_polytopes"].remove(poly) diff --git a/polytope_feature/options.py b/polytope_feature/options.py index 2e88c9679..d309c0bf3 100644 --- a/polytope_feature/options.py +++ b/polytope_feature/options.py @@ -18,7 +18,7 @@ class CyclicConfig(TransformationConfig): class MapperConfig(TransformationConfig): name: Literal["mapper"] type: str = "" - resolution: Union[int, List[int]] = 0 + resolution: Optional[Union[int, List[int]]] = 0 axes: List[str] = [""] md5_hash: Optional[str] = None local: Optional[List[float]] = None @@ -61,7 +61,7 @@ class Config(ConfigModel): axis_config: List[AxisConfig] = [] compressed_axes_config: List[str] = [""] pre_path: Optional[Dict[str, path_subclasses_union]] = {} - alternative_axes: List[GribJumpAxesConfig] = [] + alternative_axes: Optional[List[GribJumpAxesConfig]] = [] class PolytopeOptions(ABC): diff --git a/polytope_feature/polytope.py b/polytope_feature/polytope.py index f18f10c4d..c787ceb00 100644 --- a/polytope_feature/polytope.py +++ b/polytope_feature/polytope.py @@ -1,9 +1,15 @@ -import logging from typing import List +from .datacube.backends.datacube import Datacube +from .datacube.datacube_axis import UnsliceableDatacubeAxis +from .datacube.tensor_index_tree import TensorIndexTree +from .engine.hullslicer import HullSlicer +from .engine.quadtree_slicer import QuadTreeSlicer from .options import PolytopeOptions from .shapes import ConvexPolytope +from .utility.combinatorics import group, tensor_product from .utility.exceptions import AxisOverdefinedError +from .utility.list_tools import unique class Request: @@ -39,33 +45,143 @@ def __repr__(self): class Polytope: - def __init__(self, datacube, engine=None, options=None, context=None): + def __init__( + self, + datacube, + options=None, + engine_options=None, + point_cloud_options=None, + context=None, + ): from .datacube import Datacube - from .engine import Engine if options is None: options = {} + if engine_options is None: + engine_options = {} - axis_options, compressed_axes_options, config, alternative_axes = PolytopeOptions.get_polytope_options(options) - + self.compressed_axes = [] self.context = context + axis_options, compressed_axes_options, config, alternative_axes = PolytopeOptions.get_polytope_options(options) self.datacube = Datacube.create( - datacube, config, axis_options, compressed_axes_options, alternative_axes, self.context + datacube, + config, + axis_options, + compressed_axes_options, + point_cloud_options, + alternative_axes, + self.context, ) - self.engine = engine if engine is not None else Engine.default() - self.time = 0 - - def slice(self, polytopes: List[ConvexPolytope]): + if engine_options == {}: + for ax_name in self.datacube._axes.keys(): + engine_options[ax_name] = "hullslicer" + self.engine_options = engine_options + self.engines = self.create_engines() + self.ax_is_unsliceable = {} + + def create_engines(self): + engines = {} + engine_types = set(self.engine_options.values()) + if "quadtree" in engine_types: + # TODO: need to get the corresponding point cloud from the datacube + quadtree_points = self.datacube.find_point_cloud() + engines["quadtree"] = QuadTreeSlicer(quadtree_points) + if "hullslicer" in engine_types: + engines["hullslicer"] = HullSlicer() + return engines + + def _unique_continuous_points(self, p: ConvexPolytope, datacube: Datacube): + for i, ax in enumerate(p._axes): + mapper = datacube.get_mapper(ax) + if self.ax_is_unsliceable.get(ax, None) is None: + self.ax_is_unsliceable[ax] = isinstance(mapper, UnsliceableDatacubeAxis) + if self.ax_is_unsliceable[ax]: + break + for j, val in enumerate(p.points): + p.points[j][i] = mapper.to_float(mapper.parse(p.points[j][i])) + # Remove duplicate points + unique(p.points) + + def slice(self, datacube, 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) + self.find_compressed_axes(datacube, polytopes) - def retrieve(self, request: Request, method="standard"): + self.remove_compressed_axis_in_union(polytopes) + + # Convert the polytope points to float type to support triangulation and interpolation + for p in polytopes: + self._unique_continuous_points(p, datacube) + + groups, input_axes = group(polytopes) + datacube.validate(input_axes) + request = TensorIndexTree() + combinations = tensor_product(groups) + + # 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 = [] + for combi in c: + if isinstance(combi, list): + new_c.extend(combi) + else: + new_c.append(combi) + r["unsliced_polytopes"] = set(new_c) + current_nodes = [r] + for ax in datacube.axes.values(): + engine = self.find_engine(ax) + next_nodes = [] + interm_next_nodes = [] + for node in current_nodes: + engine._build_branch(ax, node, datacube, interm_next_nodes, self) + next_nodes.extend(interm_next_nodes) + interm_next_nodes = [] + current_nodes = next_nodes + + request.merge(r) + return request + + def find_engine(self, ax): + slicer_type = self.engine_options[ax.name] + return self.engines[slicer_type] + + def old_retrieve(self, request: Request, method="standard"): """Higher-level API which takes a request and uses it to slice the datacube""" - logging.info("Starting request for %s ", self.context) - self.datacube.check_branching_axes(request) + # self.datacube.check_branching_axes(request) request_tree = self.engine.extract(self.datacube, request.polytopes()) - logging.info("Created request tree for %s ", self.context) - self.datacube.get(request_tree, self.context) - logging.info("Retrieved data for %s ", self.context) + self.datacube.get(request_tree) return request_tree + + def retrieve(self, request: Request, method="standard"): + """Higher-level API which takes a request and uses it to slice the datacube""" + request_tree = self.slice(self.datacube, request.polytopes()) + self.datacube.get(request_tree) + return request_tree + + def find_compressed_axes(self, datacube, polytopes): + # First determine compressable axes from input polytopes + compressable_axes = [] + for polytope in polytopes: + if polytope.is_orthogonal: + for ax in polytope.axes(): + compressable_axes.append(ax) + # Cross check this list with list of compressable axis from datacube + # (should not include any merged or coupled axes) + for compressed_axis in compressable_axes: + if compressed_axis in datacube.compressed_axes: + self.compressed_axes.append(compressed_axis) + + k, last_value = _, datacube.axes[k] = datacube.axes.popitem() + self.compressed_axes.append(k) + + def remove_compressed_axis_in_union(self, polytopes): + for p in polytopes: + if p.is_in_union: + for axis in p.axes(): + if axis in self.compressed_axes: + if axis == self.compressed_axes[-1]: + self.compressed_axes.remove(axis) diff --git a/polytope_feature/utility/engine_tools.py b/polytope_feature/utility/engine_tools.py new file mode 100644 index 000000000..08569ff4e --- /dev/null +++ b/polytope_feature/utility/engine_tools.py @@ -0,0 +1,31 @@ +from typing import List + +from ..datacube.backends.datacube import Datacube +from ..datacube.datacube_axis import UnsliceableDatacubeAxis +from ..datacube.tensor_index_tree import TensorIndexTree +from ..shapes import ConvexPolytope +from .combinatorics import group, tensor_product, unique + + +def unique_continuous_points_in_polytope(p: ConvexPolytope, datacube: Datacube): + # TODO: should this be in utility folder since it doesn't really depend on engine...? + for i, ax in enumerate(p._axes): + mapper = datacube.get_mapper(ax) + if isinstance(mapper, UnsliceableDatacubeAxis): + break + for j, val in enumerate(p.points): + p.points[j][i] = mapper.to_float(mapper.parse(p.points[j][i])) + # Remove duplicate points + unique(p.points) + + +def find_polytope_combinations(datacube: Datacube, polytopes: List[ConvexPolytope]) -> TensorIndexTree: + # here, we find the different possible polytope combinations that cover all of the datacube axes + + for p in polytopes: + unique_continuous_points_in_polytope(p, datacube) + + groups, input_axes = group(polytopes) + datacube.validate(input_axes) + combinations = tensor_product(groups) + return combinations diff --git a/polytope_feature/utility/slicing_tools.py b/polytope_feature/utility/slicing_tools.py new file mode 100644 index 000000000..1c873d392 --- /dev/null +++ b/polytope_feature/utility/slicing_tools.py @@ -0,0 +1,136 @@ +from copy import copy +from itertools import chain + +import scipy +import scipy.spatial + +from ..shapes import ConvexPolytope +from .geometry import lerp +from .list_tools import argmax, argmin + + +def slice_in_two(polytope: ConvexPolytope, value, slice_axis_idx): + if polytope is None: + return (None, None) + else: + assert len(polytope.points[0]) == 2 + + x_lower, x_upper, _ = polytope.extents(polytope._axes[slice_axis_idx]) + + intersects = _find_intersects(polytope, slice_axis_idx, value) + + if len(intersects) == 0: + if x_upper <= value: + # The vertical slicing line does not intersect the polygon, which is on the left of the line + # So we keep the same polygon for now since it is unsliced + left_polygon = polytope + right_polygon = None + if value < x_lower: + left_polygon = None + right_polygon = polytope + else: + left_points = [p for p in polytope.points if p[slice_axis_idx] <= value] + right_points = [p for p in polytope.points if p[slice_axis_idx] >= value] + left_points.extend(intersects) + right_points.extend(intersects) + # find left polygon + try: + hull = scipy.spatial.ConvexHull(left_points) + vertices = hull.vertices + except scipy.spatial.qhull.QhullError as e: + if "less than" or "is flat" in str(e): + # NOTE: this happens when we slice a polygon that has a border which coincides with the quadrant + # line and we slice this additional border with the quadrant line again. + # This is not actually a polygon we want to consider so we ignore it + vertices = None + + if vertices is not None: + left_polygon = ConvexPolytope(polytope._axes, [left_points[i] for i in vertices]) + else: + left_polygon = None + + try: + hull = scipy.spatial.ConvexHull(right_points) + vertices = hull.vertices + except scipy.spatial.qhull.QhullError as e: + # NOTE: this happens when we slice a polygon that has a border which coincides with the quadrant + # line and we slice this additional border with the quadrant line again. + # This is not actually a polygon we want to consider so we ignore it + if "less than" or "is flat" in str(e): + vertices = None + + if vertices is not None: + right_polygon = ConvexPolytope(polytope._axes, [right_points[i] for i in vertices]) + else: + right_polygon = None + + return (left_polygon, right_polygon) + + +def _find_intersects(polytope, slice_axis_idx, value): + intersects = [] + # Find all points above and below slice axis + above_slice = [p for p in polytope.points if p[slice_axis_idx] >= value] + below_slice = [p for p in polytope.points if p[slice_axis_idx] <= value] + + # Get the intersection of every pair above and below, this will create excess interior points + for a in above_slice: + for b in below_slice: + # edge is incident with slice plane, don't need these points + if a[slice_axis_idx] == b[slice_axis_idx]: + intersects.append(b) + continue + + # Linearly interpolate all coordinates of two points (a,b) of the polytope + interp_coeff = (value - b[slice_axis_idx]) / (a[slice_axis_idx] - b[slice_axis_idx]) + intersect = lerp(a, b, interp_coeff) + intersects.append(intersect) + return intersects + + +def _reduce_dimension(intersects, slice_axis_idx): + temp_intersects = [] + for point in intersects: + point = [p for i, p in enumerate(point) if i != slice_axis_idx] + temp_intersects.append(point) + return temp_intersects + + +def slice(polytope: ConvexPolytope, axis, value, slice_axis_idx): + # TODO: maybe these functions should go in the slicing tools? + if polytope.is_flat: + if value in chain(*polytope.points): + intersects = [[value]] + else: + return None + else: + intersects = _find_intersects(polytope, slice_axis_idx, value) + + if len(intersects) == 0: + return None + + # Reduce dimension of intersection points, removing slice axis + intersects = _reduce_dimension(intersects, slice_axis_idx) + + axes = copy(polytope._axes) + axes.remove(axis) + + if len(intersects) < len(intersects[0]) + 1: + return ConvexPolytope(axes, intersects) + # Compute convex hull (removing interior points) + if len(intersects[0]) == 0: + return None + elif len(intersects[0]) == 1: # qhull doesn't like 1D, do it ourselves + amin = argmin(intersects) + amax = argmax(intersects) + vertices = [amin, amax] + else: + try: + hull = scipy.spatial.ConvexHull(intersects) + vertices = hull.vertices + + except scipy.spatial.qhull.QhullError as e: + if "less than" or "flat" in str(e): + return ConvexPolytope(axes, intersects) + # Sliced result is simply the convex hull + return ConvexPolytope(axes, [intersects[i] for i in vertices]) diff --git a/pyproject.toml b/pyproject.toml index d9a2e1f4d..50a2dbbcb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,4 +5,5 @@ requires = ["setuptools>=61.0"] build-backend = "setuptools.build_meta" [tool.pytest.ini_options] markers = ["internet: downloads test data from the internet (deselect with '-m \"not internet\"')", - "fdb: test which uses fdb (deselect with '-m \"not fdb\"')",] \ No newline at end of file + "fdb: test which uses fdb (deselect with '-m \"not fdb\"')",] +log_cli_level = "INFO" \ No newline at end of file diff --git a/tests/profiled_quadtree.profile b/tests/profiled_quadtree.profile new file mode 100644 index 000000000..597957624 Binary files /dev/null and b/tests/profiled_quadtree.profile differ diff --git a/tests/quadtree_slicer_profiler.py b/tests/quadtree_slicer_profiler.py new file mode 100644 index 000000000..ecb3f67a9 --- /dev/null +++ b/tests/quadtree_slicer_profiler.py @@ -0,0 +1,59 @@ +import cProfile + +import numpy as np +import pygribjump as gj + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box + +options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "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": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, +} +fdbdatacube = gj.GribJump() + +x = np.linspace(0, 100, 1000) +y = np.linspace(0, 100, 1000) +# create the mesh based on these arrays +X, Y = np.meshgrid(x, y) +X = X.reshape((np.prod(X.shape),)) +Y = Y.reshape((np.prod(Y.shape),)) +coords = zip(X, Y) +points = [list(coord) for coord in coords] +polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] +API = Polytope( + request=Request(polytope), + datacube=fdbdatacube, + options=options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, +) +cProfile.runctx( + "API.engines['quadtree'].extract(API.datacube, [polytope])", globals(), locals(), "profiled_extract_quadtree.pstats" +) diff --git a/tests/test_bad_request_error.py b/tests/test_bad_request_error.py index b1154cf9c..dc2fc8fc2 100644 --- a/tests/test_bad_request_error.py +++ b/tests/test_bad_request_error.py @@ -60,6 +60,5 @@ def test_fdb_datacube(self): self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_cyclic_axis_over_negative_vals.py b/tests/test_cyclic_axis_over_negative_vals.py index 37bddd75f..5ba63c2e8 100644 --- a/tests/test_cyclic_axis_over_negative_vals.py +++ b/tests/test_cyclic_axis_over_negative_vals.py @@ -28,7 +28,7 @@ def setup_method(self, method): "compressed_axes_config": ["long", "level", "step", "date"], } self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_cyclic_axis_slicer_not_0.py b/tests/test_cyclic_axis_slicer_not_0.py index 84d36d9ce..b6ca580b8 100644 --- a/tests/test_cyclic_axis_slicer_not_0.py +++ b/tests/test_cyclic_axis_slicer_not_0.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["long", "level", "step", "date"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_cyclic_axis_slicing.py b/tests/test_cyclic_axis_slicing.py index 3f8d53e8d..0f84c2c94 100644 --- a/tests/test_cyclic_axis_slicing.py +++ b/tests/test_cyclic_axis_slicing.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["long", "level", "step", "date"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_cyclic_nearest.py b/tests/test_cyclic_nearest.py index 86e37c69c..7622eb915 100644 --- a/tests/test_cyclic_nearest.py +++ b/tests/test_cyclic_nearest.py @@ -3,7 +3,6 @@ from eccodes import codes_grib_find_nearest, codes_grib_new_from_file from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -85,10 +84,8 @@ def test_regular_grid(self): Point(["latitude", "longitude"], [[39, 360 - 76.45]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -109,10 +106,8 @@ def test_regular_grid(self): Point(["latitude", "longitude"], [[39, -76.45]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_cyclic_simple.py b/tests/test_cyclic_simple.py index d9974b7e7..538932a6b 100644 --- a/tests/test_cyclic_simple.py +++ b/tests/test_cyclic_simple.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["long", "level", "step", "date"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=options, ) diff --git a/tests/test_cyclic_snapping.py b/tests/test_cyclic_snapping.py index 599dc9573..5c55eaeea 100644 --- a/tests/test_cyclic_snapping.py +++ b/tests/test_cyclic_snapping.py @@ -21,7 +21,7 @@ def setup_method(self, method): "compressed_axes_config": ["long"], } self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_datacube_axes_init.py b/tests/test_datacube_axes_init.py index d92a08cf3..5dc6b329a 100644 --- a/tests/test_datacube_axes_init.py +++ b/tests/test_datacube_axes_init.py @@ -3,7 +3,6 @@ from helper_functions import download_test_data from polytope_feature.datacube.datacube_axis import FloatDatacubeAxis -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -41,10 +40,8 @@ def setup_method(self, method): "type", ], } - self.slicer = HullSlicer() self.API = Polytope( datacube=latlon_array, - engine=self.slicer, options=self.options, ) self.datacube = self.API.datacube diff --git a/tests/test_datacube_xarray.py b/tests/test_datacube_xarray.py index f15a4f8bc..8f66eb822 100644 --- a/tests/test_datacube_xarray.py +++ b/tests/test_datacube_xarray.py @@ -24,7 +24,7 @@ def test_validate(self): array = xr.Dataset(data_vars=dict(param=(["x", "y", "z"], dims)), coords={"x": [1], "y": [1], "z": [1]}) array = array.to_array() - datacube = Datacube.create(array, axis_options={}) + datacube = Datacube.create(datacube=array, axis_options={}) datacube.validate(["x", "y", "z", "variable"]) datacube.validate(["x", "z", "y", "variable"]) @@ -53,7 +53,7 @@ def test_create(self): for d, v in array.coords.variables.items(): print(v.dtype) - datacube = Datacube.create(array, axis_options={}) + datacube = Datacube.create(datacube=array, axis_options={}) # Check the factory created the correct type of datacube assert isinstance(datacube, XArrayDatacube) diff --git a/tests/test_ecmwf_oper_data_fdb.py b/tests/test_ecmwf_oper_data_fdb.py index e848716c8..381cc566f 100644 --- a/tests/test_ecmwf_oper_data_fdb.py +++ b/tests/test_ecmwf_oper_data_fdb.py @@ -1,7 +1,6 @@ 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, Point, Select, Span @@ -59,10 +58,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -87,10 +84,8 @@ def test_fdb_datacube_point(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -116,10 +111,8 @@ def test_fdb_datacube_point_v2(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -174,10 +167,8 @@ def test_fdb_datacube_point_step_not_compressed(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_engine_slicer.py b/tests/test_engine_slicer.py index 3392ccd3a..738b462eb 100644 --- a/tests/test_engine_slicer.py +++ b/tests/test_engine_slicer.py @@ -1,17 +1,17 @@ from polytope_feature.datacube.backends.mock import MockDatacube -from polytope_feature.engine.hullslicer import HullSlicer +from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box, Polygon class TestEngineSlicer: def setup_method(self, method): - self.slicer = HullSlicer() + pass def test_2D_box(self): datacube = MockDatacube({"x": 100, "y": 100}) polytopes = Box(["x", "y"], lower_corner=[3, 3], upper_corner=[6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 total_leaves = 0 for leaf in result.leaves: @@ -21,8 +21,8 @@ def test_2D_box(self): def test_3D_box(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100}) polytopes = Box(["x", "y", "z"], lower_corner=[3, 3, 3], upper_corner=[6, 6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 * 4 total_leaves = 0 for leaf in result.leaves: @@ -32,8 +32,8 @@ def test_3D_box(self): def test_4D_box(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100, "q": 100}) polytopes = Box(["x", "y", "z", "q"], lower_corner=[3, 3, 3, 3], upper_corner=[6, 6, 6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 * 4 * 4 total_leaves = 0 for leaf in result.leaves: @@ -43,44 +43,40 @@ def test_4D_box(self): def test_triangle(self): datacube = MockDatacube({"x": 100, "y": 100}) 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 - # total_leaves = 0 - # for leaf in result.leaves: - # total_leaves += len(leaf.values) - # assert total_leaves == 4 + 3 + 2 + 1 + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, triangle) + assert len(result.leaves) == 4 + 3 + 2 + 1 def test_reusable(self): datacube = MockDatacube({"x": 100, "y": 100}) 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 + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) + assert len(result.leaves) == 4 + 3 + 2 + 1 polytopes = Box(["x", "y"], lower_corner=[3, 3], upper_corner=[6, 6]).polytope() - result = self.slicer.extract(datacube, polytopes) - result.pprint() + result = API.slice(datacube, polytopes) assert len(result.leaves) == 4 def test_2D_box_get_function(self): datacube = MockDatacube({"x": 100, "y": 100}) polytopes = Box(["x", "y"], lower_corner=[2, -2], upper_corner=[4, -1]).polytope() - result = self.slicer.extract(datacube, polytopes) + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) datacube.get(result) result.pprint() def test_3D_box_get_function(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100}) polytopes = Box(["x", "y", "z"], lower_corner=[3, 2, -2], upper_corner=[6, 2, -1]).polytope() - result = self.slicer.extract(datacube, polytopes) + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) datacube.get(result) result.pprint() def test_3D_box_get_function2(self): datacube = MockDatacube({"x": 100, "y": 100, "z": 100}) polytopes = Box(["x", "y", "z"], lower_corner=[3, 2, 1], upper_corner=[6, 2, 1]).polytope() - result = self.slicer.extract(datacube, polytopes) + API = Polytope(datacube=datacube, options={}) + result = API.slice(datacube, polytopes) datacube.get(result) result.pprint() diff --git a/tests/test_fdb_datacube.py b/tests/test_fdb_datacube.py index 301ceec30..0b20ce78a 100644 --- a/tests/test_fdb_datacube.py +++ b/tests/test_fdb_datacube.py @@ -1,7 +1,6 @@ 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, Select, Span @@ -63,10 +62,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -114,10 +111,8 @@ def test_fdb_datacube_select_grid(self): Span("longitude", 0, 0.070093457944), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_fdb_unmap_tree.py b/tests/test_fdb_unmap_tree.py index a9dcf9cbe..c875433e7 100644 --- a/tests/test_fdb_unmap_tree.py +++ b/tests/test_fdb_unmap_tree.py @@ -1,7 +1,6 @@ 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, Select @@ -60,10 +59,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_float_type.py b/tests/test_float_type.py index 51df70df8..9916d6b1f 100644 --- a/tests/test_float_type.py +++ b/tests/test_float_type.py @@ -2,7 +2,6 @@ import pytest import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select, Span @@ -19,9 +18,8 @@ def setup_method(self, method): "alt": np.arange(0.0, 20.0, 0.1), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["lat", "long", "alt"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_slicing_span(self): request = Request(Span("lat", 4.1, 4.3), Select("long", [4.1]), Select("alt", [4.1])) diff --git a/tests/test_healpix_mapper.py b/tests/test_healpix_mapper.py index 8badbffc6..b2c9a7863 100644 --- a/tests/test_healpix_mapper.py +++ b/tests/test_healpix_mapper.py @@ -5,7 +5,6 @@ from polytope_feature.datacube.transformations.datacube_mappers.mapper_types.healpix import ( HealpixGridMapper, ) -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -30,10 +29,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["longitude", "latitude", "step", "time", "isobaricInhPa", "valid_time"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.latlon_array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_hullslicer_engine.py b/tests/test_hullslicer_engine.py index 198b96b25..467ba3723 100644 --- a/tests/test_hullslicer_engine.py +++ b/tests/test_hullslicer_engine.py @@ -3,7 +3,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box @@ -19,14 +18,13 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(self.array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["step", "level"]} - self.API = Polytope(datacube=self.array, engine=self.slicer, options=options) + self.API = Polytope(datacube=self.array, options=options) def test_extract(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) assert request.axis == TensorIndexTree.root assert request.parent is None assert request.values is tuple() diff --git a/tests/test_icon_grid_unstructured.py b/tests/test_icon_grid_unstructured.py new file mode 100644 index 000000000..97bf46459 --- /dev/null +++ b/tests/test_icon_grid_unstructured.py @@ -0,0 +1,123 @@ +# import geopandas as gpd +# import matplotlib.pyplot as plt +# import numpy as np +# import iris +import math + +import pandas as pd +import pytest +import xarray as xr +from earthkit import data + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + +# from helper_functions import find_nearest_latlon + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "time": "hullslicer", + "heightAboveGround": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + + ds = data.from_source("file", "../../Downloads/icon_global_icosahedral_single-level_2025011000_000_T_2M.grib2") + grid = xr.open_dataset("../../Downloads/icon_grid_0026_R03B07_G.nc", engine="netcdf4") + + self.arr = ds.to_xarray(engine="cfgrib").t2m + self.longitudes = grid.clon.values * (180 / math.pi) + self.latitudes = grid.clat.values * (180 / math.pi) + self.points = list(zip(self.latitudes, self.longitudes)) + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + {"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [-180, 180]}]}, + ], + } + + self.API = Polytope( + datacube=self.arr, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + import datetime + + request = Request( + Select("time", [pd.Timestamp("2025-01-10")]), + Select("heightAboveGround", [2.0]), + Select("step", [datetime.timedelta(days=0)]), + Box(["latitude", "longitude"], [0, 0], [1, 5]), + ) + + result = self.API.retrieve(request) + assert len(result.leaves) == 345 + result.pprint() + + # lats = [] + # lons = [] + # eccodes_lats = [] + # eccodes_lons = [] + # tol = 1e-8 + # leaves = result.leaves + # for i in range(len(leaves)): + # cubepath = leaves[i].flatten() + # lat = cubepath["latitude"][0] + # lon = cubepath["longitude"][0] + # lats.append(lat) + # lons.append(lon) + + # # # each variable in the netcdf file is a cube + # # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # # nearest_points = find_nearest_latlon( + # # "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2", lat, lon) + # # eccodes_lat = nearest_points[0][0]["lat"] + # # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # # eccodes_lats.append(eccodes_lat) + # # eccodes_lons.append(eccodes_lon) + # # assert eccodes_lat - tol <= lat + # # assert lat <= eccodes_lat + tol + # # assert eccodes_lon - tol <= lon + # # assert lon <= eccodes_lon + tol + + # 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(lons, lats, s=18, c="red", cmap="YlOrRd") + # plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + # plt.colorbar(label="Temperature") + # plt.show() + + @pytest.mark.fdb + def test_quad_tree_slicer_extract_cyclic(self): + import datetime + + request = Request( + Select("time", [pd.Timestamp("2025-01-10")]), + Select("heightAboveGround", [2.0]), + Select("step", [datetime.timedelta(days=0)]), + Box(["latitude", "longitude"], [0, 190], [1, 195]), + ) + + result = self.API.retrieve(request) + assert len(result.leaves) == 346 + result.pprint() + leaves = result.leaves + for i in range(len(leaves)): + cubepath = leaves[i].flatten() + lon = cubepath["longitude"][0] + assert -180 <= lon <= 180 diff --git a/tests/test_icon_grid_unstructured_fdb.py b/tests/test_icon_grid_unstructured_fdb.py new file mode 100644 index 000000000..327f5f45e --- /dev/null +++ b/tests/test_icon_grid_unstructured_fdb.py @@ -0,0 +1,123 @@ +import math + +# import os +import time + +# import geopandas as gpd +# import matplotlib.pyplot as plt +import pandas as pd +import pytest +import xarray as xr +from earthkit import data + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + +# from helper_functions import find_nearest_latlon + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "date": "hullslicer", + "levtype": "hullslicer", + "param": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + + ds = data.from_source("file", "../../Downloads/icon_global_icosahedral_single-level_2025011000_000_T_2M.grib2") + grid = xr.open_dataset("../../Downloads/icon_grid_0026_R03B07_G.nc", engine="netcdf4") + + self.arr = ds.to_xarray(engine="cfgrib").t2m + self.longitudes = grid.clon.values * 180 / math.pi + self.latitudes = grid.clat.values * 180 / math.pi + self.points = list(zip(self.latitudes, self.longitudes)) + self.options = { + "axis_config": [ + {"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": "irregular", + "resolution": 0, + "axes": ["latitude", "longitude"], + "md5_hash": "f68071a8ac9bae4e965822afb963c04f", + } + ], + }, + ], + "pre_path": {"date": "20250110"}, + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + import pygribjump as gj + + request = Request( + Select("date", [pd.Timestamp("20250110T0000")]), + Select("step", [0]), + Select("param", ["167"]), + Select("levtype", ["sfc"]), + Box(["latitude", "longitude"], [0, 0], [10, 10]), + ) + + self.fdbdatacube = gj.GribJump() + + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + time0 = time.time() + result = self.API.retrieve(request) + time1 = time.time() + assert len(result.leaves) == 1 + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + result.pprint() + + # lats = [] + # lons = [] + # eccodes_lats = [] + # eccodes_lons = [] + # # tol = 1e-8 + # leaves = result.leaves + # for i in range(len(leaves)): + # cubepath = leaves[i].flatten() + # lat = cubepath["latitude"][0] + # lon = cubepath["longitude"][0] + # lats.append(lat) + # lons.append(lon) + + # # # each variable in the netcdf file is a cube + # # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # # nearest_points = find_nearest_latlon( + # # "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2", lat, lon) + # # eccodes_lat = nearest_points[0][0]["lat"] + # # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # # eccodes_lats.append(eccodes_lat) + # # eccodes_lons.append(eccodes_lon) + # # assert eccodes_lat - tol <= lat + # # assert lat <= eccodes_lat + tol + # # assert eccodes_lon - tol <= lon + # # assert lon <= eccodes_lon + tol + + # 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(lons, lats, s=18, c="red", cmap="YlOrRd") + # plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + # plt.colorbar(label="Temperature") + # plt.show() diff --git a/tests/test_icon_grid_unstructured_fdb_point.py b/tests/test_icon_grid_unstructured_fdb_point.py new file mode 100644 index 000000000..6ee53ded6 --- /dev/null +++ b/tests/test_icon_grid_unstructured_fdb_point.py @@ -0,0 +1,141 @@ +import math + +# import iris +import os +import time + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Point, Select + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "date": "hullslicer", + "levtype": "hullslicer", + "param": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + print("SETTING UP THE XARRAY") + time_now = time.time() + + # ds = data.from_source( + # "file", "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2") + + # grid = xr.open_dataset("../../Downloads/icon_grid_0047_R19B07_L.nc", engine="netcdf4") + + ds = data.from_source( + "file", "../../Downloads/icon_global_icosahedral_single-level_2025011000_000_T_2M.grib2") + + grid = xr.open_dataset("../../Downloads/icon_grid_0026_R03B07_G.nc", engine="netcdf4") + + print(time.time()-time_now) + self.arr = ds.to_xarray(engine="cfgrib").t2m + + self.longitudes = grid.clon.values * 180/math.pi + self.latitudes = grid.clat.values * 180/math.pi + + self.points = list(zip(self.latitudes, self.longitudes)) + print((min(self.latitudes), max(self.latitudes), min(self.longitudes), max(self.longitudes))) + print("FINISH SETTING UP POINTS") + self.options = { + "axis_config": [ + {"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": "irregular", "resolution": 0, "axes": [ + "latitude", "longitude"], "md5_hash": "f68071a8ac9bae4e965822afb963c04f"} + ], + }, + ], + # "pre_path": {"time": "20250110", "heightAboveGround": "2"}, + "pre_path": {"date": "20250110"}, + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + import datetime + + import pygribjump as gj + request = Request( + # Select("deptht", [0.5058], method="nearest"), + Select("date", [pd.Timestamp("20250110T0000")]), + # Select("heightAboveGround", [2.0]), + # Select("step", [datetime.timedelta(days=0)]), + Select("step", [0]), + Select("param", ["167"]), + Select("levtype", ["sfc"]), + # Select("time_counter", [pd.Timestamp("1979-02-15")]), + # Box(["latitude", "longitude"], [0, 0], [80, 80]), + Point(["latitude", "longitude"], [[0, 0]], method="nearest"), + ) + + self.fdbdatacube = gj.GribJump() + + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + time0 = time.time() + result = self.API.retrieve(request) + # result = self.API.slice(self.API.datacube, request.polytopes()) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + print(len(result.leaves)) + result.pprint() + + lats = [] + lons = [] + eccodes_lats = [] + eccodes_lons = [] + tol = 1e-8 + leaves = result.leaves + for i in range(len(leaves)): + cubepath = leaves[i].flatten() + lat = cubepath["latitude"][0] + lon = cubepath["longitude"][0] + lats.append(lat) + lons.append(lon) + + # # each variable in the netcdf file is a cube + # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # nearest_points = find_nearest_latlon( + # "../../Downloads/icon-d2_germany_icosahedral_single-level_2025011000_024_2d_t_2m.grib2", lat, lon) + # eccodes_lat = nearest_points[0][0]["lat"] + # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # eccodes_lats.append(eccodes_lat) + # eccodes_lons.append(eccodes_lon) + # assert eccodes_lat - tol <= lat + # assert lat <= eccodes_lat + tol + # assert eccodes_lon - tol <= lon + # assert lon <= eccodes_lon + tol + + 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(lons, lats, s=18, c="red", cmap="YlOrRd") + plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + plt.colorbar(label="Temperature") + plt.show() diff --git a/tests/test_incomplete_tree_fdb.py b/tests/test_incomplete_tree_fdb.py index 2fbfcaf28..b388d37fd 100644 --- a/tests/test_incomplete_tree_fdb.py +++ b/tests/test_incomplete_tree_fdb.py @@ -2,7 +2,6 @@ import pytest from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -67,10 +66,8 @@ def test_incomplete_fdb_branch(self): Select("number", ["0"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -99,10 +96,8 @@ def test_incomplete_fdb_branch_2(self): Select("number", ["0"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_local_grid_cyclic.py b/tests/test_local_grid_cyclic.py index 3f4bfd36c..32911dfce 100644 --- a/tests/test_local_grid_cyclic.py +++ b/tests/test_local_grid_cyclic.py @@ -1,7 +1,6 @@ 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 Point, Select @@ -68,10 +67,8 @@ def test_fdb_datacube(self): Point(["latitude", "longitude"], [[-20, -20]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -97,10 +94,8 @@ def test_fdb_datacube_2(self): Point(["latitude", "longitude"], [[-20, 50 + 360]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_local_regular_grid.py b/tests/test_local_regular_grid.py index 3f8f2d402..ac23bf172 100644 --- a/tests/test_local_regular_grid.py +++ b/tests/test_local_regular_grid.py @@ -1,7 +1,6 @@ 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 Point, Select @@ -67,10 +66,8 @@ def test_fdb_datacube(self): Point(["latitude", "longitude"], [[0.16, 0.176]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -96,10 +93,8 @@ def test_point_outside_local_region(self): Point(["latitude", "longitude"], [[0.16, 61]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -125,10 +120,8 @@ def test_point_outside_local_region_2(self): Point(["latitude", "longitude"], [[41, 1]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -154,10 +147,8 @@ def test_point_outside_local_region_3(self): Point(["latitude", "longitude"], [[1, 61]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -182,10 +173,8 @@ def test_point_outside_local_region_4(self): Point(["latitude", "longitude"], [[41, 1]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -210,10 +199,8 @@ def test_point_outside_local_region_5(self): Point(["latitude", "longitude"], [[-41, 1]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -238,10 +225,8 @@ def test_point_outside_local_region_6(self): Point(["latitude", "longitude"], [[-30, -21]]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -266,10 +251,8 @@ def test_point_outside_local_region_7(self): Point(["latitude", "longitude"], [[-41, 1]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -295,10 +278,8 @@ def test_point_outside_local_region_8(self): Point(["latitude", "longitude"], [[-30, -21]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -324,10 +305,8 @@ def test_point_outside_local_region_9(self): Point(["latitude", "longitude"], [[-30, -21]], method="surrounding"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_merge_octahedral_one_axis.py b/tests/test_merge_octahedral_one_axis.py index 4ddf90582..cde50d1d3 100644 --- a/tests/test_merge_octahedral_one_axis.py +++ b/tests/test_merge_octahedral_one_axis.py @@ -2,7 +2,6 @@ from earthkit import data from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -28,10 +27,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["longitude", "latitude", "surface", "step", "time", "valid_time", "number"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.latlon_array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_merge_transformation.py b/tests/test_merge_transformation.py index 5d2311c35..99fc6d926 100644 --- a/tests/test_merge_transformation.py +++ b/tests/test_merge_transformation.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["date", "time"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_multiple_param_fdb.py b/tests/test_multiple_param_fdb.py index 11285af84..cd82d4731 100644 --- a/tests/test_multiple_param_fdb.py +++ b/tests/test_multiple_param_fdb.py @@ -1,7 +1,6 @@ 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, Select @@ -59,10 +58,8 @@ def test_fdb_datacube(self): Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_octahedral_grid.py b/tests/test_octahedral_grid.py index f44ccf337..429794745 100644 --- a/tests/test_octahedral_grid.py +++ b/tests/test_octahedral_grid.py @@ -2,7 +2,6 @@ from earthkit import data from helper_functions import download_test_data, find_nearest_latlon -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -27,10 +26,8 @@ def setup_method(self, method): ], "compressed_axes_config": ["longitude", "latitude", "number", "step", "time", "surface", "valid_time"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=self.latlon_array, - engine=self.slicer, options=self.options, ) diff --git a/tests/test_orca_irregular_grid.py b/tests/test_orca_irregular_grid.py new file mode 100644 index 000000000..25665040a --- /dev/null +++ b/tests/test_orca_irregular_grid.py @@ -0,0 +1,94 @@ +# import geopandas as gpd +# import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +from earthkit import data +from helper_functions import find_nearest_latlon + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "step": "hullslicer", + "time": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + "oceanModelLayer": "hullslicer", + "valid_time": "hullslicer", + } + + ds = data.from_source("file", "../../Downloads/Reference_eORCA12_U_to_HEALPix_32.grib") + self.arr = ds.to_xarray(engine="cfgrib").avg_uox + + self.latitudes = self.arr.latitude.values + self.longitudes = self.arr.longitude.values + self.points = list(zip(self.latitudes, self.longitudes)) + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + request = Request( + Select("step", [np.timedelta64(0, "ns")]), + Select("oceanModelLayer", [65.0]), + Select("time", [pd.Timestamp("2017-09-06T00:00:00.000000000")]), + Box(["latitude", "longitude"], [65, 270], [75, 300]), + ) + + self.API = Polytope( + datacube=self.arr, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + import time + + time0 = time.time() + result = self.API.retrieve(request) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + assert len(result.leaves) == 27 + result.pprint() + + lats = [] + lons = [] + eccodes_lats = [] + eccodes_lons = [] + tol = 1e-8 + for i in range(len(result.leaves)): + cubepath = result.leaves[i].flatten() + lat = cubepath["latitude"][0] + lon = cubepath["longitude"][0] - 360 + lats.append(lat) + lons.append(lon) + nearest_points = find_nearest_latlon("../../Downloads/Reference_eORCA12_U_to_HEALPix_32.grib", lat, lon) + eccodes_lat = nearest_points[0][0]["lat"] + eccodes_lon = nearest_points[0][0]["lon"] - 360 + eccodes_lats.append(eccodes_lat) + eccodes_lons.append(eccodes_lon) + assert eccodes_lat - tol <= lat + assert lat <= eccodes_lat + tol + assert eccodes_lon - tol <= lon + assert lon <= eccodes_lon + tol + + # 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(lons, lats, s=18, c="red", cmap="YlOrRd") + # plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + # plt.colorbar(label="Temperature") + # plt.show() diff --git a/tests/test_orca_irregular_grid_seas5.py b/tests/test_orca_irregular_grid_seas5.py new file mode 100644 index 000000000..e7e70caba --- /dev/null +++ b/tests/test_orca_irregular_grid_seas5.py @@ -0,0 +1,129 @@ +# from helper_functions import find_nearest_latlon +import time + +# import geopandas as gpd +# import matplotlib.pyplot as plt +import numpy as np + +# import pandas as pd +import pytest + +# from earthkit import data +import xarray as xr + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + + +class TestQuadTreeSlicer: + def setup_method(self, method): + self.engine_options = { + "deptht": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + + ds = xr.open_dataset("../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc", engine="netcdf4") + self.arr = ds.votemper + + self.latitudes = self.arr.coords["nav_lat"].values + lats = [] + for lat in self.latitudes: + lats.extend(lat) + self.latitudes = lats + self.longitudes = self.arr.nav_lon.values + lons = [] + for lon in self.longitudes: + lons.extend(lon) + self.longitudes = lons + + # # Drop the x and y dimensions + nav_lat_flat = self.arr.nav_lat.values.ravel() + deptht_flat = self.arr.deptht.values.ravel() + interm_data = self.arr.data[0] + new_interm_data = [] + for data in interm_data: + new_data = data.ravel() + new_interm_data.append(new_data) + + # Create a new dimension, for example, "grid_index" + grid_index = np.arange(nav_lat_flat.size) + + self.arr = xr.DataArray( + new_interm_data, + dims=("deptht", "values"), + coords={ + "deptht": deptht_flat, + "values": grid_index, + }, + ) + + self.points = list(zip(self.latitudes, self.longitudes)) + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + request = Request( + Select("deptht", [0.5058], method="nearest"), + Box(["latitude", "longitude"], [-55, 50], [-50, 55]), + ) + + self.API = Polytope( + datacube=self.arr, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.points, + ) + + time0 = time.time() + result = self.API.retrieve(request) + time1 = time.time() + print("TIME TAKEN TO EXTRACT") + print(time1 - time0) + assert len(result.leaves) == 1386 + result.pprint() + + # lats = [] + # lons = [] + # eccodes_lats = [] + # eccodes_lons = [] + # # tol = 1e-8 + # leaves = result.leaves + # for i in range(len(leaves)): + # cubepath = leaves[i].flatten() + # lat = cubepath["latitude"][0] + # lon = cubepath["longitude"][0] + # lats.append(lat) + # lons.append(lon) + + # # # each variable in the netcdf file is a cube + # # # cubes = iris.load('../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.nc') + # # # iris.save(cubes, '../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2') + # # nearest_points = find_nearest_latlon("../../Downloads/votemper_ORAS5_1m_197902_grid_T_02.grib2", lat, + # lon) + # # eccodes_lat = nearest_points[0][0]["lat"] + # # eccodes_lon = nearest_points[0][0]["lon"] - 360 + # # eccodes_lats.append(eccodes_lat) + # # eccodes_lons.append(eccodes_lon) + # # assert eccodes_lat - tol <= lat + # # assert lat <= eccodes_lat + tol + # # assert eccodes_lon - tol <= lon + # # assert lon <= eccodes_lon + tol + + # 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(lons, lats, s=18, c="red", cmap="YlOrRd") + # plt.scatter(eccodes_lons, eccodes_lats, s=6, c="green") + # plt.colorbar(label="Temperature") + # plt.show() diff --git a/tests/test_override_md5_hash_options.py b/tests/test_override_md5_hash_options.py index 0233290af..237dc2ede 100644 --- a/tests/test_override_md5_hash_options.py +++ b/tests/test_override_md5_hash_options.py @@ -1,7 +1,6 @@ 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, Select, Span @@ -67,10 +66,8 @@ def test_fdb_datacube(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_point_nearest.py b/tests/test_point_nearest.py index b9fdc4bed..7b2e605b3 100644 --- a/tests/test_point_nearest.py +++ b/tests/test_point_nearest.py @@ -1,7 +1,6 @@ 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 Point, Select @@ -60,10 +59,8 @@ def test_fdb_datacube(self): Point(["latitude", "longitude"], [[0.16, 0.176]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -86,10 +83,8 @@ def test_fdb_datacube_true_point(self): Point(["latitude", "longitude"], [[0.175746921078, 0.210608424337]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -113,10 +108,8 @@ def test_fdb_datacube_true_point_2(self): Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -140,10 +133,8 @@ def test_fdb_datacube_true_point_3(self): Point(["latitude", "longitude"], [[0.035149384216, -0.01]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -169,10 +160,8 @@ def test_fdb_datacube_true_point_5(self): Point(["latitude", "longitude"], [[0.035149384216, 360 - 0.01]], method="nearest"), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -198,14 +187,11 @@ def test_fdb_datacube_true_point_4(self): Point(["latitude", "longitude"], [[0.035149384216, 359.97]], method="nearest"), ) 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_2() assert len(result.leaves) == 1 assert result.leaves[0].values == (359.929906542056,) assert result.leaves[0].axis.name == "longitude" diff --git a/tests/test_point_shape.py b/tests/test_point_shape.py index d4d6eb98f..03a261217 100644 --- a/tests/test_point_shape.py +++ b/tests/test_point_shape.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Point, Select @@ -19,9 +18,8 @@ def setup_method(self, method): "level": range(1, 130), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step", "date"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_point(self): request = Request(Point(["step", "level"], [[3, 10]]), Select("date", ["2000-01-01"])) diff --git a/tests/test_polytope_extract.py b/tests/test_polytope_extract.py new file mode 100644 index 000000000..2bc7400fa --- /dev/null +++ b/tests/test_polytope_extract.py @@ -0,0 +1,51 @@ +import numpy as np +import xarray as xr + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box + + +class TestPolytopeExtract: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + self.array = xr.DataArray( + np.random.randn(6, 129, 100), + dims=("step", "level", "values"), + coords={ + "step": [0, 3, 6, 9, 12, 15], + "level": range(1, 130), + "values": range(0, 100), + }, + ) + self.engine_options = { + "step": "hullslicer", + "level": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + } + self.quadtree_points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + self.options = { + "axis_config": [ + { + "axis_name": "values", + "transformations": [ + {"name": "mapper", "type": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + } + + # Testing different shapes + + def test_2D_box(self): + request = Request(Box(["step", "level"], [3, 10], [6, 11]), Box(["latitude", "longitude"], [0, 0], [20, 20])) + self.API = Polytope( + datacube=self.array, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.quadtree_points, + ) + result = self.API.retrieve(request) + + result.pprint() + assert len(result.leaves) == 12 diff --git a/tests/test_polytope_extract_fdb.py b/tests/test_polytope_extract_fdb.py new file mode 100644 index 000000000..d466581b6 --- /dev/null +++ b/tests/test_polytope_extract_fdb.py @@ -0,0 +1,95 @@ +import pandas as pd +import pytest + +from polytope_feature.polytope import Polytope, Request +from polytope_feature.shapes import Box, Select + + +class TestPolytopeExtract: + def setup_method(self, method): + # from polytope.datacube.backends.fdb import FDBDatacube + + # Create a dataarray with 3 labelled axes using different index types + self.engine_options = { + "step": "hullslicer", + "levtype": "hullslicer", + "latitude": "quadtree", + "longitude": "quadtree", + "class": "hullslicer", + "date": "hullslicer", + "type": "hullslicer", + "stream": "hullslicer", + "param": "hullslicer", + "expver": "hullslicer", + "domain": "hullslicer", + } + self.quadtree_points = [[10, 10], [0.035149384216, 0.0], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "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": "irregular", + "axes": ["latitude", "longitude"], + "md5_hash": "158db321ae8e773681eeb40e0a3d350f", + } + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + + # Testing different shapes + @pytest.mark.fdb + def test_2D_box(self): + import pygribjump as gj + + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20230625T120000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["an"]), + Box(["latitude", "longitude"], [0, -0.1], [10, 10]), + ) + self.fdbdatacube = gj.GribJump() + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options=self.engine_options, + point_cloud_options=self.quadtree_points, + ) + result = self.API.retrieve(request) + + assert len(result.leaves) == 3 + assert result.leaves[0].flatten()["longitude"] == (0,) + assert result.leaves[0].flatten()["latitude"] == (0.035149384216,) + assert result.leaves[1].flatten()["longitude"] == (10,) + assert result.leaves[1].flatten()["latitude"] == (5,) + assert result.leaves[2].flatten()["longitude"] == (10,) + assert result.leaves[2].flatten()["latitude"] == (10,) diff --git a/tests/test_quad_tree.py b/tests/test_quad_tree.py new file mode 100644 index 000000000..23be020fe --- /dev/null +++ b/tests/test_quad_tree.py @@ -0,0 +1,164 @@ +import pytest + +from polytope_feature.datacube.quad_tree import QuadNode +from polytope_feature.engine.quadtree_slicer import QuadTreeSlicer +from polytope_feature.polytope import Polytope +from polytope_feature.shapes import Box, ConvexPolytope +from polytope_feature.utility.slicing_tools import slice_in_two + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "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": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + slicer.quad_tree.pprint() + pass + + @pytest.mark.fdb + def test_quad_tree_query_polygon(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + polytope = Box(["lat", "lon"], [1, 1], [20, 30]).polytope()[0] + results = slicer.quad_tree.query_polygon(polytope) + assert len(results) == 3 + assert (10, 10) in [node.item for node in results] + assert (5, 10) in [node.item for node in results] + assert (5, 20) in [node.item for node in results] + points = [[10, 10], [80, 10], [-5, 5], [5, 50], [5, 10], [50, 10], [2, 10], [15, 15]] + slicer = QuadTreeSlicer(points) + polytope = ConvexPolytope(["lat", "lon"], [[-10, 1], [20, 1], [5, 20]]) + results = slicer.quad_tree.query_polygon(polytope) + assert len(results) == 4 + assert (-5, 5) in [node.item for node in results] + assert (5, 10) in [node.item for node in results] + assert (10, 10) in [node.item for node in results] + assert (2, 10) in [node.item for node in results] + + @pytest.mark.fdb + def test_slice_in_two_vertically(self): + polytope = Box(["lat", "lon"], [0, 0], [2, 2]).polytope()[0] + lower, upper = slice_in_two(polytope, 1, 0) + assert lower.points == [[0, 0], [1.0, 0.0], [1.0, 2.0], [0, 2]] + assert upper.points == [[1.0, 0.0], [2, 0], [2, 2], [1.0, 2.0]] + + @pytest.mark.fdb + def test_slice_in_two_horizontally(self): + polytope = Box(["lat", "lon"], [0, 0], [2, 2]).polytope()[0] + lower, upper = slice_in_two(polytope, 1, 1) + assert lower.points == [[0, 0], [2, 0], [2.0, 1.0], [0.0, 1.0]] + assert upper.points == [[2, 2], [0, 2], [0.0, 1.0], [2.0, 1.0]] + + @pytest.mark.fdb + def test_quad_node_is_contained_in_box(self): + node = QuadNode([1, 1], 0) + polytope = Box(["lat", "lon"], [0, 0], [2, 2]).polytope()[0] + assert node.is_contained_in(polytope) + second_node = QuadNode([3, 3], 0) + assert not second_node.is_contained_in(polytope) + third_node = QuadNode([1, 0], 0) + assert third_node.is_contained_in(polytope) + + @pytest.mark.fdb + def test_quad_node_is_contained_in_triangle(self): + node = QuadNode([1, 1], 0) + polytope = ConvexPolytope(["lat", "lon"], [[0, 0], [1, 1], [2, 0]]) + assert node.is_contained_in(polytope) + node = QuadNode([1, 0.5], 0) + assert node.is_contained_in(polytope) + second_node = QuadNode([3, 3], 0) + assert not second_node.is_contained_in(polytope) + third_node = QuadNode([1, 0], 0) + assert third_node.is_contained_in(polytope) + third_node = QuadNode([0.1, 0.5], 0) + assert not third_node.is_contained_in(polytope) + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract_single(self.API.datacube, polytope) + # tree = self.API.slice(self.API.datacube, [polytope]) + # assert len(tree.leaves) == 3 + assert len(tree) == 3 + # tree.pprint() + points = [[10, 10], [80, 10], [-5, 5], [5, 50], [5, 10], [50, 10], [2, 10], [15, 15]] + polytope = ConvexPolytope(["latitude", "longitude"], [[-10, 1], [20, 1], [5, 20]]) + tree = self.API.engines["quadtree"].extract_single(self.API.datacube, polytope) + # tree = self.API.slice(self.API.datacube, [polytope]) + # assert len(tree.leaves) == 4 + assert len(tree) == 4 + # tree.pprint() + + # @pytest.mark.skip("performance test") + @pytest.mark.fdb + def test_large_scale_extraction(self): + import time + + import numpy as np + + x = np.linspace(0, 100, 1000) + y = np.linspace(0, 100, 1000) + # create the mesh based on these arrays + X, Y = np.meshgrid(x, y) + X = X.reshape((np.prod(X.shape),)) + Y = Y.reshape((np.prod(Y.shape),)) + coords = zip(X, Y) + points = [list(coord) for coord in coords] + time0 = time.time() + polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + print(time.time() - time0) + time1 = time.time() + tree = self.API.engines["quadtree"].extract_single(self.API.datacube, polytope) + print(time.time() - time1) # = 5.919436931610107 + # print(len(tree.leaves)) # = 55100 + print(len(tree)) # = 55100 + # NOTE: maybe for 2D qhull here, scipy is not the fastest + # but use shewchuk's triangle algo: https://www.cs.cmu.edu/~quake/triangle.html? diff --git a/tests/test_quadtree_edge_cases.py b/tests/test_quadtree_edge_cases.py new file mode 100644 index 000000000..8d78f5b60 --- /dev/null +++ b/tests/test_quadtree_edge_cases.py @@ -0,0 +1,56 @@ +import pytest + +from polytope_feature.polytope import Polytope +from polytope_feature.shapes import Box + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "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": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10], [0.035149384216, 0.0]] + polytope = Box(["latitude", "longitude"], [0, 0], [15, 15]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract_single(self.API.datacube, polytope) + # tree.pprint() + assert len(tree) == 3 + assert set([leaf.index for leaf in tree]) == set([0, 4, 6]) diff --git a/tests/test_quadtree_indices.py b/tests/test_quadtree_indices.py new file mode 100644 index 000000000..ff1ebad33 --- /dev/null +++ b/tests/test_quadtree_indices.py @@ -0,0 +1,81 @@ +import pytest + +from polytope_feature.polytope import Polytope +from polytope_feature.shapes import Box, ConvexPolytope + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "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": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer_extract(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + polytope = Box(["latitude", "longitude"], [1, 1], [20, 30]).polytope()[0] + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract_single(self.API.datacube, polytope) + assert len(tree) == 3 + assert set([leaf.index for leaf in tree]) == set( + [ + 0, + 3, + 4, + ] + ) + # tree.pprint() + points = [[10, 10], [80, 10], [-5, 5], [5, 50], [5, 10], [50, 10], [2, 10], [15, 15]] + polytope = ConvexPolytope(["latitude", "longitude"], [[-10, 1], [20, 1], [5, 20]]) + self.API = Polytope( + datacube=self.fdbdatacube, + options=self.options, + engine_options={"latitude": "quadtree", "longitude": "quadtree"}, + point_cloud_options=points, + ) + tree = self.API.engines["quadtree"].extract_single(self.API.datacube, polytope) + assert len(tree) == 4 + assert set([leaf.index for leaf in tree]) == set( + [ + 0, + 2, + 4, + 6, + ] + ) + # tree.pprint() diff --git a/tests/test_quadtree_optimisation.py b/tests/test_quadtree_optimisation.py new file mode 100644 index 000000000..c1662a621 --- /dev/null +++ b/tests/test_quadtree_optimisation.py @@ -0,0 +1,61 @@ +import pytest + +from polytope_feature.engine.quadtree_slicer import QuadTreeSlicer +from polytope_feature.shapes import Box + + +class TestQuadTreeSlicer: + def setup_method(self, method): + import pygribjump as gj + + self.options = { + "axis_config": [ + {"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}, + {"axis_name": "number", "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": "irregular", "resolution": 1280, "axes": ["latitude", "longitude"]} + ], + }, + ], + "compressed_axes_config": [ + "longitude", + "latitude", + "levtype", + "step", + "date", + "domain", + "expver", + "param", + "class", + "stream", + "type", + ], + "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, + } + self.fdbdatacube = gj.GribJump() + + @pytest.mark.fdb + def test_quad_tree_slicer(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + slicer.quad_tree.pprint() + pass + + @pytest.mark.fdb + def test_quad_tree_query_polygon(self): + points = [[10, 10], [80, 10], [-5, 5], [5, 20], [5, 10], [50, 10]] + slicer = QuadTreeSlicer(points) + polytope = Box(["lat", "lon"], [0, 0], [90, 45]).polytope()[0] + results = slicer.quad_tree.query_polygon(polytope) + assert len(results) == 5 + assert (10, 10) in [node.item for node in results] + assert (5, 10) in [node.item for node in results] + assert (5, 20) in [node.item for node in results] + assert (80, 10) in [node.item for node in results] + assert (50, 10) in [node.item for node in results] diff --git a/tests/test_regular_grid.py b/tests/test_regular_grid.py index 7eaaeccb5..bea22f265 100644 --- a/tests/test_regular_grid.py +++ b/tests/test_regular_grid.py @@ -2,7 +2,6 @@ import pytest from helper_functions import download_test_data, find_nearest_latlon -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Disk, Select @@ -76,10 +75,8 @@ def test_regular_grid(self): Select("number", ["0", "1"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_request_tree.py b/tests/test_request_tree.py index 038599724..d0976b941 100644 --- a/tests/test_request_tree.py +++ b/tests/test_request_tree.py @@ -274,7 +274,8 @@ def test_get_ancestors(self): child1.add_child(grandchild1) root_node1 = TensorIndexTree() root_node1.add_child(child1) - assert greatgrandchild1.get_ancestors() == SortedList([greatgrandchild1, grandchild1, child1]) + # assert greatgrandchild1.get_ancestors() == SortedList([greatgrandchild1, grandchild1, child1]) + assert greatgrandchild1.get_ancestors() == SortedList([child1, grandchild1, greatgrandchild1]) def test_add_or_get_child(self): axis1 = IntDatacubeAxis() @@ -290,7 +291,7 @@ def test_add_or_get_child(self): root_node = TensorIndexTree() root_node.add_child(child1) assert root_node.create_child(axis1, 0, [])[0] == child1 - assert root_node.create_child(axis2, (), [])[0].parent == root_node + assert root_node.create_child(axis2, 0, [])[0].parent == root_node def test_eq(self): axis1 = IntDatacubeAxis() diff --git a/tests/test_request_trees_after_slicing.py b/tests/test_request_trees_after_slicing.py index 2aa28b4a7..af13ab352 100644 --- a/tests/test_request_trees_after_slicing.py +++ b/tests/test_request_trees_after_slicing.py @@ -3,7 +3,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.datacube_axis import IntDatacubeAxis -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box @@ -19,14 +18,13 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_path_values(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) datacube_path = request.leaves[0].flatten() request.pprint() assert datacube_path.values() == tuple([tuple([3.0]), tuple([1.0, 2, 3])]) @@ -35,7 +33,7 @@ def test_path_values(self): def test_path_keys(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) datacube_path = request.leaves[0].flatten() assert datacube_path.keys()[0] == "step" assert datacube_path.keys()[1] == "level" @@ -43,14 +41,14 @@ def test_path_keys(self): def test_path_pprint(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) datacube_path = request.leaves[0].flatten() datacube_path.pprint() def test_flatten(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) path = request.leaves[0].flatten() request.pprint() assert path["step"] == tuple([3.0]) @@ -59,7 +57,7 @@ def test_flatten(self): def test_add_child(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) request1 = request.leaves[0] request2 = request.leaves[0] # Test adding child @@ -77,13 +75,13 @@ def test_add_child(self): def test_pprint(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) request.pprint() def test_remove_branch(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) prev_request_size = len(request.leaves) request.leaves[0].remove_branch() new_request_size = len(request.leaves) diff --git a/tests/test_reverse_transformation.py b/tests/test_reverse_transformation.py index f1d9b24ed..1b752ae3a 100644 --- a/tests/test_reverse_transformation.py +++ b/tests/test_reverse_transformation.py @@ -21,7 +21,7 @@ def setup_method(self, method): "compressed_axes_config": ["lat"], } self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_reverse_transformation(self): request = Request(Select("lat", [1, 2, 3])) diff --git a/tests/test_shapes.py b/tests/test_shapes.py index ac730905a..b6408ccc8 100644 --- a/tests/test_shapes.py +++ b/tests/test_shapes.py @@ -3,7 +3,6 @@ import pytest import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import All, Select, Span @@ -25,10 +24,8 @@ def setup_method(self, method): "axis_config": [{"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}], "compressed_axes_config": ["date", "step", "level", "longitude"], } - self.slicer = HullSlicer() self.API = Polytope( datacube=array, - engine=self.slicer, options=self.options, ) @@ -76,7 +73,6 @@ def test_all_mapper_cyclic(self): "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, } self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() request = Request( Select("step", [11]), @@ -91,7 +87,7 @@ def test_all_mapper_cyclic(self): Span("latitude", 89.9, 90), All("longitude"), ) - self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, options=self.options) + self.API = Polytope(datacube=self.fdbdatacube, options=self.options) result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 1 diff --git a/tests/test_slice_date_range_fdb.py b/tests/test_slice_date_range_fdb.py index dae81cb22..ed90b903c 100644 --- a/tests/test_slice_date_range_fdb.py +++ b/tests/test_slice_date_range_fdb.py @@ -1,7 +1,6 @@ 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 @@ -61,10 +60,8 @@ def test_fdb_datacube(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -91,10 +88,8 @@ def test_fdb_datacube_select_non_existing_last(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -121,10 +116,8 @@ def test_fdb_datacube_select_non_existing_first(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -134,37 +127,37 @@ def test_fdb_datacube_select_non_existing_first(self): assert len(result.leaves[i].result) == 3 @pytest.mark.fdb - def test_fdb_datacube_select_completely_non_existing(self): + def test_fdb_datacube_disk(self): import pygribjump as gj request = Request( Select("step", [0]), Select("levtype", ["sfc"]), - Select("date", [pd.Timestamp("20230624T120000"), pd.Timestamp("20230626T120000")]), + 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"]), - Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), + Disk(["latitude", "longitude"], [0, 0], [0.1, 0.1]), ) 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() - assert len(result.leaves) == 1 - for i in range(len(result.leaves)): - assert len(result.leaves[i].result) == 0 + assert len(result.leaves) == 2 + assert len(result.leaves[0].result) == 3 + assert len(result.leaves[1].result) == 3 + assert len(result.leaves[0].values) == 3 + assert len(result.leaves[1].values) == 3 @pytest.mark.fdb - def test_fdb_datacube_disk(self): + def test_fdb_datacube_disk_2(self): import pygribjump as gj request = Request( @@ -177,54 +170,48 @@ def test_fdb_datacube_disk(self): Select("class", ["od"]), Select("stream", ["oper"]), Select("type", ["an"]), - Disk(["latitude", "longitude"], [0, 0], [0.1, 0.1]), + Disk(["latitude", "longitude"], [0.05, 0.070148090413], [0.1, 0.15]), ) 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() - assert len(result.leaves) == 2 + assert len(result.leaves) == 3 assert len(result.leaves[0].result) == 3 - assert len(result.leaves[1].result) == 3 + assert len(result.leaves[1].result) == 5 + assert len(result.leaves[2].result) == 3 assert len(result.leaves[0].values) == 3 - assert len(result.leaves[1].values) == 3 + assert len(result.leaves[1].values) == 5 + assert len(result.leaves[2].values) == 3 @pytest.mark.fdb - def test_fdb_datacube_disk_2(self): + def test_fdb_datacube_select_completely_non_existing(self): import pygribjump as gj request = Request( Select("step", [0]), Select("levtype", ["sfc"]), - Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")), + Select("date", [pd.Timestamp("20230624T120000"), pd.Timestamp("20230626T120000")]), Select("domain", ["g"]), Select("expver", ["0001"]), Select("param", ["167"]), Select("class", ["od"]), Select("stream", ["oper"]), Select("type", ["an"]), - Disk(["latitude", "longitude"], [0.05, 0.070148090413], [0.1, 0.15]), + Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]), ) 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() - assert len(result.leaves) == 3 - assert len(result.leaves[0].result) == 3 - assert len(result.leaves[1].result) == 5 - assert len(result.leaves[2].result) == 3 - assert len(result.leaves[0].values) == 3 - assert len(result.leaves[1].values) == 5 - assert len(result.leaves[2].values) == 3 + assert len(result.leaves) == 1 + for i in range(len(result.leaves)): + assert len(result.leaves[i].result) == 0 diff --git a/tests/test_slice_date_range_fdb_v2.py b/tests/test_slice_date_range_fdb_v2.py index fcacf6b7c..cb76be722 100644 --- a/tests/test_slice_date_range_fdb_v2.py +++ b/tests/test_slice_date_range_fdb_v2.py @@ -2,7 +2,6 @@ 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 Select, Span @@ -67,10 +66,8 @@ def test_fdb_datacube(self): Select("number", ["0"]), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_slicer_engine.py b/tests/test_slicer_engine.py index 289eeab6b..ba04b2672 100644 --- a/tests/test_slicer_engine.py +++ b/tests/test_slicer_engine.py @@ -3,7 +3,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope from polytope_feature.shapes import Box @@ -19,14 +18,13 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_extract(self): box = Box(["step", "level"], [3.0, 1.0], [6.0, 3.0]) polytope = box.polytope() - request = self.slicer.extract(self.xarraydatacube, polytope) + request = self.API.slice(self.xarraydatacube, polytope) assert request.axis == TensorIndexTree.root assert request.parent is None assert request.values is tuple() diff --git a/tests/test_slicer_era5.py b/tests/test_slicer_era5.py index 9f11bcef8..c3e450c07 100644 --- a/tests/test_slicer_era5.py +++ b/tests/test_slicer_era5.py @@ -3,7 +3,6 @@ from earthkit import data from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -15,14 +14,12 @@ def setup_method(self, method): ds = data.from_source("file", "./tests/data/era5-levels-members.grib") array = ds.to_xarray(engine="cfgrib").isel(step=0).t - self.slicer = HullSlicer() options = { "axis_config": [{"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}], "compressed_axes_config": ["number", "time", "latitude", "longitude", "step", "isobaricInhPa"], } self.API = Polytope( datacube=array, - engine=self.slicer, options=options, ) diff --git a/tests/test_slicer_xarray.py b/tests/test_slicer_xarray.py index b36cd051c..3aa8d3f7c 100644 --- a/tests/test_slicer_xarray.py +++ b/tests/test_slicer_xarray.py @@ -2,7 +2,6 @@ import pandas as pd import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select, Span @@ -19,9 +18,8 @@ def setup_method(self, method): "level": range(1, 130), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "step", "level"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_2D_box(self): request = Request(Box(["step", "level"], [3, 10], [6, 11]), Select("date", ["2000-01-01"])) diff --git a/tests/test_slicing_unsliceable_axis.py b/tests/test_slicing_unsliceable_axis.py index 286252065..f5a73aab0 100644 --- a/tests/test_slicing_unsliceable_axis.py +++ b/tests/test_slicing_unsliceable_axis.py @@ -3,7 +3,6 @@ import pytest import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select from polytope_feature.utility.exceptions import UnsliceableShapeError @@ -17,9 +16,8 @@ def setup_method(self, method): dims=("date", "variable", "level"), coords={"date": pd.date_range("2000-01-01", "2000-01-03", 3), "variable": ["a"], "level": range(1, 130)}, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "variable", "level"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_slicing_xarray_3D.py b/tests/test_slicing_xarray_3D.py index 72f68d0ff..b1740810e 100644 --- a/tests/test_slicing_xarray_3D.py +++ b/tests/test_slicing_xarray_3D.py @@ -7,7 +7,6 @@ from polytope_feature.datacube.backends.xarray import XArrayDatacube from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import ( Box, @@ -34,9 +33,8 @@ def setup_method(self, method): }, ) self.xarraydatacube = XArrayDatacube(array) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "step", "level"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_slicing_xarray_4D.py b/tests/test_slicing_xarray_4D.py index 9345e8db5..67db9292d 100644 --- a/tests/test_slicing_xarray_4D.py +++ b/tests/test_slicing_xarray_4D.py @@ -4,7 +4,6 @@ import xarray as xr from polytope_feature.datacube.tensor_index_tree import TensorIndexTree -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import ( Box, @@ -36,9 +35,8 @@ def setup_method(self, method): "lat": np.around(np.arange(0.0, 10.0, 0.1), 15), }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["date", "step", "level", "lat"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_snapping.py b/tests/test_snapping.py index 529f1d69a..240ef9729 100644 --- a/tests/test_snapping.py +++ b/tests/test_snapping.py @@ -1,7 +1,6 @@ import numpy as np import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -17,9 +16,8 @@ def setup_method(self, method): "step": [1, 3, 5], }, ) - self.slicer = HullSlicer() options = {"compressed_axes_config": ["level", "step"]} - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) # Testing different shapes diff --git a/tests/test_snapping_real_data.py b/tests/test_snapping_real_data.py index 9b5950742..c2431d5c0 100644 --- a/tests/test_snapping_real_data.py +++ b/tests/test_snapping_real_data.py @@ -6,7 +6,6 @@ from earthkit import data from helper_functions import download_test_data -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -17,20 +16,14 @@ def setup_method(self, method): download_test_data(nexus_url, "era5-levels-members.grib") ds = data.from_source("file", "./tests/data/era5-levels-members.grib") - array = ds.to_xarray(engine="cfgrib").isel(step=0).t - self.slicer = HullSlicer() - options = { + self.array = ds.to_xarray(engine="cfgrib").isel(step=0).t + self.options = { "axis_config": [ {"axis_name": "latitude", "transformations": [{"name": "reverse", "is_reverse": True}]}, {"axis_name": "longitude", "transformations": [{"name": "cyclic", "range": [0, 360]}]}, ], "compressed_axes_config": ["longitude", "latitude", "step", "time", "number", "isobaricInhPa"], } - self.API = Polytope( - datacube=array, - engine=self.slicer, - options=options, - ) @pytest.mark.internet def test_surrounding_on_grid_point(self): @@ -43,6 +36,11 @@ def test_surrounding_on_grid_point(self): Select("longitude", [requested_lon], method="surrounding"), Select("step", [np.timedelta64(0, "s")]), ) + + self.API = Polytope( + datacube=self.array, + options=self.options, + ) result = self.API.retrieve(request) result.pprint() country_points_plotting = [] diff --git a/tests/test_tree_protobuf_encoding.py b/tests/test_tree_protobuf_encoding.py index b177cbe2b..710c01694 100644 --- a/tests/test_tree_protobuf_encoding.py +++ b/tests/test_tree_protobuf_encoding.py @@ -49,17 +49,14 @@ def setup_method(self): def test_encoding(self): import pygribjump as gj - from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope self.options = { "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, } self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) fdb_datacube = self.API.datacube diff --git a/tests/test_tree_protobuf_encoding_fdb.py b/tests/test_tree_protobuf_encoding_fdb.py index 2e288f0d6..89a9075a2 100644 --- a/tests/test_tree_protobuf_encoding_fdb.py +++ b/tests/test_tree_protobuf_encoding_fdb.py @@ -12,7 +12,6 @@ def setup_method(self): def test_encoding(self): import pygribjump as gj - from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Box, Select @@ -61,10 +60,8 @@ def test_encoding(self): "pre_path": {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"}, } self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_type_change_transformation.py b/tests/test_type_change_transformation.py index 09f800610..2dd033f0d 100644 --- a/tests/test_type_change_transformation.py +++ b/tests/test_type_change_transformation.py @@ -1,7 +1,6 @@ import numpy as np import xarray as xr -from polytope_feature.engine.hullslicer import HullSlicer from polytope_feature.polytope import Polytope, Request from polytope_feature.shapes import Select @@ -21,8 +20,7 @@ def setup_method(self, method): "axis_config": [{"axis_name": "step", "transformations": [{"name": "type_change", "type": "int"}]}], "compressed_axes_config": ["step"], } - self.slicer = HullSlicer() - self.API = Polytope(datacube=array, engine=self.slicer, options=options) + self.API = Polytope(datacube=array, options=options) def test_merge_axis(self): request = Request(Select("step", [0])) diff --git a/tests/test_union_gj.py b/tests/test_union_gj.py index 9182fc360..13a99552d 100644 --- a/tests/test_union_gj.py +++ b/tests/test_union_gj.py @@ -1,7 +1,6 @@ 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, Select, Span, Union @@ -67,10 +66,8 @@ def test_fdb_datacube(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -108,10 +105,8 @@ def test_fdb_datacube_complete_overlap(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) @@ -149,10 +144,8 @@ def test_fdb_datacube_complete_overlap_v2(self): ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( datacube=self.fdbdatacube, - engine=self.slicer, options=self.options, ) result = self.API.retrieve(request) diff --git a/tests/test_union_point_box.py b/tests/test_union_point_box.py index ba280bc02..9132f14dc 100644 --- a/tests/test_union_point_box.py +++ b/tests/test_union_point_box.py @@ -1,7 +1,6 @@ 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, Point, Select, Union @@ -64,9 +63,7 @@ def test_fdb_datacube(self): Union(["latitude", "longitude"], box, point), ) self.fdbdatacube = gj.GribJump() - self.slicer = HullSlicer() self.API = Polytope( - request=request, datacube=self.fdbdatacube, engine=self.slicer, options=self.options,