A new collection provider XdggsCollectionProvider for xdggs datasets.

pydggsapi/dependencies/collections_providers/xdggs_collection_provider.py

@@ -0,0 +1,113 @@
+from pydggsapi.dependencies.collections_providers.zarr_collection_provider import ZarrCollectionProvider, ZarrDatasourceInfo
+from pydggsapi.schemas.api.collection_providers import CollectionProviderGetDataReturn, CollectionProviderGetDataDictReturn
+
+import xarray as xr
+import xdggs
+import logging
+from pygeofilter.ast import AstType
+from datetime import datetime
+from typing import List, Any, Dict
+from dataclasses import dataclass, field
+
+logger = logging.getLogger()
+
+
+def _legacy_convertor_translator(dggs_attribute: Dict[str, str]) -> (str, str, Dict[str, str]):
+    dggs = dggs_attribute.copy()
+    if (dggs.get("name")):
+        dggs.update({"grid_name": dggs.pop("name")})
+    if (dggs.get("refinement_level")):
+        dggs.update({"level": dggs.pop("refinement_level")})
+    compression = dggs.pop("compression")
+    id_col = dggs.pop("coordinate")
+    dggs.pop("spatial_dimension")
+    return id_col, compression, dggs
+
+
+@dataclass
+class XdggsDatasourceInfo(ZarrDatasourceInfo):
+    index_kind: str = None
+    compression: str = "none"
+    dggs_metadata: Dict[str, str] = field(default_factory=lambda: {})
+
+
+# xdggs provider inhert from ZarrCollectionProvider
+class XdggsCollectionProvider(ZarrCollectionProvider):
+
+    def __init__(self, datasources):
+        self.datasources = {}
+        try:
+            for k, v in datasources.items():
+                zones_group = {}
+                datasource = XdggsDatasourceInfo(**v)
+                dt = xr.open_datatree(datasource.filepath, engine="zarr", chunks="auto",
+                                      storage_options=datasource.storage_options, consolidated=True)
+                # we need a new datatree to hold the xdggs dataset
+                new_dt = xr.DataTree(name=k)
+                skip = [True] * len(dt.groups)
+                # id_col, there are multiple gropus and the column name of cell id can be different.
+                # the coordinate of the last group will be used as the id_col
+                for i, group in enumerate(dt.groups):
+                    if (group != "/"):
+                        ds = dt[group].to_dataset()
+                        group_name = group[1:]
+                        # xdggs infer a convention from the dataset.
+                        try:
+                            ds = xdggs.decode(ds, index_options={"index_kind": datasource.index_kind,
+                                                                 "compression": datasource.compression})
+                            datasource.id_col = ds.dggs.coord.name
+                            # save the dggs attribute for future use
+                            datasource.dggs_metadata = ds.dggs.grid_info.to_dict()
+                            new_dt[group_name] = xr.DataTree(dataset=ds)
+                            zones_group[str(ds.dggs.grid_info.level)] = group_name
+                            skip[i] = False
+                            # new_dt = new_dt.assign({group_name: xr.DataTree(dataset=ds)})
+                        except Exception as e:
+                            logger.warning(f'{__name__} error when decoding {group} in {k} to xdggs: {e} , skipped.')
+                            continue
+                # skip adding the datasource if all groups failed from decode as xdggs dataset
+                if (not all(skip)):
+                    datasource.zone_groups = zones_group
+                    datasource.filehandle = new_dt
+                    self.datasources[k] = datasource
+                    logger.info(f"{__name__} added datasource: {k}")
+        except Exception as e:
+            logger.error(f'{__name__} create datasource failed: {e}')
+            raise Exception(f'{__name__} create datasource failed: {e}')
+
+    def get_data(self, zoneIds: List[Any], res: int, datasource_id: str,
+                 cql_filter: AstType = None, include_datetime: bool = False,
+                 include_properties: List[str] = None,
+                 exclude_properties: List[str] = None,
+                 input_zoneIds_padding: bool = True,
+                 collection_timestamp: datetime = None,
+                 check_if_exists_only: bool = False) -> CollectionProviderGetDataReturn:
+        # since xarray-sql doesn't work with xdggs xindexes,
+        # I have to change it back to a normal dataset when the cql_filter is not None
+        result = CollectionProviderGetDataReturn(zoneIds=[], cols_meta={}, data=[])
+        datasource = self.datasources[datasource_id]
+        datatree = datasource.filehandle
+        try:
+            group_name = datasource.zone_groups[str(res)]
+        except KeyError as e:
+            logger.error(f'{__name__} get zone_grp for resolution {res} failed: {e}')
+            return result
+        target_ds = datatree[group_name].to_dataset()
+        index_type = type(target_ds.xindexes.get(datasource.id_col))
+        restore = False
+        if (cql_filter is not None and (not isinstance(index_type, xr.core.indexes.PandasIndex))):
+            dim = target_ds.dggs.coord.dims[0]
+            temp = target_ds.drop_indexes(datasource.id_col)
+            temp = temp.swap_dims({dim: datasource.id_col})
+            self.datasources[datasource_id].filehandle[group_name] = xr.DataTree(dataset=temp)
+            restore = True
+
+        result = super().get_data(zoneIds, res, datasource_id, cql_filter, include_datetime,
+                                  include_properties, exclude_properties, input_zoneIds_padding,
+                                  collection_timestamp, check_if_exists_only)
+        if (restore):
+            #target_ds[datasource.id_col].attrs.update(datasource.dggs_metadata)
+            #target_ds = xdggs.decode(target_ds, index_options={"index_kind": datasource.index_kind,
+            #                                                   "compression": datasource.compression})
+            self.datasources[datasource_id].filehandle[group_name] = xr.DataTree(dataset=target_ds)
+        return result

