swisstopopy 0.2.0

swisstopo geospatial Python utilities

swisstopopy

swisstopo geospatial Python utilities.

Features

STAC API utilities

Easily filter swisstopo STAC collections based on geospatial extents, dates, file extensions or data resolutions:

import contextily as cx
import swisstopopy

region = "EPFL"
client = swisstopopy.SwissTopoClient(region)

alti3d_gdf = client.get_collection_gdf(
    swisstopopy.SWISSALTI3D_COLLECTION_ID,
)
ax = alti3d_gdf.plot(alpha=0.1)
cx.add_basemap(ax, crs=alti3d_gdf.crs)

Filter to get the latest data for each tile:

latest_alti3d_gdf = swisstopopy.get_latest(alti3d_gdf)
latest_alti3d_gdf.head()

	id	collection	...	geometry
0	swissalti3d_2021_2532-1151	ch.swisstopo.swissalti3d	...	POLYGON ((6.56572 46.50684, 6.56572 46.51594, ...
1	swissalti3d_2021_2532-1151	ch.swisstopo.swissalti3d	...	POLYGON ((6.56572 46.50684, 6.56572 46.51594, ...
2	swissalti3d_2021_2532-1151	ch.swisstopo.swissalti3d	...	POLYGON ((6.56572 46.50684, 6.56572 46.51594, ...
3	swissalti3d_2021_2532-1151	ch.swisstopo.swissalti3d	...	POLYGON ((6.56572 46.50684, 6.56572 46.51594, ...
4	swissalti3d_2021_2532-1152	ch.swisstopo.swissalti3d	...	POLYGON ((6.56558 46.51584, 6.56558 46.52493, ...

or filter by other metadata attributes such as ground resolution and/or file extensions:

alti3d_gdf[
    (alti3d_gdf["assets.eo:gsd"] == 0.5)
    & alti3d_gdf["assets.href"].str.endswith(".tif")
]

	id	collection	...	geometry
0	swissalti3d_2019_2532-1151	ch.swisstopo.swissalti3d	...	POLYGON ((6.56572 46.50684, 6.56572 46.51594, ...
4	swissalti3d_2019_2532-1152	ch.swisstopo.swissalti3d	...	POLYGON ((6.56558 46.51584, 6.56558 46.52493, ...
8	swissalti3d_2019_2533-1152	ch.swisstopo.swissalti3d	...	POLYGON ((6.57861 46.51594, 6.57861 46.52503, ...
12	swissalti3d_2021_2532-1151	ch.swisstopo.swissalti3d	...	POLYGON ((6.56572 46.50684, 6.56572 46.51594, ...
16	swissalti3d_2021_2532-1152	ch.swisstopo.swissalti3d	...	POLYGON ((6.56558 46.51584, 6.56558 46.52493, ...
20	swissalti3d_2021_2533-1152	ch.swisstopo.swissalti3d	...	POLYGON ((6.57861 46.51594, 6.57861 46.52503, ...

STAC data processing

Automated generation of geospatial datasets: building footprints with estimated heights, DEM and tree canopy. For example, a tree canopy raster for any given part of Switzerland can be obtained as in:

import rasterio as rio
from rasterio import plot

dst_filepath = "tree-canopy.tif"
swisstopopy.get_tree_canopy_raster(region, dst_filepath)

with rio.open(dst_filepath) as src:
    plot.show(src)

See the overview notebook and the API documentation for more details on the geospatial dataset generation functions.

Installation

Like many other geospatial Python packages, swisstopopy requires many base C libraries that cannot be installed with pip. Accordingly, the best way to install swisstopopy is to use conda/mamba, i.e., in a given conda environment, run:

# or mamba install -c conda-forge geopandas
conda install -c conda-forge geopandas

Within the same conda environment, you can then install swisstopopy using pip:

pip install swisstopopy

Note that the get_tree_canopy_raster requires PDAL and its Python bindings, which are not installed by default with swisstopopy. Like with geopandas, the easiest way to install such requirements is using conda/mamba, e.g.: conda install -c conda-forge python-pdal.

Notes

The SwissTopoClient class can be used to process any collection of the swisstopo STAC API, and basic features succh as geospatial and datetime filtering should work out of the box. However, filtering based on further metadata such as the resolution is only fully supported for the following collections:

• "ch.swisstopo.swissalti3d", namely swissALTI3D
• "ch.swisstopo.swissimage-dop10", namely SWISSIMAGE 10 cm
• "ch.swisstopo.swisssurface3d", namely swissSURFACE3D
• "ch.swisstopo.swisssurface3d-raster", namely swissSURFACE3D Raster.

Acknowledgements

• This package was created with the martibosch/cookiecutter-geopy-package project template.