Load ODIM_H5 Volume data from German Weather Service

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

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

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

[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().magic("matplotlib inline")
except:
    pl.ion()
from wradlib.io.xarray import CfRadial, OdimH5
[2]:
import urllib
import os
import io
import glob

Download Current Filelist from server

Import it into pandas DataFrame, export to list and sort it descending in time.

Helper Class to parse response

[3]:
from html.parser import HTMLParser

class DWDHTMLParser(HTMLParser):
    def handle_starttag(self, tag, attrs):
        if tag != 'a':
            return
        self.links.append(attrs[0][1])

parser = DWDHTMLParser()
[4]:
radar = 'ESS'
DBZH = 'sweep_vol_z'
VRADH = 'sweep_vol_v'

opendata_url1 = (f"https://opendata.dwd.de/weather/radar/sites/{DBZH}/{radar.lower()}/hdf5/")
with urllib.request.urlopen(opendata_url1) as url_request:
    response = url_request.read().decode("utf-8")

parser.links = []
parser.feed(response)
filelist1 = parser.links[1:]

filelist1.sort(key=lambda x: x.split('-')[2])
filelist1.reverse()

opendata_url2 = (f"https://opendata.dwd.de/weather/radar/sites/{VRADH}/{radar.lower()}/hdf5/")
with urllib.request.urlopen(opendata_url2) as url_request:
    response = url_request.read().decode("utf-8")

parser.links = []
parser.feed(response)
filelist2 = parser.links[1:]

filelist2.sort(key=lambda x: x.split('-')[2])
filelist2.reverse()
[5]:
for f in filelist1[:10]:
    print(f)
ras07-vol5minng01_sweeph5onem_dbzh_01-2019091611112000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_00-2019091611105700-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_09-2019091611090200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_08-2019091611084900-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_07-2019091611083500-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_06-2019091611082200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_05-2019091611080000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_04-2019091611073000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_03-2019091611070600-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_02-2019091611064300-ess-10410-hd5
[6]:
for f in filelist2[:10]:
    print(f)
ras07-vol5minng01_sweeph5onem_vradh_01-2019091611112000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_00-2019091611105700-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_09-2019091611090200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_08-2019091611084900-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_07-2019091611083500-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_06-2019091611082200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_05-2019091611080000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_04-2019091611073000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_03-2019091611070600-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_02-2019091611064300-ess-10410-hd5

Clean up local folder

[7]:
flist = glob.glob('ras07*')
for f in flist:
    os.remove(f)

Download latest volume to current directory

[8]:
for f in filelist1[:10]:
    urllib.request.urlretrieve(os.path.join(opendata_url1, f), f)

for f in filelist2[:10]:
    urllib.request.urlretrieve(os.path.join(opendata_url2, f), f)

Read the files into xarray powered structure

[9]:
flist = glob.glob('ras07*')
vol = wrl.io.OdimH5(flist, dim0='azimuth', georef=True)

Inspect structure

Root Group

[10]:
vol.root
[10]:
<xarray.Dataset>
Dimensions:              (sweep: 10)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int64 0
    platform_type        <U5 'fixed'
    instrument_type      <U5 'radar'
    primary_axis         <U6 'axis_z'
    altitude_agl         float64 nan
    frequency            float64 nan
    status_xml           <U4 'None'
    latitude             (sweep) float64 51.41 51.41 51.41 ... 51.41 51.41 51.41
    longitude            (sweep) float64 6.967 6.967 6.967 ... 6.967 6.967 6.967
    altitude             (sweep) float64 185.1 185.1 185.1 ... 185.0 185.1 185.1
    time_coverage_start  <U20 '2019-09-16T11:06:21Z'
    time_coverage_end    <U20 '2019-09-16T11:11:20Z'
    sweep_group_name     (sweep) <U8 'sweep_1' 'sweep_2' ... 'sweep_10'
    sweep_fixed_angle    (sweep) float64 0.4999 2.499 3.499 ... 5.499 1.5 25.0
Attributes:
    Conventions:          Cf/Radial
    version:              H5rad 2.2
    title:                None
    institution:          WMO:10410,NOD:deess
    references:           None
    source:               None
    history:              None
    comment:              im/exported using wradlib
    instrument_name:      None
    site_name:            name of site where data were gathered
    scan_name:            name of scan strategy used, if applicable
    scan_id:              scan strategy id, if applicable. assumed 0 if missing
    platform_is_mobile:   "true" or "false", assumed "false" if missing
    ray_times_increase:   "true" or "false", assumed "true" if missing. This ...
    field_names:          array of strings of field names present in this file.
    time_coverage_start:  copy of time_coverage_start global variable
    time_coverage_end:    copy of time_coverage_end global variable
    simulated data:       "true" or "false", assumed "false" if missing. data...
    instrument:           WMO:10410,NOD:deess
[11]:
vol.root.sweep_fixed_angle
[11]:
<xarray.DataArray 'sweep_fixed_angle' (sweep: 10)>
array([ 0.499878,  2.49939 ,  3.499146,  7.998047, 12.002563, 17.001343,
        4.498901,  5.498657,  1.499634, 24.99939 ])
Dimensions without coordinates: sweep

Sweep Groups

[12]:
list(vol)
[12]:
['sweep_1',
 'sweep_2',
 'sweep_3',
 'sweep_4',
 'sweep_5',
 'sweep_6',
 'sweep_7',
 'sweep_8',
 'sweep_9',
 'sweep_10']
[13]:
vol['sweep_1']
[13]:
<xarray.Dataset>
Dimensions:       (azimuth: 360, range: 180)
Coordinates:
  * azimuth       (azimuth) float64 0.5109 1.516 2.521 ... 357.5 358.5 359.5
    elevation     (azimuth) float64 0.5054 0.5054 0.5054 ... 0.5054 0.5054
  * range         (range) float32 500.0 1500.0 2500.0 ... 178500.0 179500.0
    sweep_mode    <U20 'azimuth_surveillance'
    longitude     float64 6.967
    latitude      float64 51.41
    altitude      float64 185.1
    x             (azimuth, range) float64 4.458 13.37 ... -1.531e+03 -1.54e+03
    y             (azimuth, range) float64 499.9 1.5e+03 ... 1.784e+05 1.794e+05
    z             (azimuth, range) float64 189.4 198.2 ... 3.636e+03 3.666e+03
    gr            (azimuth, range) float64 500.0 1.5e+03 ... 1.784e+05 1.794e+05
    rays          (azimuth, range) float64 0.5109 0.5109 0.5109 ... 359.5 359.5
    bins          (azimuth, range) float32 500.0 1500.0 ... 178500.0 179500.0
    time          (azimuth) datetime64[ns] 2019-09-16T11:07:52.192000 ... 2019-09-16T11:07:52.107999744
Data variables:
    VRADH         (azimuth, range) float32 ...
    fixed_angle   float64 0.4999
    follow_mode   <U4 'none'
    prt_mode      <U5 'fixed'
    sweep_number  int64 0
    DBZH          (azimuth, range) float32 ...

plot sweeps

DBZH

[14]:
fig = pl.figure(figsize=(20, 30))
import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(4, 3, wspace=0.4, hspace=0.4)
for i, swp in enumerate(vol.values()):
    swp.DBZH.wradlib.plot(ax=gs[i], fig=fig)
    ax = pl.gca()
    ax.set_title(vol.root.sweep_fixed_angle[i])
../../_images/notebooks_fileio_wradlib_load_DWD_opendata_volumes_23_0.png

VRADH

[15]:
fig = pl.figure(figsize=(20, 30))
import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(4, 3, wspace=0.4, hspace=0.4)
for i, swp in enumerate(vol.values()):
    swp.VRADH.wradlib.plot(ax=gs[i], fig=fig)
    ax = pl.gca()
    ax.set_title(vol.root.sweep_fixed_angle[i])
../../_images/notebooks_fileio_wradlib_load_DWD_opendata_volumes_25_0.png

Plot single sweep using cartopy

[16]:
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

map_trans = ccrs.AzimuthalEquidistant(central_latitude=vol['sweep_10'].latitude.values,
                                      central_longitude=vol['sweep_10'].longitude.values)
[17]:
map_proj = ccrs.AzimuthalEquidistant(central_latitude=vol['sweep_10'].latitude.values,
                                      central_longitude=vol['sweep_10'].longitude.values)
pm = vol['sweep_10'].DBZH.wradlib.plot_ppi(proj=map_proj)
ax = pl.gca()
ax.gridlines(crs=map_proj)
print(ax)
< GeoAxes: <cartopy.crs.AzimuthalEquidistant object at 0x7f184cac2200> >
../../_images/notebooks_fileio_wradlib_load_DWD_opendata_volumes_28_1.png
[18]:
map_proj = ccrs.Mercator(central_longitude=vol['sweep_10'].longitude.values)
fig = pl.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection=map_proj)
pm = vol['sweep_10'].DBZH.wradlib.plot_ppi(ax=ax)
ax.gridlines(draw_labels=True)
[18]:
<cartopy.mpl.gridliner.Gridliner at 0x7f184b9d6668>
../../_images/notebooks_fileio_wradlib_load_DWD_opendata_volumes_29_1.png
[19]:
fig = pl.figure(figsize=(10, 8))
proj=ccrs.AzimuthalEquidistant(central_latitude=vol['sweep_10'].latitude.values,
                               central_longitude=vol['sweep_10'].longitude.values)
