Load ODIM_H5 Volume data from German Weather Service#

In this example, we obtain and read the latest 30 minutes of available volumetric radar data from German Weather Service available at opendata.dwd.de. Finally we do some plotting.

This retrieves 6 timesteps of the 10 sweeps (moments DBZH and VRADH) of the DWD volume scan of a distinct radar. This amounts to 120 data files which are combined into one volumetric Cf/Radial2 like xarray powered structure.

Exports to single file Odim_H5 and Cf/Radial2 format are shown at the end of this tutorial.

Note

The following code is based on xarray, xarray-datatree and xradar. It claims multiple data files and presents them in a DataTree.

[1]:
import wradlib as wrl
import warnings

warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()
[2]:
import urllib3
import os
import io
import glob
import shutil
import datetime

Download radar volumes of latest 30 minutes from server using wetterdienst#

wetterdienst is a neat package for easy retrieval of data primarily from DWD. For further information have a look at their documentation.

[3]:
from wetterdienst.provider.dwd.radar import (
    DwdRadarDataFormat,
    DwdRadarDataSubset,
    DwdRadarParameter,
    DwdRadarValues,
)
from wetterdienst.provider.dwd.radar.sites import DwdRadarSite
[4]:
elevations = range(10)

end_date = datetime.datetime.utcnow()
start_date = end_date - datetime.timedelta(minutes=30)

results_velocity = []
results_reflectivity = []

for el in elevations:
    # Horizontal Doppler Velocity
    request_velocity = DwdRadarValues(
        parameter=DwdRadarParameter.SWEEP_VOL_VELOCITY_H,
        start_date=start_date,
        end_date=end_date,
        site=DwdRadarSite.ESS,
        elevation=el,
        fmt=DwdRadarDataFormat.HDF5,
        subset=DwdRadarDataSubset.POLARIMETRIC,
    )

    # Horizontal Reflectivity
    request_reflectivity = DwdRadarValues(
        parameter=DwdRadarParameter.SWEEP_VOL_REFLECTIVITY_H,
        start_date=start_date,
        end_date=end_date,
        elevation=el,
        site=DwdRadarSite.ESS,
        fmt=DwdRadarDataFormat.HDF5,
        subset=DwdRadarDataSubset.POLARIMETRIC,
    )

    # Submit requests.
    results_velocity.append(request_velocity.query())
    results_reflectivity.append(request_reflectivity.query())
[5]:
import wetterdienst

wetterdienst.__version__
[5]:
'0.49.0'

Acquire data as memory buffer#

[6]:
%%time
volume_velocity = []
for item1 in results_velocity:
    files = []
    for item2 in item1:
        files.append(item2.data)
    volume_velocity.append(files)
  3%|▎         | 7/240 [00:03<02:03,  1.88it/s]
  3%|▎         | 7/240 [00:01<00:55,  4.20it/s]
  3%|▎         | 7/240 [00:02<01:10,  3.30it/s]
  3%|▎         | 7/240 [00:02<01:10,  3.29it/s]
  3%|▎         | 7/240 [00:02<01:07,  3.45it/s]
  2%|▎         | 6/240 [00:01<01:07,  3.45it/s]
  2%|▎         | 6/240 [00:01<01:05,  3.55it/s]
  2%|▎         | 6/240 [00:02<01:25,  2.75it/s]
  2%|▎         | 6/240 [00:02<01:31,  2.55it/s]
  2%|▎         | 6/240 [00:01<01:16,  3.08it/s]
CPU times: user 3.67 s, sys: 56 ms, total: 3.72 s
Wall time: 21.6 s

[7]:
volume_velocity = [v[-6:] for v in volume_velocity]
volume_velocity = np.array(volume_velocity).T.tolist()
[8]:
%%time
volume_reflectivity = []
for item1 in results_reflectivity:
    files = []
    for item2 in item1:
        files.append(item2.data)
    volume_reflectivity.append(files)
  3%|▎         | 7/240 [00:02<01:20,  2.89it/s]
  3%|▎         | 7/240 [00:02<01:21,  2.85it/s]
  3%|▎         | 7/240 [00:02<01:19,  2.94it/s]
  3%|▎         | 7/240 [00:02<01:17,  3.02it/s]
  3%|▎         | 7/240 [00:02<01:19,  2.95it/s]
  3%|▎         | 7/240 [00:03<01:50,  2.11it/s]
  3%|▎         | 7/240 [00:01<00:56,  4.12it/s]
  2%|▎         | 6/240 [00:01<00:59,  3.91it/s]
  2%|▎         | 6/240 [00:01<01:03,  3.69it/s]
  2%|▎         | 6/240 [00:03<02:05,  1.87it/s]
CPU times: user 3.65 s, sys: 95 ms, total: 3.74 s
Wall time: 23.4 s

