Digital elevation models (DEM) are representation of surface topography measured as elevation. There are processed a 3D raster representaion of a land surface, typically build as a 2D array, each cell having an associated elevation. A DEM allows to derive quantities such as slope,aspect and hillshade, which we will explore in this chapter.

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
import rasterio
from rasterio.plot import plotting_extent

In this section we will use three packages in particular to handle SRTM data:

  • EarthPy
  • elevation
  • richDEM

We will have to load these packages as well:

In [3]:
import earthpy as et
import earthpy.plot as ep
import earthpy.spatial as es

import elevation
import richdem as rd

Here we will work with a DEM covering Germany. 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 and can be accessed here.

In [4]:
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")

You may want to check for the file size prior to loading the file. Processing large raster images can result in very long computing times.

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

The file size is 285.254776 MB

richDEM provides the LoadGDAL function for loading DEM data.

In [6]:
srtm_germ = rd.LoadGDAL(filepath, no_data=np.nan)
In [7]:
plt.imshow(srtm_germ, interpolation='none', cmap ="terrain", vmin = 0, vmax = 2000)
plt.colorbar(label = "reflectance *1000")

plt.show()

Check meta data of srtm_germ, since it is slightly different compared to loading it with rasterio.

In [8]:
srtm_germ.projection
Out[8]:
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
In [9]:
srtm_germ.geotransform
Out[9]:
(4.999583333333334,
 0.0008333333333333301,
 0.0,
 56.000416666666666,
 0.0,
 -0.0008333333333333301)

Note: We are also able to directly download elevation data from from the NASA SRTM mission through the elevation package. You can checkout this tutorial from Oakley & Joseph 2020 if you want to learn how.

Crop the SRTM of Germany

As we saw in the previous section, we are able to crop the SRTM raster data with the help of the EarthPy library. Consequently, we will crop the rasterdata set prior to working with it as follows.

First, we will load the required shapefile data, convert to desired crs ("epsg:32633").

In [10]:
# Retrieve Federal States

import zipfile, requests, io

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")
germany_admin = gpd.read_file(
    "../data/vg5000_01-01.utm32s.shape.ebenen/vg5000_ebenen_0101/VG5000_LAN.shp"
)
In [11]:
srtm_germ.projection
Out[11]:
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
In [12]:
# load shapefile of Berlin
berlin = germany_admin[germany_admin["GEN"] == "Berlin"]
berlin = berlin.to_crs("epsg:4326")

Now, we are ready to crop the srtm to the size of berlin.

In [13]:
with rasterio.open(filepath) as src:
    srtm_berlin, srtm_berlin_meta = es.crop_image(src, berlin)

extent = plotting_extent(srtm_berlin[0], srtm_berlin_meta["transform"])
In [14]:
# mask the nodata and plot the newly cropped raster layer
srtm_berlin_ma = np.ma.masked_equal(srtm_berlin[0], -9999.0) 


### plot ###
fig, ax = plt.subplots(figsize=(8, 8))

berlin.boundary.plot(ax=ax, color="red", zorder=10)

ep.plot_bands(srtm_berlin_ma, title = "DEM SRTM Berlin", 
              cmap='terrain', vmin = 0, vmax = 120,
              cbar=True, extent = extent,
              ax = ax,)

plt.show()

Computing slope and aspect¶

The rd.TerrainAttribute function is very easy to use, when computing the slope and aspect for each pixel. The function provides the atrrib argument which can be set to 'slope_riserun' and 'aspect', respectively. To visualization use the rd.rdShow function.

Note: The richDEM documentation is a helpful ressource for various DEM derivates and calculations.

Slope

Slope is the steepness or the degree of incline of a surface and an be expressed in degrees or as a percentage. The slope of a cell is calculated by using a central difference estimation of a surface fitted to the cell and its neighbours. The slope chosen is the maximum of this surface (Horn 1981).

In [15]:
slope = rd.TerrainAttribute(srtm_germ, attrib='slope_riserun')
rd.rdShow(slope, axes=False, cmap='RdYlGn_r', figsize=(8, 5.5))
plt.show()

