Luftmesswerte mit Python holen und darstellen (Software)

Martin Vogel , Dortmund / Bochum, Mon, 11.03.2019, 20:52 (vor 2147 Tagen)

Das Landesamt für Natur, Umwelt und Verbraucherschutz NRW (LANUV) stellt die Daten seiner 63 Luftschadstoffmessstationen allen interessierten Bürgerinnen und Bürgern über das Internet zur Verfügung. Gemessen werden Feinstaub bis 10µm Partikelgröße (PM10), Stickstoffdioxid (NO₂), Stickstoffmonoxid (NO), Ozon (O₃) sowie gelegentlich einige andere Werte.

Von aktuellem Interesse sind diese Werte, weil davon Fahrverbote insbesondere für Dieselfahrzeuge abhängig gemacht werden sollen. Die Annahme der Umweltpolitik ist, dass dadurch die gefährliche Feinstaub- und NO₂-Belastung gesenkt werden kann.

Belastungen aus Fahrzeugverkehr erkennt man im Tagesverlauf an einer typischen „Kamelhöckerkurve” mit Höchstwerten zu den Stoßzeiten zwischen 7 und 8 Uhr morgens und 17 bis 19 Uhr nachmittags.
[image]

Eine spannende Entdeckung war es für mich, dass der abendliche Höchstwert im Jahresmittel beim Stickstoffdioxid deutlich über dem morgendlichen Wert liegt, obwohl der Stickstoffmonoxidwert abends meistens unter den morgens gemessenen Werten liegt. Noch besser: das abendliche NO₂-Maximum liegt im Jahresmittel deutlich später als zur Zeit des heimkehrenden Berufsverkehrs.

Sind nach Feierabend etwa besonders viele Trecker und Lastwagen unterwegs? Aus persönlicher Erfahrung wissen wir, dass diese Erklärung eher abseitig wäre. Zum Glück können wir auf Millionen von Messwerten zugreifen, in denen wir nach einer Antwort auf unsere Fragen suchen können.

Beim Geodatenportal NRW kann man die stündlich erfassten Messwerte aller Stationen für die letzten 365 Tage und teilweise sogar für die letzten Jahre herunterladen und auswerten.

Optimalerweise lässt man das von einem Python-Skript erledigen. Das kann die Stundenwerte eines ganzen Jahres dann auch gleich als 3D-Diagramm darstellen. Die x-Achse (unten links) zeigt die Stunden jedes Tages und die y-Achse (unten rechts) zeigt den jeweiligen Tag des Jahres (1 für den 1. Januar und 365 für den 31. Dezember).
[image]

Hier fallen zwei Bereiche im Frühling und Herbst ganz besonders auf, die extreme NO₂-Belastungen in den späten Abendstunden von etwa 19 bis 23 Uhr zeigen. Der Verdacht liegt nahe, dass es sich hierbei gar nicht um Verkehrsbelastungen handelt, sondern um Holzöfen und Kamine, die als extreme Feinstaub- und NO₂-Schleudern bekannt sind. Leider werden die Feinstaubmesswerte in NRW durch Mittelung über 24 Stunden unbrauchbar gemacht, sodass keine Zuordnung zu einer Emissionsquelle mehr möglich wird, die Stickstoffdioxidkurven sprechen jedoch für sich.

An Messpunkten, die berüchtigt für ihren Verkehrsdreck sind, wie dem Graf-von-Galen-Ring in der Straßenschlucht hinter dem Finanzamt in Hagen, sind die Messwerte tatsächlich tagsüber am höchsten.
[image]

An den meisten anderen Orten jedoch sind die typischen Abend- und Nachtspitzen im Frühjahr und Herbst deutlich zu erkennen. Eine komplette Sammlung von Grafiken habe ich als Google-Photo-Album abgelegt.

Wer eigene Experimente mit den Daten durchführen will, darf sich gerne mein Python-Programm zur grafischen Auswertung herunterladen. Die Grafiken werden in der aktuellen Version nicht auf dem Bildschirm dargestellt, sondern direkt als einzelne 1000 × 500 Pixel große PNG-Grafiken auf die Platte geschrieben.

#!/usr/bin/env python3
 
"""
Auswertung der Luftschadstoffmessungen der NRW-Stationen anhand der
stündlichen Messwerte der letzten 365 Tage.
 
Das Programm benötigt Python 3 ab Version 3.6 und Matplotlib.
"""
 
from urllib.request import urlopen
from html import unescape
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
 
# Liste der Stationen holen
url="https://www.lanuv.nrw.de/fileadmin/lanuv/luft/immissionen/aktluftqual/"\
     "eu_luft_akt.htm"
htmldata = urlopen(url).readlines()
 
Station=[]
# Daten aus der HTML-Tabelle auf der Seite zusammensuchen
for Zeile in htmldata:
    data = unescape(Zeile.decode("windows-1252")).split("<td")
    if len(data)>2:
        Name = data[1].split("<")[0][1:]
        Kürzel = data[2].split("<")[0][1:]
        Station.append((Kürzel, Name))
 
# Zahl der Tage seit Jahresbeginn (ohne Schaltjahr)
def tag(t,m):
    return t+sum([0,31,28,31,30,31,30,31,31,30,31,30,31][:m])
 
