xarray ODIM backend

Contents

xarray ODIM backend#

In this example, we read ODIM_H5 (HDF5) data files using the xradar odim backend.

[1]:

import glob
import os
import wradlib as wrl
import warnings

warnings.filterwarnings("ignore")
import matplotlib.pyplot as pl
import numpy as np
import xradar as xd
import datatree as xt
import xarray as xr

try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    pl.ion()

Load ODIM_H5 Volume Data#

[2]:

fpath = "hdf5/knmi_polar_volume.h5"
f = wrl.util.get_wradlib_data_file(fpath)
vol = xd.io.open_odim_datatree(f)

Downloading file 'hdf5/knmi_polar_volume.h5' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/knmi_polar_volume.h5' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

Goereferencing#

[6]:

swp = vol["sweep_0"].ds
swp = swp.assign_coords(sweep_mode=swp.sweep_mode)
swp = swp.pipe(wrl.georef.georeference_dataset)

Plotting#

[7]:

swp.DBZH.plot.pcolormesh(x="x", y="y")
pl.gca().set_aspect("equal")

../../_images/notebooks_fileio_wradlib_odim_backend_14_0.png

[8]:

fig = pl.figure(figsize=(10, 10))
swp.DBZH.wrl.vis.plot(proj="cg", fig=fig)

[8]:

<matplotlib.collections.QuadMesh at 0x7f0b7008c8d0>

../../_images/notebooks_fileio_wradlib_odim_backend_15_1.png

[9]:

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

map_trans = ccrs.AzimuthalEquidistant(
    central_latitude=swp.latitude.values, central_longitude=swp.longitude.values
)

[10]:

map_proj = ccrs.AzimuthalEquidistant(
    central_latitude=swp.latitude.values, central_longitude=swp.longitude.values
)
pm = swp.DBZH.wrl.vis.plot(proj=map_proj)
ax = pl.gca()
ax.gridlines(crs=map_proj)
print(ax)

< GeoAxes: +proj=aeqd +ellps=WGS84 +lon_0=4.7899699211120605 +lat_0=52.953338623046875 +x_0=0.0 +y_0=0.0 +no_defs +type=crs >

../../_images/notebooks_fileio_wradlib_odim_backend_17_1.png

[11]:

map_proj = ccrs.Mercator(central_longitude=swp.longitude.values)
fig = pl.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=map_proj)
pm = swp.DBZH.wrl.vis.plot(ax=ax)
ax.gridlines(draw_labels=True)

[11]:

<cartopy.mpl.gridliner.Gridliner at 0x7f0b726b79d0>

../../_images/notebooks_fileio_wradlib_odim_backend_18_1.png

[12]:

import cartopy.feature as cfeature


def plot_borders(ax):
    borders = cfeature.NaturalEarthFeature(
        category="physical", name="coastline", scale="10m", facecolor="none"
    )
    ax.add_feature(borders, edgecolor="black", lw=2, zorder=4)


map_proj = ccrs.Mercator(central_longitude=swp.longitude.values)
fig = pl.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=map_proj)

DBZH = swp.DBZH
pm = DBZH.where(DBZH > 0).wrl.vis.plot(ax=ax)
plot_borders(ax)
ax.gridlines(draw_labels=True)

[12]:

<cartopy.mpl.gridliner.Gridliner at 0x7f0b62e27cd0>

../../_images/notebooks_fileio_wradlib_odim_backend_19_1.png

[13]:

import matplotlib.path as mpath

theta = np.linspace(0, 2 * np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)

map_proj = ccrs.AzimuthalEquidistant(
    central_latitude=swp.latitude.values,
    central_longitude=swp.longitude.values,
)
fig = pl.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=map_proj)
ax.set_boundary(circle, transform=ax.transAxes)

pm = swp.DBZH.wrl.vis.plot(proj=map_proj, ax=ax)
ax = pl.gca()
ax.gridlines(crs=map_proj)

[13]:

<cartopy.mpl.gridliner.Gridliner at 0x7f0b727934d0>

../../_images/notebooks_fileio_wradlib_odim_backend_20_1.png

[14]:

fig = pl.figure(figsize=(10, 8))
proj = ccrs.AzimuthalEquidistant(
    central_latitude=swp.latitude.values, central_longitude=swp.longitude.values
)
ax = fig.add_subplot(111, projection=proj)
pm = swp.DBZH.wrl.vis.plot(ax=ax)
ax.gridlines()

[14]:

<cartopy.mpl.gridliner.Gridliner at 0x7f0b72df5510>

../../_images/notebooks_fileio_wradlib_odim_backend_21_1.png

[15]:

swp.DBZH.wrl.vis.plot()

[15]:

<matplotlib.collections.QuadMesh at 0x7f0b702e9210>

../../_images/notebooks_fileio_wradlib_odim_backend_22_1.png

Create simple plot#

Using xarray features a simple plot can be created like this. Note the sortby('time') method, which sorts the radials by time.

[17]:

swp.DBZH.sortby("time").plot(x="range", y="time", add_labels=False)

[17]:

<matplotlib.collections.QuadMesh at 0x7f0b702cbe50>

../../_images/notebooks_fileio_wradlib_odim_backend_26_1.png

[18]:

fig = pl.figure(figsize=(5, 5))
pm = swp.DBZH.wrl.vis.plot(proj={"latmin": 33e3}, fig=fig)

../../_images/notebooks_fileio_wradlib_odim_backend_27_0.png

Mask some values#

[19]:

swp["DBZH"] = swp["DBZH"].where(swp["DBZH"] >= 0)
swp["DBZH"].plot()

[19]:

<matplotlib.collections.QuadMesh at 0x7f0b62d27f10>

../../_images/notebooks_fileio_wradlib_odim_backend_29_1.png