Die Mündung des Syrdarja und den Aralsee mit Sentinel-2-Bildern überwachen und mit dem NDWI auf CODE-DE-Lab analysieren
Vergleichen Sie mit Sentinel-2A-Level-2A-Bildern von CODE-DE-Lab die wasserbezogenen spektralen Muster an der Mündung des Syrdarja und am Aralsee in Kasachstan von Mai bis Juli 2023. Das Beispiel durchsucht den OData-Katalog und erstellt Raster des Normalized Difference Water Index (NDWI) aus der Oberflächenreflexion im grünen und nahinfraroten Spektralbereich. Der NDWI misst weder Wassertiefe noch Wasservolumen. Die verfügbaren Aufnahmen bilden keine tägliche Zeitreihe.
Was ist NDWI?
Der NDWI ist ein normalisierter Differenzindex aus der Reflexion in zwei Spektralbändern: Grün (G) und Nahinfrarot (NIR). Er wird mit folgender Formel berechnet:
Dabei gilt:
G – Reflexion im grünen Band,
NIR – Reflexion im nahinfraroten Band.
NDWI-Werte liegen zwischen -1 und 1. Höhere Werte weisen auf Wasser hin. Niedrigere Werte können andere Landbedeckungen wie Vegetation oder bebaute Flächen anzeigen.
Anwendungsbereiche des NDWI
Raster mit NDWI-Werten lassen sich in verschiedenen Bereichen verwenden:
- Überwachung von Wasserressourcen
Vergleich der Wasserflächen von Seen, Flüssen und Stauseen, auch in Gebieten mit Dürre oder übermäßigen Niederschlägen.
- Katastrophenmanagement
Schnelle Erkennung überfluteter Gebiete und Bewertung von Schäden.
- Umweltschutz
Überwachung von Feuchtgebieten und anderen aquatischen Ökosystemen zum Schutz ihrer Artenvielfalt.
- Stadtplanung
Analyse der Auswirkungen städtischer Entwicklung auf nahe gelegene Wasserressourcen.
- Landwirtschaft
Bewertung der Bewässerung landwirtschaftlicher Flächen und Überwachung der Pflanzengesundheit.
Ablauf
Variablen für die Daten definieren
Mit der OData-API arbeiten
Den OData-Endpunkt konfigurieren
Die OData-Abfrage erstellen
Produktpfade ermitteln
Den NDWI berechnen
Funktionen zum Lesen, Berechnen und Speichern von Rasterdaten definieren
NDWI-Raster erzeugen
Ergebnisse analysieren
Voraussetzungen
Nr. 1 Hosting und Zugriff auf EODATA
Sie benötigen ein Hosting-Konto bei CODE-DE-Lab mit Zugriff auf https://identity.cloudferro.com/auth/realms/CODE-DE3/account/. Führen Sie das Notebook in einer Umgebung aus, in der die Produktpfade aus dem Katalog unter /eodata lesbar sind. Katalogabfragen liefern Metadaten; sie binden die Bilddaten weder ein noch laden sie diese herunter. Für eine virtuelle Linux-Maschine bereiten Sie die Einbindung anhand von EODATA unter Linux mit S3FS auf CODE-DE-Lab einbinden vor.
Nr. 2 Jupyter Notebook und Python-Pakete
Richten Sie Jupyter Notebook anhand der offiziellen Installationsanleitung ein. Installieren Sie die folgenden Pakete in der Python-Umgebung des Notebook-Kernels. Rasterio benötigt JPEG2000-Unterstützung, um die Sentinel-2-Bänder zu lesen.
python -m pip install requests numpy rasterio matplotlib
Laden Sie sentinel2_aral_lake.ipynb herunter und öffnen Sie das Notebook in dieser Umgebung. Das Notebook erstellt das Ausgabeverzeichnis aral_ndwi in seinem Arbeitsverzeichnis. Dieses muss beschreibbar sein.
Nr. 3 Katalog-Endpunkt und Produktmetadaten
Die OData-Abfrageoptionen sind in Anleitung zur EODATA-Katalog-API auf CODE-DE-Lab beschrieben. Verwenden Sie den unten angegebenen Endpunkt für CODE-DE-Lab. Prüfen Sie, ob der gewählte Katalog das Gebiet am Aralsee und den Zeitraum Mai–Juli 2023 abdeckt. Wenn die Abfrage keine Produkte liefert, wählen Sie ein Gebiet und einen Zeitraum, die auf Ihrer Plattform verfügbar sind. Allein der Wechsel des Katalog-Endpunkts ermöglicht keinen Zugriff auf die Bilddaten. Die Verarbeitung benötigt das vollständige Level-2A-SAFE-Produkt einschließlich MTD_MSIL2A.xml, der Bänder B03 und B08 mit 10 m Auflösung und der Szenenklassifikation (SCL) mit 20 m Auflösung. Die Dokumentation zur Sentinel-2-Verarbeitung erläutert die Reflexions-Offsets und Quantifizierungswerte aus diesen Metadaten.
import requests
Variablen für die Daten definieren
Wählen Sie den Aufnahmezeitraum und das Gebiet an der Mündung des Syrdarja und am Aralsee. Die BBOX (Bounding Box) enthält die Koordinaten Xmin, Ymin, Xmax und Ymax im Koordinatenreferenzsystem EPSG:4326.
Definieren Sie die drei Variablen in Python:
# Timerange of data that we want to recieve
start_date = '2023-05-01'
end_date = '2023-07-31'
# Geometry in form of a BBOX
bbox = [60.2347,46.1731,61.1160,46.6541]
So werden die Variablen in Jupyter Notebook dargestellt:
Mit der OData-API arbeiten
Die OData-API durchsucht EODATA-Produkte mit räumlichen und zeitlichen Filtern. Die Antwort enthält den Produktpfad im Feld S3Path. Für die folgenden Verarbeitungsschritte müssen diese Pfade in Ihrer Umgebung unter /eodata zugänglich sein.
Den OData-Endpunkt konfigurieren
Setzen Sie im heruntergeladenen Notebook ODATA_CATALOGUE_URL auf https://catalogue.code-de.org/odata/v1.
ODATA_CATALOGUE_URL = "https://catalogue.code-de.org/odata/v1"
API_URL = ODATA_CATALOGUE_URL.rstrip("/") + "/Products"
Die OData-Abfrage erstellen
Die Abfrage wählt Sentinel-2A-Level-2A-Produkte aus, die die Bounding Box schneiden, zwischen dem 1. Mai und dem 31. Juli 2023 aufgenommen wurden und höchstens 10 Prozent Wolkenbedeckung auf Produktebene aufweisen. Die spätere Verarbeitung beschneidet jede Granule auf die Bounding Box und wendet eine Qualitätsmaske auf Pixelebene an. requests kodiert die Abfrageparameter für die URL.
minx, miny, maxx, maxy = bbox
polygon = (
f"POLYGON(({minx} {miny},{maxx} {miny},{maxx} {maxy},"
f"{minx} {maxy},{minx} {miny}))"
)
query_params = {
"$filter": (
"Collection/Name eq 'SENTINEL-2' and startswith(Name,'S2A') and "
"Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' "
"and att/OData.CSC.StringAttribute/Value eq 'S2MSI2A') and "
"Attributes/OData.CSC.DoubleAttribute/any(att:att/Name eq 'cloudCover' "
"and att/OData.CSC.DoubleAttribute/Value le 10) and "
f"ContentDate/Start ge {start_date}T00:00:00.000Z and "
f"ContentDate/Start le {end_date}T23:59:59.999Z and "
f"OData.CSC.Intersects(area=geography'SRID=4326;{polygon}')"
),
"$orderby": "ContentDate/Start asc",
"$top": 1000,
"$select": "Id,Name,S3Path",
}
specified_url = requests.Request("GET", API_URL, params=query_params).prepare().url
specified_url
Senden Sie die Anfrage und folgen Sie @odata.nextLink, solange weitere Ergebnisseiten verfügbar sind:
products = []
next_url = specified_url
while next_url:
response = requests.get(next_url, timeout=90)
response.raise_for_status()
page = response.json()
products.extend(page["value"])
next_url = page.get("@odata.nextLink")
if not products:
raise RuntimeError("No Sentinel-2 Level-2A products matched this region and date range. Check catalogue coverage before changing the query parameters.")
Zeigen Sie das erste Produkt an:
products[0]
Beispiel für die zurückgegebenen Felder:
{
"Id": "49ec07df-edc5-4410-80ee-28312a787128",
"Name": "S2A_MSIL2A_20230503T064621_N0510_R020_T40TGS_20240903T161656.SAFE",
"S3Path": "/eodata/Sentinel-2/MSI/L2A_N0500/2023/05/03/S2A_MSIL2A_20230503T064621_N0510_R020_T40TGS_20240903T161656.SAFE"
}
Produktpfade ermitteln
Lesen Sie die Produktpfade aus dem Feld S3Path:
dataset = [product["S3Path"] for product in products]
dataset
Den NDWI berechnen
Rechnen Sie die digitalen Werte von B03 und B08 vor der NDWI-Berechnung in Oberflächenreflexion um. Die mit Processing Baseline 04.00 eingeführten bandspezifischen Offsets heben sich im Nenner des Index nicht auf. Gleitkommaarithmetik verhindert außerdem einen Ganzzahlüberlauf beim Addieren oder Subtrahieren der Bandwerte.
from pathlib import Path
import xml.etree.ElementTree as ET
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.mask import mask
from rasterio.vrt import WarpedVRT
from rasterio.warp import transform_geom
Funktionen zum Lesen, Berechnen und Speichern von Rasterdaten
Die folgenden Funktionen lesen den Quantifizierungswert und die B03/B08-Offsets aus MTD_MSIL2A.xml. Sie beschneiden jede Granule auf die Bounding Box und passen die SCL-Ebene mit 20 m Auflösung durch Nearest-Neighbour-Resampling an die Bänder mit 10 m Auflösung an. Nur die SCL-Klassen 4, 5 und 6 (Vegetation, unbewachsener Boden und Wasser) bleiben erhalten. NoData-Pixel, negative Reflexionswerte und Nenner kleiner oder gleich null werden ausgeschlossen. Fehlende Kalibrierungsmetadaten führen zum Abbruch der Verarbeitung, damit keine falschen Werte unbemerkt verwendet werden.
def read_scaling(product_dir):
metadata = product_dir / "MTD_MSIL2A.xml"
root = ET.parse(metadata).getroot()
elements = list(root.iter())
def value(tag):
return next(e.text for e in elements if e.tag.split("}")[-1] == tag)
quantification = float(value("BOA_QUANTIFICATION_VALUE"))
baseline = float(value("PROCESSING_BASELINE"))
offsets = {
e.attrib["band_id"]: float(e.text)
for e in elements if e.tag.split("}")[-1] == "BOA_ADD_OFFSET"
}
# Sentinel-2 metadata numbers B03 and B08 as band IDs 2 and 7.
if baseline < 4.0 and not offsets:
offsets = {"2": 0.0, "7": 0.0}
if not np.isfinite(quantification) or quantification <= 0:
raise ValueError(f"Invalid reflectance quantification in {metadata}")
if any(key not in offsets or not np.isfinite(offsets[key]) for key in ("2", "7")):
raise ValueError(f"Missing or invalid B03/B08 offsets in {metadata}")
return quantification, offsets
def find_band(granule, band, resolution):
matches = sorted(granule.glob(f"IMG_DATA/R{resolution}m/*_{band}_{resolution}m.jp2"))
if len(matches) != 1:
raise ValueError(f"Expected one {band} file at {resolution} m in {granule}")
return matches[0]
def process_granule(product_dir, granule, bbox, output_dir):
quantification, offsets = read_scaling(product_dir)
minx, miny, maxx, maxy = bbox
geometry = {
"type": "Polygon",
"coordinates": [[[minx, miny], [maxx, miny], [maxx, maxy],
[minx, maxy], [minx, miny]]],
}
with (
rasterio.open(find_band(granule, "B03", 10)) as green_src,
rasterio.open(find_band(granule, "B08", 10)) as nir_src,
rasterio.open(find_band(granule, "SCL", 20)) as scl_src,
):
if green_src.crs is None or scl_src.crs is None:
raise ValueError(f"Missing raster CRS in {granule}")
if (green_src.crs != nir_src.crs or green_src.transform != nir_src.transform
or green_src.shape != nir_src.shape):
raise ValueError(f"B03 and B08 grids differ in {granule}")
shape = transform_geom("EPSG:4326", green_src.crs, geometry)
try:
green_dn, transform = mask(green_src, [shape], crop=True, filled=False, indexes=1)
except ValueError as exc:
if "do not overlap" in str(exc):
return None
raise
nir_dn, nir_transform = mask(nir_src, [shape], crop=True, filled=False, indexes=1)
if transform != nir_transform or green_dn.shape != nir_dn.shape:
raise ValueError(f"Cropped B03 and B08 grids differ in {granule}")
# Resample the categorical 20 m SCL grid onto the 10 m band grid.
with WarpedVRT(scl_src, crs=green_src.crs, transform=green_src.transform,
width=green_src.width, height=green_src.height,
resampling=Resampling.nearest) as scl_vrt:
scl, _ = mask(scl_vrt, [shape], crop=True, filled=False, indexes=1)
green = (green_dn.data.astype("float32") + offsets["2"]) / quantification
nir = (nir_dn.data.astype("float32") + offsets["7"]) / quantification
denominator = green + nir
valid = (
~np.ma.getmaskarray(green_dn) & ~np.ma.getmaskarray(nir_dn)
& ~np.ma.getmaskarray(scl)
& (green_dn.data != 0) & (nir_dn.data != 0)
& np.isin(scl.data, [4, 5, 6]) # vegetation, bare soil, water
& np.isfinite(green) & np.isfinite(nir)
& (green >= 0) & (nir >= 0) & (denominator > 0)
)
if not valid.any():
return None
ndwi = np.full(green.shape, -9999.0, dtype="float32")
np.divide(green - nir, denominator, out=ndwi, where=valid)
output = output_dir / f"{product_dir.stem}_{granule.name}_NDWI.tif"
with rasterio.open(output, "w", driver="GTiff", height=ndwi.shape[0],
width=ndwi.shape[1], count=1, dtype="float32",
crs=green_src.crs, transform=transform,
nodata=-9999.0, compress="deflate") as dst:
dst.write(ndwi, 1)
return output
NDWI-Raster erzeugen
Erstellen Sie das Ausgabeverzeichnis und verarbeiten Sie jedes zurückgegebene Produkt. Granulen außerhalb der Bounding Box oder ohne gültige Pixel werden übersprungen. Jede GeoTIFF-Datei behält das ursprüngliche projizierte Raster ihrer Granule bei und verwendet -9999 als NoData-Wert. Überlappende Kacheln werden nicht zu einem Mosaik zusammengefügt.
output_dir = Path.cwd() / "aral_ndwi"
output_dir.mkdir(parents=True, exist_ok=True)
outputs = []
for item in dict.fromkeys(dataset):
product_dir = Path(item)
if not product_dir.is_dir():
raise FileNotFoundError(
f"Product is not accessible: {product_dir}. "
"Mount EODATA at /eodata in the environment running this notebook."
)
granules = sorted(p for p in (product_dir / "GRANULE").iterdir() if p.is_dir())
if not granules:
raise ValueError(f"No granules found in {product_dir}")
for granule in granules:
output = process_granule(product_dir, granule, bbox, output_dir)
if output is not None:
outputs.append(output)
print(output)
if not outputs:
raise RuntimeError("No valid pixels remain inside the bounding box after masking.")
Zeigen Sie das erste erzeugte Raster mit einer festen NDWI-Farbskala an. Maskierte Pixel sind transparent. Wählen Sie einen anderen Eintrag in outputs, um ein anderes Raster anzuzeigen.
import matplotlib.pyplot as plt
from rasterio.plot import plotting_extent
with rasterio.open(outputs[0]) as src:
ndwi = src.read(1, masked=True)
fig, ax = plt.subplots(figsize=(10, 6))
image = ax.imshow(ndwi, extent=plotting_extent(src), cmap="RdYlBu",
vmin=-1, vmax=1)
ax.set_title(outputs[0].name)
ax.set_xlabel(f"Easting ({src.crs})")
ax.set_ylabel(f"Northing ({src.crs})")
fig.colorbar(image, ax=ax, label="NDWI")
plt.show()
Analyse
Der folgende Vergleich für Mai, Juni und Juli 2023 stammt aus dem früheren Verarbeitungsablauf. Er wurde nicht mit der korrigierten Reflexionsumrechnung und Qualitätsmaske neu erzeugt.
Historische Abbildung aus dem früheren Verarbeitungsablauf.
Verwenden Sie neu erzeugte Raster, um die gewählten Aufnahmen zu bewerten. Höhere NDWI-Werte können auf Wasser hinweisen. Ein Schwellenwert allein bestimmt jedoch weder Wassertiefe noch Wasservolumen oder die Ursache einer Veränderung. Wolkenmaskierung und Abdeckung unterscheiden sich zwischen den Aufnahmen. Die Ergebnisse sind einzelne Granulen, keine täglichen oder monatlichen Komposite. Vergleichen Sie übereinstimmende Gebiete mit gültigen Beobachtungen, bevor Sie Rückschlüsse auf die Ausdehnung der Wasserflächen ziehen.
Weiterführende Schritte
Weitere Beispiele für die Untersuchung der Umwelt mit Satellitenbildern:
Städtische Ausbreitung:
Überschwemmungen:
Analyse und Überwachung von Überschwemmungen mit Python und Sentinel-2 Daten auf CODE-DE-Lab
Luftverschmutzung:
Verarbeitung von Sentinel-5P-Daten zur Luftverschmutzung mit Jupyter Notebook auf CODE-DE-Lab