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

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

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

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

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

../../_images/notebooks_fileio_wradlib_iris_backend_21_1.png

[15]:

swp.DBZH.wradlib.plot_ppi()

[15]:

<matplotlib.collections.QuadMesh at 0x7f3830ccff10>

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

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

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