Analyse von modellierten und gemessenen Luftqualitätsdaten
Inhalt
In diesem Jupyter-Notebook wird gezeigt, wie COSMO-MUSCAT-Simulationsdaten geladen und plottet werden und wie Messdaten zur weiteren Auswertung genutzt werden können.
Jupyter Notebook ist eine Browser-basierte Programmierlösung, um wissenschaftliche Daten mit der Programmiersprache
Pythonzu analysieren.Python ist eine Programmiersprache, die man leicht erlernen kann
COSMO-MUSCAT ist ein Model für Meteorologie und Luftqualität.
Importieren von Modulen
Python arbeitet mit Modulen, die importiert werden, um bestimmte Aufgaben einfacher zu machen
[ ]:
import numpy as np # Standard Bibliothek für mathematische Operationen
import pandas as pd # Standard Bibliothek für Datenverarbeitung
import cartopy.crs as ccrs # Bibliothek für Karten / Kartographie
import cartopy.feature as cfeature
import xarray as xr # Einlesen von Daten (speziell gut für wissenschaftliche Daten im netCDF-Format )
from matplotlib import pyplot as plt # Standard Bibliothek für Plots
from netCDF4 import Dataset
Daten
Daten einlesen
Die Modelldaten werden aus einer Datei mit der Endung .nc eingelesen. Dieses Dateiformat heißt NetCDF und wird häufig für wissenschaftliche Daten verwendet – zum Beispiel in der Wetter- und Klimaforschung. Für diese Übung nutzen wir die Datei ‘model_data.nc’.
Mit der Python-Bibliothek xarray können solche Dateien einfach geöffnet und die enthaltenen Daten weiterverarbeitet werden.
[ ]:
# Hier legen wir den Pfad zur NetCDF-Datei fest.
# Das './' am Anfang bedeutet, dass sich die Datei im gleichen Ordner wie dieses Notebook befindet.
datei_pfad = './model_data.nc'
[ ]:
# Nun lesen wir die NetCDF-Datei mit der xarray-Bibliothek (abgekürzt als 'xr') ein.
# Die eingelesenen Daten werden in der Variable 'model_daten' gespeichert.
model_daten = xr.open_dataset(datei_pfad)
Daten anschauen
[ ]:
# Einen ersten Überblick, was wir nun als Variable 'model_daten' eingelesen haben:
model_daten
Datenstruktur
Der Datensatz besteht aus mehreren Koordinaten und Messgrößen. Jede Messgröße wird jeweils für jede Kombination der Koordinaten angegeben (Zeit- und Raumbezug).
Koordinaten:
time (1128 Zeitpunkte)
rlon, rlat (313 × 200 Gitterpunkte in relativen Längen- und Breitengraden)
Datenvariablen:
TEMP: Temperatur in K (Form: [time, rlat, rlon])
RAIN: Niederschlagsmenge
SO2 : Schwefeldioxid
NO2 : Stickstoffdioxid
PM25: Feinstaub PM2.5-Konzentration
PM10: Feinstaub PM10-Konzentration
Du kannst die jeweiligen Variablen aufklappen und die Informationen dazu ansehen.
Koordinaten
Wir können nun einzelne Koordinaten aufrufen. Dazu rufen wir wieder die Variable ‘model_daten’ auf und geben die gewünschte Koordinate nach dem ‘.’ an.
[ ]:
model_daten.time
Die Modeldaten bestehen aus 833 Zeitschritten mit jeweils stündlichen Daten für März 2021.
Die Koordinaten lat und lon geben die Längen- und Breitengrade jeder Modelgitterzelle an.
Karte erstellen
[ ]:
# Wir wählen einen Zeitpunkt aus, den wir uns anschauen möchten und die Variable (Datenvariablen: z.B.: PM25, PM10).
tag = '2021-03-15'
variable = 'PM25'
# Alle Stunden dieses Tages auswählen
tagesdaten = model_daten[variable].sel(time=tag)
# Median über den Tag berechnen
tagesmedian = tagesdaten.median(dim='time')
[ ]:
# Koordinaten auslesen
lon = model_daten['lon'].values
lat = model_daten['lat'].values
lon2d, lat2d = np.meshgrid(lon, lat)
# Festlegen der Projektion:
# Karten und Globusdarstellungen müssen immer eine Projektion nutzen, um die Erdoberfläche (3D) in eine 2D-Karte zu überführen.
# Unsere Daten liegen in echten geografischen Koordinaten (Breiten- und Längengrad), deshalb verwenden wir hier die sogenannte "PlateCarree"-Projektion.
# Diese stellt geografische Koordinaten direkt und einfach auf einer flachen Karte dar.
geo_crs = ccrs.PlateCarree()
# Nun können wir die eigentliche Karte zeichnen
plt.figure(figsize=(10, 6))
ax = plt.axes(projection=geo_crs)
ax.set_title(f'{variable} – Tagesmedian am {tag}')
# Wir fügen noch Küsten und Ländergrenzen ein
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':', edgecolor='gray')
# Und die Daten in die Karte plotten
mesh = ax.pcolormesh(lon2d, lat2d, tagesmedian, transform=geo_crs,
cmap='BuPu', shading='auto')
# Jetzt noch die Farbskala festlegen
cbar = plt.colorbar(mesh, ax=ax)
einheit = 'μg m⁻³'
cbar.set_label(f'{variable} [{einheit}]')
plt.tight_layout()
# Grafik speichern (optional)
# plt.savefig(f'./Karte_{variable}_{tag}.png', format="png", dpi=300, bbox_inches='tight', transparent=False)
plt.show()
Zeitserie erstellen
Wir können uns auch einen Ort aussuchen und eine Zeitserie der Daten erstellen
[ ]:
# Dazu müssen wir Längen- und Breitengrad des gewünschten Ortes angeben:
ort = 'Leipzig'
lat_ort = 51.34
lon_ort = 12.37
# und die Variable, die wir uns ansehen wollen:
variable = 'PM25'
[ ]:
# Den Gitterpunkt finden, der dem gewünschten Ort am nächsten liegt
# Die Methode `.sel(..., method="nearest")` sucht den nächsten Wert im Gitter
daten_ort = model_daten[variable].sel(lat=lat_ort, lon=lon_ort, method="nearest")
# Zeitachse extrahieren (für den Plot)
zeit = model_daten['time'].values
# Plot erstellen
plt.figure(figsize=(10, 5))
plt.plot(zeit, daten_ort, label=f'{variable} in {ort}')
plt.title(f'Zeitreihe von {variable} am Ort {ort}')
plt.xlabel('Zeit')
plt.ylabel(f'{variable} [μg m⁻³]')
plt.grid(True)
plt.tight_layout()
# Grafik speichern (optional)
# plt.savefig(f'./Zeitreihe_{variable}_{ort}.png', format="png", dpi=300)
plt.show()
Vergleich mit Messdaten
Um Einzuschätzen, wie gut unser Model funktioniert, vergleichen wir es mit Messdaten. Dazu können wir Daten des EEA Messnetzes nehmen, das ihr in der Vorbereitung bereits recherchiert habt.
Ihr könnte euch eine Messstation aussuchen, die euch Interessiert und die Daten für die Zeit unseres Modellaufs herunterladen. https://eeadmz1-downloads-webapp.azurewebsites.net/
Wähle als Dataset: ‘Primary validated data (E1a)’ aus und lade die Daten als Parquet files herunter. Speichere sie dann unter deinen eigenen Dokumenten ab. Unten auf der Seite können auch die Metadaten abgerufen werden. Metadata (deutsch: Metadaten) sind „Daten über Daten“. Sie beschreiben die Eigenschaften, den Kontext und die Herkunft von eigentlichen Messdaten – also z. B. wann, wo, wie und womit eine Messung durchgeführt wurde. Hier könnt ihr zum Beispiel die Koordinaten eurer ausgewählten Messstation nachschauen.
Einlesen der Messdaten
[ ]:
# Die heruntergeladenen EEA Daten liegen im Dateiformat .parquet vor
# Wir legen zuerst fest, wo ihr die Daten gespeichert habt, hier müsste ihr den Pfad zu euren Daten anpassen:
eea_daten_pfad = './SPO.DE_DESN059_PM2_dataGroup1.parquet'
[ ]:
# Die Daten werden nun eingelesen
eea_daten = pd.read_parquet(eea_daten_pfad)
# Wir können uns die ersten Zeilen der eingelesen Daten zur Kontrolle ansehen
print(eea_daten.head())
Messstation
[ ]:
# Name und Koordinaten der ausgewählten Messstation müssen hier händisch eingetragen werden.
# Somit kann aus den Modeldaten die richtigen Werte zum Vergleich ausgewählt werden.
ort = 'Leipzig'
lat_ort = 51.34
lon_ort = 12.37
# Variable für die EEA Daten vorliegen:
variable = 'PM25'
Zeitserie erstellen
[ ]:
# Nur gültige Daten auswählen (Validity == 1)
eea_daten = eea_daten[eea_daten['Validity'] == 1]
# Zeit und Messwerte extrahieren
zeit_messung = pd.to_datetime(eea_daten['Start'])
messung = eea_daten['Value']
# Modell-Daten am nächstgelegenen Punkt extrahieren
model = model_daten[variable].sel(lat=lat_ort, lon=lon_ort, method='nearest')
zeit_model = model_daten['time'].values
# Plot: Model vs. Messung
plt.figure(figsize=(12, 5))
plt.plot(zeit_model, model, label=f'Model – {variable}', color='blue')
plt.plot(zeit_messung, messung, label=f'EEA-Messung – {variable}', color='orange', alpha=0.7)
plt.title(f'{variable} – Vergleich Model vs. EEA-Messung ({ort})')
plt.xlabel('Zeit')
plt.ylabel(f'{variable} [μg m⁻³]')
plt.legend()
plt.grid(True)
plt.tight_layout()
# Grafik speichern (optional)
# plt.savefig(f'PM25_Vergleich_{ort}.png', dpi=300)
plt.show()