Städtische Ausbreitung mit Sentinel-1-SAR-Daten und Schwellenwerten für Pixelwerte und Polarisationen auf CODE-DE-Lab überwachen
In diesem Artikel verwenden Sie Schwellenwertverfahren und vergleichen die Werte der VV- und VH-Polarisationen, um städtische Gebiete zu analysieren. Als Beispiel laden Sie Sentinel-1-Aufnahmen des Stadtgebiets von Warschau in Polen herunter.
Die Analyse von Aufnahmen mit Synthetic Aperture Radar (SAR) macht Unterschiede in der Reflexion von Radarwellen sichtbar. Diese Unterschiede sind entscheidend, um bebaute Gebiete zu erkennen. Die charakteristischen Reflexionsmuster in SAR-Aufnahmen ermöglichen es, Ballungsräume und kleinere städtische Siedlungen abzugrenzen. Gebäude, Straßen und andere künstliche Infrastrukturelemente reflektieren Radarwellen besonders stark. Dadurch lassen sich bebaute Gebiete deutlich von natürlichem Gelände unterscheiden. Eine solche Analyse, beispielsweise des Ballungsraums Warschau, kann folgenden Zwecken dienen:
Überwachung des Stadtwachstums und der Stadtentwicklung
Stadtplanung
Bewertung natürlicher und menschengemachter Gefahren
Umweltmanagement
Beispielergebnis des ursprünglichen Workflows:
Das Bild zeigt deutlich den Ballungsraum Warschau und das angrenzende Gebiet zwischen den Ballungsräumen Warschau und Łódź. Für dieses Gebiet ist derzeit eine Entwicklung geplant, die beide Ballungsräume miteinander verbinden soll. Im zentralen Teil des Bildausschnitts ist ein großer Flughafenkomplex geplant.
Was wir behandeln werden
Auf Daten zugreifen
Städtische Gebiete mithilfe von Schwellenwerten abgrenzen
Einzelne Pixel quantitativ analysieren
Voraussetzungen
Nr. 1 EODATA-Zugriff
Melden Sie sich im CODE-DE-Lab Tenant Manager an. Erstellen Sie unter Configuration → EODATA Access persönliche S3-Zugangsdaten und halten Sie den Zugriffsschlüssel und den geheimen Schlüssel bereit. Folgen Sie der CODE-DE-Lab-Migrationsanleitung für die Kontoaktivierung und den EODATA-Zugriff. Die Katalogabfrage ist öffentlich; zum Herunterladen der Messdaten benötigen Sie diese EODATA-Zugangsdaten.
Nr. 2 Verbindungseinstellungen für Katalog und S3
Die Abfrage benötigt einen OData-Katalog mit den Feldern Collection, Attributes, ContentDate und S3Path sowie dem räumlichen Filter OData.CSC.Intersects. Der S3Path des ausgewählten Produkts muss auf den Bucket verweisen, der über Ihren S3-Endpunkt erreichbar ist. Eine STAC-URL kann in diesem OData-Workflow nicht verwendet werden.
Geben Sie bei den Eingabeaufforderungen diese CODE-DE-Lab-Einstellungen an oder setzen Sie die Umgebungsvariablen vor dem Start von Jupyter:
Umgebungsvariable |
Wert |
|---|---|
EODATA_CATALOGUE_URL |
|
EODATA_S3_ENDPOINT |
|
EODATA_S3_BUCKET |
eodata |
EODATA_S3_REGION |
s3.fra1-3.cloudferro.com |
Dieser Workflow verwendet den OData-Katalog. Für neue CODE-DE-Lab-Katalogworkflows wird der STAC-Katalog empfohlen.
Nr. 3 Jupyter Notebook oder JupyterLab mit Python 3.9 oder neuer
Verwenden Sie eine Umgebung mit Netzwerkzugriff auf den ausgewählten Katalog und S3-Endpunkt. Sie benötigt genügend Arbeitsspeicher für zwei vollständige Sentinel-1-Messdatenarrays und die daraus abgeleiteten Arrays sowie Speicherplatz für die heruntergeladenen TIFF-Dateien und Ergebnisse. Falls Jupyter noch nicht installiert ist, folgen Sie der Jupyter-Installationsanleitung.
Nr. 4 Python-Pakete
Installieren Sie folgende Pakete in der Python-Umgebung Ihres Notebooks:
python -m pip install requests boto3 numpy rasterio matplotlib
Laden Sie urban_sprawl.ipynb herunter und führen Sie die Zellen der Reihe nach aus. Das Notebook und der folgende Code verwenden dieselben Verbindungseinstellungen. Die S3-Konfiguration liest EODATA_ACCESS_KEY und EODATA_SECRET_KEY aus der Umgebung oder fragt die Werte verdeckt ab, falls sie nicht gesetzt sind. Speichern Sie keine echten Schlüssel in den Zellen.
Die gespeicherten Abbildungen und numerischen Ausgaben veranschaulichen einen früheren Durchlauf; das ausgewählte Produkt und Ihre Ergebnisse können abweichen.
Schritt 1: Auf Daten in der Cloud zugreifen
1. Bibliotheken importieren und Variablen deklarieren
Die Abfrage wählt ein IW-GRDH-Produkt mit VV/VH-Polarisation im angegebenen Zeitraum und im Untersuchungsgebiet Warschau aus.
import os
from pathlib import Path
from getpass import getpass
import boto3
from botocore.config import Config
import requests
import numpy as np
import rasterio
import matplotlib.pyplot as plt
# Use the OData Products URL documented by your data provider.
CATALOGUE_URL = (
os.environ.get("EODATA_CATALOGUE_URL")
or input("OData Products URL (ending in /odata/v1/Products): ")
).strip().rstrip("/")
if not CATALOGUE_URL.startswith("https://") or not CATALOGUE_URL.endswith("/odata/v1/Products"):
raise ValueError("Provide an HTTPS OData Products URL ending in /odata/v1/Products.")
# Warsaw area of interest.
collection = "SENTINEL-1"
product_type = "GRD"
start_date = "2024-05-15T00:00:00.000Z"
end_date = "2024-05-17T00:05:00.000Z"
polygon_coords = (
"20.851215 52.261859, 21.271151 52.261859, "
"21.271151 52.017798, 20.851215 52.017798, 20.851215 52.261859"
)
query_params = {
"$filter": (
f"Collection/Name eq '{collection}' and "
"contains(Name,'_IW_GRDH_1SDV_') and "
"Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' "
f"and att/OData.CSC.StringAttribute/Value eq '{product_type}') and "
f"ContentDate/Start gt {start_date} and ContentDate/End lt {end_date} and "
f"OData.CSC.Intersects(area=geography'SRID=4326;POLYGON(({polygon_coords}))')"
),
"$orderby": "ContentDate/Start asc",
"$top": 1,
}
2. S3 konfigurieren und eine Funktion zum Auflisten von TIFF-Dateien definieren
S3_ENDPOINT = (
os.environ.get("EODATA_S3_ENDPOINT") or input("EODATA S3 endpoint (https://...): ")
).strip().rstrip("/")
S3_BUCKET = (os.environ.get("EODATA_S3_BUCKET") or input("EODATA bucket [eodata]: ").strip() or "eodata").strip()
S3_REGION = (os.environ.get("EODATA_S3_REGION") or input("S3 signing region (blank for SDK default): ")).strip() or None
if not S3_ENDPOINT.startswith("https://"):
raise ValueError("Provide the HTTPS endpoint for your EODATA S3 service.")
# Use credentials authorized to list and download the selected EODATA objects.
# Environment variables are optional; otherwise enter the keys at hidden prompts.
access_key = os.environ.get("EODATA_ACCESS_KEY") or getpass("EODATA access key: ")
secret_key = os.environ.get("EODATA_SECRET_KEY") or getpass("EODATA secret key: ")
if not access_key or not secret_key:
raise ValueError("Provide both EODATA access keys authorized for the configured S3 service.")
s3 = boto3.client(
"s3",
aws_access_key_id=access_key,
aws_secret_access_key=secret_key,
endpoint_url=S3_ENDPOINT,
region_name=S3_REGION,
config=Config(s3={"addressing_style": "path"}),
)
def list_tiff_files_in_s3(prefix):
tiff_files = []
paginator = s3.get_paginator("list_objects_v2")
for page in paginator.paginate(Bucket=S3_BUCKET, Prefix=prefix):
for obj in page.get("Contents", []):
if obj["Key"].lower().endswith((".tif", ".tiff")):
tiff_files.append(obj["Key"])
return sorted(tiff_files)
3. Den Produktpfad aus der API-Antwort abrufen
response = requests.get(CATALOGUE_URL, params=query_params, timeout=90)
response.raise_for_status()
products = response.json().get("value", [])
if not products:
raise RuntimeError("No Sentinel-1 products matched the date range and area of interest.")
product = products[0]
product_name = product["Name"]
s3_path = product["S3Path"]
bucket_prefix = f"/{S3_BUCKET}/"
if not s3_path.startswith(bucket_prefix):
raise ValueError(f"The catalogue path {s3_path!r} does not match bucket {S3_BUCKET!r}.")
product_path = s3_path.removeprefix(bucket_prefix).rstrip("/") + "/"
print(product_name)
4. TIFF-Dateien aus S3 herunterladen und lokal speichern
# Keep each selected product in its own directory to avoid mixing acquisitions.
output_dir = Path("sen_1") / product_name
output_dir.mkdir(parents=True, exist_ok=True)
tiff_files = list_tiff_files_in_s3(product_path + "measurement/")
if not tiff_files:
raise RuntimeError("No measurement TIFF files were found for the selected product.")
for key in tiff_files:
file_name = Path(key).name
local_file_path = output_dir / file_name
s3.download_file(S3_BUCKET, key, str(local_file_path))
print(f"Downloaded: {file_name}")
Downloaded: s1a-iw-grd-vh-20240516t162754-20240516t162819-053899-068d12-002.tiff
Downloaded: s1a-iw-grd-vv-20240516t162754-20240516t162819-053899-068d12-001.tiff
Schritt 2: Städtische Gebiete mithilfe von Schwellenwerten abgrenzen
1. Daten einlesen und vorverarbeiten
Lesen Sie die Messdaten mit rasterio aus dem Produktverzeichnis des vorherigen Abschnitts ein. Die ursprünglichen Profile und gegebenenfalls vorhandene Bodenpasspunkte (GCPs) bleiben für den Export erhalten.
def find_files_with_names(folder_path):
files = sorted(Path(folder_path).glob("*"))
vv_files = [p for p in files if "-vv-" in p.name.lower() and p.suffix.lower() in (".tif", ".tiff")]
vh_files = [p for p in files if "-vh-" in p.name.lower() and p.suffix.lower() in (".tif", ".tiff")]
if len(vv_files) != 1 or len(vh_files) != 1:
raise ValueError("The product directory must contain exactly one VV and one VH measurement TIFF.")
return vv_files[0], vh_files[0]
input_file_vv, input_file_vh = find_files_with_names(output_dir)
print(input_file_vv)
with rasterio.open(input_file_vv) as src_vv:
vv_band = src_vv.read(1)
profile_vv = src_vv.profile
gcps_vv, gcps_crs_vv = src_vv.gcps
with rasterio.open(input_file_vh) as src_vh:
vh_band = src_vh.read(1)
profile_vh = src_vh.profile
gcps_vh, gcps_crs_vh = src_vh.gcps
if vv_band.shape != vh_band.shape:
raise ValueError("VV and VH measurements must have matching dimensions.")
# Retain the source profiles and ground control points; do not invent a CRS or transform.
2. Pixelwerte in Dezibel umrechnen
Nach dem Einlesen können Sie die Pixelwerte in Dezibel umrechnen. Die Dezibelskala (dB) macht Unterschiede in der Intensität der Radarreflexion in Sentinel-1-Daten deutlicher sichtbar. Das erleichtert die Analyse städtischer Gebiete mit ihren charakteristischen Reflexionseigenschaften.
Das größte Maximum im Histogramm der VV-Polarisation liegt im Bereich von 15–20 dB. Dies deutet darauf hin, dass die meisten Oberflächen im Untersuchungsgebiet eine mäßige Radarrückstreuung aufweisen. Der Schwellenwert von 23,5 dB hebt Pixel mit den höchsten Rückstreuwerten hervor. Diese können Gebieten mit dichter städtischer Infrastruktur oder anderen stark reflektierenden Strukturen entsprechen.
Das größte Maximum im Histogramm der VH-Polarisation liegt im Bereich von 13–18 dB. Dies deutet darauf hin, dass dieser Kanal empfindlicher auf Oberflächen mit geringerer Radarrückstreuung reagiert. Mit dem Schwellenwert von 21,5 dB für den VH-Kanal lassen sich Pixel mit höherer Rückstreuung erkennen. Dadurch können Geländestrukturen mit etwas geringerer, aber weiterhin deutlicher Rückstreuung besser unterschieden werden.
def convert_to_db(image):
# Create a copy of the image to avoid modifying the original array
image = np.array(image, dtype=np.float32, copy=True)
# Replace non-positive values with a small positive value to avoid taking log of zero or negative
image[image <= 0] = np.finfo(float).eps
# Calculate dB values
with np.errstate(divide='ignore', invalid='ignore'):
image_db = 10 * np.log10(image)
# Replace infinite values with the maximum finite value
max_value = np.nanmax(image_db[np.isfinite(image_db)])
image_db[np.isinf(image_db)] = max_value
return image_db
vv_band_db = convert_to_db(vv_band)
vh_band_db = convert_to_db(vh_band)
3. Daten analysieren und Schwellenwerte anwenden
Erstellen Sie als Nächstes Histogramme der Pixelwerte in Dezibel für die VV- und VH-Bänder. Die Histogramme zeigen die Verteilung der Radarreflexionsintensitäten und helfen dabei, die Eigenschaften des untersuchten Gebiets zu verstehen. Die Schwellenwerte von 23,5 dB für VV und 21,5 dB für VH beruhen auf Vorwissen oder auf den Anforderungen an die Analyse städtischer Gebiete. Sie dienen dazu, die relevanten Bereiche in den Bildern abzugrenzen. Das hier verwendete Verfahren setzt Pixelwerte oberhalb des Schwellenwerts auf 1 und darunterliegende Werte auf 0. So lassen sich Strukturen oder Gebiete mit den für städtische Umgebungen typischen Radarreflexionseigenschaften abgrenzen.
fig, axs = plt.subplots(1, 2, figsize=(18, 5))
# Histogram of VV Image in dB
axs[0].hist(vv_band_db.flatten(), bins=50, color='blue', alpha=0.7, label='VV dB')
axs[0].set_title('Histogram of VV Image in dB')
axs[0].set_xlabel('Radar Reflection Value (dB)')
axs[0].set_ylabel('Number of Pixels')
axs[0].legend()
axs[0].grid(True, linestyle='--', alpha=0.5)
# Histogram of VH Image in dB
axs[1].hist(vh_band_db.flatten(), bins=50, color='red', alpha=0.7, label='VH dB')
axs[1].set_title('Histogram of VH Image in dB')
axs[1].set_xlabel('Radar Reflection Value (dB)')
axs[1].set_ylabel('Number of Pixels')
axs[1].legend()
axs[1].grid(True, linestyle='--', alpha=0.5)
plt.show()
# Function to apply a threshold to an image
# Pixels greater than the threshold are set to 1, others to 0
def thresholding(image, threshold):
return (image > threshold).astype(np.uint8)
# Define threshold values for VV and VH bands
threshold_value_vv = 23.5
threshold_value_vh = 21.5
# Apply thresholding to the VV and VH band
vv_band_threshold = thresholding(vv_band_db, threshold_value_vv)
vh_band_threshold = thresholding(vh_band_db, threshold_value_vh)
# Combine the thresholded images using a logical AND operation
# The result is 1 where both VV and VH are above their respective thresholds
combined_threshold = np.logical_and(vv_band_threshold, vh_band_threshold).astype(np.uint8)
4. Bilder darstellen
In den Teilbildern können Sie die Unterschiede der Radarsignalreflexion zwischen den Polarisationen sichtbar machen. Die Originalbilder und die in Dezibel umgerechneten Bilder der VV- und VH-Polarisationen zeigen unterschiedliche Strukturen und Intensitäten in städtischen Gebieten. In den Originalbildern sind die Unterschiede in der Signalreflexion kaum erkennbar. Dies verdeutlicht den Nutzen der Dezibelskala. Wenn Sie die Informationen beider Polarisationen mithilfe von Schwellenwerten kombinieren, lassen sich bebaute Gebiete besser maskieren. Dadurch werden die städtischen Strukturen des Untersuchungsgebiets deutlicher sichtbar.
# Create a figure with a 3x2 grid of subplots, setting the overall figure size
fig, ax = plt.subplots(3, 2, figsize=(18, 18))
# Display the original VV image
ax[0, 1].imshow(vv_band, cmap='viridis')
ax[0, 1].set_title('Original VV Image')
# Display the VV image in decibels
ax[1, 1].imshow(vv_band_db, cmap='viridis')
ax[1, 1].set_title('VV Image in dB')
# Display the VV image after thresholding
ax[2, 1].imshow(vv_band_threshold, cmap='inferno')
ax[2, 1].set_title('VV Image after Thresholding')
# Display the original VH image
ax[0, 0].imshow(vh_band, cmap='viridis')
ax[0, 0].set_title('Original VH Image')
# Display the VH image in decibels
ax[1, 0].imshow(vh_band_db, cmap='viridis')
ax[1, 0].set_title('VH Image in dB')
# Display the VH image after thresholding
ax[2, 0].imshow(vh_band_threshold, cmap='inferno')
ax[2, 0].set_title('VH Image after Thresholding')
# Turn off the axes for all subplots to improve the visual presentation
for i in range(2):
for j in range(3):
ax[j, i].axis('off')
# Show the complete figure
plt.show()
# Display the combined threshold image using a specific colormap
plt.figure(figsize=(18, 9))
im = plt.imshow(combined_threshold, cmap='twilight_shifted')
plt.title('Combined Threshold Image')
cbar = plt.colorbar(im)
cbar.set_label('Signal Strength')
plt.axis('off')
plt.show()
Mit dieser Methode lassen sich bebaute Gebiete einfach maskieren. Im verarbeiteten Bild ist der Ballungsraum Warschau auf der rechten Seite deutlich zu erkennen. Entlang wichtiger Verkehrsachsen wie Bahnstrecken und Schnellstraßen zeigt sich eine sternförmige Konzentration der Bevölkerung. Auch kleinere städtische Zentren sind sichtbar. Zum Stadtzentrum hin nimmt die Zahl der maskierten Pixel zu. Dies weist auf eine dichtere Bebauung im Zentrum als in den Vororten hin.
5. Quantitative Analyse
Nach der Verarbeitung der Satellitenbilder können Sie eine quantitative Analyse durchführen. Die Eigenschaften der Satellitendaten ermöglichen es, einzelne Pixelwerte zu untersuchen. Dies ist für detaillierte und statistische Auswertungen entscheidend. Die Flächenberechnung nimmt Pixel von 10 × 10 Metern an und verwendet die Schwellenwertmaske als Näherung für bebaute Flächen. Das Ergebnis ist eine beispielhafte Schätzung und keine validierte Landbedeckungsklassifikation.
# Calculate the number of urban pixels by summing up the binary values in the combined threshold image
urban_pixels = np.sum(combined_threshold)
# Define the area covered by a single pixel (in square meters)
pixel_area = 10 * 10 # Assuming each pixel represents a 10m x 10m area
# Calculate the total urban area in square meters
urban_area = urban_pixels * pixel_area
# Calculate the total number of pixels in the combined threshold image
total_pixels = combined_threshold.size
# Calculate the percentage of urbanized area in the image
urban_percentage = (urban_pixels / total_pixels) * 100
# Convert urban area from square meters to hectares
urban_area_hectares = urban_area / 10_000
# Convert total area from square meters to hectares
total_area_hectares = (total_pixels * pixel_area) / 10_000
print(f"Percentage of Urbanized Area: {urban_percentage:.2f}%")
print(f"Urbanized Area: {urban_area_hectares:.2f} ha")
print(f"Total Image Area: {total_area_hectares:.2f} ha")
Percentage of Urbanized Area: 9.08%
Urbanized Area: 401728.18 ha
Total Image Area: 4426150.50 ha
6. Ergebnisse speichern
Die Ausgabedateien behalten die ursprünglichen Profile und Bodenpasspunkte. Diese Schritte führen keine Geländekorrektur oder Reprojektion der Messdaten durch.
output_file_db_vv = "vv_db_urban_warsaw.tif"
output_file_db_vh = "vh_db_urban_warsaw.tif"
output_file_threshold_combined = "warsaw_urban_threshold_combined.tif"
def write_result(path, data, profile, gcps, gcps_crs):
output_profile = profile.copy()
output_profile.update(dtype=data.dtype.name, count=1, nodata=None)
with rasterio.open(path, "w", **output_profile) as dst:
dst.write(data, 1)
if gcps:
dst.gcps = (gcps, gcps_crs)
write_result(output_file_db_vv, vv_band_db.astype(np.float32), profile_vv, gcps_vv, gcps_crs_vv)
write_result(output_file_db_vh, vh_band_db.astype(np.float32), profile_vh, gcps_vh, gcps_crs_vh)
write_result(output_file_threshold_combined, combined_threshold, profile_vv, gcps_vv, gcps_crs_vv)
Was als Nächstes zu tun ist
Der folgende Artikel verwendet ähnliche Verfahren zur Verarbeitung von Sentinel-5P-Daten zur Luftverschmutzung:
Verarbeitung von Sentinel-5P-Daten zur Luftverschmutzung mit Jupyter Notebook auf CODE-DE-Lab.
Weitere Informationen finden Sie im folgenden Artikel:
Analyse und Überwachung von Überschwemmungen mit Python und Sentinel-2 Daten auf CODE-DE-Lab.