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