Ziel

Im Tutorial 2 von https://pv4ing.ch sind drei Übungen vorgesehen. Die erste Übung erstellt eine Zeitreihe mit fiktiven Daten. Die zweite Übung verwendet Messdaten als Zeitreihe und die dritte Übung führt eine Berechnung der Jahresstrahlungsenergie durch, bei unterschiedlichen Ausrichtungen.

Tutorial 2 - Übung 1: Einstieg

Zuerst machen wir eine einfache Zeitreihenrechnung.

Schritt 1 - Synthetische Messwerte für 1 Tag

Im folgenden Beispiel, nehmen wir keine realen Messwerte, sondern "bauen" uns selbst die Messwerte über die Sinusfunktion. Zuerst erstellen wir einen Array (Zahlenreihe) für die Zeit. Wir nehmen an, dass die Sonne um 6:00 aufgeht und um 18:00 untergeht. Über Mittag erre icht die Strahlung einen Wert von 650 W/m2.

In [1]:
# Bibliotheken Import
import numpy as np
import matplotlib.pyplot as plt
import datetime as dt

# Funktionsdefinition (benötigen wir später)
cos  = lambda arg : np.cos(np.deg2rad(arg))
sin  = lambda arg : np.sin(np.deg2rad(arg))
acos = lambda arg : np.rad2deg(np.arccos(arg))
asin = lambda arg : np.rad2deg(np.arcsin(arg))

t = np.linspace(6,18,1000)
H = 650 * sin((t-6)/(18-6)*180)

Strahlungsarray (Zahlenreihe für die Strahlung in W/m2). Beachte das die Funktion sin() oben mit dem lambda-Operator definiert wurde und in Grad angegeben wird. Dies ist nicht bei der Funktion np.sin() der Fall.

Überprüfe ob der Zeitarray aus lauter Spaltenwerte (Zeile mit vielen Zahlen) oder aus lauter Zeilenwerte (Spalte mit vielen Zahlen) besteht. Spalte = column. Zeile = row.

In [2]:
np.shape(t)
Out[2]:
(1000,)

Stelle den Strahlungsverlauf als figure mit der Linienfarbe grün (g=green) dar.

In [3]:
plt.figure(2, figsize=(8,4))
plt.plot(t,H, 'g')
plt.xlabel('Zeit [h]')
plt.ylabel('Strahlung [W/m2]')
Out[3]:
Text(0, 0.5, 'Strahlung [W/m2]')

Nun können wir die Strahlungsenergie $W_H$ für diesen einen Tag berechnen. Beachte das die Berechung der "Energie = Leistung mal Zeit" ist. Hier ist der Zeitschritt "deltaT".

Wir geben einen String aus "Strahlungs...." indem wir das Rechenergebnis darstellen. Das Ergebnis wird im Text mit {} angegeben. Dieser {}-Bereich wird ersetzt mit Zahlenwerten welche in .format() angegeben wird. In unserem Fall sind dies die beiden Werte A und wH. In der ersten {1:3.0f} geben wir an, dass wir das Element "wH" mit Index=1 wollen. Das erste Element "A" hat den Index=0. Dann folgt ein ":" und die Formatierung:

3.0f heisst eine Float-Zahl mit drei Stellen und 0 Nachkommastellen.

In [4]:
A = 1 # Fläche 1 m2
deltaT = t[1] - t[0]
wH = np.sum(H)*deltaT*A/1000
print('Strahlungsenergie {1:3.0f} kWh/m2/Tag auf {0:1.2f} m2 Fläche'.format(A,wH))
Strahlungsenergie   5 kWh/m2/Tag auf 1.00 m2 Fläche

Nun haben wir folgende Punkte gelernt:

  1. Erstellen eines Arrays.
  2. Darstellen von Arrays, inklusiv Beschriftung.
  3. Berechnen eines Flächenintegrals aus diskreten Werten, d.h. Aufsummieren von Werten.

Schritt 2 - Synthetische Messwerte für 7 Tag

