Die besten (Python-)Tools für die Fernerkundung

Emilius Richter

2. August 2022

Satellite that is used for remote sensing

Schätzungsweise 906 Erdbeobachtungssatelliten befinden sich derzeit im Erdorbit und stellen der Wissenschaft und Industrie täglich mehrere Terabyte an Daten zur Verfügung. Die Satelliten arbeiten sowohl mit Radar als auch optischen Sensoren und decken dabei verschiedene Spektralbereiche mit unterschiedlicher spektraler, räumlicher und zeitlicher Auflösung ab. Durch dieses breite Spektrum an geographischen Daten, ist es möglich, dass Fernerkundungsmethoden in vielen Industriebranchen und staatlichen Einrichtungen neue Anwendungsbereiche finden. Auf unserer Webseite finden Sie einige Projekte, in denen wir erfolgreich Satellitendaten eingesetzt haben, und mögliche Anwendungsfälle von Fernerkundungsmethoden für verschiedene Industrien.

Bekannte Satellitensysteme und -programme sind z.B. Sentinel-1 (Radar) und Sentinel-2 (optisch) von der ESA, Landsat (optisch) von der NASA, TerraSAR-X und TanDEM-X (beide Radar) von der DLR und PlanetScope (optisch) von Planet.

Es gibt im Wesentliche zwei Arten an geographischen Daten: Rasterdaten und Vektordaten.

Rasterdaten

Rasterdaten sind ein Gitter von regelmäßig angeordneten Pixeln, wobei jeder Pixel mit einem geographischen Standort verbunden ist, und werden als Matrix dargestellt. Die Pixelwerte hängen von der Art der Informationen ab, die gespeichert werden, z.B. Helligkeitswerte bei digitalen Bildern oder Temperaturwerte bei Wärmebildern. Die Größe der Pixel bestimmen außerdem die räumliche Auflösung des Rasters. Geographischen Rasterdaten werden also dazu verwendet, Satellitenbilder zu repräsentieren.

Rasterbilder enthalten in der Regel mehrere Bänder bzw. Kanäle, z.B. einen roten, grünen und blauen Kanal. Bei Satellitendaten gibt es zudem oft infrarote und/oder ultraviolette Bänder.

Vektordaten

Vektordaten repräsentieren geographische Eigenschaften auf der Erdoberfläche, wie z.B. Städte, Ländergrenzen, Straßen, Gewässer, Besitzrechte etc.. Solche Eigenschaften werden durch ein oder mehrere miteinander verbundene Vertices repräsentiert, wobei ein Vertex durch x-, y- und z-Werte eine Position im Raum festlegt. Ein einzelner Vertex ist ein Punkt, mehrere verbundene Vertices sind eine Linie und mehrere (>3) verbundene und geschlossene Vertices werden als Polygon bezeichnet. Die x-, y- und z-Werte sind dabei immer auf das entsprechende Koordinatenreferenzsystem (CRS) bezogen, das in Vektordateien als Metainformation gespeichert ist. Die gebräuchlichsten Dateiformate für Vektordaten sind GeoJSON, KML und SHAPEFILE. Um diese Daten prozessieren und analysieren zu können, werden verschiedene Tools benötigt.

Im Folgenden stelle ich die Tools vor, mit denen wir bei dida die besten Erfahrungen gemacht haben und die in unseren Fernerkundungsprojekten regelmäßig zum Einsatz kommen. Ich stelle ein Tool nach dem anderen vor, in folgende Kategorien gruppiert:

Abrufen von Satellitendaten
- EOBrowser
- Sentinelsat
- Sentinelhub
Verarbeitung von Rasterdaten
- Rasterio
- Pyproj
- SNAP (new)
- pyroSAR
- Rioxarray (new)
Verarbeitung von Vektordaten
- Shapely
- Python-geojson
- Geojson.io
- Geopandas
- Fiona
Bereitstellung geographischer Daten
- QGIS
- GeoServer
- Leafmap (new)
Verarbeitung meteorologischer Satellitendaten
- Wetterdienst
- Wradlib

Abrufen von Satellitendaten

EOBrowser