pydggsapi/dependencies/collections_providers/zarr_collection_provider.py

@@ -22,6 +22,8 @@
 from dataclasses import dataclass
 import logging
 
+import time
+
 logger = logging.getLogger()
 
 
@@ -95,7 +97,6 @@ def get_data(self, zoneIds: List[Any], res: int, datasource_id: str,
                 table = read_xarray_table(ds)
                 ctx = SessionContext()
                 ctx.register_table('ds', table)
-
                 if ("*" in datasource.data_cols):
                     incl = ",".join(include_properties) if include_properties else "*"
                     excl = datasource.exclude_data_cols or []
@@ -109,13 +110,23 @@ def get_data(self, zoneIds: List[Any], res: int, datasource_id: str,
                         incl -= set(exclude_properties)
                     cols = f"{','.join(incl)}, {id_col}"
                 sql = f"""select {cols} from ds where ("{id_col}" in ({', '.join(f"'{z}'" for z in zoneIds)})) and ({cql_sql}) """
+                start = time.time()
                 zarr_result = xr.Dataset.from_dataframe(ctx.sql(sql).to_pandas().set_index(id_col))
+                logger.debug(f"query time: {time.time() - start}")
             else:
                 cols = OrderedSet(ds.data_vars) if ("*" in datasource.data_cols) else OrderedSet(datasource.data_cols)
                 cols = list(cols - OrderedSet(datasource.exclude_data_cols))
+                # for xdggs: when dimension (cells) != cell_ids to_dataframe()
+                # will result in an extra column for the index.
                 dimension = ds[id_col].dims[0]
-                zarr_result = ds.query({dimension: f'{id_col} in {zoneIds}'})
+                start = time.time()
+                pos = np.where(ds[id_col].isin(zoneIds))[0]
+                zarr_result = ds.isel({dimension: pos})
+                # zarr_result = ds.query({dimension: f'{id_col} in {zoneIds}'})
+                logger.debug(f"query time: {time.time() - start}")
                 zarr_result = zarr_result[cols]
+                if ("xdggs" in str(zarr_result.xindexes.get(id_col).__class__)):
+                    zarr_result = zarr_result.drop_indexes(id_col).swap_dims({dimension: id_col})
         except Exception as e:
             # Zarr will raise exception if nothing matched
             logger.error(f'{__name__} {datasource_id} sel failed: {e}')

pydggsapi/routers/tiles_api.py

@@ -63,7 +63,6 @@ async def query_mvt_tiles(
         id_col = "zone_id"
     bbox, tile = mercator.getWGS84bbox(tilesreq.z, tilesreq.x, tilesreq.y)
     res_info = mercator.get(tile.z)
-    print(tile.z)
     tile_width_km = float(res_info["Tile width deg lons"]) / 0.01 * 0.4  # in tile_width_km
     zone_level = dggrs_provider.get_zone_level_by_cls(tile_width_km)
     if (tilesreq.relative_depth != 0):
@@ -96,8 +95,16 @@ async def query_mvt_tiles(
             convertedlist = gpd.GeoDataFrame({'zone_id': converted.target_zoneIds,
                                               'target_id': converted.zoneIds}).set_index('target_id')
             zoneslist = convertedlist.join(zoneslist).reset_index().set_index('zone_id')
-    zones_data = collection_provider.get_data(zoneslist.index.to_list(), zone_level,
-                                              collection.collection_provider.datasource_id, input_zoneIds_padding=False)
+    try:
+        zones_data = collection_provider.get_data(zoneslist.index.to_list(), zone_level,
+                                                  collection.collection_provider.datasource_id, input_zoneIds_padding=False)
+    except Exception as e:
+        logger.error(f"{__name__} get_data error: {e}")
+        content = mapbox_vector_tile.encode({"name": tilesreq.collectionId, "features": []},
+                                            quantize_bounds=bbox,
+                                            default_options={"transformer": transformer.transform})
+        return Response(bytes(content), media_type="application/x-protobuf")
+
     if (len(zones_data.zoneIds) == 0):
         content = mapbox_vector_tile.encode({"name": tilesreq.collectionId, "features": []},
                                             quantize_bounds=bbox,

pyproject.toml

@@ -50,6 +50,7 @@ zarr = "^3.2.1"
 s3fs = ">=2026.2.0,<2027.0.0"
 gcsfs = ">=2026.2.0,<2027.0.0"
 healpix-geo = "0.4.1"
+xdggs = "0.7.0"
 
 
 [tool.poetry.group.dev.dependencies]

tests/dggs_api_config_testing.json

@@ -67,7 +67,7 @@
         "description": "A HEALPix dataset converted from a clipped TWI near Elva with 10M spatial resolution. The dataset",
 	"extent": {"spatial": {"bbox": [[26.2718,  58.2267, 26.5815,  58.3295]]}},
         "collection_provider": {
-          "providerId": "zarr",
+          "providerId": "xdggs",
           "dggrsId": "healpixgeo_nested",
           "dggrs_zoneid_repr": "int",
           "min_refinement_level": 5,
@@ -227,6 +227,17 @@
 	      }
         }
      }
+    },
+    "3": {
+      "xdggs": {
+        "classname": "xdggs_collection_provider.XdggsCollectionProvider",
+        "datasources": {
+          "est_topo_twi_10m_elva_healpix_nested": {
+		  	  "filepath": "testing_data/est_topo_twi_clipped_10m_healpix_nested_datatree.zarr",
+			  "index_kind": "moc"
+	      }
+        }
+     }
     }
   }
 }