Nun wollen wir die Strahlungsenergie von 7 Tagen berechnen, von Stunde 0 bis Stunde 7 x 24. Die Sinusfunktion liefert uns für die Nacht (18:00 bis 6:00) negative Werte. Wir korrigieren die negativen Strahlungswerte zu 0.

In [5]:
H[H<0] = 0

Der Ausdruck "H<0" generiert ein Selektor-Array der Länge von H, bei dem die Positionen mit "1" oder "True" gekennzeichnet sind, bei welchen der Wert < 0 ist. Der Ausdruck H[Selektor-Array] = 0 weisst allen selektierten Positionen eine 0 zu.

Denkbar ist auch "M[N<1] = 2". M und N müssen gleich lange Arrays sein. Im Array M werden alle Stellen auf 2 gesetzt, bei denen der Wert im Array N < 1 ist.

Tutorial 2 - Übung 2: Berechnung Strahlungsenergie von mehreren Tagen mit Messdaten

Nun verwenden wir nicht "künstliche" Strahlungswerte (generiert durch die Sinusfunktion), sondern reale Messwert. Hierzu benötigen wir eine Datendatei. Wir erstellen einen Plot über 7 Tage mit der Globalstrahlung "hGlob", der Diffusstrahlung "hDif" und der Direktstrahlung "hDir".

Bei Energiesimulationen ist es empfehlenswert als Kommentar die Einheit anzugeben, sodass verständlich ist um was es sich handelt, z.B. [grad] oder [radiant] oder [W] oder [kWh].

Verwende für diese Übung ein neues py-File und importiere wieder die benötigten Bibliotheken (numpy, matplotlib, datetime).

Datenimport mit Pandas

Zum Importieren verwenden wir die Bibliothek pandas. Diese zeichnet sich durch gute Verarbeitung von Daten aus, vorallem Zeitreihen. Das Einlesen geht mittels read_csv, dabei erhalten wir ein DataFrame (df). Dies ist vergleichbar mit einem Excel-Tabellenblatt mit einem Header, d.h. Überschriftzeile (Header) und Links einem laufenden Index.

Als Übung verwenden wir Daten von PVGIS über das Tool "Solar Radiation". Wir wählen einen Ort, geben die Ausrichtung an und laden die Daten als csv-Datei. pvgis

Sieh dir die csv-Datei mit dem Texteditor an. Sie enthält einige Information zu Beginn und am Ende. Bei Einlesen sollen dies Informationen ignoriert werden. Dies wird mit skiprows und skipfooter angegeben. Diese Funktionen sind mit engine = 'python' verfügbar. Standardmässig verwendet read_csv die 'c'-engine, da diese schneller ist, jedoch nicht so mächtig.

In [2]:
import pandas as pd

df = pd.read_csv('pvgis.csv', delimiter=',', decimal='.', skiprows = 8, skipfooter = 10, header = 0, engine = 'python')
df['tutc'] = pd.to_datetime(df['time'], format='%Y%m%d:%H%M' )

Pandas erkennt automatisch, welche Daten "strings" sind und welche Daten "Zahlen", d.h. float sind. Pandas erkennt aber nicht das Datumsformat. Hier helfen wir, indem wir das Format angeben und die Zeitwerte als maschinenlesbares Format in eine neue Spalte tutc speichern.

Nach dem Import, schau dir die importen Daten als DataFrame (df) an, mit

  • print(df) oder
  • list(df) oder
  • describe(df).

Zeitreihenrechnung

Für die nachfolgenden Systemberechnungen verwenden wir die csv-Datei auf github https://github.com/markstaler/pv4ing. Hier ist neben Strahlung und Temperatur, noch ein Lastprofil und Weiteres enthalten. Mit "save_as" kann die Datei gespeicher werden. Achtung, wenn eine csv-Datei wird standardmässig mit Excel geöffnet, dabei formatiert Excel die Zahlen und das Datum und die Uhrzeit mit den Einstellungen des Betriebssystems. Dies geht in 9 von 10 Fällen gut, aber nicht immer. Deshalb die Daten nicht speichern, beim Öffnen mit Excel. Besser die Datei mit dem Editor/Notepad öffnen.

