Daten plotten und fitten

Inhalt

Daten plotten und fitten#

Zunächst werden die für dieses Jupyter Notebook benötigten Libraries geladen:

#Benötigte Libraries:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


import plotly.offline as py
py.init_notebook_mode(connected=True)
import plotly.graph_objs as go
import plotly.tools as tls
import seaborn as sns
import time
import warnings
warnings.filterwarnings('ignore')

# MatplotLib Settings:
plt.style.use('default') # Matplotlib Style wählen
plt.rcParams['font.size'] = 10; # Schriftgröße

.csv-Datei als DataFrame einlesen#

Im Folgenden Nutzen wir globale Klimadaten, die auf der Webseite der NASA zu finden sind: https://data.giss.nasa.gov/gistemp/. Hierbei handelt es sich um Temperaturdaten, die Anomalien gegenüber dem Mittelwert in den Jahren 1951-1980 aufweisen. Es werden Daten von Dateien (online oder offline) eingelesen mit der Python Bilbiothek pandas. Die Daten werden in sogenannten DataFrames hier mit dem Namen global_mean abgespeichert.

#link = "https://data.giss.nasa.gov/gistemp/graphs_v4/graph_data/Global_Mean_Estimates_based_on_Land_and_Ocean_Data/graph.csv"
link = 'data/graph.csv'
global_mean = pd.read_csv(link, header = 1) #DataFrame erstellen

Wir geben das DataFrame aus um uns die Messdaten einmal anzusehen:

global_mean.head(10) # Ausgabe der ersten 5 Spalten
#global_mean.tail(5) # Ausgabe der letzten 5 Spalten
#global_mean # Ausgabe des DataFrames 

	Year	No_Smoothing	Lowess(5)
0	1880	-0.17	-0.09
1	1881	-0.08	-0.13
2	1882	-0.11	-0.17
3	1883	-0.18	-0.20
4	1884	-0.28	-0.24
5	1885	-0.33	-0.26
6	1886	-0.32	-0.27
7	1887	-0.36	-0.27
8	1888	-0.17	-0.27
9	1889	-0.11	-0.26

In der ersten Spalte befinden sich lediglich die Indizes der Messungen. Die zweite Spalte beinhaltet das Jahr und die dritte Spalte zeigt den gemessenen globalen Temperaturunterschied im Vergleich zur gemittelten Temeratur der Jahre 1951-1980. Die letzte Spalte zeigt die gleichen Messwerte, jedoch gefiltert.

Einzelne Spalten kann man sich anzeigen lassen, indem den Spaltel-Namen des zugehörigen DataFrames nutzt:

global_mean['Year']

    1880
    1881
    1882
    1883
    1884
       ... 
  2017
  2018
  2019
  2020
  2021
Name: Year, Length: 142, dtype: int64

Daten plotten und Diagramm sichern mit ‘matplotlib’#

Als Beispiel für eine gelungene grafische Darstellung wollen wir die beiden Spalten, No_Smoothing and Lowess(5) gegenüber der Zeitachse Year plotten. Hierfür benützen wir die Python Library matplotlib. Einmal geplottet kann das zuletzt angezeigte Diagramm in verschiedenen Formaten mit plt.savefig('klima_plot1.png') abgespeichert werden. Wenn nicht anders angegeben, wird die Datei im gleichen Ordner angelegt.

import matplotlib.pyplot as plt

#plt.style.use('default')
plt.figure(figsize=(10,5))
plt.rcParams['font.size'] = 10; # Schriftgröße
plt.plot(global_mean["Year"],global_mean["No_Smoothing"], ls="--", lw=1, marker="s", ms=3, color="tab:gray", alpha=0.5, label="Werte");
plt.plot(global_mean["Year"],global_mean["Lowess(5)"], lw=4,  color="tab:blue", label="Glättung (NASA)");
plt.xlabel('Jahr')
plt.ylabel("Jahresmitteltemperaturabweichung [°C]")
plt.legend();
plt.grid();
plt.savefig('klima_plot1.png')
plt.savefig('klima_plot1.pdf')

../_images/6afb4b024b9c6d8ffb38dfe686701eb34f2015bc1770674cc38a416b0f853a91.png

Daten verarbeiten#

Die Bibliothek pandas ist sehr umfangreich und wird viel zur Datenverarbeitung genutzt. Im folgenden dazu einige Beispiele:

Statistische Größen: Mittelwert, Standardabweichung, Min, Max#

Für jede Spalte lassen sich statistische Größen wie z.B. die Anzahl der Einträge pro Spalte, deren Mittelwert, Standardabweichung, Minimal- und Maximalwert bestimmen:

global_mean.describe() #Dataframe Statistik

	Year	No_Smoothing	Lowess(5)