tests/test_cql/test_cql_data-retrieval_query_zone_data.py

@@ -2,6 +2,7 @@
 from pydggsapi.schemas.ogc_dggs.dggrs_zones_data import ZonesDataDggsJsonResponse, ZonesDataGeoJson
 from pydggsapi.dependencies.collections_providers.parquet_collection_provider import ParquetCollectionProvider
 from pydggsapi.dependencies.collections_providers.zarr_collection_provider import ZarrCollectionProvider
+from pydggsapi.dependencies.collections_providers.xdggs_collection_provider import XdggsCollectionProvider
 from fastapi.testclient import TestClient
 import pytest
 import importlib
@@ -116,7 +117,14 @@
         # because zone data retrieval separates the handling for temporal/non-temporal collections
         # and we are not testing the datetime parameter here, so it is default to non-temporal collections
         continue
-    if (isinstance(cp, ZarrCollectionProvider)):
+    if (isinstance(cp, XdggsCollectionProvider)):
+        # It is a data tree
+        dt = cp.datasources[datasource_id].filehandle
+        zone_groups = cp.datasources[datasource_id].zone_groups
+        id_col = cp.datasources[datasource_id].id_col
+        ds1 = dt[zone_groups[str(rf1)]].to_dataset().drop_indexes(id_col).reset_coords().to_dataframe().reset_index(drop=True)
+        ds1 = ds1.drop('spatial_ref', axis=1).reset_index() if ("saptial_ref" in ds1.columns) else ds1
+    elif (isinstance(cp, ZarrCollectionProvider)):
         # It is a data tree
         dt = cp.datasources[datasource_id].filehandle
         zone_groups = cp.datasources[datasource_id].zone_groups

