xarray CfRadial2 backend

Contents

xarray CfRadial2 backend#

In this example, we read CfRadial2 data files using the xarray cfradial2 backend.

[1]:

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 CfRadial2 Volume Data#

[2]:

fpath = "netcdf/cfrad.20080604_002217_000_SPOL_v36_SUR_cfradial2.nc"
f = wrl.util.get_wradlib_data_file(fpath)
vol = xt.open_datatree(f)

Downloading file 'netcdf/cfrad.20080604_002217_000_SPOL_v36_SUR_cfradial2.nc' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/netcdf/cfrad.20080604_002217_000_SPOL_v36_SUR_cfradial2.nc' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

[3]:

# fix: remove when available in xradar
for k in vol.groups[1:]:
    vol[k].ds = (
        vol[k]
        .ds.assign(sweep_fixed_angle=vol[k].ds.attrs["fixed_angle"])
        .swap_dims(time="azimuth")
        .sortby("azimuth")
    )

Goereferencing#

[7]:

swp = vol["sweep_0"].ds.copy()
swp = swp.assign_coords(sweep_mode=swp.sweep_mode)
swp = swp.pipe(wrl.georef.georeference_dataset)

Plotting#

[8]:

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_15_0.png

[9]:

fig = pl.figure(figsize=(10, 10))
swp.DBZ.wrl.vis.plot(proj="cg", fig=fig)

[9]:

<matplotlib.collections.QuadMesh at 0x7f04009bba90>

../../_images/notebooks_fileio_wradlib_cfradial2_backend_16_1.png

[10]:

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
)

[11]:

map_proj = ccrs.AzimuthalEquidistant(
    central_latitude=swp.latitude.values, central_longitude=swp.longitude.values
)
pm = swp.DBZ.wrl.vis.plot(proj=map_proj)
ax = pl.gca()
ax.gridlines(crs=map_proj)
print(ax)

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_18_1.png

[12]:

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.DBZ.wrl.vis.plot(ax=ax)
ax.gridlines(draw_labels=True)

[12]:

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_19_1.png

[13]:

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)

DBZ = swp.DBZ
pm = DBZ.where(DBZ > 0).wrl.vis.plot(ax=ax)
plot_borders(ax)
ax.gridlines(draw_labels=True)

[13]:

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_20_1.png

[14]:

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.DBZ.wrl.vis.plot(proj=map_proj, ax=ax)
ax = pl.gca()
ax.gridlines(crs=map_proj)

[14]:

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_21_1.png

[15]:

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.DBZ.wrl.vis.plot(ax=ax)
ax.gridlines()

[15]:

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_22_1.png

[16]:

swp.DBZ.wrl.vis.plot()

[16]:

<matplotlib.collections.QuadMesh at 0x7f04006fde90>

../../_images/notebooks_fileio_wradlib_cfradial2_backend_23_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.

[18]:

swp.DBZ.sortby("time").plot(x="range", y="time", add_labels=False)

[18]:

<matplotlib.collections.QuadMesh at 0x7f04005b1e90>

../../_images/notebooks_fileio_wradlib_cfradial2_backend_27_1.png

[19]:

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

../../_images/notebooks_fileio_wradlib_cfradial2_backend_28_0.png