feat: added deps, tests, fix tox/gdal/pdal setup with gh actions

Author: martibosch

.github/workflows/tests.yml

@@ -32,14 +32,21 @@ jobs:
           environment-name: test-env
           create-args: >-
             python=${{ matrix.python-version }}
+            gdal>=3.3
+            pdal
             pip
+            python-pdal
+            tox
+            tox-current-env
+            tox-gh
 
-      - name: install dependencies
-        run: pip install tox tox-gh-actions
+      - name: install current package with dev/doc/test extras
+        run: pip install ".[dev,doc,test]"
 
       - name: test with tox
-        run: tox
+        run: tox --current-env
         env:
+          TOX_GH_MAJOR_MINOR: ${{ matrix.python-version }}
           CONDA_EXE: mamba
 
       - name: upload coverage reports to Codecov

docs/environment.yml

@@ -6,6 +6,7 @@ dependencies:
   - gdal=3.10.2
   - myst-parser=4.0.1
   - sphinx=8.2.3
+  - pdal=2.8.4
   - pip=25.0.1
   - pydata-sphinx-theme=0.16.1
   - setuptools=75.8.2

pyproject.toml

@@ -19,6 +19,18 @@ classifiers = [
   "Programming Language :: Python :: 3.12",
   "Programming Language :: Python :: 3.13"
 ]
+dependencies = [
+  "geopandas>=1.0.0",
+  "numpy>=2.0.0",
+  "osmnx>=2.0.0",
+  "pooch",
+  "pyregeon",
+  "pystac-client",
+  "rasterio",
+  "rasterstats",
+  "shapely>=2.0.0",
+  "tqdm"
+]
 
 [project.license]
 text = "GPL-3.0"
@@ -40,6 +52,9 @@ doc = [
   "sphinx",
   "sphinxemoji"
 ]