tests/test_cql/test_cql_zone-query_dggrs_zones.py

@@ -2,6 +2,7 @@
 from pydggsapi.schemas.ogc_dggs.dggrs_zones import ZonesResponse, ZonesGeoJson
 from pydggsapi.dependencies.collections_providers.parquet_collection_provider import ParquetCollectionProvider
 from pydggsapi.dependencies.collections_providers.zarr_collection_provider import ZarrCollectionProvider
+from pydggsapi.dependencies.collections_providers.xdggs_collection_provider import XdggsCollectionProvider
 from fastapi.testclient import TestClient
 import pytest
 import importlib
@@ -103,7 +104,20 @@
     rf2 = max_rf - 1
     ds1, ds2 = None, None
     cp = collection_providers_dict[collection.collection_provider.providerId]
-    if (isinstance(cp, ZarrCollectionProvider)):
+    if (isinstance(cp, XdggsCollectionProvider)):
+        # It is a data tree
+        dt = cp.datasources[datasource_id].filehandle
+        zone_groups = cp.datasources[datasource_id].zone_groups
+        id_col = cp.datasources[datasource_id].id_col
+        ds1 = dt[zone_groups[str(rf1)]].to_dataset().drop_indexes(id_col).reset_coords().to_dataframe().reset_index(drop=True)
+        ds2 = dt[zone_groups[str(rf2)]].to_dataset().drop_indexes(id_col).reset_coords().to_dataframe().reset_index(drop=True)
+        if ("cell_ids" in ds1.columns or "cell_ids" in ds2.columns):
+            ds1 = ds1.rename(columns={"cell_ids": "zone_id"})
+            ds2 = ds2.rename(columns={"cell_ids": "zone_id"})
+        ds1 = ds1.drop('spatial_ref', axis=1).reset_index() if ("saptial_ref" in ds1.columns) else ds1
+        ds2 = ds2.drop('spatial_ref', axis=1).reset_index() if ("saptial_ref" in ds2.columns) else ds2
+        test_rf = {rf1: [ds1], rf2: [ds2]}
+    elif (isinstance(cp, ZarrCollectionProvider)):
         # It is a data tree
         dt = cp.datasources[datasource_id].filehandle
         zone_groups = cp.datasources[datasource_id].zone_groups
@@ -178,8 +192,9 @@ def test_cql_zone_query_dggrs_zones():
                 if (validation_data.empty):
                     continue
                 print(f"CQL to test: {cql_string}")
