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 pl
import numpy as np
import xarray as xr

try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    pl.ion()
from wradlib.io import open_odim_mfdataset

/home/runner/micromamba-root/envs/wradlib-notebooks/lib/python3.11/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

[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.20.3'

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)

  2%|▎         | 6/240 [00:03<02:29,  1.56it/s]
  2%|▎         | 6/240 [00:01<01:03,  3.68it/s]
  2%|▎         | 6/240 [00:01<01:03,  3.67it/s]
  2%|▎         | 6/240 [00:02<01:26,  2.72it/s]
  2%|▎         | 6/240 [00:01<01:15,  3.10it/s]
  2%|▎         | 6/240 [00:02<01:18,  2.97it/s]
  2%|▎         | 6/240 [00:01<00:57,  4.09it/s]
  2%|▎         | 6/240 [00:01<01:00,  3.88it/s]
  2%|▎         | 6/240 [00:01<00:59,  3.95it/s]
  2%|▎         | 6/240 [00:01<01:16,  3.05it/s]

CPU times: user 4.68 s, sys: 169 ms, total: 4.85 s
Wall time: 19.8 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)

  2%|▎         | 6/240 [00:02<01:25,  2.75it/s]
  2%|▎         | 6/240 [00:02<01:20,  2.89it/s]
  2%|▎         | 6/240 [00:02<01:20,  2.91it/s]
  2%|▎         | 6/240 [00:01<01:14,  3.13it/s]
  2%|▎         | 6/240 [00:02<01:28,  2.64it/s]
  2%|▎         | 6/240 [00:01<01:16,  3.05it/s]
  2%|▎         | 6/240 [00:02<01:23,  2.81it/s]
  2%|▎         | 6/240 [00:01<00:58,  4.00it/s]
  2%|▎         | 6/240 [00:01<01:00,  3.88it/s]
  2%|▎         | 6/240 [00:01<01:04,  3.61it/s]

CPU times: user 4.98 s, sys: 88.2 ms, total: 5.07 s
Wall time: 19.3 s

[9]:

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

Export to Cf/Radial2¶

This exports the radar volume at given timestep including all moments into one Cf/Radial2 compliant data file.

[30]:

xradar.io.to_cfradial2(vol0, "dwd_cfradial2.nc")

Import again and check equality¶

[31]:

vol1 = xradar.io.open_odim_datatree("dwd_odim.h5")
vol2 = open_datatree("dwd_cfradial2.nc")

[32]:

xr.testing.assert_equal(vol1.root.ds, vol2.root.ds)
for grp in vol1.groups[1:]:
    xr.testing.assert_equal(
        vol1[grp]
        .to_dataset()
        .reset_coords(["latitude", "longitude", "altitude"], drop=True),
        vol2[grp].to_dataset().swap_dims(time="azimuth").sortby("azimuth"),
    )

Load ODIM_H5 Volume data from German Weather Service¶

Download radar volumes of latest 30 minutes from server using `wetterdienst`¶

Acquire data as memory buffer¶

Read the data into xarray powered structure¶

Inspect structure¶

Root Group¶

Sweep Groups¶

plot sweeps¶

DBZH¶

VRADH¶

Plot single sweep using cartopy¶

Inspect radar moments¶

Plot Quasi Vertical Profile¶

Export to OdimH5¶

Export to Cf/Radial2¶

Import again and check equality¶

Load ODIM_H5 Volume data from German Weather Service¶

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

Acquire data as memory buffer¶

Read the data into xarray powered structure¶

Inspect structure¶

Root Group¶

Sweep Groups¶

plot sweeps¶

DBZH¶

VRADH¶

Plot single sweep using cartopy¶

Inspect radar moments¶

Plot Quasi Vertical Profile¶

Export to OdimH5¶

Export to Cf/Radial2¶

Import again and check equality¶

Download radar volumes of latest 30 minutes from server using `wetterdienst`¶