EOBrowser ist eine Web-Anwendung von sentinelhub, mit der man sich verschiedene Satellitendaten, wie z.B. Sentinel-1, Sentinel-2, Landsat, MODIS etc., anzeigen lassen kann. Die gewünschte Region, für die die Daten angezeigt werden sollen, kann man über die Suchleiste, durch Hochladen von KML oder GeoJSON Dateien oder durch Einzeichnen von Polygonen auswählen. Über den Schalter “Advanced search” können außerdem bestimmte Parameter angepasst werden, die die Suche verfeinern. Die Suche liefert sowohl einen Link, um die entsprechenden Rohdaten herunterzuladen, als auch verschiedene Visualisierungen. Für Sentinel-2 Daten kann man sich dabei z.B. das RGB-Bild aber auch den visualisierten Normalized Difference Water Index (NDWI) oder Normalized Difference Vegetation Index (NDWI) anzeigen lassen.

Sentinelsat

Sentinelsat ist ein Pythonpaket, mit dem man die Sentineldaten vom Copernicus Open Access Hub anfragen kann. Die Implementierung ist dabei sehr simpel und die Anfrage basiert auf Parametern, wie beispielsweise Region (z.B. in Form einer GeoJSON Datei), Datum, Platform (z.B. Sentinel-2), Produkttyp (z.B. L1C), maximale Wolkenbedeckung etc.

Sentinelsat liefert allerdings keine Tools, mit der die heruntergeladenen Rohdaten prozessiert und in geeignete Formate umgewandelt werden können. Das kann aber mit Hilfe von rasterio (siehe unten) relativ einfach implementiert werden.

Sentinelhub

Das Pythonpaket sentinelhub ist die Python-Schnittstelle, um sentinelhub Dienste nutzen zu können. Die Methoden und Funktionen erlauben es, die verschiedenen unterstützten Satellitendaten (wie beim EOBrowser) herunterzuladen und zu prozessieren. Beispielsweise gibt es Methoden, um sich Bilder mit einer Cloud Mask oder Mosaike der am wenigsten bewölkten Aufnahmen anzeigen zu lassen. Die vielen Funktionen sind sowohl in der Python- als auch in der API-Dokumentation ausführlich und verständlich erklärt. Sentinelhub bietet darüber hinaus an, die Satellitendaten und Metainformationen über die OGC Dienste WMS, WMTS, WCS und WFS abzufragen.

Leider stehen viele Funktionen von sentinelhub nicht kostenlos zur Verfügung. Je nach Plan können einige Features nicht genutzt werden und bestehen Limits bei der Anfrage von Daten.

Verarbeitung von Rasterdaten

Rasterio

Rasterio ist ein Pythonpaket, mit dem man geographische Rasterdaten, wie z.B. GeoTIFFs, verarbeiten kann. Dafür stehen viele Module und Funktionen zur Verfügung, um beispielsweise Rohdaten von Satelliten zu prozessieren, Rasterdaten einzulesen, deren geographischen Metadaten abzurufen, deren Koordinaten zu transformieren, die Bilder zu beschneiden, mehrere Bilder zusammenzufügen und die Daten in andere Formate abzuspeichern. Diese große Anzahl an Funktionen und die einfache Implementierung machen Rasterio zu einem Standardtool bei der Analyse von Satellitendaten.

Das untere Code-Beispiel zeigt, wie man mit Rasterio eine Sentinel-2 Szene einliest und als RGB-Bild darstellt und abspeichert:

import numpy as np
import rasterio
from rasterio import plot

file = 'S2BMSIL2A20210903T144719N0301R139T19LDG20210903T184347.tif'
# opening the GeoTIFF file and reading the R, G and B bands
with rasterio.open(file) as src:
	b, g, r = src.read(1), src.read(2), src.read(3)
	rgb = np.array([r, g, b])
	meta = src.meta.copy()

# updating the metadata to 3 bands
meta.update({'count': 3})

# saving the data as new GeoTIFF
with rasterio.open(file[:-4] + '_rgb.tif', 'w', **meta) as dest:
	dest.write(rgb)

# adjusting the RGB image for creating a nice plot
rgb[rgb > 3000] = 3000
img = ((rgb / rgb.max()) * 255).astype(np.uint8)
plot.show(img)

Das resultierende Bild sollte wie folgt aussehen:

Pyproj

Die Pythonbibliothek pyproj ist ein Python-Interface für PROJ, eine Softwarebibliothek, um geographische Koordinaten von einem Koordinatenreferenzsystem (CRS) in ein anderes zu transformieren. Sowohl geodätische Transformationen als auch kartographische Projektionen werden dabei unterstützt. Pyproj ist ein sehr simples Tool und unerlässlich bei Fernerkundungsprojekten, um die Koordinatenreferenzsysteme, die sich möglicherweise durch verschiedene genutzte Datenquellen ergeben, zu vereinheitlichen.