-                print(f"test data: {validation_data}")
-                print(f"Success test case with dggs zones query ({dggrsid}, bbox: {aoi.bounds}, zone_level={rf}, compact=False")
+                print(f" - test data: {validation_data}")
+                print(f" - Success test case with dggs zones query ({dggrsid}, bbox: {aoi.bounds}, zone_level={rf}, \
+                        compact=False, cql: {cql_string})")
                 response = client.get(f'/dggs-api/dggs/{dggrsid}/zones', params={"bbox": ",".join(bounds),
                                                                                  "zone-level": rf,
                                                                                  "compact-zones": False,
@@ -196,18 +211,21 @@ def test_cql_zone_query_dggrs_zones():
             test_values = dataset.iloc[iloc_pos].values
             for cql_string in cql_composition:
                 cql_string = cql_string.format(var_name=column_names[0], value1=test_values[0][0], value2=test_values[1][0])
-                print(f"Success test case with dggs zones query ({dggrsid}, bbox: {aoi.bounds}, zone_level={rf}, \
-                        compact=False, cql: {cql_string})")
-                response = client.get(f'/dggs-api/dggs/{dggrsid}/zones', params={"bbox": ",".join(bounds),
-                                                                                        "zone-level": rf,
-                                                                                        "compact-zones": False,
-                                                                                        "filter": cql_string})
                 if ("AND" in cql_string):
                     validation_data = dataset[(dataset[column_names[0]] >= test_values[0][0]) &
                                               (dataset[column_names[0]] <= test_values[1][0])]
                 elif ("OR" in cql_string):
                     validation_data = dataset[(dataset[column_names[0]] == test_values[0][0]) |
                                               (dataset[column_names[0]] == test_values[1][0])]
+                print(f"CQL to test: {cql_string}")
+                print(f" - test data: {validation_data}")
+                print(f"Success test case with dggs zones query ({dggrsid}, bbox: {aoi.bounds}, zone_level={rf}, \
+                        compact=False, cql: {cql_string})")
+                response = client.get(f'/dggs-api/dggs/{dggrsid}/zones', params={"bbox": ",".join(bounds),
+                                                                                 "zone-level": rf,
+                                                                                 "compact-zones": False,
+                                                                                 "filter": cql_string})
+
                 if (validation_data.shape[0] > 0):
                     assert response.status_code == 200
                     zones = ZonesResponse(**response.json())

tests/test_temporal_collections/test_temporal_zone-query_dggrs_zones.py

@@ -12,6 +12,7 @@
 import tempfile
 import shapely
 import numpy as np
+import xarray as xr
 import pandas as pd
 from dggal import Application, pydggal_setup, CRS, ogc, epsg, GeoExtent, Array, GeoPoint
 from dggal import IVEA7H, ISEA7H_Z7, rHEALPix, HEALPix
@@ -107,8 +108,15 @@
         # It is a data tree
         dt = cp.datasources[datasource_id].filehandle
         zone_groups = cp.datasources[datasource_id].zone_groups
-        ds1 = dt[zone_groups[rf1]].to_dataset().to_dataframe()
-        ds2 = dt[zone_groups[rf2]].to_dataset().to_dataframe()
+        id_col = cp.datasources[datasource_id].id_col
+        ds1 = dt[zone_groups[rf1]].to_dataset()
+        ds2 = dt[zone_groups[rf2]].to_dataset()
+        if (not isinstance(ds1.xindexes.get(id_col), xr.core.indexes.PandasIndex)):
+            ds1 = ds1.drop_indexes(id_col).reset_coords()
+        if (not isinstance(ds2.xindexes.get(id_col), xr.core.indexes.PandasIndex)):
+            ds2 = ds2.drop_indexes(id_col).reset_coords()
+        ds1 = ds1.to_dataframe()
+        ds2 = ds2.to_dataframe()
         ds1 = ds1.drop('spatial_ref', axis=1).reset_index() if ("saptial_ref" in ds1.columns) else ds1
         ds2 = ds2.drop('spatial_ref', axis=1).reset_index() if ("saptial_ref" in ds2.columns) else ds2
     elif (isinstance(cp, ParquetCollectionProvider)):