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), construct 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.

import datetime as dt
import warnings

import numpy as np
import wradlib as wrl
import wradlib_data
import xarray as xr
import xradar as xd
from IPython.display import display
from osgeo import osr

warnings.filterwarnings("ignore")

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

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")

swp_list[0]

<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]],
      shape=(115200, 3))
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:
    axis:           X
    long_name:      Easting
    standard_name:  projection_x_coordinate
    units:          metre

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")

# 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
)

# 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.601973

# 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/5644f25c0995c8f841f68f5d55fdebaf9b4212df16351fa920eb066751ee6ff0.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.