Das folgende Beispiel zeigt, wie man eine einfache Koordinatentransformation der Begrenzungen des obigen Sentinel-2 Bildes implementiert:

import rasterio
import pyproj

file = 'S2B_MSIL2A_20210903T144719_N0301_R139_T19LDG_20210903T184347.tif'
# reading the GeoTIFF file and returning the bounds and the coordinate reference
# system (CRS) which is EPS:32719 for Sentinel-2
with rasterio.open(file) as src:
    bounds = src.bounds
    crs_from = src.crs

# defining the CRS to be transformed into
epsg_id = 4326
crs_to = pyproj.CRS.from_epsg(epsg_id)

# defining the transformer which transforms from one CRS to another
proj = pyproj.Transformer.from_crs(crs_from=crs_from, crs_to=crs_to)
bounds_4326 = proj.transform_bounds(*bounds)

print('bounds in EPSG 32719:', bounds)
print('bounds in EPSG 4326:', bounds_4326)

Der Output wäre:

bounds in EPSG 32719: BoundingBox(left=399960.0, bottom=8590240.0, right=509760.0, top=8700040.0)
bounds in EPSG 4326: (-12.752446747147742, -69.92157939546101, -11.7580241630838, -68.91008576435853)

SNAP

SNAP ist die Abkürzung für Sentinel Application Platform und ist eine von der ESA entwickelte umfassende Toolbox für die wissenschaftliche Analyse von Satellitendaten. SNAP bietet eine intuitive Plattform für die Verarbeitung, Modellierung und Visualisierung von Satellitenbildern, insbesondere von den Sentinel-Missionen. Die Software ist für die Verarbeitung großer Mengen von Satellitendaten optimiert. Sie erleichtert auch die Verarbeitung von SAR-Daten, z. B. von Sentinel-1.

Die SNAP-Tools und -Methoden kann man auch in Python verwenden. Eine hilfreiche Anleitung, wie man die SNAP-Module für Python installiert und konfiguriert, finden Sie hier.

pyroSAR

Mit der Pythonbibliothek pyroSAR ist es möglich, Synthetic Aperture Radar (SAR) Satellitendaten, wie z.B. Sentinel-1, TerraSAR-X, TanDEM-X oder ALOS-2/PALSAR-2, zu prozessieren. Dafür stehen unter anderem Methoden, die auf SNAP und GAMMA Remote Sensing Tools zugreifen, und speziell für Sentinel-1 entwickelte Funktionen zur Verfügung.

Rioxarray

xarray ist ein Python-Paket, das Arrays mit Labels wie Dimension, Koordinaten und anderen spezifischen Attributen versieht. Damit macht es die Arbeit mit hochdimensionalen Arrays intuitiver. Rioxarray kombiniert die Funktionalitäten von rasterio mit allen Vorteilen von xarray.

Rioxarray vereinfacht die Arbeit mit hochdimensionalen Rasterdaten. Es gibt einen guten Überblick über die Attribute einer Rasterdatei, wie z.B Grenzen, Dimensionen, Koordinatenreferenzsystem und mehr. Leider ist das Paket nicht besonders ausführlich dokumentiert. Und für anspruchsvollere Operationen mit Satellitendaten sollte man weiterhin rasterio verwenden.

Verarbeitung von Vektordaten

Shapely

Shapely ist eine Pythonbibliothek zur Analyse und Manipulation von planaren geographischen Features. Die verschiedenen geographischen Objekte werden von einzelnen Klassen repräsentiert. Ein Punkt wird durch die Point Klasse definiert, eine Kurve durch die LineString und LinearRing Klassen und eine Oberfläche durch die Polygon Klasse. Eine Sammlung von Objekten werden durch die Klassen MultiPoint, MultiLineString und MultiPolygon implementiert. Diese Klassen sind der Standard bei der Verarbeitung und Manipulation von Vektordaten und werden auch von anderen Softwarebibliotheken und Pythonpaketen verwendet.

Shapely besitzt eine Vielzahl ein eingebauten Methoden, die sowohl die geometrischen Eigenschaften eines Objektes wie z.B. Fläche, Begrenzung und Länge bestimmen, aber auch die Beziehungen zu anderen Objekten wie z.B. ob ein Objekt ein anderes enthält, ob ein Objekt sich mit einem anderen überschneidet, ob ein Objekt ein anderes berührt u.v.m.

Shapely selbst besitzt keine Funktionen, um Vektordaten einzulesen, wofür andere Bibliotheken wie python-geojson, fiona oder geopandas (siehe unten) benötigt werden. Die Transformation der Features in shapley Objekte ist aber simpel.

