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.
# First, let's import the needed libraries.
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
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.
import rasterio
## 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.
srtm_germany.count
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
).
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
).
srtm_germany.bounds
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.
srtm_germany.transform
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.
srtm_germany.crs
CRS.from_epsg(4326)
You may access all meta data at once using .meta
# 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.
print(f"The file size is {os.path.getsize('../data/srtm_germany_dtm.tif')/1e6} MB")
The file size is 285.254776 MB
# 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"
<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.
# 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")
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()