Landbedeckung (ESA WorldCover)

In diesem Notebook lernst du, wie du mit der Bibliothek cdec eine Landbedeckungskarte für ein bestimmtes Gebiet herunterladen, auf einer Karte darstellen und anschließend als Rasterdatensatz untersuchen kannst.

Die ESA-WorldCover-Karte beschreibt die Erdoberfläche mit elf Klassen, beispielsweise Baumbedeckung, Grasland, Ackerland, bebaute Flächen und Gewässer. Der Datensatz basiert auf Sentinel-1- und Sentinel-2-Daten aus dem Jahr 2021 und besitzt eine räumliche Auflösung von 10 Metern.

Im Gegensatz zu Höhen- oder Temperaturdaten enthält dieses Raster keine kontinuierlichen Messwerte. Jeder Pixelwert ist eine Nummer, die für eine bestimmte Landbedeckungsklasse steht.

1 Räumliche Ausdehnung festlegen

Zunächst wird festgelegt, für welches Gebiet die Landbedeckung 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?

ESA WorldCover 2021 ist eine fertige, statische Landbedeckungskarte. Für den Download muss deshalb kein eigener Zeitraum ausgewählt werden.

Die Funktion download_land_cover() verwendet automatisch den Zeitraum des WorldCover-Produkts von 2021. Die Karte zeigt somit die für dieses Jahr abgeleitete Landbedeckung und keine tagesaktuelle Situation.

3 Landbedeckung 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 Landbedeckungsklassen als auch die räumlichen Informationen gespeichert.

Anschließend wird die Bibliothek cdec importiert. Mit der Funktion download_land_cover() wird die WorldCover-Karte für das zuvor definierte Gebiet heruntergeladen.

Dabei werden zwei Angaben übergeben:

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

Standardmäßig wird die Datei als Cloud Optimized GeoTIFF (COG) gespeichert. Für die Übersichtsebenen wird die Methode nearest verwendet, damit beim Verkleinern keine neuen oder gemischten Klassenwerte entstehen.

file_path = "data/land_cover.tif"

import cdec

cdec.download_land_cover(
    file_path,
    spatial_extent=spatial_extent,
)

Authenticated using refresh token.

PosixPath('/Users/soeren/PycharmProjects/cdec/snippets/data_access/notebooks/data/land_cover.tif')

4 Rasterdaten auf einer Karte anzeigen

Die heruntergeladene Landbedeckung 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 WorldCover-Farbpalette ordnet jeder Klasse eine feste Farbe zu. Mit add_legend() wird zusätzlich eine Legende eingeblendet, über die sich die Farben den Landbedeckungsklassen zuordnen lassen.

from geolibre import Map

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

layer_id = m.add_cog(
    file_path,
    name="ESA WorldCover 2021",
)

m.add_legend(
    builtin="esa_worldcover",
    title="Landbedeckung",
)

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 Klassenwerte 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)

WorldCover besitzt ein Rasterband vom Typ uint8. Dieser ganzzahlige Datentyp ist für die Speicherung der elf Klassen ausreichend.

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: ('uint8',)
Shape: (1080, 2193)
Coordinate Reference System: EPSG:4326

7 Landbedeckungsklassen kennenlernen

Die Zahlen im Raster sind Codes für unterschiedliche Arten der Landbedeckung. Beispielsweise steht der Wert 10 für Baumbedeckung, 40 für Ackerland und 50 für bebaute Flächen.

Die Bibliothek cdec stellt die offiziellen Klassennamen und Farben als WORLDCOVER_CLASSES und WORLDCOVER_COLORS bereit. Daraus kann eine übersichtliche Tabelle erzeugt werden.

import pandas as pd

class_table = pd.DataFrame(
    {
        "Klassenwert": cdec.WORLDCOVER_CLASSES.keys(),
        "Klasse": cdec.WORLDCOVER_CLASSES.values(),
        "Farbe": cdec.WORLDCOVER_COLORS.values(),
    }
)

class_table

	Klassenwert	Klasse	Farbe
0	10	Tree cover	#006400
1	20	Shrubland	#ffbb22
2	30	Grassland	#ffff4c
3	40	Cropland	#f096ff
4	50	Built-up	#fa0000
5	60	Bare / sparse vegetation	#b4b4b4
6	70	Snow and ice	#f0f0f0
7	80	Permanent water bodies	#0064c8
8	90	Herbaceous wetland	#0096a0
9	95	Mangroves	#00cf75
10	100	Moss and lichen	#fae6a0

8 Klassenanteile im Untersuchungsgebiet bestimmen

Mit src.read(1, masked=True) wird das Landbedeckungsband eingelesen. Pixel ohne gültige Daten werden durch masked=True von der Analyse ausgeschlossen.

Anschließend zählt numpy.unique(), wie häufig jeder Klassenwert vorkommt. Aus den Pixelzahlen werden prozentuale Anteile berechnet. So lässt sich beispielsweise untersuchen, welcher Anteil des Gebiets bebaut, bewaldet oder landwirtschaftlich genutzt ist.

import numpy as np

data = src.read(1, masked=True)
values, counts = np.unique(data.compressed(), return_counts=True)

class_distribution = pd.DataFrame(
    {
        "Klassenwert": values,
        "Klasse": [cdec.WORLDCOVER_CLASSES.get(int(value), "Unknown") for value in values],
        "Pixel": counts,
        "Anteil (%)": counts / counts.sum() * 100,
    }
)

class_distribution

	Klassenwert	Klasse	Pixel	Anteil (%)
0	10	Tree cover	892161	37.668719
1	30	Grassland	422907	17.855930
2	40	Cropland	35634	1.504535
3	50	Built-up	871015	36.775895
4	60	Bare / sparse vegetation	6971	0.294329
5	80	Permanent water bodies	134752	5.689483
6	90	Herbaceous wetland	5000	0.211109

9 Rasterdaten visualisieren

Zum Abschluss werden die eingelesenen Rasterdaten direkt im Notebook dargestellt.

from rasterio.plot import show

show(src)