đŸ Archived View for bacaliu.de âș 20220727-regenradar.gmi captured on 2023-07-10 at 13:27:40. Gemini links have been rewritten to link to archived content
âĄïž Next capture (2023-09-08)
-=-=-=-=-=-=-
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]