Digitales Höhenmodell (DEM)

In diesem Notebook lernst du, wie du mit der Bibliothek cdec ein digitales Höhenmodell (Digital Elevation Model, kurz DEM) für ein bestimmtes Gebiet herunterladen, auf einer Karte darstellen und anschließend als Rasterdatensatz untersuchen kannst.

Ein digitales Höhenmodell beschreibt die Höhe der Erdoberfläche. Damit lassen sich beispielsweise Täler, Hügel, Hanglagen oder Höhenunterschiede innerhalb eines Untersuchungsgebiets erkennen.

Verwendet wird das Copernicus DEM GLO-30 mit einer räumlichen Auflösung von ungefähr 30 Metern. Es handelt sich genauer um ein digitales Oberflächenmodell: Gebäude und Vegetation können daher Teil der gemessenen Oberfläche sein.

1 Räumliche Ausdehnung festlegen

Zunächst wird festgelegt, für welches Gebiet das Höhenmodell heruntergeladen werden soll.

Die räumliche Ausdehnung wird durch eine sogenannte Bounding Box beschrieben. Sie besteht aus vier Koordinaten, die die Grenzen des Untersuchungsgebiets angeben:

west: westliche Grenze
south: südliche Grenze
east: östliche Grenze
north: nördliche Grenze

In diesem Beispiel wird ein Gebiet rund um Bremen ausgewählt.

spatial_extent = {'west': 8.721771, 'south': 53.025027, 'east': 8.904419, 'north': 53.114965}

2 Warum ist keine zeitliche Ausdehnung erforderlich?

Das Copernicus DEM ist ein fertiges, statisches Höhenmodell. Anders als bei einer Satellitenaufnahme muss deshalb kein Aufnahmezeitpunkt ausgewählt werden.

Die Funktion download_dem() verwendet automatisch den bekannten Erfassungszeitraum der zugrunde liegenden Daten von Dezember 2010 bis Januar 2015. Für gewöhnliche Anwendungen muss dieser Zeitraum nicht selbst angegeben werden.

3 Höhenmodell herunterladen

Als Nächstes wird festgelegt, wo die heruntergeladene Rasterdatei gespeichert werden soll.

Der Dateipfad wird in der Variable file_path gespeichert. Die Endung .tif kennzeichnet eine GeoTIFF-Datei. In ihr werden sowohl die Höhenwerte als auch die räumlichen Informationen gespeichert.

Anschließend wird die Bibliothek cdec importiert. Mit der Funktion download_dem() wird das Höhenmodell für das zuvor definierte Gebiet heruntergeladen.

Dabei werden nur zwei Angaben benötigt:

file_path legt den Speicherort der Datei fest.
spatial_extent bestimmt das Untersuchungsgebiet.

Standardmäßig wird die Datei als Cloud Optimized GeoTIFF (COG) gespeichert.

file_path = "data/dem.tiff"

import cdec

cdec.download_dem(
    file_path,
    spatial_extent=spatial_extent,
)

Authenticated using refresh token.

PosixPath('/Users/soeren/PycharmProjects/cdec/snippets/data_access/notebooks/data/dem.tiff')

4 Rasterdaten auf einer Karte anzeigen

Das heruntergeladene Höhenmodell kann nun auf einer interaktiven Karte dargestellt werden.

Dazu wird zunächst mit Map() eine Karte erzeugt. Mit add_cog() wird anschließend die zuvor heruntergeladene Rasterdatei als neue Ebene hinzugefügt.

Die Farbskala terrain stellt niedrige und hohe Bereiche mit unterschiedlichen Farben dar. Dadurch werden räumliche Höhenunterschiede leichter erkennbar.

from geolibre import Map

m = Map(
    center=(8.81, 53.07),
    zoom=11
)

layer_id = m.add_cog(
    file_path,
    name="Digitales Höhenmodell",
    colormap="terrain",
)

m

5 Rasterdaten mit Rasterio einlesen

Neben der Darstellung auf einer Karte können die eigentlichen Rasterdaten auch direkt mit Python untersucht werden.