Im DataFrame ist das Datum und die Uhrzeit angegeben, jedoch auch wieder als Text-String und für pandas nicht verständlich. Dieser Text-String wandeln wir in ein maschinenlesbares Format um, hier bei geben wir das Format des Datums und der Uhrzeit an. Anschliessend legen wir diese Spalte mit Datum und Uhrzeit als Index fest (set_index) und die ursprüngliche Spalte "Time" wird gelöscht.

Mit der Festlegung von datetime als Index weiss pandas dass es sich um eine Zeitreihe handelt und kann diese, sofern gewünscht in eine beliebige Auflösung umrechnen mit resample.

Im untenstehenden Beispiel wird die Auflösung (=Frequenz) mit 30 Minuten definiert durch "30T" und die Mittelwerte in diesem Zeotintervall berechnet. Die Formatierung ist definiert durch DateOffset. Anstatt des Mittelwerts könnte auch der Minimalwert, Maximalwert, Summenwert, u.s.w. berechnet werden - resample ist ein sehr mächtiges Instrument.

Zum Rechnen verwenden wir Numpy-Array, d.h. wir nehmen die Daten aus dem DataFrame (df) in dem wir über den Spaltennamen, in eckiger Klammer, auf die Spalte zugreifen und die Werte (values) auslesen:

In [6]:
## Daten importieren
df = pd.read_csv('2017_DataExport15min.csv', delimiter=';', decimal='.', header = 0)
df['Time'] = pd.to_datetime(df['Time'], format='%Y-%m-%d %H:%M:%S' ) 
df = df.set_index(['Time']) # Datetime als Index definieren
df = df.resample('30T').mean() # Die Daten werden "resampled" zu 30 min Werte

hGlo  = df['hGlo'].values    # [W/m2] 
hDif  = df['hDif'].values    # [W/m2]
tAmb  = df['Tamb'].values    # [°C]
pVer  = df['Pload'].values   # [W] Verbrauchsprofil
zapf  = df['Zapfung'].values # [l/15min] Profil für Warmwasserbezug

Für die Sonnenstandsberechnung benötigen wir einen Array mit laufenden Stunden "lfStd", beginnend mit Stunde 0 bei Jahresanfang um die Zeitreihenberechnung umzusetzten.

Nun kann es sein, dass die einglesenen Daten zu einem beliebigem Zeitpunkt beginnen, die laufenden Stunden sind aber mit 1.Januar definiert, deshalb berechnen wir die Zeitdifferenz (zDif) zum 1.Januar des Jahres vom ersten Messwert.

Die Zeitdifferenz (zDif) ist im Format Timedelta, welches wir in einen Float umwandeln, mit Sekunden und von dort in laufende Stunden umrechnen.

In [7]:
zDif = df.index - dt.datetime(df.index[0].year, 1, 1, 0)
lfStd = zDif.total_seconds().values/3600 

deltaT = lfStd[1] - lfStd[0] # [h] Auflösung
tutc = df.index

Schritt 1

Es soll die Globalstrahlung in grün, Direktstrahlung in rot und Diffusstrahlung in blau dargestellt werden. Weiter soll in der Plot-Figure eine Legende mit den Kurvennamen dargestellt werden. Als Ausschnitt wollen wir die ersten 7 Tage betrachten, d.h. in der Plotfunktion grenzen wir den Array ein durch tutc[a:e], wobei a für den Anfangsindex und e für den Endindex steht.

Schritt 2

Wenn wir den ersten Tag im Plot ansehen, so fällt auf, dass negative Strahlungswerte vorkommen. Dies entspricht den Messergebnissen aber nicht der Realität. Deshalb sollen die Strahlungsdaten bereinigt werden.  

In [8]:
hDir = hGlo - hDif

hDir[hDir<0] = 0
hDif[hDif<0] = 0
hGlo[hGlo<0] = 0

a=int(0)           # Anfangstag
e=int(7*24/deltaT) # Endtag
t = tutc[a:e]

