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

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 0x7f66a0f5ae30>

../../_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 0x7f66a0c28a90>

../../_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 0x7f66a0c962f0>

../../_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 0x7f66a091c310>

../../_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 0x7f66a0de0eb0>

../../_images/notebooks_fileio_wradlib_iris_backend_21_1.png

[15]:

swp.DBZH.wradlib.plot_ppi()

[15]:

<matplotlib.collections.QuadMesh at 0x7f66a0e72230>

../../_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 0x7f66a100a560>

../../_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 0x7f669cd8c370>

../../_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¶