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

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

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

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

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

../../_images/notebooks_fileio_wradlib_iris_backend_21_1.png

[15]:

swp.DBZH.wradlib.plot_ppi()

[15]:

<matplotlib.collections.QuadMesh at 0x7f6ab3333d00>

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

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

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