Get the just the slope for the SRTM of Berlin. Therefore, we will transform the data loaded with rasterio into a format richDEM can process. We will have to set the right projection and geotransform of the srtm_berlin prior to calculating the slope.

In [16]:
rd_srtm_berlin = rd.rdarray(srtm_berlin[0], no_data=-9999)
rd_srtm_berlin.projection = srtm_germ.projection
rd_srtm_berlin.geotransform = srtm_germ.geotransform
In [17]:
# calculate the slope
slope = rd.TerrainAttribute(rd_srtm_berlin, attrib='slope_riserun')
# plot
rd.rdShow(slope, axes=False, cmap='RdYlGn_r', figsize=(8, 5.5))
plt.show()

Aspect

The Aspect is the orientation of slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, and so on. The aspect is calculated as the direction of the maximum slope of the cell (Horn 1981).

In [18]:
## aspect
aspect = rd.TerrainAttribute(rd_srtm_berlin, attrib='aspect')
rd.rdShow(aspect, axes=False, cmap='gnuplot', figsize=(8, 5.5))
plt.show()

Computing the hillshade¶

A hillshade is a 3D visualization of a surface and is usually rendered in greyscale. Hillshades are used to make a map look more appealing. For computing the hillshade, we will use the EarthPy library which provides the hillshade() function within the earthpy.spatial module.

The scope of this function is explained very nicely in this tutorial.

In [19]:
hillshade = es.hillshade(srtm_berlin_ma)

ep.plot_bands(
    hillshade,
    cbar=True,
    title="Hillshade made from DEM Berlin",
    figsize=(10, 6),
)
plt.show()

We can also overlay the hillshade on top of the DEM for better visualization.

In [20]:
fig, ax = plt.subplots(figsize=(10, 6))
ep.plot_bands(
    srtm_berlin_ma,
    ax=ax,
    cmap="terrain",
    title="SRTM DEM Berlin\n overlayed with hillshade",
)
ax.imshow(hillshade, cmap="Greys", alpha=0.5)
plt.show()

Computing the Flow Accumulation¶

Each cell in a DEM can be modeled as generating a certain amount of flow. When computing the flow accumulation, we cumulatively count the number of pixels that naturally drain into outlets. Flow accumulation can be used to find the drainage pattern of a terrain.

Flow accumulation calculations are always based on a so-called flow metric. A flow metric is a rule that divides the flow passing through a cell among one or more of its neighboring cells. There are several ways to do this. The methods vary in accuracy, complexity and thus computation time. If you need more metrics and methods than provided here, check out the documentation of richDEM, since it gives a nice overview about all implemented flow metrics.

Here, we will only concetrate on the D8 metric after O'Callaghan & Mark 1984, which is sufficient for the scope of this introduction. The D8 metric assigns the flow from a target cell to one of its eight neighboring cells. The chosen neighbour is the one with the steepest slope. If there is no steepest neighbour, no flow direction is assigned (O'Callaghan & Mark 1984).

In [21]:
#Get flow accumulation with no explicit weighting, default = 1
accum_d8 = rd.FlowAccumulation(rd_srtm_berlin, method='D8')
rd.rdShow(accum_d8, axes=False, cmap='jet', figsize=(8, 5.5))
plt.show()

If you want to see the flow pattern for a specific area, please note, that we would have to use a smaller area of interest. Since it is very difficult to see specific flow patterns for larger areas. Please note that there are also other flow metrics, such as the $D_{4}$ (O'Callaghan & Mark 1984) or $D_{\infty}$ (Tarboton 1997). If you want to learn more about flow accumulation and metrics in Python, check out the documentation of richDEM.

Citation

The E-Learning project SOGA-Py was developed at the Department of Earth Sciences by Annette Rudolph, Joachim Krois and Kai Hartmann. You can reach us via mail by soga[at]zedat.fu-berlin.de.

Creative Commons License
You may use this project freely under the Creative Commons Attribution-ShareAlike 4.0 International License.

Please cite as follow: Rudolph, A., Krois, J., Hartmann, K. (2023): Statistics and Geodata Analysis using Python (SOGA-Py). Department of Earth Sciences, Freie Universitaet Berlin.