In [ ]:
%pylab inline
import pandas as pd
import re
In [ ]:
from ftplib import FTP
import zipfile
from cStringIO import StringIO # um nichts auf die Festplatte zu schreiben, benutzen wir StringIO, mit dem wir ein file-artiges Objekt bekommen, das alles als String speichert.
In [ ]:
# Adresse des FTP-Servers sowie des Verzeichnisses der Wetterdaten
server = "ftp-cdc.dwd.de"
filepath = "/pub/CDC/observations_germany/climate/daily/kl/historical/"
In [ ]:
# Wir loggen uns auf dem FTP-Server ein und wechseln in das Verzeichnis
ftp = FTP(server)
ftp.login()
ftp.cwd(filepath)
In [ ]:
# Wir durchsuchen das Verzeichnis nach der Datei mit dem Code 01957 für Halle-Kröllwitz
files = []
ftp.dir(files.append)
for line in files:
if "tageswerte_01957" in line:
halle_filename = line.split()[-1]
In [ ]:
# Wir erzeugen ein Pseudofile mit StringIO und speichern das Zip-File dort. Dann öffnen wir es mithilfe des zipfile Moduls
pseudofile = StringIO()
ftp.retrbinary("RETR {}".format(halle_filename), pseudofile.write)
zf = zipfile.ZipFile(pseudofile)
In [ ]:
# namelist zeigt uns, welche Dateien im Zip-Archiv gespeichert sind.
# Wir sind nur an produkt_klima_Tageswerte... interessiert
for f in zf.namelist():
if "produkt_klima_Tageswerte" in f:
data_filename = f
In [ ]:
# Unser erster Versuch, die Daten in Pandas einzulesen. Klappt so weit ganz gut, allerdings fällt auf, dass die Spaltenbezeichner
# noch Leerzeichen am Anfang enthalten
# Das resultierende Objekt ist übrigens ein sogenannter "DataFrame"
with zf.open(data_filename, "r") as f:
data = pd.read_csv(f, sep=";")
print data.columns
print type(data)
In [ ]:
# Indem wir den Delimiter mittels einer sogenannten Regular Expression beschreiben (delimiter=\s*;\s* bedeutet, dass alles was null bis beliebig
# viele Leerzeichen am Anfang, dann ein Semikolon, und am Ende null bis beliebig viele Leerzeichen hat als Delimiter betrachtet werden soll)
# verschwinden die Leerzeichen aus den Spaltennamen.
with zf.open(data_filename, "r") as f:
data = pd.read_csv(f, delimiter="\s*;\s*")
print data.columns
In [ ]:
# Um jetzt schön plotten zu können, tun wir zwei Dinge:
# 1) wir wandeln die Integers in der Spalte "MESS_DATUM", die in der Form 19510101 (stellvertretend für Jahr, Monat, Tag) auftreten
# in echte Datumsobjekte um
# 2) wir sagen Pandas, dass wir das Datum als unseren Default-Index benutzen möchten
data.MESS_DATUM = pd.to_datetime(data.MESS_DATUM, format="%Y%m%d")
data = data.set_index("MESS_DATUM")
In [ ]:
# Wenn wir uns nun die Lufttemperatur ausgeben lassen, wird als Index automatisch das Datum der Messung angezeigt
print data["LUFTTEMPERATUR"]
In [ ]:
# Nun sortieren wir die Daten etwas aus. Oft gibt es fehlerhafte Daten, bei denen die Temperatur auf -999 gesetzt ist
# Ähnlich wie bei Numpyarrays gibt ein boolscher Ausdruck bei Pandas eine Serie von boolschen Werten zurück
good_values = data.LUFTTEMPERATUR > -999
print type(good_values)
In [ ]:
# Die gesäuberten Daten können wir nun leicht plotten
# Da wir das Messdatum als Index ausgewählt haben, benutzt Pandas es automatisch als x-Achse
data_good = data[good_values] # Wählt nur die Zeilen aus, bei denen good_values True ist.
data_good.LUFTTEMPERATUR.plot()
plt.show()
In [ ]:
# Die Bestimmung der kältesten gemessenen Temperatur ist einfach
data_good.LUFTTEMPERATUR.min()
In [ ]:
# Um den Temperaturtrend zu bestimmen, fitten wir eine Gerade durch unsere Temperaturwerte
from scipy.optimize import curve_fit
def f(x, y0, m):
return y0+x*m
(y0, m), _ = curve_fit(f, range(data_good.shape[0]), data_good.LUFTTEMPERATUR.values)
In [ ]:
# Da der kleinste Zeitunterschied ein Tag ist, hat auch die Steigung die Einheit Kelvin pro Tag
# Um nun die Temperaturänderung pro Jahr zu bekommen, multiplizieren wir mit 365
print m*365
In [ ]:
snow = data_good[["SCHNEEHOEHE"]] # hier die extra Klammer beachten!
In [ ]:
snow["month"] = snow.index.month # hier fügen wir zusätzliche Spalten ein
snow["year"] = snow.index.year
snow_good = snow[snow.SCHNEEHOEHE >= -10] # und nehmen nur die Tage, an denen vernünftige Werte gemessen wurden
In [ ]:
# Wir gruppieren unsere Daten! Das sorgt dafür, dass wir für einen bestimmten Monat alle Tage und alle Jahre angezeigt bekommen
month_groups = snow_good.groupby("month")
In [ ]:
# Z.B. für den Januar:
januaries = month_groups.get_group(1)
print type(januaries)
# Wie wir sehen ist januaries wieder ein DataFrame
print januaries
In [ ]:
# Aber eigentlich sind wir ja an den Februar Daten interessiert
allthefebruaries = month_groups.get_group(2)
In [ ]:
# Die Februare gruppieren wir jetzt nochmal nach Jahr. So bekommen wir Jahresgruppen, die jeweils nur ihren Februar enthalten
for year, month in allthefebruaries.groupby("year"):
print "this is year", year
month.SCHNEEHOEHE.plot()
plt.show()
In [ ]:
data_good["weekday"] = data_good.index.weekday
weekdaygroup = data_good[data_good.NIEDERSCHLAGSHOEHE >= 0].groupby("weekday")
#data_good[["NIEDERSCHLAGSHOEHE"]].groupby("day")
In [ ]:
weekdays = weekdaygroup.aggregate(np.mean)
weekdays_var = weekdaygroup.aggregate(np.var)
In [ ]:
niederschlag = weekdays[["NIEDERSCHLAGSHOEHE"]]
niederschlag_var = weekdays_var[["NIEDERSCHLAGSHOEHE"]]
wd_labels = ["Mo", "Tu", "We", "Th", "Fr", "Sa", "Su"]
#plt.bar(niederschlag.index, niederschlag.NIEDERSCHLAGSHOEHE)
plt.errorbar(niederschlag.index, niederschlag.NIEDERSCHLAGSHOEHE, np.sqrt(niederschlag_var.NIEDERSCHLAGSHOEHE))
plt.xticks(niederschlag.index, wd_labels)
plt.show()
In [ ]:
data_good["month"] = data_good.index.month
weekdaygroup = data_good[data_good.SONNENSCHEINDAUER >= 0].groupby("month")
weekdays = weekdaygroup.aggregate(np.mean)
weekdays_var = weekdaygroup.aggregate(np.var)
niederschlag = weekdays[["SONNENSCHEINDAUER"]]
niederschlag_var = weekdays_var[["SONNENSCHEINDAUER"]]
#wd_labels = ["Mo", "Tu", "We", "Th", "Fr", "Sa", "Su"]
#plt.bar(niederschlag.index, niederschlag.NIEDERSCHLAGSHOEHE)
plt.errorbar(niederschlag.index.values, niederschlag.SONNENSCHEINDAUER, np.sqrt(niederschlag_var.SONNENSCHEINDAUER.values))
#plt.errorbar(niederschlag.index, niederschlag.SONNENSCHEINDAUER, niederschlag_var.SONNENSCHEINDAUER)
#plt.xticks(niederschlag.index, wd_labels)
plt.show()
np.sqrt(niederschlag_var.SONNENSCHEINDAUER)
In [ ]: