adds support for ISEA7H and ISEA9R (others supported by DGGAL incoming)

Author: alpha-beta-soup

.gitignore

@@ -4,3 +4,4 @@ conda*
 .conda*
 .vscode/
 .spyproject/
+tests/data/input

README.md

@@ -265,7 +265,7 @@ In brief, to get started:
     - If you're on Windows, `pip install gdal` may be necessary before running the subsequent commands.
     - On Linux, install GDAL 3.6+ according to your platform-specific instructions, including development headers, i.e. `libgdal-dev`.
 - Create the virtual environment with `poetry init`. This will install necessary dependencies.
-- Subsequently, the virtual environment can be re-activated with `poetry shell`.
+- Subsequently, the virtual environment can be re-activated with `poetry env activate`.
 
 If you run `poetry install`, the CLI tool will be aliased so you can simply use `raster2dggs` rather than `poetry run raster2dggs`, which is the alternative if you do not `poetry install`.
 

poetry.lock

@@ -38,6 +38,7 @@ maidenhead = "^1.8"
 s2sphere = "^0.2"
 pya5 = "^0.5"
 geoarrow-pyarrow = "^0.2.0"
+dggal = "^0.0.6"
 
 [tool.poetry.group.dev.dependencies]
 pytest = "^7.2.2"

pyproject.toml

@@ -7,6 +7,8 @@
 from raster2dggs.maidenhead import maidenhead
 from raster2dggs.s2 import s2
 from raster2dggs.a5 import a5
+from raster2dggs.isea9r import isea9r
+from raster2dggs.isea7h import isea7h
 
 
 @click.group()
@@ -21,6 +23,8 @@ def cli():
 cli.add_command(maidenhead)
 cli.add_command(s2)
 cli.add_command(a5)
+cli.add_command(isea9r)
+cli.add_command(isea7h)
 
 
 def main():

raster2dggs/cli.py

@@ -261,7 +261,7 @@ def address_boundary_issues(
         }
     )
     out_meta.index = pd.Index([], name=index_col, dtype="object")
-
+    print(ddf.compute(), band_cols, out_meta)
     with TqdmCallback(desc=f"Aggregating{'/compacting' if kwargs['compact'] else ''}"):
         ddf = ddf.map_partitions(
             indexer.parent_groupby,
@@ -289,7 +289,11 @@ def address_boundary_issues(
 
             write_tasks = [
                 dask.delayed(write_partition_as_geoparquet)(
-                    part, geo_serialisation_method, output, partition_col, kwargs["compression"]
+                    part,
+                    geo_serialisation_method,
+                    output,
+                    partition_col,
+                    kwargs["compression"],
                 )
                 for part in delayed_parts
             ]

raster2dggs/common.py

@@ -7,6 +7,14 @@
 MIN_MAIDENHEAD, MAX_MAIDENHEAD = 1, 6
 MIN_S2, MAX_S2 = 0, 30
 MIN_A5, MAX_A5 = 0, 30
+MIN_ISEA9R, MAX_ISEA9R = (
+    0,
+    18,
+)  # NB correct for packing into 64 bit integer? see https://dggal.org/docs/html/dggal/Classes/DGGRS/Virtual%20Methods/getMaxDGGRSZoneLevel.html
+MIN_ISEA7H, MAX_ISEA7H = (
+    0,
+    18,
+) # NB correct?
 
 DEFAULT_NAME: str = "value"
 
@@ -33,6 +41,12 @@
     "maidenhead": lambda resolution: MIN_MAIDENHEAD,
     "s2": lambda resolution: max(MIN_S2, (resolution - DEFAULT_PARENT_OFFSET)),
     "a5": lambda resolution: max(MIN_A5, (resolution - DEFAULT_PARENT_OFFSET)),
+    "isea9r": lambda resolution: max(
+        MIN_ISEA9R, (resolution - 5)
+    ),  # NB use 5 for IS/VEA9R, and 10 for IS/VEA3H, and 8 for GNOSIS --- corresponds to ~64K sub-zones
+    "isea7h": lambda resolution: max(
+        MIN_ISEA7H, (resolution - DEFAULT_PARENT_OFFSET)
+    )
 }
 
 
@@ -44,6 +58,8 @@ def zero_padding(dggs: str) -> int:
         "maidenhead": MAX_MAIDENHEAD,
         "s2": MAX_S2,
         "a5": MAX_A5,
+        "isea9r": MAX_ISEA9R,
+        "isea7h": MAX_ISEA7H,
     }
     max_res = max_res_lookup.get(dggs)
     if max_res is None:

raster2dggs/constants.py

@@ -10,6 +10,7 @@
 import raster2dggs.indexers.maidenheadrasterindexer as maidenheadrasterindexer
 import raster2dggs.indexers.s2rasterindexer as s2rasterindexer
 import raster2dggs.indexers.a5rasterindexer as a5rasterindexer
