Luftmesswerte mit Python holen und darstellen (Software)

Martin Vogel ⌂ @, Dortmund / Bochum, Mon, 11.03.2019, 20:52 (vor 1873 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.

[image]

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.

  1. #!/usr/bin/env python3
  2.  
  3. """
  4. Auswertung der Luftschadstoffmessungen der NRW-Stationen anhand der
  5. stündlichen Messwerte der letzten 365 Tage.
  6.  
  7. Das Programm benötigt Python 3 ab Version 3.6 und Matplotlib.
  8. """
  9.  
  10. from urllib.request import urlopen
  11. from html import unescape
  12. import matplotlib.pyplot as plt
  13. from mpl_toolkits.mplot3d import Axes3D
  14. from matplotlib import cm
  15.  
  16. # Liste der Stationen holen
  17. url="https://www.lanuv.nrw.de/fileadmin/lanuv/luft/immissionen/aktluftqual/"\
  18. "eu_luft_akt.htm"
  19. htmldata = urlopen(url).readlines()
  20.  
  21. Station=[]
  22. # Daten aus der HTML-Tabelle auf der Seite zusammensuchen
  23. for Zeile in htmldata:
  24. data = unescape(Zeile.decode("windows-1252")).split("<td")
  25. if len(data)>2:
  26. Name = data[1].split("<")[0][1:]
  27. Kürzel = data[2].split("<")[0][1:]
  28. Station.append((Kürzel, Name))
  29.  
  30. # Zahl der Tage seit Jahresbeginn (ohne Schaltjahr)
  31. def tag(t,m):
  32. return t+sum([0,31,28,31,30,31,30,31,31,30,31,30,31][:m])
  33.  
  34. # Kontrollausgabe der gelesenen Stationsliste
  35. print("Messwerte folgender Stationen sind verfügbar:\n")
  36. print("Nr. Kürzel Stationsname")
  37. for i,nk in enumerate(Station):
  38. Kürzel, Stationsname = nk
  39. print(f"{i:2}: {Kürzel:^6} {Stationsname}")
  40.  
  41. # CSV-Datei der Station lesen
  42. url = "https://www.opengeodata.nrw.de/produkte/umwelt_klima/"\
  43. "luftqualitaet/luqs/konti_nach_station/"\
  44. f"OpenKontiLUQS_{Kürzel}_aktuell.csv"
  45.  
  46. # laufende Summe über alle Tage eines Jahres
  47. pm10Summe = [0 for i in range(24)]
  48. noSumme = [0 for i in range(24)]
  49. no2Summe = [0 for i in range(24)]
  50.  
  51. # 3D-Daten
  52. x = []
  53. y = []
  54. z = []
  55.  
  56. with urlopen(url) as csv:
  57. zeilen = csv.readlines()
  58. # Der Aufbau der CSV-Dateien der einzelnen Stationen ist nicht
  59. # identisch. Immerhin gibt es in allen Dateien eine Kopfzeile,
  60. # die uns sagt, in welcher Spalte unsere Messwerte zu finden sind.
  61. posno = posno2 = pospm10 = None
  62. no = no2 = pm10 = None
  63. for i,j in enumerate(zeilen[1].decode("windows-1252").split(";")[2:]):
  64. if "NO2" in j:
  65. posno2 = i
  66. elif " NO " in j:
  67. posno = i
  68. elif "PM10" in j:
  69. pospm10 = i
  70.  
  71. # Schleife über alle stündlichen Messungen der letzten 365 Tage:
  72. for zeile in zeilen[2:]:
  73. # Datum, Uhrzeit, Messwerte (z. B. für NO, NO₂ und PM10)
  74. d,u,*mw = zeile.decode("windows-1252").strip().split(";")
  75.  
  76. # Stunde, Tag und Monat dieser Messung
  77. h = int(u.split(":")[0])
  78. t = int(d.split(".")[0])
  79. m = int(d.split(".")[1])
  80.  
  81. # In welchen Spalten sind Werte, die uns interessieren?
  82. if posno is not None: no = mw[posno]
  83. if posno2 is not None: no2 = mw[posno2]
  84. if pospm10 is not None: pm10 = mw[pospm10]
  85.  
  86. # NO₂-Messwert für 3D-Diagramm eintragen (falls vorhanden)
  87. try:
  88. z.append(int(no2)) # kann schiefgehen
  89. x.append(h)
  90. y.append(tag(t,m))
  91. except (TypeError, ValueError):
  92. # Kein Messwert
  93. pass
  94.  
  95. # Stundensummen für 2D-Diagramm aufaddieren:
  96.  
  97. # Feinstaub (derzeit überall eine horizontale Linie, denn bei
  98. # allen Stationen in NRW wurden die Messwerte durch
  99. # 24-Stunden-Mittelwertbildung unbrauchbar gemacht!)
  100. try:
  101. pm10Summe[h-1] += int(pm10)
  102. except (TypeError, ValueError):
  103. pass
  104.  
  105. # Stickstoffmonoxid, NO
  106. try:
  107. noSumme[h-1] += int(no)
  108. except (TypeError, ValueError):
  109. pass
  110.  
  111. # Stickstoffdioxid, NO₂
  112. try:
  113. no2Summe[h-1] += int(no2)
  114. except (TypeError, ValueError):
  115. pass
  116.  
  117. # Jahressummen auf Tageswerte zurückrechnen
  118. for i in range(24):
  119. pm10Summe[i] /= 365
  120. no2Summe[i] /= 365
  121. noSumme[i] /= 365
  122.  
  123. # Diagramme
  124.  
  125. # figsize 10,5 ergibt 1000 × 500 Pixel in der PNG-Datei
  126. fig = plt.figure(figsize=(10,5))
  127.  
  128. # 2D-Plot der Stundenmittelwerte
  129. ax = fig.add_subplot(1, 2, 1) # 1 Zeile, 2 Spalten: Spalte 1
  130.  
  131. ax.plot(range(1,25),noSumme,"r",label="NO")
  132. ax.plot(range(1,25),no2Summe,"b",label="NO2")
  133. ax.plot(range(1,25),pm10Summe,"g",label="PM10")
  134. plt.legend()
  135. plt.xlim(1,24)
  136. plt.xticks(range(0,25,1))
  137. plt.grid(color="lightgrey")
  138. plt.xlabel("Uhrzeit")
  139. plt.ylabel("µg/m³")
  140. plt.title(f"Stundenmittelwerte Station {Kürzel}\n{Stationsname}")
  141.  
  142. #3D-Plot der Stundenwerte übers ganze Jahr
  143. if no2:
  144. ax = fig.add_subplot(1, 2, 2, projection='3d') # Spalte 2
  145.  
  146. surf = ax.plot_trisurf(x,y,z, cmap=cm.jet, linewidth=1)
  147. fig.colorbar(surf, shrink=0.5, aspect=5)
  148. plt.xlim(1,24)
  149. plt.xticks(range(0,25,3))
  150. plt.yticks(range(1,367,61))
  151. ax.set_xlabel("Uhrzeit")
  152. ax.set_ylabel("Tag des Jahres")
  153. ax.set_zlabel("µg/m³")
  154. plt.title(f"NO₂-Werte Station {Kürzel}\n{Stationsname}")
  155.  
  156.  
  157. # Teilplots ohne unnötigen Rand anordnen
  158. plt.subplots_adjust(left=0.001, bottom=0.001,
  159. right=0.999, top=0.999,
  160. wspace=0.001, hspace=0.001)
  161. plt.tight_layout()
  162. # plt.show()
  163. plt.savefig(f"{Kürzel}.png")
  164. plt.close()
  165.  
  166. print("Fertig.")

Lizenzinformationen: Die Daten des Geodatenportals NRW stehen unter der Datenlizenz Deutschland Namensnennung 2.0.

--
Dipl.-Ing. Martin Vogel
Leiter des Bauforums

Bücher:
CAD mit BricsCAD
Bauinformatik mit Python

Tags:
Python, Feinstaub, NOx

RSS-Feed dieser Diskussion
powered by my little forum