diff --git a/xdggs/accessor.py b/xdggs/accessor.py index f8b67cc..97fbc12 100644 --- a/xdggs/accessor.py +++ b/xdggs/accessor.py @@ -104,3 +104,23 @@ 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. + """ + 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") diff --git a/xdggs/h3.py b/xdggs/h3.py index 17db4ed..56266a9 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): diff --git a/xdggs/healpix.py b/xdggs/healpix.py index bb76e32..b5e6579 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): diff --git a/xdggs/index.py b/xdggs/index.py index 2c3a7c2..c6a94bc 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