+import raster2dggs.indexers.dggalrasterindexer as dggalrasterindexer
 
 """
 Match DGGS name to indexer class name
@@ -21,6 +22,8 @@
     "maidenhead": maidenheadrasterindexer.MaidenheadRasterIndexer,
     "s2": s2rasterindexer.S2RasterIndexer,
     "a5": a5rasterindexer.A5RasterIndexer,
+    "isea9r": dggalrasterindexer.ISEA9RRasterIndexer,
+    "isea7h": dggalrasterindexer.ISEA7HRasterIndexer,
 }
 
 

raster2dggs/indexerfactory.py

@@ -0,0 +1,187 @@
+"""
+@author: ndemaio
+"""
+
+from abc import abstractmethod
+from functools import reduce
+from numbers import Number
+from typing import Tuple
+
+import h3pandas  # Necessary import despite lack of explicit use
+
+import dggal
+import pandas as pd
+import pyarrow as pa
+import xarray as xr
+import numpy as np
+import shapely
+
+import raster2dggs.constants as const
+
+from raster2dggs.indexers.rasterindexer import RasterIndexer
+
+
+# TODO memoisation; needs to use hashable inputs and outputs (int)
+def apply_n_reduce(f, n, x):
+    def step(v):
+        return f(v)[0]
+
+    return reduce(lambda v, _: step(v), range(n), x)
+
+
+class DGGALRasterIndexer(RasterIndexer):
+    """
+    Provides integration for Uber's H3 DGGS.
+    """
+
+    @property
+    @abstractmethod
+    def dggrs(self) -> dggal.DGGRS:
+        raise NotImplementedError
+
+    @property
+    @abstractmethod
+    def refinementRatio(self) -> int:
+        raise NotImplementedError
+
+    def index_func(
+        self,
+        sdf: xr.DataArray,
+        resolution: int,
+        parent_res: int,
+        nodata: Number = np.nan,
+        band_labels: Tuple[str] = None,
+    ) -> pa.Table:
+        """
+        Index a raster window to a DGGRS.
+        Subsequent steps are necessary to resolve issues at the boundaries of windows.
+        If windows are very small, or in strips rather than blocks, processing may be slower
+        than necessary and the recommendation is to write different windows in the source raster.
+
+        Implementation of interface function.
+        """
+        sdf: pd.DataFrame = (
+            sdf.to_dataframe().drop(columns=["spatial_ref"]).reset_index()
+        )
+        subset: pd.DataFrame = sdf.dropna()
+        subset = subset[subset.value != nodata]
+        subset = pd.pivot_table(
+            subset, values=const.DEFAULT_NAME, index=["x", "y"], columns=["band"]
+        ).reset_index()
+        # Primary DGGSRS index
+        cells = [
+            self.dggrs.getZoneFromWGS84Centroid(resolution, dggal.GeoPoint(lon, lat))
+            for lon, lat in zip(subset["y"], subset["x"])
+        ]  # Vectorised
+        dggrs_parent = [
+            apply_n_reduce(self.dggrs.getZoneParents, resolution - parent_res, zone)
+            for zone in cells
+        ]
+        subset = subset.drop(columns=["x", "y"])
+        index_col = self.index_col(resolution)
+        subset[index_col] = pd.Series(
+            map(self.dggrs.getZoneTextID, cells), index=subset.index
+        )
+        partition_col = self.partition_col(parent_res)
+        subset[partition_col] = pd.Series(
+            map(self.dggrs.getZoneTextID, dggrs_parent), index=subset.index
+        )
+        # Rename bands
+        bands = sdf["band"].unique()
+        columns = dict(zip(bands, band_labels))
+        subset = subset.rename(columns=columns)
+        return pa.Table.from_pandas(subset)
+
+    def cell_to_children_size(self, cell: str, desired_resolution: int) -> int:
+        """
+        Determine total number of children at some offset resolution
+
+        Implementation of interface function.
+        """
+        current_resolution = self.dggrs.getZoneLevel(self.dggrs.getZoneFromTextID(cell))
+        n = desired_resolution - current_resolution
+        return self.refinementRatio**n
+
+    @staticmethod
+    def valid_set(cells: set) -> set[str]:
+        """
+        Implementation of interface function.
+        """
+        return set(filter(lambda c: (not pd.isna(c)), cells))
+
+    def parent_cells(self, cells: set, resolution) -> map:
+        """
+        Implementation of interface function.
+        """
+        # TODO appropriately handle potential for multiple parentage in dggal (e.g. ISEAH3)
+        child_resolution = self.dggrs.getZoneLevel(self.dggrs.getZoneFromTextID(next(iter(cells))))
+        return map(
+            lambda zone: self.dggrs.getZoneTextID(apply_n_reduce(
+                self.dggrs.getZoneParents, child_resolution - resolution, self.dggrs.getZoneFromTextID(zone)
+            )), cells
+        )
+
+    def expected_count(self, parent: str, resolution: int):
+        """
+        Implementation of interface function.
+        """
+        return self.cell_to_children_size(parent, resolution)
+
+    def cell_to_point(self, cell: str) -> shapely.geometry.Point:
+        geo_point : dggal.GeoPoint = self.dggrs.getZoneWGS84Centroid(self.dggrs.getZoneFromTextID(cell))
+        return shapely.Point(
+            geo_point.lon, geo_point.lat
+        )
+
+    def cell_to_polygon(self, cell: str) -> shapely.geometry.Polygon:
+        geo_points : List[dggal.GeoPoint] = self.dggrs.getZoneWGS84Vertices(self.dggrs.getZoneFromTextID(cell))
+        return shapely.Polygon(
+            tuple([(p.lon, p.lat) for p in geo_points])
+        )
+
+
+# def ISEA7HRasterIndexer(DGGALRasterIndexer):
+#     def __init__(self, dggs: str):
+#         super().__init__(dggs)
+#         self.dggrs = dggal.ISEA7H()
+
+# NB  All zones of GNOSIS Global Grid and ISEA9R have single parents, whereas ISEA3H zones have one parent if they are a centroid child, and three parents otherwise if they are a vertex child.  See dggrs.getMaxParents()
+
+
+class ISEA9RRasterIndexer(DGGALRasterIndexer):
+    """
+    A raster indexer for the ISEA9R DGGS, an axis-aligned and equal-area DGGH based on the Icosahedral Snyder Equal-Area (ISEA) planar projection using rhombuses with a refinement ratio of 9.
+    """
+
+    def __init__(self, dggs: str):
+        super().__init__(dggs)
+        dggal.pydggal_setup(dggal.Application(appGlobals=globals()))
+        self._dggrs = dggal.ISEA9R()
+        self._refinementRatio = self.dggrs.getRefinementRatio()
+
+    @property
+    def dggrs(self) -> dggal.DGGRS:
+        return self._dggrs
+
+    @property
+    def refinementRatio(self) -> int:
+        return self._refinementRatio
+
+class ISEA7HRasterIndexer(DGGALRasterIndexer):
+    """
+    A raster indexer for the ISEA7H DGGS, an equal-area hexagonal grid with a refinement ratio of 7 defined in the ISEA projection.
+    """
+
+    def __init__(self, dggs: str):
+        super().__init__(dggs)
+        dggal.pydggal_setup(dggal.Application(appGlobals=globals()))
+        self._dggrs = dggal.ISEA7H()
+        self._refinementRatio = self.dggrs.getRefinementRatio()
+
+    @property
+    def dggrs(self) -> dggal.DGGRS:
+        return self._dggrs
+
+    @property
+    def refinementRatio(self) -> int:
+        return self._refinementRatio

raster2dggs/indexers/dggalrasterindexer.py

@@ -11,7 +11,7 @@
 import numpy as np
 
 from .. import constants as const
-from .. interfaces import IRasterIndexer
+from ..interfaces import IRasterIndexer
 
 
 class RasterIndexer(IRasterIndexer):
@@ -31,6 +31,12 @@ def __init__(self, dggs: str):
         """
         self.dggs = dggs
 
