xarray IRIS backend¶

In this example, we read IRIS (sigmet) data files using the wradlib iris xarray backend.

[1]:

import glob
import gzip
import io
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()

/home/runner/micromamba-root/envs/wradlib-notebooks/lib/python3.11/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Load IRIS Volume Data¶

[2]:

fpath = "sigmet/SUR210819000227.RAWKPJV"
f = wrl.util.get_wradlib_data_file(fpath)
vol = wrl.io.open_iris_dataset(f, reindex_angle=False)

Inspect RadarVolume¶

[3]:

display(vol)

<wradlib.RadarVolume>
Dimension(s): (sweep: 1)
Elevation(s): (0.5)

Goereferencing¶

[6]:

swp = vol[0].copy().pipe(wrl.georef.georeference_dataset)

Plotting¶

[7]:

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

../../_images/notebooks_fileio_wradlib_iris_backend_14_0.png

[8]:

fig = pl.figure(figsize=(10, 10))
swp.DBZH.wradlib.plot_ppi(proj="cg", fig=fig)

[8]:

<matplotlib.collections.QuadMesh at 0x7fdc770caf90>

../../_images/notebooks_fileio_wradlib_iris_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.wradlib.plot_ppi(proj=map_proj)
ax = pl.gca()
ax.gridlines(crs=map_proj)
print(ax)

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

../../_images/notebooks_fileio_wradlib_iris_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.wradlib.plot_ppi(ax=ax)
ax.gridlines(draw_labels=True)

[11]:

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

../../_images/notebooks_fileio_wradlib_iris_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).wradlib.plot_ppi(ax=ax)
plot_borders(ax)
ax.gridlines(draw_labels=True)

[12]:

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

../../_images/notebooks_fileio_wradlib_iris_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.wradlib.plot_ppi(proj=map_proj, ax=ax)
ax = pl.gca()
ax.gridlines(crs=map_proj)

[13]:

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

../../_images/notebooks_fileio_wradlib_iris_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.wradlib.plot_ppi(ax=ax)
ax.gridlines()

[14]:

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

../../_images/notebooks_fileio_wradlib_iris_backend_21_1.png

[15]:

swp.DBZH.wradlib.plot_ppi()

[15]:

<matplotlib.collections.QuadMesh at 0x7fdc76e1cd50>

../../_images/notebooks_fileio_wradlib_iris_backend_22_1.png

Create simple plot¶

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

[17]:

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

[17]:

<matplotlib.collections.QuadMesh at 0x7fdc76d10c90>

../../_images/notebooks_fileio_wradlib_iris_backend_26_1.png

[18]:

fig = pl.figure(figsize=(5, 5))
pm = swp.DBZH.wradlib.plot_ppi(proj={"latmin": 3e3}, fig=fig)

../../_images/notebooks_fileio_wradlib_iris_backend_27_0.png

Export to ODIM and CfRadial2¶

Need to remove DB_XHDR since it can’t be represented as ODIM/CfRadial2 moment.

[21]:

vol[0] = vol[0].drop("DB_XHDR", errors="ignore")
vol[0].DBZH.sortby("rtime").plot(y="rtime")

[21]:

<matplotlib.collections.QuadMesh at 0x7fdc769e7290>

../../_images/notebooks_fileio_wradlib_iris_backend_32_1.png

[22]:

vol.to_odim("iris_as_odim.h5")
vol.to_cfradial2("iris_as_cfradial2.nc")

Check equality¶

We have to drop the time variable when checking equality since IRIS has millisecond resolution.

[25]:

xr.testing.assert_allclose(vol.root.drop("time"), vola.root.drop("time"))
xr.testing.assert_allclose(
    vol[0].drop(["rtime", "time", "DB_HCLASS2"]), vola[0].drop(["rtime", "time"])
)
xr.testing.assert_allclose(vol.root.drop("time"), volb.root.drop("time"))
xr.testing.assert_allclose(vol[0].drop("time"), volb[0].drop("time"))
xr.testing.assert_allclose(vola.root, volb.root)
xr.testing.assert_allclose(vola[0].drop("rtime"), volb[0].drop(["rtime", "DB_HCLASS2"]))

xarray IRIS backend¶

Load IRIS Volume Data¶

Inspect RadarVolume¶

Inspect root group¶

Inspect sweep group(s)¶

Goereferencing¶

Plotting¶

Inspect radar moments¶

Create simple plot¶

Mask some values¶

Export to ODIM and CfRadial2¶

Import again¶

Check equality¶

More Iris loading mechanisms¶

Use `xr.open_dataset` to retrieve explicit group¶

xarray IRIS backend¶

Load IRIS Volume Data¶

Inspect RadarVolume¶

Inspect root group¶

Inspect sweep group(s)¶

Goereferencing¶

Plotting¶

Inspect radar moments¶

Create simple plot¶

Mask some values¶

Export to ODIM and CfRadial2¶

Import again¶

Check equality¶

More Iris loading mechanisms¶

Use xr.open_dataset to retrieve explicit group¶

Use `xr.open_dataset` to retrieve explicit group¶