plt.figure(3, figsize=(8,4))# Grösse des Plots (figsize) in Zoll
plt.plot(t,hGlo[a:e],'g', label ='Globalstrahlung')
plt.plot(t,hDir[a:e],'r', label ='Direktstrahlung')
plt.plot(t,hDif[a:e],'b', label ='Diffusstrahlung')
plt.xlabel('Zeit [h]');
plt.ylabel('Strahlung [W/m2]');
plt.legend(loc="upper left")
Out[8]:
<matplotlib.legend.Legend at 0x213aa67e9a0>

Schritt 3 - Dynamische Diagramme mit Bokeh

Für eine Visalisierung eines Strahlungsverlauf übers Jahr, wäre es nützlich wenn in die Daten "gezoomt" werden könnte. Dies ist möglich mit dem Paket "bokeh". Es funktionert ähnlich wie Matplotlib. Es wird eine html-Datei erzeugt (mittels "output_file('filename.html')"), welche das Diagramm enthält und Werkzeuge zum Verschieben, Zommen, usw.

In [9]:
from bokeh.io import show  
from bokeh.plotting import figure 

x = np.linspace(0,7,100)
y = np.sin(x)

p1 = figure()
p1.line(x,y)

show(p1)

Nun zur Maximalvariante einer Darstellung:

Die Übergabe der Daten erfolgt direkt als Dataframe (kann auch als Array erfolgen, wie bei Matplotlib). Im Dataframe erfänzen wir noch die Direktstrahlung und die Uhrzeit als String ("tutcStr") in einem, durch un, definierbaren Format.

Es werden die einzelnen Werkzeuge parametriert. "pan, box_zoom, reset, save" verstehen sich von selbst. Mit "hover" kann ein Tooltip ein/ausgeschalten werden. Tooltip heisst das die Daten der Kennlinie ausgegeben werden an der Position wo die Maus über der Kennlinie ist. Die Daten können die Messwerte sein mit Zugriff "@hGlo" oder der Wert aus der Kennlinie mit $y", welcher interpoliert wird.

Bei der Erzeugung des Diagramm-Objekt kann die Grösse eingestellt werden und die Parameter für Werkzeuge (tools) und Datenausgabe (tooltips).

Sehr schöne Funktion ist p1.legend.click_policy="hide" dabei können die Kurven ausgeblendet werden, wenn die jeweilige Legende angeklickt wird.

In [10]:
from bokeh.io import output_file, output_notebook, show  # Bokeh is for making charts in a browser
from bokeh.plotting import figure 
from bokeh.embed import components


# zum Dateframe ergänzen wir die Uhrzeit als Text
tmez = tutc + pd.Timedelta('1 hours') # Ursprüngliche Zeit in UTC wird in MEZ umgewandelt (+1)
tmezStr = tmez.strftime("%d.%m.%Y %H:%M")
 
df['tmez'] = tmez    
df['tmezStr'] = tmezStr
df['hDir'] = df['hGlo'].values - df['hDif'].values

output_file('myHTML.html') # das Diagramm wird in einem separaten html-File exportiert
output_notebook()   # das Diagramm wirm unten im JupyterNotebook dargestellt

tt = [
      ('Zeit MEZ:', '@tmezStr'),
      ('Globalstrahlung Diagramm:', '$y'),
      ('Globalstragkzbg Messwert:', '@hGlo'),
      ]

p1 = figure(title='Strahlungsverlauf', plot_width=800, plot_height=600, 
            tooltips = tt,  
            x_axis_type='datetime')
p1.line(x='tmez',y='hGlo', source=df, legend_label = 'Global', color="green")
p1.xaxis.axis_label = 'Zeit [MEZ]'
p1.yaxis.axis_label = 'Strahlung [W/m²]'
p1.legend.location = 'top_right'
p1.legend.click_policy="hide" # Kurve ein/ausschaltbar
p1.toolbar.logo = None # deaktivieren des Bokeh-Logo
show(p1)
Loading BokehJS ...