In dem unten aufgeführten Code-Beispiel werden ein paar der oben vorgestellten Standardmethoden anhand einer GeoJSON Datei demonstriert:

import shapely.geometry
import geojson

file = 'natural_areas.geojson'
# reading GeoJSON file with geojson package (see below)
with open('natural_areas.geojson') as f:
    natural_areas = geojson.load(f)

natural_areas_shapely = []
# converting every polygon into a shapely polygon object
for polygon in natural_areas['features']:
    p = shapely.geometry.asShape(polygon['geometry'])
    # since the GeoJSON were downloaded from another source,
    # errors in the polygon can occur (e.g. it crosses itself)
    if p.is_valid:
      natural_areas_shapely.append(p)
    else:
      natural_areas_shapely.append(p.buffer(0))

print('geom_type:', natural_areas_shapely[1].geom_type)
print('area in km^2:', natural_areas_shapely[1].area)
print('bounds:', natural_areas_shapely[1].bounds)
print('object contains point:', natural_areas_shapely[1].contains(shapely.geometry.Point(-76.86, -6.81)))
print('object contains another object:', natural_areas_shapely[1].intersects(natural_areas_shapely[2]))

Der Output wäre:

geom_type: Polygon
area in km^2: 0.010795427253500077
bounds: (-76.952331, -6.932597, -76.794716, -6.783479)
object contains point: True
object contains another object: False

Python-geojson

GeoJSON ist ein Standardformat, um geographische Features wie Punkte, Linien oder Polygone und deren Koordinaten abzuspeichern. Das gleichnamige Pythonpaket liefert Funktionen, um diese Daten einzulesen, zu bearbeiten und als GeoJSON abzuspeichern. Durch die einfache Bedienbarkeit, ist geojson ein sehr nützliches Tool bei der Arbeit mit Vektordaten, die als GeoJSON vorliegen.

Im obigen Code-Beispiel ist demonstriert, wie man eine GeoJSON Datei mit dem geojson Pekt einliest.

Geojson.io

Geojson.io ist ein Geobrowser, mit dem man sich geographische Features auf einer Karte darstellen lassen kann. Dabei kann zwischen Karten von OpenStreetMap, Mapbox oder Satellitenbildern von Maxar auswählen. Die Satellitenbilder haben zudem eine hohe räumliche Auflösung, sind aber oftmals nicht ganz aktuell. Geojson.io bietet außerdem die Möglichkeit, GeoJSON Dateien hochzuladen und anzeigen zu lassen sowie Punkte, Linien und Polygone zu zeichnen. Die Features kann man wiederum als GeoJSON, Shapefile, KML, WKT und in anderen Formaten abspeichern. Dadurch ist geojson.io ein sehr hilfreiches Tool, um geographische Vektordaten zu erstellen und zu visualisieren ohne dabei Programmierkenntnisse zu besitzen.

Geopandas

Geopandas ist eine Erweiterung von pandas, einer Standardbibliothek für die Datenanalyse und -manipulation. Die aus pandas bekannten Datentypen und Methoden werden durch geographische und geometrische Objekte und Methoden erweitert, die auf shapely basieren. Geopandas besitzt außerdem Funktionen, um Vektordaten zu lesen und zu speichern. Diese Funktionen wiederum hängen von fiona ab (s.u.).

In geopandas lassen sich shapely-Objekte in einer GeoSeries zusammenfassen, auf die man die verschiedenen geometrischen Methoden aber auch klassische Operationen zur Datenanalyse und -manipulation anwenden kann. Zusätzlich besitzt geopandas praktische Tools zur Geocodierung.

Fiona

Mit fiona können alle von GDAL unterstützten Raster- und Vektordateiformate eingelesen und Daten auch als solche gespeichert werden.

Darüber hinaus besitzt fiona ein paar grundlegende Funktionen für die geometrischen Analyse, z.B. die Bestimmung von Grenzen eines Objekts, und Koordinatentransformationen, die auf shapely und pyproj basieren.

Bereitstellung geographischer Daten

QGIS

QGIS ist ein kostenloses geographisches Informationssystem für den Desktop. Es ist eine Software, mit der man geographische Informationen (Raster- und Vektordaten) darstellen, bearbeiten, erstellen, analysieren und veröffentlichen kann.

