xarray GAMIC backend

Contents

xarray GAMIC backend#

In this example, we read GAMIC (HDF5) data files using the xradar gamic backend.

[1]:

import glob
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 ODIM_H5 Volume Data#

[2]:

fpath = "hdf5/DWD-Vol-2_99999_20180601054047_00.h5"
f = wrl.util.get_wradlib_data_file(fpath)
vol = xd.io.open_gamic_datatree(f)

Downloading file 'hdf5/DWD-Vol-2_99999_20180601054047_00.h5' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/DWD-Vol-2_99999_20180601054047_00.h5' 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_gamic_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 0x7f40dda9f3d0>

../../_images/notebooks_fileio_wradlib_gamic_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.4569489 +lat_0=50.9287272 +x_0=0.0 +y_0=0.0 +no_defs +type=crs >

../../_images/notebooks_fileio_wradlib_gamic_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 0x7f40dd884750>

../../_images/notebooks_fileio_wradlib_gamic_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 0x7f40dd9ed0d0>

../../_images/notebooks_fileio_wradlib_gamic_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 0x7f40dd8d2a50>

../../_images/notebooks_fileio_wradlib_gamic_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 0x7f40ddba9510>

../../_images/notebooks_fileio_wradlib_gamic_backend_21_1.png

[15]:

swp.DBZH.wrl.vis.plot()

[15]:

<matplotlib.collections.QuadMesh at 0x7f40dd584790>

../../_images/notebooks_fileio_wradlib_gamic_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 0x7f40ddb4a8d0>

../../_images/notebooks_fileio_wradlib_gamic_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_gamic_backend_27_0.png

Mask some values#

[19]:

swp["DBZH"] = swp["DBZH"].where(swp["DBZH"] >= 0)
swp["DBZH"].plot()

[19]:

<matplotlib.collections.QuadMesh at 0x7f40dd512990>

../../_images/notebooks_fileio_wradlib_gamic_backend_29_1.png

Export to ODIM and CfRadial2#

[20]:

xd.io.to_odim(vol, "gamic_as_odim.h5")
xd.io.to_cfradial2(vol, "gamic_as_cfradial2.nc")