[9]:
volume_reflectivity = [v[-6:] for v in volume_reflectivity]
volume_reflectivity = np.array(volume_reflectivity).T.tolist()

Read the data into xarray powered structure#

[10]:
from datatree import DataTree, open_datatree
import xradar


def concat_radar_datatree(objs, dim="volume_time"):
    root_ds = [obj["/"].ds for obj in objs]
    root = xr.concat(root_ds, dim=dim)
    dtree = DataTree(data=root, name="root")
    for grp in objs[0].groups[1:]:
        ngrps = [obj[grp[1:]].ds for obj in objs]
        ngrp = xr.concat(ngrps, dim=dim)
        DataTree(ngrp, name=grp[1:], parent=dtree)
    return dtree
[11]:
dsl = []
reindex_angle = dict(
    tolerance=1.0, start_angle=0, stop_angle=360, angle_res=1.0, direction=1
)
for r, v in zip(volume_reflectivity, volume_velocity):
    ds0 = [
        xr.open_dataset(r0, engine="odim", group="sweep_0", reindex_angle=reindex_angle)
        for r0 in r
    ]
    ds1 = [
        xr.open_dataset(v0, engine="odim", group="sweep_0", reindex_angle=reindex_angle)
        for v0 in v
    ]
    ds = [xr.merge([r0, v0], compat="override") for r0, v0 in zip(ds0, ds1)]
    ds.insert(0, xr.open_dataset(r[0], group="/"))
    dsl.append(ds)
# this takes some private functions from xradar, take care here
trees = [
    DataTree(data=xradar.io.backends.common._assign_root(ds), name="root") for ds in dsl
]
trees = [
    xradar.io.backends.common._attach_sweep_groups(tree, ds[1:])
    for tree, ds in zip(trees, dsl)
]
vol = concat_radar_datatree(trees, dim="volume_time")
# align sweep_numbers to cover for single sweep single moment layout of DWD
for i, swp in enumerate(vol.groups[1:]):
    vol[swp]["sweep_number"] = i
[12]:
vol
[12]:
<xarray.DatasetView>
Dimensions:              (volume_time: 6)
Dimensions without coordinates: volume_time
Data variables:
    volume_number        (volume_time) int64 0 0 0 0 0 0
    platform_type        (volume_time) <U5 'fixed' 'fixed' ... 'fixed' 'fixed'
    instrument_type      (volume_time) <U5 'radar' 'radar' ... 'radar' 'radar'
    time_coverage_start  (volume_time) <U20 '2023-09-28T16:13:24Z' ... '2023-...
    time_coverage_end    (volume_time) <U20 '2023-09-28T16:18:22Z' ... '2023-...
    longitude            (volume_time) float64 6.967 6.967 6.967 ... 6.967 6.967
    altitude             (volume_time) float64 185.1 185.1 185.1 ... 185.1 185.1
    latitude             (volume_time) float64 51.41 51.41 51.41 ... 51.41 51.41
Attributes:
    Conventions:      ODIM_H5/V2_2
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using xradar
    instrument_name:  None
[13]:
vol["sweep_9"]
[13]:
<xarray.DatasetView>
Dimensions:            (range: 240, azimuth: 360, volume_time: 6)
Coordinates:
  * range              (range) float32 125.0 375.0 625.0 ... 5.962e+04 5.988e+04
  * azimuth            (azimuth) float64 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (volume_time, azimuth) float64 24.98 24.98 ... 24.98
    time               (volume_time, azimuth) datetime64[ns] 2023-09-28T16:13...
    longitude          float64 6.967
    latitude           float64 51.41
    altitude           float64 185.1
Dimensions without coordinates: volume_time
Data variables:
    DBZH               (volume_time, azimuth, range) float32 -64.0 ... -64.0
    sweep_mode         (volume_time) <U20 'azimuth_surveillance' ... 'azimuth...
    sweep_number       int64 9
    prt_mode           (volume_time) <U7 'not_set' 'not_set' ... 'not_set'
    follow_mode        (volume_time) <U7 'not_set' 'not_set' ... 'not_set'
    sweep_fixed_angle  (volume_time) float64 25.0 25.0 25.0 25.0 25.0 25.0
    VRADH              (volume_time, azimuth, range) float64 -128.0 ... -128.0

Inspect structure#

Root Group#

