The SRTM digital elevation model (DEM) for Germany was retrieved from OPENDEM on June 25, 2017. Please note that the zipped file is about 86.4 MB in size.

In [2]:
# First, let's import the needed libraries.
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
In [3]:
import requests, zipfile, io, os

url = "http://www.userpage.fu-berlin.de/~soga/300/30100_data_sets/spatial/srtm_germany_dtm.zip"
r = requests.get(url)
z = zipfile.ZipFile(io.BytesIO(r.content))
z.extractall("../data")

To load the raster data set we apply the rasterio package providing open() functionality. The open() function takes a path string or path-like object pointing to a file of a supported raster format and returns an opened dataset object. Rasterio will open it using the proper GDAL format driver. Dataset objects have some of the same attributes as Python file objects.

In [4]:
import rasterio
In [5]:
## open file
srtm_germany = rasterio.open("../data/srtm_germany_dtm.tif", driver="GTiff")

Properties of the raster data stored in the example GeoTIFF can be accessed through attributes of the opened dataset object. Dataset objects have bands, which can be counted using .count. Our SRTM data has a band count of 1.

In [6]:
srtm_germany.count
Out[6]:
1

A dataset band is an array of values containing the distribution of a variable (here elevation) in rows and columns. Check the dimensions of the array of DN values (.width and .height).

In [7]:
print(
    f"Number of columns =  {srtm_germany.width}, Number of rows = {srtm_germany.height}"
)

Number of columns =  13201, Number of rows = 10801

The georeferencing of the dataset is stored in the meta data. Every pixels of a dataset is contained within a spatial bounding box (.bounds).

In [8]:
srtm_germany.bounds
Out[8]:
BoundingBox(left=4.999583333333334, bottom=46.99958333333337, right=16.000416666666624, top=56.000416666666666)

The value of DatasetReader.bounds attribute is derived from a more fundamental attribute: the dataset’s geospatial transform, which describes the origin coordinates of the raster.

In [9]:
srtm_germany.transform
Out[9]:
Affine(0.0008333333333333301, 0.0, 4.999583333333334,
       0.0, -0.0008333333333333301, 56.000416666666666)

These coordinate values are relative to the origin of the dataset’s coordinate reference system (CRS). EPSG:4326 identifies a particular coordinate reference system.

In [10]:
srtm_germany.crs
Out[10]:
CRS.from_epsg(4326)

You may access all meta data at once using .meta

In [11]:
# metadata functions from rasterio
print(srtm_germany.meta)

{'driver': 'GTiff', 'dtype': 'int16', 'nodata': None, 'width': 13201, 'height': 10801, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.0008333333333333301, 0.0, 4.999583333333334,
       0.0, -0.0008333333333333301, 56.000416666666666)}

When working with large data sets it is recommendable to check the size of the file in advance. The processing time increases with the size of the data set.

In [12]:
print(f"The file size is {os.path.getsize('../data/srtm_germany_dtm.tif')/1e6} MB")

The file size is 285.254776 MB
In [13]:
# plot with rasterio.plot, which provides matplotlib functionality
from rasterio.plot import show

show(srtm_germany, title="Digital Elevation Model - SRTM", cmap="RdYlGn_r")  # "terrain"
Out[13]:
<AxesSubplot: title={'center': 'Digital Elevation Model - SRTM'}>

An useful function is the mask() function from from rasterio.mask, which creates a new raster object that has the same values as the input raster object, except for the masked cells, which are set to 0 by default. We apply the mask() function to mask out the area outside the borders of Germany.

For the shapefiles of Germany we rely on shapefiles provided by the German Federal Agency for Cartography and Geodesy (under this licence). You may directly download the shapefiles here. For the purpose of this tutorial, we also provided the data here.

In [14]:
# Retrieve Federal States

import zipfile

url = "https://daten.gdz.bkg.bund.de/produkte/vg/vg5000_0101/aktuell/vg5000_01-01.utm32s.shape.ebenen.zip"

r = requests.get(url)
z = zipfile.ZipFile(io.BytesIO(r.content))
z.extractall(path="../data")
boundary_germany = gpd.read_file(
    "../data//vg5000_01-01.utm32s.shape.ebenen/vg5000_ebenen_0101/VG5000_STA.shp"
)

boundary_germany = boundary_germany.to_crs("epsg:4326")
In [15]:
fig, ax = plt.subplots(figsize=(5, 15))

srtm_plot = show(
    srtm_germany,
    ax=ax,
    cmap="RdYlGn_r",
    title="German Boundary and \n Digital Elevation Model - SRTM\n",
)
boundary_germany.plot(ax=ax, facecolor="none", edgecolor="lightgrey")

## add colorbar
dummy_colorbar = srtm_plot.get_images()[0]
fig.colorbar(dummy_colorbar, ax=ax, fraction=0.046, pad=0.04, label="Elevation [m]")


plt.show()