ax = fig.add_subplot(111, projection=proj)
pm = vol['sweep_10'].DBZH.wradlib.plot_ppi(ax=ax)
ax.gridlines()
[19]:
<cartopy.mpl.gridliner.Gridliner at 0x7f184bb0ad68>
../../_images/notebooks_fileio_wradlib_load_DWD_opendata_volumes_30_1.png

Inspect radar moments

The dataarrays can be accessed by key or by attribute. Each dataarray has the datasets dimensions and coordinates of it’s parent dataset. There are attributes connected which are defined by Cf/Radial and/or ODIM_H5 standard.

[20]:
vol['sweep_10'].DBZH
[20]:
<xarray.DataArray 'DBZH' (azimuth: 360, range: 60)>
array([[ 21.74871 ,  18.915619,  17.87555 , ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [ 20.945953,  20.951813,  18.605057, ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [ 21.236   ,  19.914673,  18.731041, ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       ...,
       [ 24.67263 ,  21.253578,  18.941986, ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [ 24.584732,  20.140266,  18.31208 , ..., -64.00293 , -64.00293 ,
        -64.00293 ],
       [ 24.795677,  20.773094,  17.740776, ..., -64.00293 , -64.00293 ,
        -64.00293 ]], dtype=float32)
Coordinates:
  * azimuth     (azimuth) float64 0.5109 1.522 2.521 3.521 ... 357.5 358.5 359.5
    elevation   (azimuth) float64 25.02 25.02 25.02 25.02 ... 25.02 25.02 25.02
  * range       (range) float32 500.0 1500.0 2500.0 ... 57500.0 58500.0 59500.0
    sweep_mode  <U20 'azimuth_surveillance'
    longitude   float64 6.967
    latitude    float64 51.41
    altitude    float64 185.1
    x           (azimuth, range) float64 4.04 12.12 20.2 ... -458.6 -466.4
    y           (azimuth, range) float64 453.1 1.359e+03 ... 5.285e+04 5.376e+04
    z           (azimuth, range) float64 396.4 819.2 ... 2.509e+04 2.552e+04
    gr          (azimuth, range) float64 453.1 1.359e+03 ... 5.286e+04 5.376e+04
    rays        (azimuth, range) float64 0.5109 0.5109 0.5109 ... 359.5 359.5
    bins        (azimuth, range) float32 500.0 1500.0 2500.0 ... 58500.0 59500.0
    time        (azimuth) datetime64[ns] 2019-09-16T11:08:54.356999936 ... 2019-09-16T11:08:54.322500096
Attributes:
    standard_name:  radar_equivalent_reflectivity_factor_h
    long_name:      Equivalent reflectivity factor H
    units:          dBZ
[21]:
vol['sweep_10'].sweep_mode
[21]:
<xarray.DataArray 'sweep_mode' ()>
array('azimuth_surveillance', dtype='<U20')
Coordinates:
    sweep_mode  <U20 'azimuth_surveillance'
    longitude   float64 6.967
    latitude    float64 51.41
    altitude    float64 185.1
[22]:
vol.root
[22]:
<xarray.Dataset>
Dimensions:              (sweep: 10)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int64 0
    platform_type        <U5 'fixed'
    instrument_type      <U5 'radar'
    primary_axis         <U6 'axis_z'
    altitude_agl         float64 nan
    frequency            float64 nan
    status_xml           <U4 'None'
    latitude             (sweep) float64 51.41 51.41 51.41 ... 51.41 51.41 51.41
    longitude            (sweep) float64 6.967 6.967 6.967 ... 6.967 6.967 6.967
    altitude             (sweep) float64 185.1 185.1 185.1 ... 185.0 185.1 185.1
    time_coverage_start  <U20 '2019-09-16T11:06:21Z'
    time_coverage_end    <U20 '2019-09-16T11:11:20Z'
    sweep_group_name     (sweep) <U8 'sweep_1' 'sweep_2' ... 'sweep_10'
    sweep_fixed_angle    (sweep) float64 0.4999 2.499 3.499 ... 5.499 1.5 25.0
Attributes:
    Conventions:          Cf/Radial
    version:              H5rad 2.2
    title:                None
    institution:          WMO:10410,NOD:deess
    references:           None
    source:               None
    history:              None
    comment:              im/exported using wradlib
    instrument_name:      None
    site_name:            name of site where data were gathered
    scan_name:            name of scan strategy used, if applicable
    scan_id:              scan strategy id, if applicable. assumed 0 if missing
    platform_is_mobile:   "true" or "false", assumed "false" if missing
    ray_times_increase:   "true" or "false", assumed "true" if missing. This ...
    field_names:          array of strings of field names present in this file.
    time_coverage_start:  copy of time_coverage_start global variable
    time_coverage_end:    copy of time_coverage_end global variable
    simulated data:       "true" or "false", assumed "false" if missing. data...
    instrument:           WMO:10410,NOD:deess
[23]:
#vol.root = vol.root.assign(longitude=vol.root.longitude[0])
#vol.root = vol.root.assign(latitude=vol.root.latitude[0])
#vol.root = vol.root.assign(altitude=vol.root.altitude[0])

Export to OdimH5

This exports the radar volume including all moments into one ODIM_H5 compliant data file.

[24]:
vol.to_odim('testodim.h5')

Export to Cf/Radial2

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

[25]:
vol.to_cfradial2('testcfradial2.nc')

Import again and check equality

Small time differences are possible so drop times before comparison

[26]:
try:
    vol1 = OdimH5('testodim.h5', georef=True)
    vol2 = CfRadial('testcfradial2.nc', georef=True)
    xr.testing.assert_equal(vol1.root, vol2.root)
    xr.testing.assert_equal(vol1['sweep_1'].drop('time'), vol2['sweep_1'].drop('time'))
except:
    pass