💾 Archived View for bacaliu.de › 20220727-regenradar.gmi captured on 2023-09-08 at 16:12:06. Gemini links have been rewritten to link to archived content
⬅️ Previous capture (2023-07-10)
-=-=-=-=-=-=-
Der [DWD] stellt auf 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:
Zwecks dessen schreibe ich ein Skript, welches dieses darstellt. Folgendermaßen stelle ich *schematisch* dar, wie das funktioniert.
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.
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].
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
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.
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
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.
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.
“Nominatim,”. 2022. July 12, 2022, URL: , retrieved on July 27, 2022.
Van de Ven, B. 2022. “Bokeh,” July 7, 2022, URL: , retrieved on July 27, 2022.
License: CC BY-4.0 [Impressum und Datenschutz]