count	142.000000	142.000000	142.000000
mean	1950.500000	0.053169	0.053169
std	41.135953	0.364161	0.352324
min	1880.000000	-0.490000	-0.420000
25%	1915.250000	-0.200000	-0.227500
50%	1950.500000	-0.065000	-0.040000
75%	1985.750000	0.252500	0.235000
max	2021.000000	1.010000	0.930000

Maximum und zugehöriges Jahr ausgeben#

Dies können wir benutzen, um beispielsweise die maximale Temperaturdifferenz zu ermitteln:

global_mean["No_Smoothing"].std()

0.36416107752007404

Wenn wir wissen wollen, wann diese maximale Temperaturdifferenz auftrat, muss der zugehörige Index dieses Events gespeichert werden:

index_max = global_mean["No_Smoothing"].idxmax()
print(index_max)

Diesen Index können wir nun benutzen, um mittels .loc Befehl den Eintrag zu diesem Index auszugeben:

global_mean.loc[index_max]

Year            2016.00
No_Smoothing       1.01
Lowess(5)          0.87
Name: 136, dtype: float64

Daten aufsteigend/absteigend sortieren#

Mit .sort_values("Spaltenname") können wir auch die Tabelle nach dem definierten Spaltennamen sortieren (standardmäßig in ansteigender Reihenfolge):

global_mean.sort_values("No_Smoothing", ascending=False).head(10)

	Year	No_Smoothing	Lowess(5)
140	2020	1.01	0.93
136	2016	1.01	0.87
139	2019	0.98	0.93
137	2017	0.92	0.91
135	2015	0.90	0.83
138	2018	0.85	0.92
141	2021	0.84	0.93
134	2014	0.75	0.79
130	2010	0.72	0.65
133	2013	0.68	0.74

Mit ascending=False wird absteigend sortiert:

global_mean.sort_values("No_Smoothing", ascending = False)

	Year	No_Smoothing	Lowess(5)
140	2020	1.01	0.93
136	2016	1.01	0.87
139	2019	0.98	0.93
137	2017	0.92	0.91
135	2015	0.90	0.83
...	...	...	...
30	1910	-0.44	-0.42
31	1911	-0.45	-0.40
24	1904	-0.47	-0.31
37	1917	-0.47	-0.30
29	1909	-0.49	-0.41

142 rows × 3 columns

WICHTIG: Der Dataframe global_mean wird dadurch nicht verändert, es handelt sich nur um eine Anzeige!

Daten glätten#

Die von der NASA verwendete Glättung ist die LOcally WEighted Scatter-plot Smoother (LOWESS). Dabei wird in einem lokal zu definierenden Bereich eine lineare Regression durchgeführt. Eine genauere Erklärung zur Methode findet ihr auf Youtube.

Es gibt natürlich viele Methoden und Filter, um Daten zu glätten. Wir wollen nun versuchen, die Methode der NASA zu rekonstruieren. Hierfür benutzen wir die Python Library statsmodels und erstellen eine weitere Spalte Lowess(own) in unserem DataFrame global_mean. In diese Spalte schreiben wir die geglätteten Werte von den Rohdaten global_mean["No_Smoothing"] indem wir die Funktion lowess aufrufen. Details zu Nutzung der Funktion findet ihr https://www.statsmodels.org:

an erster Stelle in der Funktion werden die Y-Werte eingegeben, hier global_mean["No_Smoothing"]
an zweiter Stelle in der Funktion werden die X-Werte eingegeben, hier global_mean["Year"]
die Option frac ist eine Zahl zwischen 0 und 1. Dies ist der Anteil der Daten, der bei der Schätzung der einzelnen y-Werte verwendet wird.
Ausgegeben wird zweidimensionalas Array. Die erste Spalte enthält die sortierten x-Werte und die zweite Spalte die zugehörigen geschätzten y-Werte. Um die zweite Spalte in den DataFrame zu speichern, wählen wir diese mit [:,1] aus.

from statsmodels.nonparametric.smoothers_lowess import lowess

global_mean["Lowess(eigene)"] = lowess(global_mean["No_Smoothing"],global_mean["Year"], frac=1/14)[:,1]

global_mean

	Year	No_Smoothing	Lowess(5)	Lowess(eigene)
0	1880	-0.17	-0.09	-0.092698
1	1881	-0.08	-0.13	-0.129873
2	1882	-0.11	-0.17	-0.167172
3	1883	-0.18	-0.20	-0.203172
4	1884	-0.28	-0.24	-0.239029
...	...	...	...	...
137	2017	0.92	0.91	0.911976
138	2018	0.85	0.92	0.921410
139	2019	0.98	0.93	0.926072
140	2020	1.01	0.93	0.929365
141	2021	0.84	0.93	0.931468

142 rows × 4 columns

