xarray Rainbow5 backend

Contents

xarray Rainbow5 backend#

In this example, we read Rainbow5 data files using the xradar rainbow 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 xradar as xd
import datatree as xt
import xarray as xr

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

Load Rainbow5 Volume Data#

[2]:

fpath = "rainbow/2013051000000600dBZ.vol"
f = wrl.util.get_wradlib_data_file(fpath)
vol = xd.io.open_rainbow_datatree(f, reindex_angle=False)

Downloading file 'rainbow/2013051000000600dBZ.vol' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/rainbow/2013051000000600dBZ.vol' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

Goereferencing#

[6]:

swp = vol["sweep_0"].ds.copy()
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_rainbow_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 0x7f589ed297d0>

../../_images/notebooks_fileio_wradlib_rainbow_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=6.379967 +lat_0=50.856633 +x_0=0.0 +y_0=0.0 +no_defs +type=crs >

../../_images/notebooks_fileio_wradlib_rainbow_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 0x7f589ec6cd90>

../../_images/notebooks_fileio_wradlib_rainbow_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 0x7f589ebefe10>

../../_images/notebooks_fileio_wradlib_rainbow_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 0x7f589ebedf50>

../../_images/notebooks_fileio_wradlib_rainbow_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 0x7f589ea14690>

../../_images/notebooks_fileio_wradlib_rainbow_backend_21_1.png

[15]:

swp.DBZH.wrl.vis.plot()

[15]:

<matplotlib.collections.QuadMesh at 0x7f589ec150d0>

../../_images/notebooks_fileio_wradlib_rainbow_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 0x7f58a6f02c50>

../../_images/notebooks_fileio_wradlib_rainbow_backend_26_1.png

[18]:

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

../../_images/notebooks_fileio_wradlib_rainbow_backend_27_0.png

Mask some values#

[19]:

dbzh = swp["DBZH"].where(swp["DBZH"] >= 0)
dbzh.plot(x="x", y="y")

[19]:

<matplotlib.collections.QuadMesh at 0x7f589e9a1550>

../../_images/notebooks_fileio_wradlib_rainbow_backend_29_1.png

Export to ODIM and CfRadial2#

[20]:

xd.io.to_odim(vol, "rainbow_as_odim.h5")
xd.io.to_cfradial2(vol, "rainbow_as_cfradial2.nc")