Recipe #2: Reading and visualizing an ODIM_H5 polar volume

Recipe #2: Reading and visualizing an ODIM_H5 polar volume#

This recipe shows how extract the polar volume data from an ODIM_H5 hdf5 file (KNMI example file from OPERA), contruct a 3-dimensional Cartesian volume and produce a diagnostic plot. The challenge for this file is that for each elevation angle, the scan strategy is different.

[2]:

import datetime as dt
from osgeo import osr

# read the data (sample file in WRADLIB_DATA)
filename = wrl.util.get_wradlib_data_file("hdf5/knmi_polar_volume.h5")

raw_dt = xd.io.open_odim_datatree(filename)
display(raw_dt)

<xarray.DatasetView> Size: 388B
Dimensions:              (sweep: 14)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int64 8B 0
    platform_type        <U5 20B 'fixed'
    instrument_type      <U5 20B 'radar'
    time_coverage_start  <U20 80B '2011-06-10T11:40:02Z'
    time_coverage_end    <U20 80B '2011-06-10T11:43:54Z'
    longitude            float32 4B 4.79
    altitude             float32 4B 50.0
    latitude             float32 4B 52.95
    sweep_fixed_angle    (sweep) float32 56B 0.3 0.4 0.8 1.1 ... 15.0 20.0 25.0
    sweep_group_name     (sweep) int64 112B 0 1 2 3 4 5 6 7 8 9 10 11 12 13
Attributes:
    Conventions:      ODIM_H5/V2_2
    instrument_name:  None
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using xradar
[3]:

proj = osr.SpatialReference()
proj.ImportFromEPSG(32632)
for key in list(raw_dt.children):
    if "sweep" in key:
        raw_dt[key].ds = raw_dt[key].ds.wrl.georef.georeference(crs=proj)

[4]:

swp_list = []
for key in list(raw_dt.children):
    if "sweep" in key:
        ds = raw_dt[key].ds
        xyz = (
            xr.concat(
                [
                    ds.coords["x"].reset_coords(drop=True),
                    ds.coords["y"].reset_coords(drop=True),
                    ds.coords["z"].reset_coords(drop=True),
                ],
                "xyz",
            )
            .stack(npoints=("azimuth", "range"))
            .transpose(..., "xyz")
        )
        swp_list.append(xyz)
xyz = xr.concat(swp_list, "npoints")

[5]:

swp_list[0]

[5]:

<xarray.DataArray 'x' (npoints: 115200, xyz: 3)> Size: 3MB
array([[2.17262259e+05, 5.87587733e+06, 5.26327115e+01],
       [2.17329665e+05, 5.87687562e+06, 5.79865115e+01],
       [2.17397078e+05, 5.87787390e+06, 6.34581474e+01],
       ...,
       [2.33423522e+05, 6.19241546e+06, 7.64849169e+03],
       [2.33475527e+05, 6.19341275e+06, 7.69116241e+03],
       [2.33527539e+05, 6.19441005e+06, 7.73395058e+03]])
Coordinates:
  * npoints  (npoints) object 922kB MultiIndex
  * azimuth  (npoints) float32 461kB 0.5 0.5 0.5 0.5 ... 359.5 359.5 359.5 359.5
  * range    (npoints) float32 461kB 500.0 1.5e+03 ... 3.185e+05 3.195e+05
Dimensions without coordinates: xyz
Attributes:
    standard_name:  east_west_distance_from_radar
    units:          meters
[6]:

data_list = []
for key in list(raw_dt.children):
    if "sweep" in key:
        ds = raw_dt[key].ds
        data = ds.DBZH.stack(npoints=("azimuth", "range"))
        data_list.append(data)
data = xr.concat(data_list, "npoints")

[7]:

# generate 3-D Cartesian target grid coordinates
sitecoords = (raw_dt.longitude.values, raw_dt.latitude.values, raw_dt.altitude.values)
maxrange = 200000.0
minelev = 0.1
maxelev = 25.0
maxalt = 5000.0
horiz_res = 2000.0
vert_res = 250.0
trgxyz, trgshape = wrl.vpr.make_3d_grid(
    sitecoords, proj, maxrange, maxalt, horiz_res, vert_res
)

[8]:

# interpolate to Cartesian 3-D volume grid
tstart = dt.datetime.now()
gridder = wrl.vpr.CAPPI(
    xyz.values,
    trgxyz,
    # gridshape=trgshape,
    maxrange=maxrange,
    minelev=minelev,
    maxelev=maxelev,
)
vol = np.ma.masked_invalid(gridder(data.values).reshape(trgshape))
print("3-D interpolation took:", dt.datetime.now() - tstart)

3-D interpolation took: 0:00:00.701437

[9]:

# diagnostic plot
trgx = trgxyz[:, 0].reshape(trgshape)[0, 0, :]
trgy = trgxyz[:, 1].reshape(trgshape)[0, :, 0]
trgz = trgxyz[:, 2].reshape(trgshape)[:, 0, 0]
wrl.vis.plot_max_plan_and_vert(
    trgx,
    trgy,
    trgz,
    vol,
    unit="dBZH",
    levels=range(-32, 60),
    cmap="turbo",
)
../../_images/notebooks_workflow_recipe2_11_0.png

Note

In order to run the recipe code, you need to extract the sample data into a directory pointed to by environment variable WRADLIB_DATA.