feat: get_collection_gdf, several fixes, improvements and docs

Author: martibosch

README.md

@@ -1,6 +1,31 @@
 # API
 
+## STAC client
+
+```{eval-rst}
+.. automodule:: swisstopopy.stac
+   :members:
+```
+
+## Geospatial utilities
+
+### Buildings
+
+```{eval-rst}
+.. automodule:: swisstopopy.buildings
+   :members:
+```
+
+### Digital elevation model (DEM)
+
+```{eval-rst}
+.. automodule:: swisstopopy.dem
+   :members:
+```
+
+### Tree canopy raster
+
 ```{eval-rst}
-.. automodule:: swisstopopy
+.. automodule:: swisstopopy.tree_canopy
    :members:
 ```

docs/api.md

@@ -4,6 +4,7 @@ hidden:
 maxdepth: 1
 ---
 
+overview
 api
 contributing
 changelog

docs/index.md

@@ -64,6 +64,9 @@ test = [
 [project.urls]
 Repository = "https://github.com/martibosch/swisstopopy"
 
+[tool.codespell]
+skip = "docs/overview.ipynb"
+
 [tool.commitizen]
 major_version_zero = true
 name = "cz_conventional_commits"

docs/overview.ipynb

@@ -4,6 +4,8 @@
 Raster and swissALTI3D products provided by swisstopo's STAC API.
 """
 
+import warnings
+
 import geopandas as gpd
 import osmnx as ox
 import pandas as pd
@@ -73,9 +75,7 @@ def get_bldg_gdf(
     region: RegionType,
     *,
     region_crs: CRSType = None,
-    surface3d_datetime: DatetimeLike | None = None,
-    alti3d_datetime: DatetimeLike | None = None,
-    alti3d_res: float = 0.5,
+    item_datetime: DatetimeLike | None = None,
     **pooch_retrieve_kwargs: utils.KwargsType,
 ) -> gpd.GeoDataFrame:
     """Get buildings geo-data frame with height information.
@@ -88,12 +88,10 @@ def get_bldg_gdf(
         Coordinate reference system (CRS) of the region. Required if `region` is a naive
         geometry or a list of bounding box coordinates. Ignored if `region` already has
         a CRS.
-    surface3d_datetime, alti3d_datetime : datetime-like, optional
-        Datetime to filter swissSURFACE3D and swissALTI3D data respectively, forwarded
-        to `pystac_client.Client.search`. If None, the latest data for each tile is
-        used.
-    alti3d_res : {0.5, 2}, default 2
-        Resolution of the swissALTI3D data to get, can be 0.5 or 2 (meters).
+    item_datetime : datetime-like, optional
+        Datetime to filter swissSURFACE3D Raster and swissALTI3D data to use (must be
+        the same for both collections), forwarded to `pystac_client.Client.search`.
+        If None, the latest data for each tile is used.
     pooch_retrieve_kwargs : mapping, optional
         Additional keyword arguments to pass to `pooch.retrieve`.
 
@@ -106,42 +104,59 @@ def get_bldg_gdf(
     # 1. we first need to project the region to OSM CRS (EPSG:4326) to query via osmnx
     # 2. we drop the "node" column to keep only the "way" and "relation" columns that
     # correspond to polygon geometries
-    region_gser = RegionMixin._process_region_arg(region, crs=region_crs)
+    region_gser = RegionMixin._process_region_arg(region, crs=region_crs)["geometry"]
     bldg_gdf = (
         ox.features_from_polygon(
             region_gser.to_crs(ox.settings.default_crs).iloc[0],
             tags=OSMNX_TAGS,
         )
         .to_crs(stac.CH_CRS)
-        .drop("node")
+        .drop("node", errors="ignore")
     )
 
     # use the STAC API to get building heights from swissSURFACE3D and swissALTI3D
     client = stac.SwissTopoClient(region_gser)
 
+    def _get_and_process_gdf(collection_id):
+        gdf = client.get_collection_gdf(
+            collection_id,
+            datetime=item_datetime,
+        )
+        if gdf.empty:
+            # warn and return gdf without heights
+            warnings.warn(
+                f"No '{collection_id}' data available for the specified region and "
+                "datetime."
+            )
+            return gdf
+        # filter to get tiff images only
+        gdf = gdf[gdf["assets.href"].str.endswith(".tif")]
+        # filter to get the resolution data at the specified resolution
+        # AFAIK, there is no swissSURFACE3D Raster collection with 2 m resolution, so we
+        # can only set the resolution to 0.5 m
+        gdf = gdf[gdf["assets.eo:gsd"] == 0.5]
+        # if no datetime specified, get the latest data for each tile (location)
+        if item_datetime is None:
+            gdf = stac.get_latest(gdf)
+        return gdf
+
     # surface3d-raster (raster dsm)
-    surface3d_gdf = client.gdf_from_collection(
-        stac.SWISSSURFACE3D_RASTER_COLLECTION_ID,
-        datetime=surface3d_datetime,
+    surface3d_raster_gdf = _get_and_process_gdf(
+        stac.SWISSSURFACE3D_RASTER_COLLECTION_ID
     )
-    # filter to get tiff images only
-    surface3d_gdf = surface3d_gdf[surface3d_gdf["assets.href"].str.endswith(".tif")]
-    # if no datetime specified, get the latest data for each tile (location)
-    if surface3d_datetime is None:
-        surface3d_gdf = stac.get_latest(surface3d_gdf)
 
     # alti3d (raster dem)
-    alti3d_gdf = client.gdf_from_collection(
+    alti3d_gdf = _get_and_process_gdf(
         stac.SWISSALTI3D_COLLECTION_ID,
-        datetime=alti3d_datetime,
     )
-    # filter to get tiff images only
-    alti3d_gdf = alti3d_gdf[alti3d_gdf["assets.href"].str.endswith(".tif")]
-    # filter to get the resolution data at the specified resolution
-    alti3d_gdf = alti3d_gdf[alti3d_gdf["assets.eo:gsd"] == alti3d_res]
-    # if no datetime specified, get the latest data for each tile (location)
-    if alti3d_datetime is None:
-        alti3d_gdf = stac.get_latest(alti3d_gdf)
+
+    # both collections need to have items for the selected datetimes
+    if surface3d_raster_gdf.empty or alti3d_gdf.empty:
+        # warn and return gdf without heights
+        # note that warnings will have already been raised above, so this follows from
+        # there
+        warnings.warn("Returning building footprints without height information.")
+        return bldg_gdf
 
     # compute the building heights as zonal statistics. To that end, we first compute
     # a "building height raster" as the difference between the swissSURFACE3D (surface
@@ -155,7 +170,7 @@ def get_bldg_gdf(
 
     # we need to project the gdf of tiles to the same CRS as the actual swissSURFACE3D
     # and swissALTI3D products (again, EPSG:2056)
-    tile_gdf = surface3d_gdf.sjoin(
+    tile_gdf = surface3d_raster_gdf.sjoin(
         alti3d_gdf, how="inner", predicate="contains"
     ).to_crs(stac.CH_CRS)
 

figures/tiles.png

@@ -20,7 +20,7 @@ def get_dem_raster(
     *,
     region_crs: CRSType = None,
     alti3d_datetime: DatetimeLike | None = None,
-    alti3d_res: float = 2,
+    alti3d_res: float = 0.5,
     pooch_retrieve_kwargs: utils.KwargsType = None,
     rio_merge_kwargs: utils.KwargsType = None,
 ) -> None:
@@ -39,7 +39,7 @@ def get_dem_raster(
     alti3d_datetime : datetime-like, optional
         Datetime to filter swissALTI3D data, forwarded to `pystac_client.Client.search`.
         If None, the latest data for each tile is used.
-    alti3d_res : {0.5, 2}, default 2
+    alti3d_res : {0.5, 2}, default 0.5
         Resolution of the swissALTI3D data to get, can be 0.5 or 2 (meters).
     pooch_retrieve_kwargs, rio_merge_kwargs : mapping, optional
         Additional keyword arguments to respectively pass to `pooch.retrieve` and
@@ -48,10 +48,15 @@ def get_dem_raster(
     """
     # use the STAC API to get the DEM from swissALTI3D
     client = stac.SwissTopoClient(region, region_crs=region_crs)
-    alti3d_gdf = client.gdf_from_collection(
+    alti3d_gdf = client.get_collection_gdf(
         stac.SWISSALTI3D_COLLECTION_ID,
         datetime=alti3d_datetime,
     )
+    if alti3d_gdf.empty:
+        raise ValueError(
+            "Cannot compute DEM raster: no data available for the specified region and "
+            "datetime."
+        )
 
     # filter to get tiff images only
     alti3d_gdf = alti3d_gdf[alti3d_gdf["assets.href"].str.endswith(".tif")]

figures/tree-canopy.png

@@ -2,14 +2,14 @@
 
 import warnings
 from collections.abc import Iterator
-from copy import deepcopy
 
 import geopandas as gpd
+import numpy as np
 import pandas as pd
 import pystac_client
 from pyregeon import CRSType, RegionMixin, RegionType
 from pystac_client.item_search import DatetimeLike
-from shapely.geometry import shape
+from shapely import box
 
 # import all constants too
 # __all__ = ["get_latest", "SwissTopoClient"]
@@ -45,17 +45,16 @@
 # see pystac-client.readthedocs.io/en/stable/tutorials/stac-metadata-viz.html#GeoPandas
 def _items_to_gdf(items: Iterator) -> gpd.GeoDataFrame:
     """Convert a list of STAC Items into a geo-data frame."""
-    _items = []
-    for i in items:
-        _i = deepcopy(i)
-        _i["geometry"] = shape(_i["geometry"])
-        _items.append(_i)
-    gdf = gpd.GeoDataFrame(pd.json_normalize(_items))
-    for field in ["properties.datetime", "properties.created", "properties.updated"]:
-        if field in gdf:
-            gdf[field] = pd.to_datetime(gdf[field])
-    # gdf.set_index("properties.datetime", inplace=True)
-    return gdf
+    # TODO: use polars or stacrs to improve performance
+    gdf = pd.json_normalize(list(items))
+    if gdf.empty:
+        # there is no data for the specified datetime, warn and return empty gdf
+        warnings.warn(
+            "No data available for the specified datetime and collection. Returning an "
+            "empty data frame."
+        )
+        return gdf
+    return gpd.GeoDataFrame(gdf, geometry=box(*np.array(gdf["bbox"].tolist()).T))
 
 
 def _postprocess_items_gdf(items_gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
@@ -96,28 +95,56 @@ def expand_row(row):
             axis="columns",
         )
 
-    return pd.concat(
+    gdf = pd.concat(
         [expand_row(row) for _, row in items_gdf.iterrows()], ignore_index=True
     )
 
+    # set datetime columns to datetime
+    for field in ["properties.datetime", "properties.created", "properties.updated"]:
+        if field in gdf:
+            gdf[field] = pd.to_datetime(gdf[field])
+
+    return gdf
+
 
 def get_latest(
-    gdf: gpd.GeoDataFrame,
+    collection_gdf: gpd.GeoDataFrame,
     *,
-    id_col: str = "id",
+    tile_id_col: str = "id",
     datetime_col: str = "properties.datetime",
 ) -> gpd.GeoDataFrame:
-    """Get the latest item for each tile."""
+    """Get the latest item for each tile and file metadata.
+
+    Parameters
+    ----------
+    collection_gdf : geopandas.GeoDataFrame
+        Collection geo-data frame.
+    tile_id_col : str, default "id"
+        Column name for the tile ID.
+    datetime_col : str, default "properties.datetime"
+        Column name for the datetime.
+
+    Returns
+    -------
+    collection_gdf : geopandas.GeoDataFrame
+        Collection geo-data frame with the latest item for each tile and file metadata.
+    """
+    by = [
+        collection_gdf[tile_id_col].str.split("_").str[-1],
+        collection_gdf["assets.type"],
+    ]
+    if "assets.eo:gsd" in collection_gdf:
+        by.append(collection_gdf["assets.eo:gsd"])
     return (
-        gdf.sort_values(
+        collection_gdf.sort_values(
             datetime_col,
             ascending=False,
         )
-        .groupby(gdf[id_col].str.split("_").str[-1])
+        .groupby(by)
         .first()
-        .rename_axis(index={id_col: "tile_id"})
-        .reset_index()
-        .set_crs(gdf.crs)
+        .rename_axis(index={tile_id_col: "tile_id"})
+        .reset_index(drop=True)
+        .set_crs(collection_gdf.crs)
     )
 
 
@@ -162,7 +189,7 @@ def __init__(
             # keyword argument in `pystac_client.client.Search`.
             self.region = None
 
-    def gdf_from_collection(
+    def get_collection_gdf(
         self,
         collection_id: str,
         *,
@@ -182,12 +209,21 @@ def gdf_from_collection(
             Coordinate reference system (CRS) of the returned geo-data frame. If None,
             the CRS of the collection tiles will be used - note that this is not
             (necessarily) the same as the CRS of the tile data itself.
+
+        Returns
+        -------
+        collection_gdf : geopandas.GeoDataFrame
+            Geo-data frame of the collection tiles.
         """
         if dst_crs is None:
             dst_crs = self._client.get_collection(collection_id).extra_fields["crs"][0]
         search = self._client.search(
             collections=[collection_id], intersects=self.region, datetime=datetime
         )
-        return gpd.GeoDataFrame(
-            _postprocess_items_gdf(_items_to_gdf(search.items_as_dicts()))
-        ).set_crs(dst_crs)
+        gdf = _items_to_gdf(search.items_as_dicts())
+        if gdf.empty:
+            # the warning has already been raised in `_items_to_gdf`, do not raise it
+            # again, just return an empty data frame
+            # TODO: best approach to handle this? i.e., where to raise the warning?
+            return gdf
+        return gpd.GeoDataFrame(_postprocess_items_gdf(gdf)).set_crs(dst_crs)