💾 Archived View for bacaliu.de › 20220727-regenradar.txt captured on 2023-07-10 at 13:51:47.
-=-=-=-=-=-=-
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
REGENRADAR
DWD-Daten zu Diagramm mit Python-Bokeh
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
2022-07-27
Inhaltsverzeichnis
──────────────────
1. Der Plan
2. Das Vorgehen
.. 1. Laden der Daten
.. 2. Speichern als .pickle
.. 3. Laden der Arrays und Timestamps
.. 4. Stereographische Umwandlung
.. 5. Fehlerprüfung: richtige Koordinaten?
.. 6. Endlich: das Plotten
3. Bibliography
Nav
Footer
1 Der Plan
══════════
Der [DWD] stellt auf <https://opendata.dwd.de> auch einige Radar-Daten
zum Download bereit. Diese möchte ich nutzen, um eine kurzfristige
Niederschlagsprognose zu erstelen. Insbesondere für das Planen von
Radtouren oder Wanderungen ist es von Vorteil 5-minütige Daten zu
haben, statt nur 1-Stündige des MOSMIX, die zudem nur alle 6 Stunden
aktualisiert werden.
Bei der Recherche stoße ich auf Bokeh (Van de Ven). Vorteile gegenüber
Matplotlib sehe ich mehrere:
⁃ Geschwindigkeit: Um einige Matplotlib-Figuren abzuspeichern werden
mehrere Sekunden benötigt, Bokeh schafft das Erstellen in einem
Bruchteil der Zeit, weil erst am Endgerät gerendert wird
• daher habe ich auch `mosmix.py' auf Bokeh umgestellt
⁃ Interaktivität: Der User kann zoomen, verschieben etc. und die
Achsbeschriftungen „wandern mit“, im gegensatz dazu ein langes
`.png' in einem x-scrollable-DIV zu stecken.
[DWD] <https://www.dwd.de>
2 Das Vorgehen
══════════════
Zwecks dessen schreibe ich ein Skript, welches dieses
darstellt. Folgendermaßen stelle ich *schematisch* dar, wie das
funktioniert.
2.1 Laden der Daten
───────────────────
┌────
│ def reload_data():
│ link = "https://opendata.dwd.de/weather/radar/composit/rv/DE1200_RV_LATEST.tar.bz2"
│
│ os.system(f"rm {PATH}/DE1200*")
│ os.system(f"wget {link} -O {os.path.join(PATH, 'data.bz2')}")
│ os.system(f"tar -xf {os.path.join(PATH, 'data.bz2')} -C {PATH}")
└────
Programmlisting 1: Herunterladen
Die Daten liegen hinter dem angegebenen link. Man lösche zunächst
vorhandene Daten, die alle mit "DE1200" beginnen, nicht aber die
später im gleichen Ordner erzeugten `.pkl'-Dateien, damit auch
/während/ des neuladens Diagramme erstellt werden können. Dann lade
man den link herunter und extrahiere die `.tar.gz'. Es entstehen für
jede 5 Minuten eine Datei, die einfach `DE1200{TIMESTAMP}' heißt.
2.2 Speichern als .pickle
─────────────────────────
┌────
│ def bin_to_array(f, shape):
│ start = len(f) - (shape[0]*shape[1])*2
│ array = np.zeros((shape[0], shape[1]))
│ index = start
│ for i in reversed(range(shape[0])):
│ for j in range(shape[1]):
│ num = (f[index+1] % 16)*256+f[index]
│ des = f[index]//16
│ if des >= 8:
│ array[i, j] = None
│ else:
│ num = num/100*12
│ if des >= 4 and num:
│ num *= -1
│ array[i, j] = num
│ index += 2
│ return array
│
│
│ def load_arrays():
│ files = os.listdir(PATH)
│ files = filter(lambda f: 'DE1200' in f, files)
│ files = list(files)
│ files.sort()
│
│ arrays = []
│ timestamps = []
│ for n, file in enumerate(tqdm(files, desc="lese Binärdaten")):
│ num = int(file.split("_")[2])
│ tstr = re.findall("[0-9]{10}", file)[0]
│ dt = datetime.strptime("20"+tstr, "%Y%m%d%H%M")
│ dt += timedelta(minutes=num)
│ timestamps.append(dt)
│ with open(os.path.join(PATH, file), "rb") as f:
│ f = f.read()
│ array = bin_to_array(f, (1200, 1100))
│ arrays.append(array)
│
│ arrays = np.array(arrays)
│
│ with open(os.path.join(PATH, "timestamps.pkl"), "wb") as f:
│ pickle.dump(timestamps, f)
│
│ with open(os.path.join(PATH, "arrays.pkl"), "wb") as f:
│ pickle.dump(arrays, f)
└────
Programmlisting 2: Umwandeln und Abspeichern
`bin_to_array' stellt hier eine Hifsfunktion dar um die Radardaten
aufzubereiten. Diese liest die Binärdaten ein und führt einfache
Kalkulationen durch um z.B. das Vorzeichen umzukehren wenn ein
bestimmtes Bit gesetzt ist. Näheres dazu findet man im [PDF zu den
Radarprodukten des DWD].
Diese Funktion wird von `load_arrays' für jede der Dateien aufgerufen
wird um ein Array zu erhalten. Es werden dann jeweils Array und
Zeitstempel (aus dem Dateinamen) in Listen geschrieben und letztlich
in `timestamp.pkl' und `arrays.pkl' gepickelt.
[PDF zu den Radarprodukten des DWD]
<https://www.dwd.de/DE/leistungen/radolan/radolan_info/radolan_radvor_op_komposit_format_pdf.pdf?__blob=publicationFile&v=5>
2.3 Laden der Arrays und Timestamps
───────────────────────────────────
Wie es sich für ein gutes Programm gehört, werden getter und setter
genutzt.
┌────
│ def get_arrays():
│ with open(os.path.join(PATH, "arrays.pkl"), "rb") as f:
│ return pickle.load(f)
└────
Programmlisting 3: get_arrays
┌────
│ def get_timestamps():
│ with open(os.path.join(PATH, "timestamps.pkl"), "rb") as f:
│ timestamps = pickle.load(f)
│ timestamps = [UTC.localize(t).astimezone(CET).replace(tzinfo=None) for t in timestamps]
│ return timestamps
└────
Programmlisting 4: get_timestamps
In `get_timestamps' wandle ich zwar die Zeiten in CET um, aber ersetze
die `tzinfo' des `datetime'-Objekts jeweils durch `None'. Grund dafür
ist, dass `Bokeh' nicht besonders gut mit Zeitzonen umgeht, und ich
mir den Stress erspare. Das Radar ist nur für Abdeckung in
Deutschland, also reicht es auch nur dessen Zeitzone zu verwenden.
2.4 Stereographische Umwandlung
───────────────────────────────
Der Benutzer wird seine Koordinaten in lat/lon angeben oder durch
eingabe eines Suchtermes über (“Nominatim” 2022) welche
generieren. Diese müssen nun aber umgewandelt werden in die Pixel, die
das Rader bietet. Auch hier half o.g. PDF-Dokument des DWD, um eine
Funktion zu definieren.
┌────
│ def latlon2xy(lat, lon):
│ R = 6370.04
│ M = lambda p: (1+np.sin(60/180*np.pi))/(1+np.sin(p/180*np.pi))
│
│ p = lat
│ l = lon
│
│ x = R*M(p)*np.cos(p/180*np.pi)*np.sin((l-10)/180*np.pi)
│ y = -R*M(p)*np.cos(p/180*np.pi)*np.cos((l-10)/180*np.pi)
│ return x - (-543.4628665604781), 1200 - y -4808.644645330335
└────
Programmlisting 5: latlon2xy
2.5 Fehlerprüfung: richtige Koordinaten?
────────────────────────────────────────
Um Fehler aufgrund unerlaubter Koordinaten zu vermeiden rufe ich bei
jeder Nutzereingabe `is_valid_koordinate' auf.
┌────
│ def is_valid_koordinate(lat, lon):
│ x, y = latlon2xy(lat, lon)
│ arrays = get_arrays()
│ x_max, y_max = arrays[0].shape
│
│ return all([x>=0, x<x_max, y>=0, y<y_max])
└────
Programmlisting 6: valide Koordinaten?
Die Arrays werden geladen und deren Shape gibt Auskunft darüber ob die
gewählten Koordinaten im Bereich liegen.
2.6 Endlich: das Plotten
────────────────────────
┌────
│ def bok_rain(lat=51.85, lon=7.49):
│ arrays = get_arrays()
│ timestamps = get_timestamps()
│
│ x, y = latlon2xy(lat, lon)
│ radius=[1, 2, 4]
│
│ max_intense = [
│ a[
│ max(int(y-max(radius)), 0):min(int(y+max(radius)), 1200),
│ max(int(x-max(radius)), 0):min(int(x+max(radius)), 1100)
│ ].max()
│ for a in arrays
│ ]
│
│ plot = figure(
│ title=f"Niederschlagsradar für {lat}ºN/{lon}ºE",
│ x_range=(
│ timestamps[0],
│ timestamps[-1]
│ ),
│ x_axis_type="datetime",
│ tools=["save", "xpan", "xwheel_zoom", "reset", "xzoom_in", "xzoom_out"],
│ active_scroll="xwheel_zoom",
│ y_range=[0, max(max(max_intense)*1.01, 1)],
│ )
│
│ try:
│ for radius, alpha in zip(
│ [1, 2, 4],
│ [1, 0.7, 0.4]
│ ):
│ intense = [
│ a[
│ max(int(y-radius), 0):min(int(y+radius), 1200),
│ max(int(x-radius), 0):min(int(x+radius), 1100)
│ ].max()
│ for a in arrays
│ ]
│
│ plot.line(
│ np.array(timestamps),
│ intense,
│ legend_label=f"Radius von {radius}km",
│ line_color="#458588",
│ line_alpha=alpha
│ )
│
│ plot.line(
│ np.array(timestamps),
│ [
│ sum(intense[:i])/12
│ for i, _ in enumerate(intense)
│ ],
│ line_dash=[4, 4],
│ line_color="#458588",
│ line_alpha=alpha/2
│ )
│ except ():
│ print("AAAAAAAAAAAAAAAAAAAA no Error-Handling")
│
│ plot.width = 1024
│ plot.height = 384
│ plot.background_fill_alpha = 0
│ plot.border_fill_alpha = 0
│ plot.sizing_mode = "stretch_width"
│
│ plot.outline_line_color = "gray"
│ plot.title.text_color = "gray"
│
│ plot.legend.label_text_color = "gray"
│ plot.legend.background_fill_color = "gray"
│ plot.legend.background_fill_alpha = 0.125
│
│ for axis in [plot.xaxis, plot.yaxis]:
│ axis.axis_label_text_color = "gray"
│ axis.axis_line_color = "gray"
│ axis.major_label_text_color = "gray"
│ axis.major_tick_line_color = "gray"
│ axis.minor_tick_line_color = "gray"
│
│ for grid in [plot.xgrid, plot.ygrid]:
│ grid.grid_line_color = "gray"
│ grid.grid_line_alpha = 0.5
│
│ script, divs = components(plot, CDN)
│ return script, divs
└────
Programmlisting 7: plotten
Die tolle Funktion mit dem Namen `bok_rain' (im gegensatz zu
`plot_rain' eine Implementation mit `Bokeh', daher der Name) ruft nun
die Daten auf und berechnet zunächst die maximale intensität im
maximalen Radius über die gesamte Zeit.
Die Radien dienen dazu, auch Regen in der Umgebung berücksichtigen zu
können; zieht etwa in 80 Minuten ein sehr Starker und kleiner Schauer
in 2km vorbei, so sollte man dennoch mit Regen rechnen, denn er könnte
wachsen.
Dieser Radius wird dann als obere Schranke für das Diagramm
eingestellt, minimal aber 1 (in Liter pro Quadratmeter und Stunde).
Bokeh geben wir als Tools nur die x-zoom varianten, denn ein Zoom in
alle Richtungen lässt einen schnell den Plot verlieren. Dann wird
versucht (und bei Fehlern passiert nichts) für jeden Radius mit
unterschiedlichen Transparenzwerten eine Linie zu plotten.
Nach einigen Einstellungen bezüglich der Farben, viel "gray", damit
das Diagramm im hellen und dunklen Modus gleich gut funktioniert, wird
es mit `components' als skript und div exportiert und auf der Seite
eingebettet.
3 Bibliography
══════════════
“Nominatim,”. 2022. July 12, 2022, URL:
<https://nominatim.openstreetmap.org/ui/search.html>, retrieved on
July 27, 2022.
Van de Ven, B. 2022. “Bokeh,” July 7, 2022, URL: <https://bokeh.org>,
retrieved on July 27, 2022.
Nav
═══
⁃ Tags: [Meta] - [Wetter] - [Python]
⁃ Formats: [md] - [txt] - [html] - [gmi]
[Meta] <./tags/Meta.org>
[Wetter] <./tags/Wetter.org>
[Python] <./tags/Python.org>
[md] <./20220727-regenradar.md>
[txt] <./20220727-regenradar.txt>
[html] <./20220727-regenradar.html>
[gmi] <./20220727-regenradar.gmi>
Footer
══════
License: CC BY-4.0
[Impressum und Datenschutz]
[Impressum und Datenschutz] <./impressum-datenschutz.gmi>