+pdal = [
+  "pdal"
+]
 test = [
   "coverage[toml]",
   "pytest",
@@ -122,9 +137,9 @@ requires = [
 
 [tool.tox.env.lint]
 commands = [
-  "python -m build",
-  "sphinx-build docs docs/_build",
-  "twine check dist/*"
+  ["python", "-m", "build"],
+  ["sphinx-build", "docs", "docs/_build"],
+  ["twine", "check", "dist/*"]
 ]
 extras = [
   "dev",
@@ -141,9 +156,6 @@ whitelist_externals = [
 commands = [
   ["pytest", "-s", "--cov=swisstopopy", "--cov-report=xml", "tests"]
 ]
-conda_deps = [
-  "gdal>=3.3"
-]
 extras = [
   "test"
 ]

swisstopopy/buildings.py

@@ -112,7 +112,7 @@ def get_bldg_gdf(
             region_gser.to_crs(ox.settings.default_crs).iloc[0],
             tags=OSMNX_TAGS,
         )
-        .to_crs(stac.SWISSALTI3D_CRS)
+        .to_crs(stac.CH_CRS)
         .drop("node")
     )
 
@@ -157,7 +157,7 @@ def get_bldg_gdf(
     # and swissALTI3D products (again, EPSG:2056)
     tile_gdf = surface3d_gdf.sjoin(
         alti3d_gdf, how="inner", predicate="contains"
-    ).to_crs(stac.SWISSALTI3D_CRS)
+    ).to_crs(stac.CH_CRS)
 
     # we could do a data frame apply approach returning a series of of building heights
     # that correspond to a single zonal statistic (e.g., "mean"). However, we use

swisstopopy/dem.py

@@ -4,17 +4,15 @@
 """
 
 import pooch
-from osgeo import gdal
 from pyregeon import CRSType, RegionType
 from pystac_client.item_search import DatetimeLike
+from rasterio import merge
 from tqdm import tqdm
 
-from swisstopopy import stac, utils
+from swisstopopy import settings, stac, utils
 
 __all__ = ["get_dem_raster"]
 
-DST_OPTIONS = ["TILED:YES"]
-
 
 def get_dem_raster(
     region: RegionType,
@@ -24,7 +22,7 @@ def get_dem_raster(
     alti3d_datetime: DatetimeLike | None = None,
     alti3d_res: float = 2,
     pooch_retrieve_kwargs: utils.KwargsType = None,
-    gdal_warp_kwargs: utils.KwargsType = None,
+    rio_merge_kwargs: utils.KwargsType = None,
 ) -> None:
     """Get digital elevation model (DEM) raster.
 
@@ -43,9 +41,10 @@ def get_dem_raster(
         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).
-    pooch_retrieve_kwargs, gdal_warp_kwargs : mapping, optional
+    pooch_retrieve_kwargs, rio_merge_kwargs : mapping, optional
         Additional keyword arguments to respectively pass to `pooch.retrieve` and
-        `gdal.Warp`.
+        `rasterio.merge.merge`.  If the latter is None, the default values from
+        `settings.RIO_MERGE_DST_KWARGS` are used.
     """
     # use the STAC API to get the DEM from swissALTI3D
     client = stac.SwissTopoClient(region=region, region_crs=region_crs)
@@ -65,16 +64,21 @@ def get_dem_raster(
     if pooch_retrieve_kwargs is None:
         pooch_retrieve_kwargs = {}
 
-    if gdal_warp_kwargs is None:
-        _gdal_warp_kwargs = {}
+    if rio_merge_kwargs is None:
+        _rio_merge_kwargs = {}
     else:
-        _gdal_warp_kwargs = gdal_warp_kwargs.copy()
-    _gdal_warp_kwargs.update(creationOptions=DST_OPTIONS)
+        _rio_merge_kwargs = rio_merge_kwargs.copy()
+    _rio_merge_kwargs.update(dst_kwds=settings.RIO_MERGE_DST_KWARGS)
 
     img_filepaths = []
     for url in tqdm(
         alti3d_gdf["assets.href"],
     ):
         img_filepath = pooch.retrieve(url, known_hash=None, **pooch_retrieve_kwargs)
         img_filepaths.append(img_filepath)
-    _ = gdal.Warp(dst_filepath, img_filepaths, format="GTiff", **_gdal_warp_kwargs)
+    # merge the images into the final raster
+    merge.merge(
+        img_filepaths,
+        dst_path=dst_filepath,
+        **_rio_merge_kwargs,
+    )

swisstopopy/settings.py

@@ -0,0 +1,3 @@
+"""Settings."""
+
+RIO_MERGE_DST_KWARGS = {"tiled": True, "blockxsize": 512, "blockysize": 512}

swisstopopy/stac.py

@@ -154,8 +154,8 @@ def __init__(
             # rather than inheriting from `RegionMixin`, we just use the
             # `_process_region_arg` static method
             self.region = (
-                RegionMixin._process_region_arg(region, region_crs=region_crs)
-                .to_crs(CLIENT_CRS)
+                RegionMixin._process_region_arg(region, crs=region_crs)
+                .to_crs(CLIENT_CRS)["geometry"]
                 .iloc[0]
             )
         else:

swisstopopy/tree_canopy.py

@@ -10,17 +10,15 @@
 import pdal
 import pooch
 import rasterio as rio
-from osgeo import gdal
 from pyregeon import CRSType, RegionType
 from pystac_client.item_search import DatetimeLike
+from rasterio import merge
 from tqdm import tqdm
 
-from swisstopopy import stac, utils
+from swisstopopy import settings, stac, utils
 
 __all__ = ["get_tree_canopy_raster"]
 
-DST_OPTIONS = ["TILED:YES"]
-
 
 def rasterize_lidar(
     lidar_filepath: utils.PathType,
@@ -61,10 +59,10 @@ def get_tree_canopy_raster(
     dst_tree_val: int = 1,
     dst_nodata: int = 0,
     dst_dtype: npt.DTypeLike = "uint32",
-    lidar_tree_values: Sequence[int] | None = None,
+    lidar_tree_values: int | Sequence[int] | None = 3,
     rasterize_lidar_kwargs: utils.KwargsType = None,
     pooch_retrieve_kwargs: utils.KwargsType = None,
-    gdal_warp_kwargs: utils.KwargsType = None,
+    rio_merge_kwargs: utils.KwargsType = None,
 ) -> None:
     """Get tree canopy raster.
 
@@ -96,9 +94,11 @@ def get_tree_canopy_raster(
     lidar_tree_values : int or sequence of int, default 3.
         LiDAR classification values to use for tree canopy. If None, defaults to
         the "Vegetation" class value of swissSURFACE3D, i.e., 3.
-    rasterize_lidar_kwargs, pooch_retrieve_kwargs, gdal_warp_kwargs : mapping, optional
+    rasterize_lidar_kwargs, pooch_retrieve_kwargs, rio_merge_kwargs : mapping, optional
         Additional keyword arguments to respectively pass to
-        `swisstopopy.tree_canopy.rasterize_lidar`, `pooch.retrieve` and `gdal.Warp`.
+        `swisstopopy.tree_canopy.rasterize_lidar`, `pooch.retrieve` and
+        `rasterio.merge.merge`. If the latter is None, the default values from
+        `settings.RIO_MERGE_DST_KWARGS` are used.
     """
     # use the STAC API to get the tree canopy from swissSURFACE3D
     # TODO: dry with `dem.get_dem_raster`?
@@ -126,18 +126,19 @@ def get_tree_canopy_raster(
         output_type="count",
         data_type="uint32",
         nodata=dst_nodata,
-        default_srs=stac.SWISSALTI3D_CRS,
+        default_srs=stac.CH_CRS,
     )
     if pooch_retrieve_kwargs is None:
         pooch_retrieve_kwargs = {}
 
     if isinstance(lidar_tree_values, int):
         lidar_tree_values = [lidar_tree_values]
-    if gdal_warp_kwargs is None:
-        _gdal_warp_kwargs = {}
+
+    if rio_merge_kwargs is None:
+        _rio_merge_kwargs = {}
     else:
-        _gdal_warp_kwargs = gdal_warp_kwargs.copy()
-    _gdal_warp_kwargs.update(creationOptions=DST_OPTIONS)
+        _rio_merge_kwargs = rio_merge_kwargs.copy()
+    _rio_merge_kwargs.update(dst_kwds=settings.RIO_MERGE_DST_KWARGS)
 
     img_filepaths = []
     with tempfile.TemporaryDirectory() as tmp_dir:
@@ -165,6 +166,13 @@ def get_tree_canopy_raster(
                     tmp_dir,
                     f"{path.splitext(path.basename(img_filepath))[0]}-counts.tif",
                 )
+                _ = rasterize_lidar(
+                    las_filepath,
+                    counts_filepath,
+                    lidar_tree_values,
+                    **_rasterize_lidar_kwargs,
+                )
+
                 try:
                     _ = rasterize_lidar(
                         las_filepath,
@@ -193,5 +201,9 @@ def get_tree_canopy_raster(
             # add path to list
             img_filepaths.append(img_filepath)
 
-        # creationOptions=dst_options
-        _ = gdal.Warp(dst_filepath, img_filepaths, format="GTiff", **_gdal_warp_kwargs)
+        # merge tiles into the final raster
+        merge.merge(
+            img_filepaths,
+            dst_path=dst_filepath,
+            **_rio_merge_kwargs,
+        )

tests/test_swisstopopy.py

@@ -21,16 +21,19 @@ def test_region(self):
         # test without region (all collection items)
         client = swisstopopy.SwissTopoClient()
 
-        # test init all collections
-        for collection_id in self.collection_ids:
-            gdf = client.gdf_from_collection(collection_id)
-            # test that we get a non empty geo-data frame
-            self.assertIsInstance(gdf, gpd.GeoDataFrame)
-            self.assertFalse(gdf.empty)
+        # since this is slow, test init one collection only
+        collection_id = self.collection_ids[0]
+        gdf = client.gdf_from_collection(collection_id)
 
-        # now only test one collection (the last one)
-        region_client = swisstopopy.SwissTopoClient(self.region)
+        # test with region
+        region_client = swisstopopy.SwissTopoClient(region=self.nominatim_query)
         # test that there are at most as many items as when not filtering spatially
         self.assertLessEqual(
             len(region_client.gdf_from_collection(collection_id).index), len(gdf.index)
         )
+        # test init all collections
+        for collection_id in self.collection_ids:
+            gdf = region_client.gdf_from_collection(collection_id)
+            # test that we get a non empty geo-data frame
+            self.assertIsInstance(gdf, gpd.GeoDataFrame)
+            self.assertFalse(gdf.empty)