Da QGIS mit vielen Funktionen und Features ausgestattet ist, erschließt sich das gesamte Potential der Software vor allem Personen, die viel damit arbeiten. Aber die Darstellung von Satellitendaten ist einfach umzusetzen und hilft sehr bei der Exploration von geographischen Daten. In QGIS können auch geographische Daten, die über WMS, WFS, WCS, WPS oder WMTS veröffentlich wurden, abgerufen und angeschaut werden.

GeoServer

GeoServer ist ein kostenloser Server, mit dem geographische Daten veröffentlicht werden können. Dabei werden die Open Geospatial Consortium Standardprotokolle Web Map Service (WMS), Web Feature Service (WFS), Web Coverage Service (WCS), Web Processing Service (WPS) und Web Map Tile Service (WMTS) verwendet.

Geographische Raster- und Vektordaten lassen sich über die oben genannten Protokolle veröffentlichen und anschließend in Web GIS Anwendungen (z.B. über Leaflet, Mapbox oder OpenLayers) oder in GIS-Softwares (z.B. QGIS oder ArcGIS) darstellen. Dadurch hat man die Möglichkeit, seine geographischen Daten interaktiv ohne lange Ladezeiten und große Datentransfers bereitzustellen.

Leafmap

Leafmap ist ein Python-Paket für interaktive Visualisierungen von Geodaten in Jupyter Notebooks. Mit ein paar Zeilen Code ist es einfach, interaktive Karten zu gestalten. Leafmap ist insofern sehr leistungsfähig, da es auf verschiedenen geographischen Karten- und Analysetools wie ipyleaflet, whiteboxtools und folium aufbaut und diese integriert.

Leafmap unterstützt alle gängigen Vektor- und Rasterdatenformate. Man kann aus verschiedenen Basemaps wählen, darunter Satellitenkarten, Hybridkarten, Straßen, Geländekarten und vielen mehr. Es ist außerdem möglich, eigene Polygone zu zeichnen und diese als Vektordateien zu exportieren oder für weitere Operationen zu verwenden. Das Paket bietet auch eine Upload-Funktion. Die Karten sind in hohem Maße konfigurierbar, da viele Stiloptionen für Polygone, Linien, Marker, Popups und Legenden existieren. Falls die Geodaten auf einem Tiel Server gehostet werden, z.B. über Geoserver, kann leafmap Raster- und Vektordaten von Diensten wie WMS oder WFS übertragen.

Verarbeitung meteorologischer Satellitendaten

Wetterdienst

Wetterdienst ist ein Pythonpaket, mit dem Wetterdaten des Deutschen Wetterdiensts (DWD) abgerufen werden können. Diese beinhalten historische Wetterdaten, Vorhersagen basierend auf Wettermodellen und Radardaten (u.a. RADOLAN und RADVOR). Die zeitliche Auflösung der einfachen Radardaten (wie z.B. Temperatur, atmosphärischer Niederschlag etc.) reicht dabei von monatlich bis jährlich. RADOLAN-Daten werden mit einer minütlichen bis täglichen zeitlichen Auflösung geliefert.

Derzeit gibt es leider keine eingebauten Funktionen, um Zeitreihen der RADOLAN-Daten zu erstellen, die eigentlich standardmäßig für die Analyse von Wetterdaten verwendet werden.

Wradlib

Wradlib ist eine Pythonbibliothek, mit der man verschiedene Wetterradardaten prozessieren kann, u.a. RADOLAN, DX und BUFR. Hilfreich ist vor allem, dass die Wetterdaten in geeignete Rasterformate, wie z.B. GeoTIFF, umgewandelt werden können.

Fazit

Durch die große Anzahl an Erdbeobachtungssatelliten finden Fernerkundungsmethoden in vielen Bereichen neue Anwendungsgebiete und stellen wertvolle Informationen über die Erdoberfläche und deren Beschaffenheit zur Verfügung. dida hat in einigen Projekten erfolgreich nachgewiesen, wie Fernerkundungsmethoden in Kombination mit Machine Learning neue Erkenntnisse liefern können, die durch klassische menschliche und manuelle Methoden nur schwer zu erreichen sind. Um allerdings die verschiedenen Satellitendaten prozessieren und analysieren zu können, werden unterschiedliche Tools benötigt. Einige davon, die wir bei dida regelmäßig einsetzen, wurden in diesem Artikel erklärt und demonstriert. Ich hoffe dadurch, eine technische Grundlage geschaffen und den Einstieg in dein nächstes Fernerkundungsprojekt erleichtert zu haben. Wir haben auch einige Projekte im Bereich der Fernerkundung durchgeführt, z. B. die automatische Klassifizierung von Pflanzenarten und die automatische Erkennung von Kleinbergbau.