import matplotlib.pyplot as plt
plt.style.use('default')
plt.figure(figsize=(10,5))
plt.rcParams['font.size'] = 10;
plt.plot(global_mean["Year"],global_mean["No_Smoothing"], ls="-", lw=1, marker="s", ms=3, color="tab:gray", alpha=0.5, label="Werte");
plt.plot(global_mean["Year"],global_mean["Lowess(5)"], lw=4,  color="tab:blue", label="Glättung (NASA)");
plt.plot(global_mean["Year"],global_mean["Lowess(eigene)"], lw=2,ls='-' ,  color="tab:red", label="Glättung (eigene)");
plt.xlabel('Jahr')
plt.ylabel("Jahresmitteltemperaturabweichung [°C]")
plt.legend();
plt.grid();

../_images/f15818f6066a57cfe22942bd2b8ce751964ab2f2c07a53b9c156a56174860d37.png

Messunsicherheiten als Fehlerbalken hinzufügen#

Bei diesem Datenset stehen uns leider keine Messunsicherheiten zur Verfügung. Um Sie jedoch als Fehlerbalken miteinzubeziehen, wollen wir im Folgenden annehmen, dass der Temperaturunterschied auf 0.25K genau messen werden konnte und fügen die unserem Datensatz hinzu:

global_mean["uncertainty"] = 0.25
print(global_mean)

     Year  No_Smoothing  Lowess(5)  Lowess(eigene)  uncertainty
  1880         -0.17      -0.09       -0.092698         0.25
  1881         -0.08      -0.13       -0.129873         0.25
  1882         -0.11      -0.17       -0.167172         0.25
  1883         -0.18      -0.20       -0.203172         0.25
  1884         -0.28      -0.24       -0.239029         0.25
..    ...           ...        ...             ...          ...
2017          0.92       0.91        0.911976         0.25
2018          0.85       0.92        0.921410         0.25
2019          0.98       0.93        0.926072         0.25
2020          1.01       0.93        0.929365         0.25
2021          0.84       0.93        0.931468         0.25

[142 rows x 5 columns]

Grafisch darstellen tun wir Messunsicherheiten mittels Fehlerbalken und der Matplotlib-Funktion plt.errorbar.

plt.errorbar(global_mean["Year"],global_mean["No_Smoothing"], yerr=global_mean["uncertainty"], ls="-", lw=1, marker="s", ms=3, color="tab:gray", alpha=0.5, label="Werte");
plt.plot(global_mean["Year"],global_mean["Lowess(5)"], lw=4,  color="tab:blue", label="Glättung (NASA)");
plt.xlabel('Jahr')
plt.ylabel("Jahresmitteltemperaturabweichung [°C]")
plt.legend();
plt.grid();

../_images/3fa974c71aaf5b2595fc232bba03a51809ed6fe55134df6bd95e040713044876.png

Ausgleichsgerade berechnen und plotten#

Mittels linearer Regression kann der Temperaturanstieg aus den Daten berechnet werden. Hierfür wird die Python Library numpy benutzt und die Funktion polyfit aufgerufen und in als model gespeichert. Diese Funktion benutzt die Least-Square Methode für polynomische Modelle. Weitere Informationen zu der Funktion findet ihr hier. Mit der Option cov=True wird die Kovarianz-Matrix berechnet, welche die Unsicherheiten für die Fit-Parameter beinhaltet.

import numpy as np
import pandas as pd

x=global_mean["Year"]
y=global_mean["No_Smoothing"]
y_err = global_mean["uncertainty"]
model = np.polyfit(x, y, deg=1, w=1/y_err, cov=True) # 1. Wert = Anstieg , 2. Wert = Schnittpunkt mit y-Achse
y_model = model[0][0]*x+model[0][1] # Modell einer linearen Regression

plt.ylabel("Jahresmitteltemperaturabweichung [°C]")
plt.xlabel("Jahr")
plt.errorbar(global_mean["Year"],global_mean["No_Smoothing"], yerr=global_mean["uncertainty"], ls="-", lw=1, marker="s", ms=3, color="tab:gray", alpha=0.5, label="Werte");
plt.plot(x,y_model, ls="-", lw=3, color="tab:red", label=f"lineare Regression y=({model[0][0]*1000:.3f}+-{np.sqrt(model[1][0][0]*1000):.3f})1e-3*x+({model[0][1]:.3f}+-{np.sqrt(model[1][1][1]):.3f})");
plt.legend(fontsize=12);
plt.grid();

../_images/53df56b61c34b7931ff4e77a66103fc5d5c72a5be581bce4dce3522b6662db85.png

model

(array([ 7.72740894e-03, -1.50191421e+01]),
 array([[ 1.33257719e-07, -2.59919181e-04],
        [-2.59919181e-04,  5.07196269e-01]]))

