From cb0ae833b31f9d1a3946d52356c856c18c3a1d49 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Wed, 22 Nov 2023 10:46:12 +0100 Subject: [PATCH 1/5] add repr for shapes and Point shape which supports surrounding method --- polytope/polytope.py | 6 +++ polytope/shapes.py | 66 +++++++++++++++++++++++++++++ tests/data/era5-levels-members.grib | 3 ++ tests/data/foo.grib | 3 ++ tests/data/healpix.grib | 3 ++ tests/test_point_shape.py | 36 ++++++++++++++++ 6 files changed, 117 insertions(+) create mode 100644 tests/data/era5-levels-members.grib create mode 100644 tests/data/foo.grib create mode 100644 tests/data/healpix.grib create mode 100644 tests/test_point_shape.py diff --git a/polytope/polytope.py b/polytope/polytope.py index 5e8fc7133..f6d4a723e 100644 --- a/polytope/polytope.py +++ b/polytope/polytope.py @@ -29,6 +29,12 @@ def polytopes(self): polytopes.extend(shape.polytope()) return polytopes + def __repr__(self): + return_str = "" + for shape in self.shapes: + return_str += shape.__repr__() + "\n" + return return_str + class Polytope: def __init__(self, datacube, engine=None, axis_options={}): diff --git a/polytope/shapes.py b/polytope/shapes.py index 94562abc5..5a824533a 100644 --- a/polytope/shapes.py +++ b/polytope/shapes.py @@ -38,6 +38,9 @@ def extents(self, axis): def __str__(self): return f"Polytope in {self.axes} with points {self.points}" + def __repr__(self): + return f"Polytope in {self.axes} with points {self.points}" + def axes(self): return self._axes @@ -60,6 +63,31 @@ def axes(self): def polytope(self): return [ConvexPolytope([self.axis], [[v]], self.method) for v in self.values] + def __repr__(self): + return f"Select in {self.axis} with points {self.values}" + + +class Point(Shape): + """Matches several discrete value""" + + def __init__(self, axes, values, method=None): + self._axes = axes + self.values = values + self.method = method + self.polytopes = [] + for i in range(len(axes)): + polytope_points = [v[i] for v in self.values] + self.polytopes.append(ConvexPolytope([axes[i]], [polytope_points], method)) + + def axes(self): + return self._axes + + def polytope(self): + return self.polytopes + + def __repr__(self): + return f"Point in {self._axes} with points {self.values}" + class Span(Shape): """1-D range along a single axis""" @@ -77,6 +105,9 @@ def axes(self): def polytope(self): return [ConvexPolytope([self.axis], [[self.lower], [self.upper]])] + def __repr__(self): + return f"Span in {self.axis} with range from {self.lower} to {self.upper}" + class All(Span): """Matches all indices in an axis""" @@ -84,12 +115,17 @@ class All(Span): def __init__(self, axis): super().__init__(axis) + def __repr__(self): + return f"All in {self.axis}" + class Box(Shape): """N-D axis-aligned bounding box (AABB), specified by two opposite corners""" def __init__(self, axes, lower_corner=None, upper_corner=None): dimension = len(axes) + self._lower_corner = lower_corner + self._upper_corner = upper_corner self._axes = axes assert len(lower_corner) == dimension assert len(upper_corner) == dimension @@ -121,6 +157,9 @@ def axes(self): def polytope(self): return [ConvexPolytope(self.axes(), self.vertices)] + def __repr__(self): + return f"Box in {self._axes} with with lower corner {self._lower_corner} and upper corner{self._upper_corner}" + class Disk(Shape): """2-D shape bounded by an ellipse""" @@ -159,6 +198,9 @@ def axes(self): def polytope(self): return [ConvexPolytope(self.axes(), self.points)] + def __repr__(self): + return f"Disk in {self._axes} with centred at {self.centre} and with radius {self.radius}" + class Ellipsoid(Shape): # Here we use the formula for the inscribed circle in an icosahedron @@ -211,12 +253,18 @@ def _icosahedron_edge_length_coeff(self): def polytope(self): return [ConvexPolytope(self.axes(), self.points)] + def __repr__(self): + return f"Ellipsoid in {self._axes} with centred at {self.centre} and with radius {self.radius}" + class PathSegment(Shape): """N-D polytope defined by a shape which is swept along a straight line between two points""" def __init__(self, axes, shape: Shape, start: List, end: List): self._axes = axes + self._start = start + self._end = end + self._shape = shape assert shape.axes() == self.axes() assert len(start) == len(self.axes()) @@ -240,12 +288,18 @@ def axes(self): def polytope(self): return self.polytopes + def __repr__(self): + return f"PathSegment in {self._axes} obtained by sweeping a {self._shape.__repr__()} \ + between the points {self._start} and {self._end}" + class Path(Shape): """N-D polytope defined by a shape which is swept along a polyline defined by multiple points""" def __init__(self, axes, shape, *points, closed=False): self._axes = axes + self._shape = shape + self._points = points assert shape.axes() == self.axes() for p in points: @@ -267,6 +321,10 @@ def axes(self): def polytope(self): return self.union.polytope() + def __repr__(self): + return f"Path in {self._axes} obtained by sweeping a {self._shape.__repr__()} \ + between the points {self._points}" + class Union(Shape): """N-D union of two shapes with the same axes""" @@ -277,6 +335,7 @@ def __init__(self, axes, *shapes): assert s.axes() == self.axes() self.polytopes = [] + self._shapes = shapes for s in shapes: self.polytopes.extend(s.polytope()) @@ -287,6 +346,9 @@ def axes(self): def polytope(self): return self.polytopes + def __repr__(self): + return f"Union in {self._axes} of the shapes {self._shapes}" + class Polygon(Shape): """2-D polygon defined by a set of exterior points""" @@ -297,6 +359,7 @@ def __init__(self, axes, points): for p in points: assert len(p) == 2 + self._points = points triangles = tripy.earclip(points) self.polytopes = [] @@ -313,3 +376,6 @@ def axes(self): def polytope(self): return self.polytopes + + def __repr__(self): + return f"Polygon in {self._axes} with points {self._points}" diff --git a/tests/data/era5-levels-members.grib b/tests/data/era5-levels-members.grib new file mode 100644 index 000000000..90d45deed --- /dev/null +++ b/tests/data/era5-levels-members.grib @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c144fde61ca5d53702bf6f212775ef2cc783bdd66b6865160bf597c1b35ed898 +size 2361600 diff --git a/tests/data/foo.grib b/tests/data/foo.grib new file mode 100644 index 000000000..9c5efa68b --- /dev/null +++ b/tests/data/foo.grib @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:12533bae9b3aa078a3298a82f289c7dbba6dbdb9a21e4e5e24e84b8695a75dee +size 26409360 diff --git a/tests/data/healpix.grib b/tests/data/healpix.grib new file mode 100644 index 000000000..693c98c12 --- /dev/null +++ b/tests/data/healpix.grib @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:df57790df280627b2883c01fa4b56986d9938ed69b7434fe09c4bede99d2895f +size 37030 diff --git a/tests/test_point_shape.py b/tests/test_point_shape.py new file mode 100644 index 000000000..0bc203d61 --- /dev/null +++ b/tests/test_point_shape.py @@ -0,0 +1,36 @@ +import numpy as np +import pandas as pd +import xarray as xr + +from polytope.datacube.backends.xarray import XArrayDatacube +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Point, Select + + +class TestSlicing3DXarrayDatacube: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + array = xr.DataArray( + np.random.randn(3, 6, 129), + dims=("date", "step", "level"), + coords={ + "date": pd.date_range("2000-01-01", "2000-01-03", 3), + "step": [0, 3, 6, 9, 12, 15], + "level": range(1, 130), + }, + ) + self.xarraydatacube = XArrayDatacube(array) + self.slicer = HullSlicer() + self.API = Polytope(datacube=array, engine=self.slicer) + + def test_point(self): + request = Request(Point(["step", "level"], [[3, 10]]), Select("date", ["2000-01-01"])) + result = self.API.retrieve(request) + assert len(result.leaves) == 1 + assert result.leaves[0].axis.name == "level" + + def test_point_surrounding_step(self): + request = Request(Point(["step", "level"], [[2, 10]], method="surrounding"), Select("date", ["2000-01-01"])) + result = self.API.retrieve(request) + assert len(result.leaves) == 6 From df874d4b6ce7e23b6e0908ff0e5ce7394c5f564c Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Wed, 22 Nov 2023 10:50:34 +0100 Subject: [PATCH 2/5] add grib to gitignore --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 037872e66..7b8e19f62 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,5 @@ polytope.egg-info .pytest_cache *.prof -*.idx \ No newline at end of file +*.idx +*.grib \ No newline at end of file From d6c6c4a51454e151bed5efb489fdd1ac2ef7ea32 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Wed, 22 Nov 2023 10:51:18 +0100 Subject: [PATCH 3/5] remove data from tests --- tests/data/era5-levels-members.grib | 3 --- tests/data/foo.grib | 3 --- tests/data/healpix.grib | 3 --- 3 files changed, 9 deletions(-) delete mode 100644 tests/data/era5-levels-members.grib delete mode 100644 tests/data/foo.grib delete mode 100644 tests/data/healpix.grib diff --git a/tests/data/era5-levels-members.grib b/tests/data/era5-levels-members.grib deleted file mode 100644 index 90d45deed..000000000 --- a/tests/data/era5-levels-members.grib +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:c144fde61ca5d53702bf6f212775ef2cc783bdd66b6865160bf597c1b35ed898 -size 2361600 diff --git a/tests/data/foo.grib b/tests/data/foo.grib deleted file mode 100644 index 9c5efa68b..000000000 --- a/tests/data/foo.grib +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:12533bae9b3aa078a3298a82f289c7dbba6dbdb9a21e4e5e24e84b8695a75dee -size 26409360 diff --git a/tests/data/healpix.grib b/tests/data/healpix.grib deleted file mode 100644 index 693c98c12..000000000 --- a/tests/data/healpix.grib +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:df57790df280627b2883c01fa4b56986d9938ed69b7434fe09c4bede99d2895f -size 37030 From 94c61b6d79b6fa05c2426a867b284cd774f664de Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Wed, 22 Nov 2023 10:56:53 +0100 Subject: [PATCH 4/5] remove double print for convex polytope --- polytope/shapes.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/polytope/shapes.py b/polytope/shapes.py index 5a824533a..306d1a64a 100644 --- a/polytope/shapes.py +++ b/polytope/shapes.py @@ -38,9 +38,6 @@ def extents(self, axis): def __str__(self): return f"Polytope in {self.axes} with points {self.points}" - def __repr__(self): - return f"Polytope in {self.axes} with points {self.points}" - def axes(self): return self._axes From 3d6b37cc5ec206872857a03d946c543304966c71 Mon Sep 17 00:00:00 2001 From: Mathilde Leuridan Date: Wed, 22 Nov 2023 14:49:41 +0100 Subject: [PATCH 5/5] make All shape work --- polytope/datacube/datacube_axis.py | 5 ++ polytope/shapes.py | 2 +- tests/test_shapes.py | 74 ++++++++++++++++++++++++++++++ 3 files changed, 80 insertions(+), 1 deletion(-) create mode 100644 tests/test_shapes.py diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 77a223710..625051f7c 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -1,3 +1,4 @@ +import math from abc import ABC, abstractmethod from copy import deepcopy from typing import Any, List @@ -18,6 +19,10 @@ def update_range(): def to_intervals(range): update_range() + if range[0] == -math.inf: + range[0] = cls.range[0] + if range[1] == math.inf: + range[1] = cls.range[1] axis_lower = cls.range[0] axis_upper = cls.range[1] axis_range = axis_upper - axis_lower diff --git a/polytope/shapes.py b/polytope/shapes.py index 306d1a64a..cc77273a7 100644 --- a/polytope/shapes.py +++ b/polytope/shapes.py @@ -89,7 +89,7 @@ def __repr__(self): class Span(Shape): """1-D range along a single axis""" - def __init__(self, axis, lower=None, upper=None): + def __init__(self, axis, lower=-math.inf, upper=math.inf): assert not isinstance(lower, list) assert not isinstance(upper, list) self.axis = axis diff --git a/tests/test_shapes.py b/tests/test_shapes.py new file mode 100644 index 000000000..ba355dd4c --- /dev/null +++ b/tests/test_shapes.py @@ -0,0 +1,74 @@ +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +from polytope.datacube.backends.FDB_datacube import FDBDatacube +from polytope.datacube.backends.xarray import XArrayDatacube +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import All, Select, Span + + +class TestSlicing3DXarrayDatacube: + def setup_method(self, method): + # Create a dataarray with 3 labelled axes using different index types + array = xr.DataArray( + np.random.randn(3, 6, 129, 360), + dims=("date", "step", "level", "longitude"), + coords={ + "date": pd.date_range("2000-01-01", "2000-01-03", 3), + "step": [0, 3, 6, 9, 12, 15], + "level": range(1, 130), + "longitude": range(0, 360), + }, + ) + self.xarraydatacube = XArrayDatacube(array) + self.options = {"longitude": {"transformation": {"cyclic": [0, 360]}}} + self.slicer = HullSlicer() + self.API = Polytope(datacube=array, engine=self.slicer, axis_options=self.options) + + def test_all(self): + request = Request(Select("step", [3]), Select("date", ["2000-01-01"]), All("level"), Select("longitude", [1])) + result = self.API.retrieve(request) + assert len(result.leaves) == 129 + + def test_all_cyclic(self): + request = Request(Select("step", [3]), Select("date", ["2000-01-01"]), Select("level", [1]), All("longitude")) + result = self.API.retrieve(request) + # result.pprint() + assert len(result.leaves) == 360 + + @pytest.mark.skip(reason="can't install fdb branch on CI") + def test_all_mapper_cyclic(self): + self.options = { + "values": { + "transformation": { + "mapper": {"type": "octahedral", "resolution": 1280, "axes": ["latitude", "longitude"]} + } + }, + "date": {"transformation": {"merge": {"with": "time", "linkers": ["T", "00"]}}}, + "step": {"transformation": {"type_change": "int"}}, + "longitude": {"transformation": {"cyclic": [0, 360]}}, + } + self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "step": 11} + self.fdbdatacube = FDBDatacube(self.config, axis_options=self.options) + self.slicer = HullSlicer() + self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, axis_options=self.options) + + request = Request( + Select("step", [11]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20230710T120000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["151130"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Span("latitude", 89.9, 90), + All("longitude"), + ) + result = self.API.retrieve(request) + # result.pprint() + assert len(result.leaves) == 20