Kapitel-4-Spektrum.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kapitel 4 - Spektrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ziel dieses Kapitels ist es, sich das Signal nochmal aus einer etwas anderen Persepektive anzuschauen - unter anderem auch ein bisschen so, wie man sich Audio-Signale anschaut."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 Auto-Korrelation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die erste Frage, die wir uns anschauen wollen, ist, wie \"schnell\" sich die Temperatur mit der Zeit verändert, und wie schnell das Klimasystem \"vergisst\", wie warm es an einem bestimmten Zeitpunkt war. Man spricht dann auch davon, dass ein System ein \"Gedächtnis\" hat.\n",
    "\n",
    "Dazu müssen wir zunächst wieder die Daten laden:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Programmierbibliotheken importieren\n",
    "\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from scipy import optimize\n",
    "from scipy.signal import periodogram\n",
    "from statsmodels.tsa.stattools import acf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = pd.read_csv(\n",
    "    \"DeinDateiname\",  # \"produkt_klima_tag_12345678_12345678_01234.txt\"\n",
    "    sep=\";\",\n",
    "    skipinitialspace=True,\n",
    ")\n",
    "\n",
    "data = xr.DataArray(d[\"TMK\"], coords={\"time\": d[\"MESS_DATUM\"]})\n",
    "data[\"time\"] = pd.to_datetime(data.time, format=\"%Y%m%d\")\n",
    "data = data.sel(time=slice(np.datetime64(\"1972-01-01\"), None))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Um die Frage beantworten zu können müssen wir außerdem noch 2 Dinge tun:\n",
    "\n",
    "1) den Lauf der Jahreszeiten aus den Daten entfernen. Der würde sonst unsere Ergebnisse beeinflussen. Zwei Tage im Sommer sind sich untereinander viel ähnlicher als ein Tag im Sommer und einer im Winter - das würde unsere Analyse auch als \"Gedächtnis\" erkennen. Wenn wir nicht wissen würden, dass es einen Lauf der Jahreszeiten gäbe wäre das spannend für uns - wir könnten mit unserer Analyse herausfinden, dass es Jahreszeiten gibt! Aber da wir das schon wissen, schauen wir uns nur \"Anomalien\" an. Das heißt, statt zB die Temperatur am 1. April 2023 zu nutzen, rechnen wir aus wie viel wärmer oder kälter als ein durchschnittlicher 1. April der 1. April 2023 war. Und dasselbe machen wir mit allen anderen Tagen in dem Datensatz.\n",
    "\n",
    "2) den Klimawandel würden wir auch sehen in unserer Analyse. Dass es den gibt wissen wir auch schon. Deshalb versuchen wir den Einfluss des Klimawandels zu minimieren, indem wir die \"beste Gerade\" aus [Kapitel 3](Kapitel-3-Klimawandel.ipynb) von den Daten \"abziehen\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die \"rohen\" Daten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Daten, nachdem wir den Lauf der Jahreszeiten abgezogen haben:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "climatology = data.rolling(time=15, center=True).mean().groupby(\"time.dayofyear\").mean()\n",
    "anomaly = (data.groupby(\"time.dayofyear\") - climatology).unstack()\n",
    "anomaly.plot(label=\"Data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Man sieht eine Menge \"rauschen\" und einen leichten Aufwärtstrend - größtenteils durch den Klimawandel. Den versuchen wir jetzt auch zu entfernen: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x, m, b):\n",
    "    return m * x + b\n",
    "\n",
    "\n",
    "days_since_start = (anomaly.time - anomaly.time[0]) / np.timedelta64(1, \"D\")\n",
    "popt, pcov = optimize.curve_fit(f, days_since_start, anomaly.values)\n",
    "detrended_anomaly = anomaly - f(\n",
    "    (anomaly.time - anomaly.time[0]) / np.timedelta64(1, \"D\"), *popt\n",
    ")\n",
    "detrended_anomaly.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jetzt geht es nicht mehr Berg auf!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.stattools import acf\n",
    "\n",
    "auto_correlation = acf(detrended_anomaly, missing=\"raise\", nlags=100)\n",
    "lags = np.arange(len(auto_correlation))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Damit können wir jetzt unsere Analyse durchführen: Wir rechnen die sogenannte Autokorrelation aus (\"Auto\" ist Griechisch für \"Selbst\" - daher kommen übrigens auch die Namen für das \"Auto-Mobil\" und den \"Auto-mat\"): \n",
    "\n",
    "Um die Autokorrelation zu verstehen, stelle dir vor, du suchst dir einen beliebigen Tag aus dem Datensatz aus. Dann suchst du dir einen zweiten Punkt in einem bestimmten Abstand (zum Beispiel 7 Tage in der Zukunft). Die Autokorrelation gibt dann an, wie ähnlich sich der erste und zweite Datenpunkt dann im Durchschnitt sind.\n",
    "\n",
    "Eine hohe Autokorrelation bedeutet, dass beide Werte immer ähnlich sind. So hat zum Beispiel die Helligkeit in deinem Zimmer wahrscheinlich eine hohe Autokorrelation bei einem Zeitunterschied von 24 Stunden - denn egal welchen Zeitpunkt man sich als ersten Zeitpunkt aussucht, zum zweiten Zeitpunkt 24 Stunden danach ist es wieder die selbe Uhrzeit - und dementsprechend kommt ähnlich viel Licht durch das Fenster und hat man das Licht im Zimmer angeschaltet oder nicht. Und deshalb wären sich die Datenpunkte wahrscheinlich ähnlich.\n",
    "\n",
    "Wenn die Temperaturen (Temperaturanomalien) einfach gewürfelt wären, also komplett zufällig, dann wäre der Korrelations-Koeffizient 0, egal wie weit weg der zweite Datenpunkt vom ersten ist (außer du hast zwei mal den selben Punkt gewählt:D). Ist der Korrelationskoeffizient 1, dann lässt sich der zweite Datenpunkte perfekt aus dem ersten Vorhersagen.\n",
    "\n",
    "Wie sieht das also für unsere Temperaturanomalien aus?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(auto_correlation)\n",
    "plt.title(\"Auto-Korrelation\")\n",
    "plt.ylabel(\"Korrelations-Koeffizient\")\n",
    "plt.xlabel(\"Zeitlicher Abstand zwischen Datenpunkten [Tage]\")\n",
    "plt.xlim(0, 30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Was fällt dir auf? Welchen Wert hat der Korrelations-Koeffizient für einen Zeitunterschied von 0? Macht das Sinn, wenn du den Text oben liest?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deine Antwort:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Der Korrelations-Koeffizient ist überall größer als Null, besonders in den ersten Tagen. Das heißt, das Wettersystem hat ein \"Gedächtnis\"! Die Temperaturen an unterschiedlichen Tagen sind also nicht wie \"zusammengewürfelt\" - sondern wenn man einen recht kleinen Abstand zwischen zwei Tagen aussucht, dann ist das Wetter an diesen Tagen (im Durchschnitt) ähnlich!\n",
    "\n",
    "Das macht ja auch Sinn: angenommen, es gäbe keine guten Wetterberichte im Internet. Dann wäre eine einfache Methode um das Wetter vorherzusagen einfach \"Das Wetter morgen ist so wie das Wetter heute\" - und damit würde man in vielen Fällen gar nicht so schlecht liegen - zumindest für die Temperatur.\n",
    "\n",
    "\n",
    "Woher kommt das Gedächtnis? Die Antwort dazu ist lang - damit du eine grobe Idee bekommst: \n",
    "- Wettersysteme (zB Hoch- und Tiefdruckgebiete) sind groß und bewegen sich nicht sonderlich schnell - deshalb ändert sich die Temperatur an einzelnen Orten auch nur langsam, während sich die Wettersysteme über sie hinweg bewegen. Anders sieht das bei Niederschlägen aus: Wolken sind eher klein (im Vergleich zur gesamten Erde). Unter anderem deshalb hängen Niederschläge zeitlich nicht so sehr voneinander ab, und der Korrelationskoeffizient vor Niederschlagsmengen wären kleiner als der für die Temperatur, wenn man ihn in die Abbildung oben einzeichnen würde.\n",
    "- Manche Komponenten im Klimasystem haben eine Speicherwirkung: zum Beispiel können die Ozeane Wärme speichern und wieder abgeben, und die Erde kann Wasser speichern, dass beim Verdunsten das Wetter beeinflussen kann - durch veränderte Luftfeuchtigkeit und Verdunstungskühlung. Das führt dazu, dass Hitze- oder Kältewellen nicht einfach \"vergessen\" werden.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Leistungsdichte-Spektrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Als nächstes Schauen wir uns das Leistungsdichte-Spektrum des Klimas an. Das ist ein etwas kompliziertes Wort, aber im Prinzip beschreibt diese Größe nur, wie viel sich die Temperatur auf verschiedenen Zeitskalen verändert. Sind zum Beispiel die Temperaturänderungen zwischen einzelnen Tagen größer als die zwischen Wochen? \n",
    "\n",
    "Diese Art Spektrum kennst du vielleicht auch von Musikinstrumenten. Bei denen sieht man dann wenn man zB ein \"A\" spielt, dass sehr viel \"Leistung\" im Frequenz-bereich eines \"A\" Liegt (440 Hz), und dazu ein bisschen bei den Obertönen des \"A\" - das sind Resonanzen, die bei vielen Musik-Instrumenten mitklingen, und dadurch den Klang des Instruments ausmachen. Wenn dich das interessiert gibt es [hier](https://www.youtube.com/watch?v=CJVSg6vy1A0) ein kurzes Video, ist aber für den Rest nicht so wichtig.\n",
    "\n",
    "![Spektrum](spectrum.png) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Und wie sieht das jetzt für unsere Temperaturzeitreihe aus? Wir rechnen das aus, einmal für die rohen Daten und einmal für die Daten, bei denen wir die Jahreszeiten und den Langzeittrend entfernt haben."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_spectrum = periodogram(data)\n",
    "power_spectrum_detrended_deseasonalised = periodogram(detrended_anomaly)\n",
    "\n",
    "plt.plot(*power_spectrum, label=\"Mit Klimawandel und Jahreszeiten\")\n",
    "plt.plot(\n",
    "    *power_spectrum_detrended_deseasonalised,\n",
    "    label=\"Ohne Klimawandel und Jahreszeiten\",\n",
    ")\n",
    "plt.ylim([10**-10, 10**6])\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"Frequenz [1/Tage]\")\n",
    "plt.ylabel(\"Leistungsdichte\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Auf der y-Achse sieht man, wie viel sich die Temperatur an der Mess-Station auf welcher Zeitskala verändert. Lange Zeitskalen (= niedrige Frequenzen) sind links in der Abbildung und kurze Zeitskalen (= hohe Frequenzen) sind rechts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Man sieht zwei Unterschiede zwischen den beiden Kurven:\n",
    "\n",
    "1) Für sehr geringe Frequenzen (das heißt sehr lange Zeitskalen) ist ein Signal deutlich höher als das andere.\n",
    "2) Ein Signal hat eine Spitze bei einer Frequenz von etwa ($\\frac{1}{365}$ Tage) und das andere nicht.\n",
    "\n",
    "Hast du eine Idee, was die Unterschiede bedeuten könnten?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Idee für 1:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Idee für 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wenn wir uns nicht tägliche Durchschnittsdaten anschauen würden, sondern zB stündliche, würde man auch eine sehr hohe Spitze bei einer Frequenz von $\\frac{1}{1}$ Tag sehen. Kannst du dir vorstellen wieso?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Antwort:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An der Form des Spektrums lässt sich auch einiges über das Klimasystem ablesen und wie groß die Speicher darin sind - aber das reicht zu weit für unser Tutorial... mehr dazu lernt man im Studium :) Wir haben ja schon wirklich viel geschafft heute!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wenn du möchtest kannst du dir deine Notebooks speichern. Das geht am besten, indem du oben links auf \"File\" klickst und dann auf \"Save and Export Notebook As...\", und dort HTML aussuchst. Dann kannst du die Dateien herunterladen - die lassen sich mit Firefox, Chrome oder Safari öffnen. Und dort kann man sie dann auch drucken lassen :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Danke für's durcharbeiten - Ich hoffe es hat Spaß gemacht!\n",
    "\n",
    "![Earth](earth.png)\n",
    "\n",
    "(Bild: Eyes Roger auf Wikimedia Commons)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "student_project",
   "language": "python",
   "name": "python3"
  },
  "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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}