# Kontrollausgabe der gelesenen Stationsliste
print("Messwerte folgender Stationen sind verfügbar:\n")
print("Nr. Kürzel Stationsname")
for i,nk in enumerate(Station):
    Kürzel, Stationsname = nk
    print(f"{i:2}: {Kürzel:^6} {Stationsname}")
 
    # CSV-Datei der Station lesen
    url = "https://www.opengeodata.nrw.de/produkte/umwelt_klima/"\
          "luftqualitaet/luqs/konti_nach_station/"\
          f"OpenKontiLUQS_{Kürzel}_aktuell.csv"
 
    # laufende Summe über alle Tage eines Jahres
    pm10Summe = [0 for i in range(24)]
    noSumme   = [0 for i in range(24)]
    no2Summe  = [0 for i in range(24)]
 
    # 3D-Daten
    x = []
    y = []
    z = []
 
    with urlopen(url) as csv:
        zeilen = csv.readlines()
        # Der Aufbau der CSV-Dateien der einzelnen Stationen ist nicht
        # identisch. Immerhin gibt es in allen Dateien eine Kopfzeile,
        # die uns sagt, in welcher Spalte unsere Messwerte zu finden sind.
        posno = posno2 = pospm10 = None
        no = no2 = pm10 = None
        for i,j in enumerate(zeilen[1].decode("windows-1252").split(";")[2:]):
            if "NO2" in j:
                posno2 = i
            elif " NO " in j:
                posno = i
            elif "PM10" in j:
                pospm10 = i
 
        # Schleife über alle stündlichen Messungen der letzten 365 Tage:
        for zeile in zeilen[2:]:
            # Datum, Uhrzeit, Messwerte (z. B. für NO, NO₂ und PM10)
            d,u,*mw = zeile.decode("windows-1252").strip().split(";")
 
            # Stunde, Tag und Monat dieser Messung
            h = int(u.split(":")[0])
            t = int(d.split(".")[0])
            m = int(d.split(".")[1])
 
            # In welchen Spalten sind Werte, die uns interessieren?
            if posno is not None: no = mw[posno]
            if posno2 is not None: no2 = mw[posno2]
            if pospm10 is not None: pm10 = mw[pospm10]
 
            # NO₂-Messwert für 3D-Diagramm eintragen (falls vorhanden)
            try:
                z.append(int(no2)) # kann schiefgehen
                x.append(h)
                y.append(tag(t,m))
            except (TypeError, ValueError):
                # Kein Messwert
                pass
 
            # Stundensummen für 2D-Diagramm aufaddieren:
 
            # Feinstaub (derzeit überall eine horizontale Linie, denn bei
            # allen Stationen in NRW wurden die Messwerte durch
            # 24-Stunden-Mittelwertbildung unbrauchbar gemacht!)
            try:
                pm10Summe[h-1] += int(pm10)
            except (TypeError, ValueError):
                pass
 
            # Stickstoffmonoxid, NO
            try:
                noSumme[h-1] += int(no)
            except (TypeError, ValueError):
                pass
 
            # Stickstoffdioxid, NO₂
            try:
                no2Summe[h-1] += int(no2)
            except (TypeError, ValueError):
                pass
 
    # Jahressummen auf Tageswerte zurückrechnen
    for i in range(24):
        pm10Summe[i] /= 365
        no2Summe[i] /= 365
        noSumme[i] /= 365
 
    # Diagramme
 
    # figsize 10,5 ergibt 1000 × 500 Pixel in der PNG-Datei
    fig = plt.figure(figsize=(10,5))
 
    # 2D-Plot der Stundenmittelwerte
    ax = fig.add_subplot(1, 2, 1) # 1 Zeile, 2 Spalten: Spalte 1
 
    ax.plot(range(1,25),noSumme,"r",label="NO")
    ax.plot(range(1,25),no2Summe,"b",label="NO2")
    ax.plot(range(1,25),pm10Summe,"g",label="PM10")
    plt.legend()
    plt.xlim(1,24)
    plt.xticks(range(0,25,1))
    plt.grid(color="lightgrey")
    plt.xlabel("Uhrzeit")
    plt.ylabel("µg/m³")
    plt.title(f"Stundenmittelwerte Station {Kürzel}\n{Stationsname}")
 
    #3D-Plot der Stundenwerte übers ganze Jahr
    if no2:
        ax = fig.add_subplot(1, 2, 2, projection='3d') # Spalte 2
 
        surf = ax.plot_trisurf(x,y,z, cmap=cm.jet, linewidth=1)
        fig.colorbar(surf, shrink=0.5, aspect=5)
        plt.xlim(1,24)
        plt.xticks(range(0,25,3))
        plt.yticks(range(1,367,61))
        ax.set_xlabel("Uhrzeit")
        ax.set_ylabel("Tag des Jahres")
        ax.set_zlabel("µg/m³")
        plt.title(f"NO₂-Werte Station {Kürzel}\n{Stationsname}")
 
 
    # Teilplots ohne unnötigen Rand anordnen
    plt.subplots_adjust(left=0.001, bottom=0.001,
                        right=0.999, top=0.999,
                        wspace=0.001, hspace=0.001)
    plt.tight_layout()
    # plt.show()
    plt.savefig(f"{Kürzel}.png")
    plt.close()
 
print("Fertig.")