[14]:
vol.root
[14]:
<xarray.DatasetView>
Dimensions:              (volume_time: 6)
Dimensions without coordinates: volume_time
Data variables:
    volume_number        (volume_time) int64 0 0 0 0 0 0
    platform_type        (volume_time) <U5 'fixed' 'fixed' ... 'fixed' 'fixed'
    instrument_type      (volume_time) <U5 'radar' 'radar' ... 'radar' 'radar'
    time_coverage_start  (volume_time) <U20 '2023-09-28T16:13:24Z' ... '2023-...
    time_coverage_end    (volume_time) <U20 '2023-09-28T16:18:22Z' ... '2023-...
    longitude            (volume_time) float64 6.967 6.967 6.967 ... 6.967 6.967
    altitude             (volume_time) float64 185.1 185.1 185.1 ... 185.1 185.1
    latitude             (volume_time) float64 51.41 51.41 51.41 ... 51.41 51.41
Attributes:
    Conventions:      ODIM_H5/V2_2
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using xradar
    instrument_name:  None

Sweep Groups#

[15]:
vol["sweep_0"]
[15]:
<xarray.DatasetView>
Dimensions:            (range: 720, azimuth: 360, volume_time: 6)
Coordinates:
  * range              (range) float32 125.0 375.0 625.0 ... 1.796e+05 1.799e+05
  * azimuth            (azimuth) float64 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (volume_time, azimuth) float64 5.482 5.482 ... 5.493
    time               (volume_time, azimuth) datetime64[ns] 2023-09-28T16:15...
    longitude          float64 6.967
    latitude           float64 51.41
    altitude           float64 185.1
Dimensions without coordinates: volume_time
Data variables:
    DBZH               (volume_time, azimuth, range) float32 -64.0 ... -64.0
    sweep_mode         (volume_time) <U20 'azimuth_surveillance' ... 'azimuth...
    sweep_number       int64 0
    prt_mode           (volume_time) <U7 'not_set' 'not_set' ... 'not_set'
    follow_mode        (volume_time) <U7 'not_set' 'not_set' ... 'not_set'
    sweep_fixed_angle  (volume_time) float64 5.499 5.499 5.499 5.499 5.499 5.499
    VRADH              (volume_time, azimuth, range) float64 -128.0 ... -128.0

plot sweeps#

DBZH#

[16]:
vol["sweep_0"].isel(volume_time=0)
[16]:
<xarray.DatasetView>
Dimensions:            (range: 720, azimuth: 360)
Coordinates:
  * range              (range) float32 125.0 375.0 625.0 ... 1.796e+05 1.799e+05
  * azimuth            (azimuth) float64 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (azimuth) float64 5.482 5.482 5.482 ... 5.482 5.482 5.482
    time               (azimuth) datetime64[ns] 2023-09-28T16:15:54.520000 .....
    longitude          float64 6.967
    latitude           float64 51.41
    altitude           float64 185.1
Data variables:
    DBZH               (azimuth, range) float32 -64.0 -64.0 ... -64.0 -64.0
    sweep_mode         <U20 'azimuth_surveillance'
    sweep_number       int64 0
    prt_mode           <U7 'not_set'
    follow_mode        <U7 'not_set'
    sweep_fixed_angle  float64 5.499
    VRADH              (azimuth, range) float64 -128.0 -128.0 ... -128.0 -128.0
[17]:
fig, gs = plt.subplots(
    4, 3, figsize=(20, 30), sharex=True, sharey=True, constrained_layout=True
)

for i, grp in enumerate(vol.groups[1:]):
    swp = vol[grp].isel(volume_time=0).ds
    swp = swp.assign_coords(sweep_mode=swp.sweep_mode)
    swp.DBZH.wrl.georef.georeference().wrl.vis.plot(ax=gs.flat[i], fig=fig)
    ax = plt.gca()
    ax.set_title(swp.sweep_fixed_angle.values)

fig.delaxes(gs.flat[-2])
fig.delaxes(gs.flat[-1])
../../_images/notebooks_workflow_recipe4_24_0.png

VRADH#

[18]:
fig, gs = plt.subplots(
    4, 3, figsize=(20, 30), sharex=True, sharey=True, constrained_layout=True
)

for i, grp in enumerate(vol.groups[1:]):
    swp = vol[grp].isel(volume_time=0).ds
    swp = swp.assign_coords(sweep_mode=swp.sweep_mode)
    swp.VRADH.wrl.georef.georeference().wrl.vis.plot(ax=gs.flat[i], fig=fig)
    ax = plt.gca()
    ax.set_title(swp.sweep_fixed_angle.values)

fig.delaxes(gs.flat[-2])
fig.delaxes(gs.flat[-1])
../../_images/notebooks_workflow_recipe4_26_0.png

Plot single sweep using cartopy#

[19]:
vol0 = vol.isel(volume_time=0)
swp = vol0["sweep_9"].ds
# need to assign sweep_mode as coordinate
swp = swp.assign_coords(sweep_mode=swp.sweep_mode)
[20]:
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

