From 3bdf9b0ab8dfbafb52285873410bf743e4df0102 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Fri, 26 Jul 2024 16:30:13 +0200 Subject: [PATCH 01/15] index and accessor implementations for `parents` --- xdggs/accessor.py | 15 +++++++++++++++ xdggs/index.py | 3 +++ 2 files changed, 18 insertions(+) diff --git a/xdggs/accessor.py b/xdggs/accessor.py index f8b67cc4..6ebd9fa7 100644 --- a/xdggs/accessor.py +++ b/xdggs/accessor.py @@ -104,3 +104,18 @@ def cell_centers(self): "longitude": (self.index._dim, lon_data), } ) + + def parents(self, resolution: int) -> xr.DataArray: + """determine the parent cell ids of the cells + + Parameters + ---------- + resolution : int + The parent resolution. Must be smaller than the current resolution. + + Returns + ------- + parents : DataArray + The parent cell ids, one for each input cell. + """ + return self.index.parents(resolution) diff --git a/xdggs/index.py b/xdggs/index.py index 2c3a7c2c..c6a94bcd 100644 --- a/xdggs/index.py +++ b/xdggs/index.py @@ -72,6 +72,9 @@ def _replace(self, new_pd_index: PandasIndex): def cell_centers(self) -> tuple[np.ndarray, np.ndarray]: return self._grid.cell_ids2geographic(self._pd_index.index.values) + def parents(self, resolution: int) -> np.ndarray: + return self._grid.parents(self._pd_index.index.values, resolution=resolution) + @property def grid_info(self) -> DGGSInfo: return self._grid From a387db3e59da16e273f9337062004c769061a9af Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Fri, 26 Jul 2024 16:50:50 +0200 Subject: [PATCH 02/15] return the parents as a `DataArray` --- xdggs/accessor.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/xdggs/accessor.py b/xdggs/accessor.py index 6ebd9fa7..97fbc12e 100644 --- a/xdggs/accessor.py +++ b/xdggs/accessor.py @@ -118,4 +118,9 @@ def parents(self, resolution: int) -> xr.DataArray: parents : DataArray The parent cell ids, one for each input cell. """ - return self.index.parents(resolution) + data = self.index.parents(resolution) + + params = self.grid_info.to_dict() + params["resolution"] = resolution + + return self.coord.copy(data=data).assign_attrs(**params).rename("parents") From cf066e55d917d7a69ffe37b831307f0a805f8359 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Fri, 26 Jul 2024 16:51:19 +0200 Subject: [PATCH 03/15] implement the healpix version of the `parents` computation --- xdggs/healpix.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/xdggs/healpix.py b/xdggs/healpix.py index bb76e326..b5e6579a 100644 --- a/xdggs/healpix.py +++ b/xdggs/healpix.py @@ -135,6 +135,19 @@ def cell_ids2geographic(self, cell_ids): def geographic2cell_ids(self, lon, lat): return healpy.ang2pix(self.nside, lon, lat, lonlat=True, nest=self.nest) + def parents(self, cell_ids, resolution): + if resolution >= self.resolution: + raise ValueError( + f"resolution is not a parent: {resolution} >= {self.resolution}" + ) + + offset = self.resolution - resolution + x, y, f = healpy.pix2xyf(self.nside, cell_ids, nest=self.nest) + x_ = x >> offset + y_ = y >> offset + + return healpy.xyf2pix(2**resolution, x_, y_, f, nest=self.nest) + @register_dggs("healpix") class HealpixIndex(DGGSIndex): From 6531f599b523ea27dc2bdf868f5cccb03025394b Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Fri, 26 Jul 2024 16:57:34 +0200 Subject: [PATCH 04/15] same for h3 --- xdggs/h3.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/xdggs/h3.py b/xdggs/h3.py index 17db4ed1..56266a99 100644 --- a/xdggs/h3.py +++ b/xdggs/h3.py @@ -9,6 +9,7 @@ import numpy as np import xarray as xr +from h3ronpy.arrow import change_resolution from h3ronpy.arrow.vector import cells_to_coordinates, coordinates_to_cells from xarray.indexes import PandasIndex @@ -45,6 +46,14 @@ def cell_ids2geographic( def geographic2cell_ids(self, lon, lat): return coordinates_to_cells(lat, lon, self.resolution, radians=False) + def parents(self, cell_ids, resolution): + if resolution >= self.resolution: + raise ValueError( + f"resolution is not a parent: {resolution} >= {self.resolution}" + ) + + return change_resolution(cell_ids, resolution) + @register_dggs("h3") class H3Index(DGGSIndex): From 49c0502592d3d90227b224eb96635756d390cd64 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Fri, 24 Jan 2025 14:23:26 +0100 Subject: [PATCH 05/15] switch to level --- xdggs/accessor.py | 10 +++++----- xdggs/healpix.py | 13 ++++++------- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/xdggs/accessor.py b/xdggs/accessor.py index e1ca83c3..9b7cb2f6 100644 --- a/xdggs/accessor.py +++ b/xdggs/accessor.py @@ -174,23 +174,23 @@ def cell_boundaries(self): boundaries, coords={self._name: self.cell_ids}, dims=self.cell_ids.dims ) - def parents(self, resolution: int) -> xr.DataArray: + def parents(self, level: int) -> xr.DataArray: """determine the parent cell ids of the cells Parameters ---------- - resolution : int - The parent resolution. Must be smaller than the current resolution. + level : int + The parent resolution level. Must be smaller than the current resolution. Returns ------- parents : DataArray The parent cell ids, one for each input cell. """ - data = self.index.parents(resolution) + data = self.index.parents(level) params = self.grid_info.to_dict() - params["resolution"] = resolution + params["resolution"] = level return self.coord.copy(data=data).assign_attrs(**params).rename("parents") diff --git a/xdggs/healpix.py b/xdggs/healpix.py index 36d4f99b..ef89e4ac 100644 --- a/xdggs/healpix.py +++ b/xdggs/healpix.py @@ -308,18 +308,17 @@ def cell_boundaries(self, cell_ids: Any, backend="shapely") -> np.ndarray: return backend_func(vertices) - def parents(self, cell_ids, resolution): - if resolution >= self.resolution: - raise ValueError( - f"resolution is not a parent: {resolution} >= {self.resolution}" - ) + def parents(self, cell_ids, level): + if level >= self.level: + raise ValueError(f"given level is not a parent: {level} >= {self.level}") - offset = self.resolution - resolution + # TODO: reimplement using `cdshealpix` + offset = self.level - level x, y, f = healpy.pix2xyf(self.nside, cell_ids, nest=self.nest) x_ = x >> offset y_ = y >> offset - return healpy.xyf2pix(2**resolution, x_, y_, f, nest=self.nest) + return healpy.xyf2pix(self.nside, x_, y_, f, nest=self.nest) @register_dggs("healpix") From 660ce8e97611489f5e7589523bad55ec8138ee49 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Tue, 17 Jun 2025 17:02:36 +0200 Subject: [PATCH 06/15] depend on `healpix-geo` --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 94f0f46d..bb384f28 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ requires-python = ">=3.10" dependencies = [ "xarray", "cdshealpix", + "healpix-geo>=0.0.3", "h3ronpy", "typing-extensions", "lonboard>=0.9.3", From fb9a593fba5656b07a2e824e2e612b95f8dfb043 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Tue, 17 Jun 2025 17:21:56 +0200 Subject: [PATCH 07/15] define `zoom_to` --- xdggs/accessor.py | 10 ++++++++++ xdggs/grid.py | 3 +++ xdggs/index.py | 4 ++-- 3 files changed, 15 insertions(+), 2 deletions(-) diff --git a/xdggs/accessor.py b/xdggs/accessor.py index 1852c981..34dbd2cb 100644 --- a/xdggs/accessor.py +++ b/xdggs/accessor.py @@ -174,6 +174,16 @@ def cell_boundaries(self): boundaries, coords={self._name: self.cell_ids}, dims=self.cell_ids.dims ) + def zoom_to(self, level): + zoomed = self.index.zoom_to(level=level) + + if zoomed.ndim == 1: + dims = self.cell_ids.dims + else: + dims = [*self.cell_ids.dims, "children"] + + return xr.DataArray(zoomed, coords={self._name: self.cell_ids}, dims=dims) + def explore(self, *, cmap="viridis", center=None, alpha=None, coords=None): """interactively explore the data using `lonboard` diff --git a/xdggs/grid.py b/xdggs/grid.py index eab18437..abac079e 100644 --- a/xdggs/grid.py +++ b/xdggs/grid.py @@ -49,6 +49,9 @@ def geographic2cell_ids(self, lon, lat): def cell_boundaries(self, cell_ids, backend="shapely"): raise NotImplementedError() + def zoom_to(self, cell_ids, level: int): + raise NotImplementedError() + def translate_parameters(mapping, translations): def translate(name, value): diff --git a/xdggs/index.py b/xdggs/index.py index 8a59113b..9f63d2de 100644 --- a/xdggs/index.py +++ b/xdggs/index.py @@ -95,8 +95,8 @@ def cell_centers(self) -> tuple[np.ndarray, np.ndarray]: def cell_boundaries(self) -> np.ndarray: return self.grid_info.cell_boundaries(self._pd_index.index.values) - def parents(self, resolution: int) -> np.ndarray: - return self._grid.parents(self._pd_index.index.values, resolution=resolution) + def zoom_to(self, level: int) -> np.ndarray: + return self._grid.zoom_to(self._pd_index.index.values, level=level) @property def grid_info(self) -> DGGSInfo: From 47ed097d1c976df98f075fe73c1e1e55add2cae3 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Tue, 17 Jun 2025 17:25:22 +0200 Subject: [PATCH 08/15] docstrings and typing --- xdggs/accessor.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/xdggs/accessor.py b/xdggs/accessor.py index 34dbd2cb..d0656c35 100644 --- a/xdggs/accessor.py +++ b/xdggs/accessor.py @@ -174,7 +174,22 @@ def cell_boundaries(self): boundaries, coords={self._name: self.cell_ids}, dims=self.cell_ids.dims ) - def zoom_to(self, level): + def zoom_to(self, level: int): + """Change the refinement level of the cell ids to `level`. + + Parameters + ---------- + level : int + The refinement level to change to. Can be smaller than the dataset's + level to compute parents, or bigger to fetch the children. In the + latter case, the array will have an additional `"children"` + dimension. + + Returns + ------- + zoomed : xr.DataArray + The children or parents of the current cells. + """ zoomed = self.index.zoom_to(level=level) if zoomed.ndim == 1: From f927149ee0e33de7164ebc0c2aa8755ceb331515 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Tue, 17 Jun 2025 17:25:32 +0200 Subject: [PATCH 09/15] implement `zoom_to` for healpix --- xdggs/healpix.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/xdggs/healpix.py b/xdggs/healpix.py index ef89e4ac..642fecb1 100644 --- a/xdggs/healpix.py +++ b/xdggs/healpix.py @@ -308,17 +308,16 @@ def cell_boundaries(self, cell_ids: Any, backend="shapely") -> np.ndarray: return backend_func(vertices) - def parents(self, cell_ids, level): - if level >= self.level: - raise ValueError(f"given level is not a parent: {level} >= {self.level}") + def zoom_to(self, cell_ids, level): + if self.indexing_scheme == "ring": + raise ValueError( + "Scaling does not make sense for the 'ring' scheme." + " Please convert to a nested scheme first." + ) - # TODO: reimplement using `cdshealpix` - offset = self.level - level - x, y, f = healpy.pix2xyf(self.nside, cell_ids, nest=self.nest) - x_ = x >> offset - y_ = y >> offset + from healpix_geo.nested import zoom_to - return healpy.xyf2pix(self.nside, x_, y_, f, nest=self.nest) + return zoom_to(cell_ids, self.level, level) @register_dggs("healpix") From 907f9ce8db775a85c27f862c86d21f3adb337283 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Tue, 17 Jun 2025 17:25:44 +0200 Subject: [PATCH 10/15] partially implement `zoom_to` for H3 Fetching the children will need some thought, the children are returned to a flattened list while we need a 2D array (with masks) --- xdggs/h3.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/xdggs/h3.py b/xdggs/h3.py index c3d5bfbf..55c63ce1 100644 --- a/xdggs/h3.py +++ b/xdggs/h3.py @@ -12,12 +12,14 @@ import xarray as xr try: + from h3ronpy import change_resolution from h3ronpy.vector import ( cells_to_coordinates, cells_to_wkb_polygons, coordinates_to_cells, ) except ImportError: + from h3ronpy.arrow import change_resolution from h3ronpy.arrow.vector import ( cells_to_coordinates, cells_to_wkb_polygons, @@ -201,13 +203,11 @@ def cell_boundaries(self, cell_ids, backend="shapely"): raise ValueError("invalid backend: {backend!r}") return backend_func(wkb) - def parents(self, cell_ids, resolution): - if resolution >= self.resolution: - raise ValueError( - f"resolution is not a parent: {resolution} >= {self.resolution}" - ) + def zoom_to(self, cell_ids, level): + if level > self.level: + raise ValueError("extracting children is not supported for H3, yet.") - return change_resolution(cell_ids, resolution) + return np.asarray(change_resolution(cell_ids, level)) @register_dggs("h3") From ffbdc731e322de78886b2b67e46a049393aec2c1 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Tue, 17 Jun 2025 17:33:54 +0200 Subject: [PATCH 11/15] check the implementation for healpix --- xdggs/tests/test_healpix.py | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/xdggs/tests/test_healpix.py b/xdggs/tests/test_healpix.py index bc75eef7..3e880fe6 100644 --- a/xdggs/tests/test_healpix.py +++ b/xdggs/tests/test_healpix.py @@ -371,6 +371,48 @@ def test_geographic2cell_ids( np.testing.assert_equal(actual, expected) + @pytest.mark.parametrize( + ["level", "cell_ids", "new_level", "expected"], + ( + pytest.param( + 1, + np.array([0, 4, 8, 12, 16]), + 0, + np.array([0, 1, 2, 3, 4]), + id="level1-parents", + ), + pytest.param( + 1, + np.array([0, 1, 2, 3]), + 2, + np.array( + [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [12, 13, 14, 15]] + ), + id="level1-children", + ), + pytest.param( + 1, + np.array([0, 4]), + 3, + np.stack([np.arange(16), 4 * 4**2 + np.arange(16)]), + id="level1-grandchildren", + ), + ), + ) + def test_zoom_to(self, level, cell_ids, new_level, expected): + grid = healpix.HealpixInfo(level=level, indexing_scheme="nested") + + actual = grid.zoom_to(cell_ids, level=new_level) + + np.testing.assert_equal(actual, expected) + + def test_zoom_to_ring(self): + cell_ids = np.array([1, 2, 3]) + grid = healpix.HealpixInfo(level=1, indexing_scheme="ring") + + with pytest.raises(ValueError, match="Scaling does not make sense.*'ring'.*"): + grid.zoom_to(cell_ids, level=0) + @pytest.mark.parametrize( ["mapping", "expected"], From 03ac22f3a11ac903857e5ff12a40a472c229cbc2 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Wed, 18 Jun 2025 11:28:26 +0200 Subject: [PATCH 12/15] install `healpix-geo` from PyPI --- ci/environment.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ci/environment.yml b/ci/environment.yml index d33d87d0..d45bdc07 100644 --- a/ci/environment.yml +++ b/ci/environment.yml @@ -16,3 +16,6 @@ dependencies: - arro3-core - cdshealpix - h3ronpy + - pip + - pip: + - healpix-geo From cf4f3a662cad001b0cfef8fa58ced9e982c9c5eb Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Wed, 18 Jun 2025 11:36:10 +0200 Subject: [PATCH 13/15] check `h3`'s implementation of `zoom_to` --- xdggs/tests/test_h3.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/xdggs/tests/test_h3.py b/xdggs/tests/test_h3.py index e1c2629d..5c1b6967 100644 --- a/xdggs/tests/test_h3.py +++ b/xdggs/tests/test_h3.py @@ -209,6 +209,38 @@ def test_cell_boundaries(self, level, cell_ids, backend, expected_coords): shapely.testing.assert_geometries_equal(converter(actual), expected) + @pytest.mark.parametrize( + ["level", "cell_ids", "new_level", "expected"], + ( + pytest.param( + 3, + np.array([0x832833FFFFFFFFF, 0x832834FFFFFFFFF, 0x832835FFFFFFFFF]), + 2, + np.array([]), + id="parents", + ), + pytest.param( + 3, + np.array([0x832833FFFFFFFFF, 0x832834FFFFFFFFF, 0x832835FFFFFFFFF]), + 1, + np.array([]), + id="grandparents", + ), + ), + ) + def test_zoom_to(self, level, cell_ids, new_level, expected): + grid = h3.H3Info(level=level) + + actual = grid.zoom_to(cell_ids, level=new_level) + np.testing.assert_equal(actual, expected) + + def test_zoom_to_children_not_implemented(self): + grid = h3.H3Info(level=3) + cell_ids = np.arange(3) + + with pytest.raises(NotImplementedError): + grid.zoom_to(cell_ids, level=4) + @pytest.mark.parametrize("level", levels) @pytest.mark.parametrize("dim", dims) From 1938fcd74d33f8c7e83dff79edc9dde211e33895 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Wed, 18 Jun 2025 11:50:56 +0200 Subject: [PATCH 14/15] change to `NotImplementedError` --- xdggs/h3.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/xdggs/h3.py b/xdggs/h3.py index 55c63ce1..e379ef35 100644 --- a/xdggs/h3.py +++ b/xdggs/h3.py @@ -205,7 +205,9 @@ def cell_boundaries(self, cell_ids, backend="shapely"): def zoom_to(self, cell_ids, level): if level > self.level: - raise ValueError("extracting children is not supported for H3, yet.") + raise NotImplementedError( + "extracting children is not supported for H3, yet." + ) return np.asarray(change_resolution(cell_ids, level)) From 93e7aa4d2dc89e6fe455a44088923d1ce7797c01 Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Wed, 18 Jun 2025 11:55:36 +0200 Subject: [PATCH 15/15] actually fill in the expected cell ids --- xdggs/tests/test_h3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/xdggs/tests/test_h3.py b/xdggs/tests/test_h3.py index 5c1b6967..c9802830 100644 --- a/xdggs/tests/test_h3.py +++ b/xdggs/tests/test_h3.py @@ -216,14 +216,14 @@ def test_cell_boundaries(self, level, cell_ids, backend, expected_coords): 3, np.array([0x832833FFFFFFFFF, 0x832834FFFFFFFFF, 0x832835FFFFFFFFF]), 2, - np.array([]), + np.array([0x822837FFFFFFFFF, 0x822837FFFFFFFFF, 0x822837FFFFFFFFF]), id="parents", ), pytest.param( 3, np.array([0x832833FFFFFFFFF, 0x832834FFFFFFFFF, 0x832835FFFFFFFFF]), 1, - np.array([]), + np.array([0x81283FFFFFFFFFF, 0x81283FFFFFFFFFF, 0x81283FFFFFFFFFF]), id="grandparents", ), ),