+    def __dask_tokenize__(self):
+        """
+        Only include stable, immutable fields that define behaviour
+        """
+        return (type(self).__name__, self.dggs)
+
     def index_col(self, resolution):
         pad_width = const.zero_padding(self.dggs)
         return f"{self.dggs}_{resolution:0{pad_width}d}"
@@ -134,7 +140,21 @@ def compaction(
         if not unprocessed_indices:
             return df
         band_cols = self.band_cols(df)
+        # partition_col = self.partition_col(parent_res)
         compaction_map = {}
+
+        # sub = df.loc[list(unprocessed_indices)]
+        # for parent, group in sub.groupby(partition_col, sort=False, observed=True):
+        #     if parent in compaction_map:
+        #         continue
+        #     expected = self.expected_count(parent, resolution)
+
+        #     # All band columns identical across the group?
+        #     if len(group) == expected and (group[band_cols].nunique() == 1).all():
+        #         compact_row = group.iloc[0].copy()
+        #         compact_row.name = parent  # parent becomes the new index
+        #         compaction_map[parent] = compact_row
+        #         unprocessed_indices -= set(group.index)
         for r in range(parent_res, resolution):
             parent_cells = self.parent_cells(unprocessed_indices, r)
             parent_groups = df.loc[list(unprocessed_indices)].groupby(