Dafür wird die Bibliothek rasterio importiert. Mit rio.open() wird die heruntergeladene GeoTIFF-Datei geöffnet. Das resultierende Rasterobjekt wird in der Variable src gespeichert.

Über dieses Objekt können sowohl die Höhenwerte als auch verschiedene Metadaten des Rasterdatensatzes abgerufen werden.

import rasterio as rio

src = rio.open(file_path)

6 Eigenschaften des Rasterdatensatzes untersuchen

Bevor die eigentlichen Pixelwerte betrachtet werden, lohnt sich ein Blick auf einige grundlegende Eigenschaften des Rasterdatensatzes.

Ausgegeben werden:

src.count: Anzahl der enthaltenen Rasterbänder
src.dtypes: Datentypen der gespeicherten Pixelwerte
src.shape: Höhe und Breite des Rasters
src.crs: verwendetes Koordinatenreferenzsystem (Coordinate Reference System)

Das Höhenmodell besitzt ein Rasterband. Jeder gültige Pixel enthält einen Höhenwert in Metern.

print(f"Number of bands: {src.count}")
print(f"Data type: {src.dtypes}")
print(f"Shape: {src.shape}")
print(f"Coordinate Reference System: {src.crs}")

Number of bands: 1
Data type: ('float32',)
Shape: (325, 440)
Coordinate Reference System: EPSG:4326

7 Höhenwerte auslesen

Mit src.read(1, masked=True) wird das erste Rasterband eingelesen. Der Parameter masked=True sorgt dafür, dass Pixel ohne gültigen Messwert bei späteren Berechnungen nicht berücksichtigt werden.

Das Ergebnis ist ein zweidimensionales Array. Jede Position enthält den Höhenwert des zugehörigen Pixels.

elevation = src.read(1, masked=True)
elevation

masked_array(
  data=[[0.5, 0.5, 0.5, ..., 7.164416790008545, 5.125278949737549,
         4.158882141113281],
        [3.188593864440918, 0.5, 0.5, ..., 4.772930145263672,
         4.257691383361816, 3.544872999191284],
        [3.9439425468444824, 3.7329113483428955, 3.13291072845459, ...,
         4.901821613311768, 4.782098293304443, 4.58300256729126],
        ...,
        [6.815716743469238, 6.867824077606201, 6.7665300369262695, ...,
         4.547393798828125, 4.910810947418213, 4.887389659881592],
        [7.122479438781738, 7.132204532623291, 6.9344000816345215, ...,
         4.787200927734375, 4.7809882164001465, 4.893062114715576],
        [7.009529113769531, 7.0337748527526855, 6.890302658081055, ...,
         4.669248104095459, 5.100966930389404, 4.800761699676514]],
  mask=[[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],
  fill_value=nan,
  dtype=float32)

8 Höhenunterschiede untersuchen

Aus den gültigen Pixelwerten lassen sich die minimale, maximale und mittlere Höhe des Untersuchungsgebiets berechnen.

Diese Kennwerte geben einen ersten Eindruck vom Relief. In einem kleinen, flachen Gebiet können bereits geringe Höhenunterschiede für Entwässerung, Hochwasserrisiken oder die Planung von Infrastruktur relevant sein.

print(f"Minimum elevation: {elevation.min():.1f} m")
print(f"Maximum elevation: {elevation.max():.1f} m")
print(f"Mean elevation: {elevation.mean():.1f} m")

Minimum elevation: -13.9 m
Maximum elevation: 43.3 m
Mean elevation: 6.5 m

9 Rasterdaten visualisieren

Zum Abschluss werden die Höhenwerte direkt im Notebook dargestellt.

Die Funktion imshow() aus plotly.express erzeugt eine interaktive Darstellung. Die Farbskala terrain hilft dabei, niedrige und hohe Bereiche voneinander zu unterscheiden; die Farblegende zeigt die Höhe in Metern.

import plotly.express as px

fig = px.imshow(
    elevation.filled(float("nan")),
    color_continuous_scale="gray",
    labels={"color": "Höhe (m)"},
    height=700,
)

fig.show()

Wiederverwendung

Dieses Material steht unter der Lizenz CC BY-SA 4.0. (Lizenz Anzeigen)