Übersicht
Web Feature Service (WFS)
Web Feature Services (WFS) sind eine weitere Möglichkeit um Zugriff auf Geodaten möglich zu machen. Im Gegenteil zu Web Map Services (WMS) gibt der WFS nur Objekte mit einer Geometrie und Attributen und keine ganzen Bilder zurück. Außerdem können die WFS-Objekte mit Hilfe eines räumlichen oder inhaltlichen Filters eingegrenzt werden.
Ground Motion Service
Boden- und Massenbewegungen zu erkennen, erfordert kontinuierliches Überwachung und stete Auswertung von aktuellen als auch historischen Daten. In diesem Service werden fortgeschrittene Techniken der differentiellen SAR-Interferometrie (DInSAR) auf Sentinel-1 Radar Satellitenbilder angewandt, um Boden- und Massenbewegungen abzuleiten.
In dem Service werden einerseits die Daten des European Ground Motion Service EGMS aufbereitet, andererseits kann auf den Subsidence and Landslide Monitoring Service Austria SuLaMoSA, Prototyp entwickelt in ASAP-14, zurückgegriffen werden. Beide Services liefern für auswertbare Punkte, d.h. Punkte, die über die Zeit ihre Radar-Rückstreueigenschaften nicht wesentlich verändern, zumindest Zeitreihen der Ost-West Verschiebung und der Hoch-Tief Bewegung. Die räumliche Auflösung beträgt für EGMS ~100x100 m² und ist momentan für den Zeitraum 2016-2020 vorhanden. Für SuLaMoSA können die Parameter der räumlichen und zeitlichen Abdeckung frei gewählt werden.
Anwendung
WFS in QGIS
Mit Hilfe eines Geoinformationssytem (GIS) können WFS-Daten angezeigt und analysiert werden. Das bekannteste GIS ist QGIS. Wie dieses eingerichtet wird, ist in diesem Artikel beschrieben.
Nachdem QGIS erfolgreich installiert und geöffnet wurde, wird auf der linken Seite im Browser nach WFS/OGC API - Features gesucht. Hier kann mit einem Rechtsklick -> new Connection eine neue WFS-Verbindung hinzugefügt werden.
Es öffnet sich ein neues Fenster. Hier kann ein Name für den Layer angegeben werden. Um auf den WFS des EODC zugreifen zu können, muss https://features.services.eodc.eu/ als URL eingetragen werden. Mit 'OK' wird die neue Verbindung hinzugefügt.
Nun ist die neue Verbindung unter WFS/OGC API - Features sichtbar. Per drag and drop in das Kartenfeld oder Rechtsklick -> Add Layer to Project wird der gewünschte Layer dem Projekt hinzugefügt.
WFS in Jupyter
Da WFS meist für größere Datensätze verwendet werden, ist QGIS nicht immer optimal. Dann ist es hilfreich, die Daten in ein Jupyter Notebook zu laden und sich dort anzeigen zu lassen. Mit Hilfe verschiedener Python Bibliotheken lassen sich die Daten dort genauso verarbeiten, wie in einem GIS. Zudem gibt es aber noch mehr Möglichkeiten z.B. Diagramme zu erzeugen.
Das im Folgenden verwendete Jupyter Notebook kann hier im Github des EODC abgerufen werden.
Falls der Link nicht funktioniert, kann es über https://github.com/eodcgmbh -> eodc-examples -> tutorials -> Ground_Motion_WFS.ipynb geöffnet werden.
Verwendete Bibliotheken
- os: Um betriebssystemabhängige Funktionen nutzen zu können
- geopandas: Einfachere Verarbeitung von Geodaten
- contextily: Um Karten aus dem Internet (z.B. OSM) abzurufen
- owslib.ogcapi.features: Damit kann der WFS-Layer gelesen werden
- rich.console: Hiermit können schön formatierte Texte in der Console geschrieben werden (nicht relevant)
- numpy: Bibliothek für mathematische Funktionen und multidimensionale Arrays
- matplotlib: Bibliothek zur Erstellung statischer, animierter und interaktiver Visualisierungen
- datetime: Zur Anzeige und Veränderung von Daten und Zeit
- matplotlib.dates: Hiermit können Daten und Zeit an matplotlib angepasst werden
Laden des WFS
EODC_OGCAPI_URL = 'https://features.dev.services.eodc.eu/'
eodc_ogcapi = Features(EODC_OGCAPI_URL)
feature_collections = eodc_ogcapi.feature_collections()
console.print(feature_collections)
Hier wird erst die URL des WFS eingelesen und mit Hilfe der owslib.ogcapi.features Bibliothek ein Features-Objekt des WFS erzeugt. Daraus werden die verschiedenen Layer in die feature_collection geladen und ausgegeben.
Ausgabe: 'adriatic_vessels', 'Vlbg_EGMS_L3_E', 'Vlbg_EGMS_L3_U', 'STATISTIK_AUSTRIA_DSR'
Daten nach gewünschter Ausdehnung filtern
bbox = [9.7549833,47.4537861,9.7616428,47.4581097]
field_items = eodc_ogcapi.collection_items(
"Vlbg_EGMS_L3_U",
bbox = bbox,
limit = 2000,
)
df = gpd.GeoDataFrame.from_features(field_items["features"], crs="EPSG:4326")
Es wird eine bounding box, ein Rechteck mit angegebener räumlicher Ausdehnung, erzeugt. Es werden alle Objekte aus dem feature Vlbg_EGMS_L3_U (Bodenbewegung Hoch/Tief) ausgewählt, welche innerhalb der bounding box liegen. Aus diesen Objekten wird dann ein Geodataframe erzeugt.
Sortieren der Daten
df_sorted = df[sorted(df.columns)]
df_sorted.head()
Der dataframe wird nach den Spalten (alphabetisch) sortiert und die ersten fünf Zeilen ausgegeben.
Anzeigen der Daten in einer Karte
ax = df[["geometry"]].plot(
facecolor="red", figsize=(12, 6)
)
cx.add_basemap(ax, crs=df.crs.to_string(), source=cx.providers.OpenStreetMap.Mapnik);
Hier werden die im Geodataframe geladenen Objekte in einer Karte angezeigt. Als Hintergrundkarte wurde die OpenStreetMap gewählt.
Anpassung des dataframe an NumPy
df_np = df_sorted.to_numpy()
columns_names=list(df_sorted)
acquisition_dates = columns_names[0:363]
acquisition_date_str = [f"{d[:4]}-{d[4:6]}-{d[6:8]}" for d in acquisition_dates]
acquisition_date_num = [datetime.strptime(d, "%Y-%m-%d") for d in acquisition_date_str]
Mit to_numpy wird aus dem dataframe ein NumPy array erzeugt. In acquisition_dates werden die Namen der ersten 364 Spalten gespeichert, welche die Daten der jeweiligen Aufzeichnung sind. Anschließend werden die Spaltennamen in ein Datumsformat konvertiert.
Erstellung des Diagramms
fig_EGMS_timeseries, axes = plt.subplots(nrows=1, ncols=1,facecolor='white', figsize=[10, 5])
for i in range(len(df_np)):
axes.plot(acquisition_date_num,df_np[i,0:363], ".-", markerfacecolor="w")
axes.xaxis.set_major_locator(mdates.MonthLocator(interval=6))
axes.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
plt.setp(axes.get_xticklabels(), rotation=30, ha="right")
axes.set_xlabel('Date')
axes.set_ylabel('EGMS up/down [cm]')
axes.grid(visible=True,axis='y')
Hier wird ein Diagramm erzeugt, welches die Hoch-/Tiefbewegung der einzelnen Objekte/Punkte im Laufe der Zeit zeigt. Jeder Punkt wird in einer anderen Farbe dargestellt und das Format des Datums auf der X-Achse noch einmal angepasst.