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_09-2020031209490200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_08-2020031209484900-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_07-2020031209483600-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_06-2020031209482200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_05-2020031209480000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_04-2020031209473000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_03-2020031209470700-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_02-2020031209464300-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_01-2020031209462000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_dbzh_00-2020031209455700-ess-10410-hd5
[6]:
for f in filelist2[:10]:
print(f)
ras07-vol5minng01_sweeph5onem_vradh_09-2020031209490200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_08-2020031209484900-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_07-2020031209483600-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_06-2020031209482200-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_05-2020031209480000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_04-2020031209473000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_03-2020031209470700-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_02-2020031209464300-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_01-2020031209462000-ess-10410-hd5
ras07-vol5minng01_sweeph5onem_vradh_00-2020031209455700-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, standard='cf', 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'
time_coverage_start <U20 '2020-03-12T09:45:35Z'
time_coverage_end <U20 '2020-03-12T09:49:02Z'
latitude float64 51.41
longitude float64 6.967
altitude float64 185.0
altitude_agl float64 nan
sweep_group_name (sweep) <U8 'sweep_1' 'sweep_2' ... 'sweep_10'
sweep_fixed_angle (sweep) float64 2.499 3.499 5.499 ... 17.0 12.0 25.0
frequency float64 nan
status_xml <U4 'None'
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([ 2.49938965, 3.49914551, 5.49865723, 0.49987793, 1.49963379,
7.99804688, 4.49890137, 17.00134277, 12.00256348, 24.99938965])
Dimensions without coordinates: sweepSweep 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: (range: 180, time: 360)
Coordinates:
sweep_mode <U20 'azimuth_surveillance'
latitude float64 51.41
altitude float64 185.0
longitude float64 6.967
elevation (time) float64 2.505 2.505 2.505 2.505 ... 2.505 2.505 2.505
azimuth (time) float64 0.5164 1.533 2.532 3.521 ... 357.5 358.5 359.5
* range (range) float32 500.0 1500.0 2500.0 ... 178500.0 179500.0
y (time, range) float64 499.5 1.498e+03 ... 1.781e+05 1.791e+05
z (time, range) float64 206.8 250.5 ... 9.858e+03 9.923e+03
gr (time, range) float64 499.5 1.499e+03 ... 1.781e+05 1.791e+05
rays (time, range) float64 0.5164 0.5164 0.5164 ... 359.5 359.5
bins (time, range) float32 500.0 1500.0 ... 178500.0 179500.0
x (time, range) float64 4.502 13.5 ... -1.529e+03 -1.537e+03
* time (time) datetime64[ns] 2020-03-12T09:47:01.938499840 ... 2020-03-12T09:47:01.875000064
Data variables:
DBZH (time, range) float32 ...
sweep_number int64 0
follow_mode <U4 'none'
prt_mode <U5 'fixed'
fixed_angle float64 2.499
VRADH (time, 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])
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])
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 0x7feabb136810> >
[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 0x7feaba8e3190>
[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 0x7feabdd96610>
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' (time: 360, range: 60)>
array([[-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 , -0.074223, -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:
sweep_mode <U20 'azimuth_surveillance'
latitude float64 51.41
altitude float64 185.0
longitude float64 6.967
elevation (time) float64 25.0 25.0 25.0 25.0 25.0 ... 25.0 25.0 25.0 25.0
azimuth (time) float64 0.5109 1.516 2.516 3.516 ... 357.5 358.5 359.5
* range (range) float32 500.0 1500.0 2500.0 ... 57500.0 58500.0 59500.0
y (time, range) float64 453.1 1.359e+03 ... 5.286e+04 5.376e+04
z (time, range) float64 396.3 819.0 ... 2.508e+04 2.551e+04
gr (time, range) float64 453.1 1.359e+03 ... 5.286e+04 5.376e+04
rays (time, range) float64 0.5109 0.5109 0.5109 ... 359.5 359.5 359.5
bins (time, range) float32 500.0 1500.0 2500.0 ... 58500.0 59500.0
x (time, range) float64 4.04 12.12 20.2 ... -445.8 -453.6 -461.3
* time (time) datetime64[ns] 2020-03-12T09:48:54.537499904 ... 2020-03-12T09:48:54.503500032
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'
latitude float64 51.41
altitude float64 185.0
longitude float64 6.967[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'
time_coverage_start <U20 '2020-03-12T09:45:35Z'
time_coverage_end <U20 '2020-03-12T09:49:02Z'
latitude float64 51.41
longitude float64 6.967
altitude float64 185.0
altitude_agl float64 nan
sweep_group_name (sweep) <U8 'sweep_1' 'sweep_2' ... 'sweep_10'
sweep_fixed_angle (sweep) float64 2.499 3.499 5.499 ... 17.0 12.0 25.0
frequency float64 nan
status_xml <U4 'None'
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('dwd_odim.h5')
Export to Cf/Radial2¶
This exports the radar volume including all moments into one Cf/Radial2 compliant data file.
[25]:
vol.to_cfradial2('dwd_cfradial2.nc')
Import again and check equality¶
Small time differences are possible so drop times before comparison
[26]:
try:
vol1 = OdimH5('dwd_odim.h5', standard='cf', georef=True)
vol2 = CfRadial('dwd_cfradial2.nc', dim0='azimuth', 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