map_trans = ccrs.AzimuthalEquidistant(
    central_latitude=vol0.root.ds.latitude.values,
    central_longitude=vol0.root.ds.longitude.values,
)
[21]:
map_proj = ccrs.AzimuthalEquidistant(
    central_latitude=vol0.root.ds.latitude.values,
    central_longitude=vol0.root.ds.longitude.values,
)

pm = swp.DBZH.wrl.georef.georeference().wrl.vis.plot(crs=map_proj)
ax = plt.gca()
ax.gridlines(crs=map_proj)
print(ax)
< GeoAxes: +proj=aeqd +ellps=WGS84 +lon_0=6.967111 +lat_0=51.405649 +x_0=0.0 +y_0=0.0 +no_defs +type=crs >
../../_images/notebooks_workflow_recipe4_30_1.png
[22]:
map_proj = ccrs.Mercator(central_longitude=vol0.root.ds.longitude.values)
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=map_proj)
pm = swp.DBZH.wrl.georef.georeference().wrl.vis.plot(ax=ax)
ax.gridlines(draw_labels=True)
[22]:
<cartopy.mpl.gridliner.Gridliner at 0x7f0ca8875c50>
../../_images/notebooks_workflow_recipe4_31_1.png
[23]:
fig = plt.figure(figsize=(10, 8))
proj = ccrs.AzimuthalEquidistant(
    central_latitude=vol0.root.ds.latitude.values,
    central_longitude=vol0.root.ds.longitude.values,
)
ax = fig.add_subplot(111, projection=proj)
pm = swp.DBZH.wrl.georef.georeference().wrl.vis.plot(ax=ax)
ax.gridlines()
[23]:
<cartopy.mpl.gridliner.Gridliner at 0x7f0ca25f9910>
../../_images/notebooks_workflow_recipe4_32_1.png

Inspect radar moments#

The DataArrays can be accessed by key or by attribute. Each DataArray inherits dimensions and coordinates of it’s parent dataset. There are attributes connected which are defined by Cf/Radial and/or ODIM_H5 standard.

[24]:
vol["sweep_9"].isel(volume_time=0).ds.DBZH
[24]:
<xarray.DataArray 'DBZH' (azimuth: 360, range: 240)>
array([[-64.00293 , -64.00293 , -64.00293 , ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [-64.00293 , -64.00293 ,  -6.809883, ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [-64.00293 , -64.00293 , -14.275066, ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       ...,
       [-64.00293 , -64.00293 , -64.00293 , ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [-64.00293 , -64.00293 , -64.00293 , ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [-64.00293 , -64.00293 , -64.00293 , ..., -64.00293 , -64.00293 ,
        -64.00293 ]], dtype=float32)
Coordinates:
  * range      (range) float32 125.0 375.0 625.0 ... 5.962e+04 5.988e+04
  * azimuth    (azimuth) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
    elevation  (azimuth) float64 24.98 24.98 24.98 24.98 ... 24.98 24.98 24.98
    time       (azimuth) datetime64[ns] 2023-09-28T16:13:54.625999872 ... 202...
    longitude  float64 6.967
    latitude   float64 51.41
    altitude   float64 185.1
Attributes:
    _Undetect:      0.0
    long_name:      Equivalent reflectivity factor H
    standard_name:  radar_equivalent_reflectivity_factor_h
    units:          dBZ
[25]:
vol["sweep_9"].isel(volume_time=0).ds.sweep_mode
[25]:
<xarray.DataArray 'sweep_mode' ()>
array('azimuth_surveillance', dtype='<U20')
Coordinates:
    longitude  float64 6.967
    latitude   float64 51.41
    altitude   float64 185.1

Plot Quasi Vertical Profile#

[26]:
ts = vol["sweep_9"]
ts
[26]:
<xarray.DatasetView>
Dimensions:            (range: 240, azimuth: 360, volume_time: 6)
Coordinates:
  * range              (range) float32 125.0 375.0 625.0 ... 5.962e+04 5.988e+04
  * azimuth            (azimuth) float64 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (volume_time, azimuth) float64 24.98 24.98 ... 24.98
    time               (volume_time, azimuth) datetime64[ns] 2023-09-28T16:13...
    longitude          float64 6.967
    latitude           float64 51.41
    altitude           float64 185.1
Dimensions without coordinates: volume_time
Data variables:
    DBZH               (volume_time, azimuth, range) float32 -64.0 ... -64.0
    sweep_mode         (volume_time) <U20 'azimuth_surveillance' ... 'azimuth...
    sweep_number       int64 9
    prt_mode           (volume_time) <U7 'not_set' 'not_set' ... 'not_set'
    follow_mode        (volume_time) <U7 'not_set' 'not_set' ... 'not_set'
    sweep_fixed_angle  (volume_time) float64 25.0 25.0 25.0 25.0 25.0 25.0
    VRADH              (volume_time, azimuth, range) float64 -128.0 ... -128.0