đŸ Archived View for bacaliu.de âș 20220727-regenradar.md captured on 2023-07-10 at 13:51:41.
-=-=-=-=-=-=-
# Inhaltsverzeichnis
1. [Der Plan](#der-plan)
2. [Das Vorgehen](#das-vorgehen)
1. [Laden der Daten](#laden-der-daten)
2. [Speichern als .pickle](#speichern-als-pickle)
3. [Laden der Arrays und Timestamps](#laden-der-arrays)
4. [Stereographische Umwandlung](#stereographische-umwandlung)
5. [FehlerprĂŒfung: richtige Koordinaten?](#fehlerpruefung-richtige-koordinaten)
6. [Endlich: das Plotten](#endlich-das-plotten)
3. [Bibliography](#bibliography)
- [Nav](#nav)
- [Footer](#footer)
<a id="der-plan"></a>
# Der Plan
Der [DWD](https://www.dwd.de) 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 (<a href="#citeproc_bib_item_2">Van de Ven</a>). 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.
<a id="das-vorgehen"></a>
# Das Vorgehen
Zwecks dessen schreibe ich ein Skript, welches dieses darstellt. FolgendermaĂen stelle ich **schematisch** dar, wie das funktioniert.
<a id="laden-der-daten"></a>
## 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}")
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.
<a id="speichern-als-pickle"></a>
## 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)
`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](https://www.dwd.de/DE/leistungen/radolan/radolan_info/radolan_radvor_op_komposit_format_pdf.pdf?__blob=publicationFile&v=5).
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.
<a id="laden-der-arrays"></a>
## 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)
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
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.
<a id="stereographische-umwandlung"></a>
## Stereographische Umwandlung
Der Benutzer wird seine Koordinaten in lat/lon angeben oder durch eingabe eines Suchtermes ĂŒber (<a href="#citeproc_bib_item_1">âNominatimâ 2022</a>) 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
<a id="fehlerpruefung-richtige-koordinaten"></a>
## 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])
Die Arrays werden geladen und deren Shape gibt Auskunft darĂŒber ob die gewĂ€hlten Koordinaten im Bereich liegen.
<a id="endlich-das-plotten"></a>
## 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
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.
<a id="bibliography"></a>
# Bibliography
<style>.csl-entry{text-indent: -0; margin-left: 0;}</style><div class="csl-bib-body">
<div class="csl-entry"><a id="citeproc_bib_item_1"></a>âNominatim,â. 2022. July 12, 2022, URL: <a href="https://nominatim.openstreetmap.org/ui/search.html">https://nominatim.openstreetmap.org/ui/search.html</a>, retrieved on July 27, 2022.</div>
<div class="csl-entry"><a id="citeproc_bib_item_2"></a>Van de Ven, B. 2022. âBokeh,â July 7, 2022, URL: <a href="https://bokeh.org">https://bokeh.org</a>, retrieved on July 27, 2022.</div>
</div>
<a id="nav"></a>
# Nav
- Tags: [Meta](./tags/Meta.md) - [Wetter](./tags/Wetter.md) - [Python](./tags/Python.md)
<!-- BEGIN insert Backlinks (but there are no) -->
- Formats: [md](./20220727-regenradar.md) - [txt](./20220727-regenradar.txt) - [html](./20220727-regenradar.html) - [gmi](./20220727-regenradar.gmi)
<a id="footer"></a>
# Footer
License: CC BY-4.0
[Impressum und Datenschutz](./impressum-datenschutz.gmi)