Das Model beinhaltet zwei Matrizen:

model

(array([ 7.72740894e-03, -1.50191421e+01]),
 array([[ 1.33257719e-07, -2.59919181e-04],
        [-2.59919181e-04,  5.07196269e-01]]))

Im ersten Array stehen die Fit-Parameter der linearen Ausgleichsgeraden entsprechend der obigen Deklaration: y_model = model[0][0]*x+model[0][1]. Im zweiten Array, hier eine 2x2 Matrix, sind die Unsicherheiten in Form von der Kovarianz-Matrix dargestellt. Der Temperaturanstieg kann entsprechend ausgegeben werden:

print(f"Temperaturanstieg pro Jahr (von 1981 bis 2020): {model[0][0]:.3f}°C/Jahr")
print(f"Temperaturanstieg seit Beginn der Messung: {(y_model.iloc[-1]-y_model.iloc[0]):.3f}°C")

Temperaturanstieg pro Jahr (von 1981 bis 2020): 0.008°C/Jahr
Temperaturanstieg seit Beginn der Messung: 1.090°C

Warnung

Die lineare Regression bezieht hier den ganzen Zeitraum mit ein! Im folgenden betrachten wir für den Temperaturgradienten nur die Daten von 1980 bis 2020!

x=global_mean.loc[global_mean["Year"] >= 1980,"Year"]
y=global_mean.loc[global_mean["Year"] >= 1980,"No_Smoothing"]
y_err = global_mean.loc[global_mean["Year"] >= 1980,"uncertainty"]

model = np.polyfit(x, y, deg=1, w=1/y_err, cov=True) # 1. Wert = Anstieg , 2. Wert = Schnittpunkt mit y-Achse
y_model = model[0][0]*x+model[0][1] # Modell einer linearen Regression
print(f"Temperaturanstieg pro Jahr (von 1980 bis 2020): {model[0][0]:.3f}°C/Jahr")

Temperaturanstieg pro Jahr (von 1980 bis 2020): 0.019°C/Jahr

plt.ylabel("Jahresmitteltemperaturabweichung [°C]")
plt.xlabel("Jahr")
plt.errorbar(global_mean["Year"],global_mean["No_Smoothing"], yerr=global_mean["uncertainty"], ls="-", lw=1, marker="s", ms=3, color="tab:gray", alpha=0.5, label="Werte");
plt.plot(x,y_model, ls="-", lw=3, color="tab:red", label=f"lineare Regression y=({model[0][0]*1000:.3f}+-{np.sqrt(model[1][0][0]*1000):.3f})1e-3*x+({model[0][1]:.3f}+-{np.sqrt(model[1][1][1]):.3f})");
plt.legend(fontsize=12);
plt.grid();

../_images/d9f70b3fa56eaf3262f0fd2ac9a1fd68570121fed49db16d6fe29b78b138cfec.png

Beliebige Funktion fitten#

Mit der Polyfit Funktion können in erster Linie Polynome gefittet werden. Für beliebige Modellierung, wie z.B. Exponential-Funktionen, werden die Funktionen in aller Regel selber definiert. Hierfür eignet sich die scipy-Funktion curve_fit.

import pandas as pd
from scipy.optimize import curve_fit

Als erstes wird die Funktion deklariert, hier nutzen wird wieder eine lineare Regression:

def lin_reg(x, a, b):
    return a*x + b

def quadratic_fit(x, a, b, c):
    return a*x**2 + b*x + c

Die wir anschließend an die Daten anpassen:

popt, pcov = curve_fit(quadratic_fit, x, y)
print('Fit:', popt)
print('Mit der Kovarianz-Matrix', pcov)

Fit: [ 2.41606861e-04 -9.47677075e-01  9.29408439e+02]
Mit der Kovarianz-Matrix [[ 1.14888003e-08 -4.59666876e-05  4.59764893e-02]
 [-4.59666876e-05  1.83914056e-01 -1.83954621e+02]
 [ 4.59764893e-02 -1.83954621e+02  1.83996544e+05]]

plt.ylabel("Jahresmitteltemperaturabweichung [°C]")
plt.xlabel("Jahr")
plt.errorbar(global_mean["Year"],global_mean["No_Smoothing"], yerr=global_mean["uncertainty"], ls="-", lw=1, marker="s", ms=3, color="tab:gray", alpha=0.5, label="Werte");
plt.plot(x,quadratic_fit(x,*popt), ls="-", lw=3, color="tab:red", label="Quadratischer Fit y=( %5.3f * x^2 + %5.3f * x + %5.3f)" % tuple(popt));
plt.legend(fontsize=12);
plt.grid();

../_images/ac921a82538bc793bd62e8d95a3c52aa597b4ffea4662a4104975f668f63dd36.png