Overview of swisstopopy#

[ ]:

import contextily as cx
import matplotlib.pyplot as plt
import rasterio as rio
from rasterio import plot

import swisstopopy

STAC data processing into raster products#

Besides the simplified interface to navigate swisstopo STAC collections, swisstopopy includes a few functions to automatically obtain useful raster products.

Building footprints with estimated heights#

We can get building footprints from the OpenStreetMap (OSM) (using osmnx) and then infer building heights using the difference between swissSURFACE3D Raster, a digital surface model (DSM) and swissALTI3D, a digital elevation model (DEM). This is all managed through the get_bldg_gdf function, which will return a geo-data frame of building footprints, OSM attributes and inferred building height (in the “height” column):

[ ]:

bldg_gdf = swisstopopy.get_bldg_gdf(region)

ax = bldg_gdf.plot("height", legend=True)
cx.add_basemap(ax, crs=bldg_gdf.crs)

 33%|█████████████████                                  | 1/3 [00:00<00:00,  3.77it/s]/home/martibosch/miniforge3/envs/swisstopopy/lib/python3.13/site-packages/rasterstats/io.py:335: NodataWarning: Setting nodata to -999; specify nodata explicitly
  warnings.warn(
100%|███████████████████████████████████████████████████| 3/3 [00:01<00:00,  2.73it/s]
_images/overview_13_1.png

By default, the latest available swissSURFACE3D Raster and swissALTI3D data will be used, but this can be changed using the item_datetime keyword argument. Similarly, the resolution of the data (by default 0.5 m) can be changed using the item_res keyword argument (note that the resolution should be the same for both collections to compute the raster difference).

Digital elevation model (DEM)#

We can use the get_dem_raster to easily get a DEM of any part of Switzerland:

[ ]:

dst_filepath = "dem.tif"
swisstopopy.get_dem_raster(region, dst_filepath)

fig, ax = plt.subplots()
with rio.open(dst_filepath) as src:
    retted = plot.show(src, ax=ax)
fig.colorbar(retted.get_images()[0], ax=ax)

100%|█████████████████████████████████████████████████| 3/3 [00:00<00:00, 2762.44it/s]

<matplotlib.colorbar.Colorbar at 0x730d3008a850>
_images/overview_15_2.png

By default, the latest data for each tile will be used at a resolution of 0.5 m, but this can be changed via the alti3d_datetime and alti3d_res keyword arguments respectively.

Tree canopy raster#

A tree canopy raster can be obtained by filtering high vegetation points from the swissSURFACE3D dataset and rasterizing it to a user-defined resolution. This can be done using the get_tree_canopy_raster function:

[ ]:

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

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

  0%|                                                           | 0/3 [00:00<?, ?it/s]Downloading data from 'https://data.geo.admin.ch/ch.swisstopo.swisssurface3d/swisssurface3d_2019_2532-1151/swisssurface3d_2019_2532-1151_2056_5728.las.zip' to file '/tmp/tmpwk2gcgff/1c03c446a7b6f85fe9d68bdc70c90916-swisssurface3d_2019_2532-1151_2056_5728.las.zip'.
SHA256 hash of downloaded file: c8ec0d2ae2a32c3ddf9f41711275078de7b229319de017fa6ff6797b1dd97ced
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
Unzipping contents of '/tmp/tmpwk2gcgff/1c03c446a7b6f85fe9d68bdc70c90916-swisssurface3d_2019_2532-1151_2056_5728.las.zip' to '/tmp/tmpwk2gcgff/1c03c446a7b6f85fe9d68bdc70c90916-swisssurface3d_2019_2532-1151_2056_5728.las.zip.unzip'
 33%|█████████████████                                  | 1/3 [00:12<00:25, 12.87s/it]Downloading data from 'https://data.geo.admin.ch/ch.swisstopo.swisssurface3d/swisssurface3d_2019_2532-1152/swisssurface3d_2019_2532-1152_2056_5728.las.zip' to file '/tmp/tmpwk2gcgff/86f648067f304b5e665d67014c7593f3-swisssurface3d_2019_2532-1152_2056_5728.las.zip'.
SHA256 hash of downloaded file: dc0490ca940e5e1c81a3e7cbc21054ebbe216be1c64e91b9af5e325b73df1921
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
Unzipping contents of '/tmp/tmpwk2gcgff/86f648067f304b5e665d67014c7593f3-swisssurface3d_2019_2532-1152_2056_5728.las.zip' to '/tmp/tmpwk2gcgff/86f648067f304b5e665d67014c7593f3-swisssurface3d_2019_2532-1152_2056_5728.las.zip.unzip'
 67%|██████████████████████████████████                 | 2/3 [00:26<00:13, 13.27s/it]Downloading data from 'https://data.geo.admin.ch/ch.swisstopo.swisssurface3d/swisssurface3d_2019_2533-1152/swisssurface3d_2019_2533-1152_2056_5728.las.zip' to file '/tmp/tmpwk2gcgff/b6c098c0bf5104e790d8733cb8067080-swisssurface3d_2019_2533-1152_2056_5728.las.zip'.
SHA256 hash of downloaded file: a29847b7eae9c59783222a7236957d4dd27cc156b390e0f51c8381a78b2ba21b
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
Unzipping contents of '/tmp/tmpwk2gcgff/b6c098c0bf5104e790d8733cb8067080-swisssurface3d_2019_2533-1152_2056_5728.las.zip' to '/tmp/tmpwk2gcgff/b6c098c0bf5104e790d8733cb8067080-swisssurface3d_2019_2533-1152_2056_5728.las.zip.unzip'
100%|███████████████████████████████████████████████████| 3/3 [00:42<00:00, 14.23s/it]
_images/overview_17_1.png

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