wradlib time series data and quasi vertical profiles

Overview

Within this notebook, we will cover:

Reading radar sweep timeseries data into xarray based RadarVolume
Examination of RadarVolume and Sweep
Calculation of Quasivertical Profiles and Plotting

Prerequisites

Concepts	Importance	Notes
Matplotlib Basics	Helpful	Basic Plotting
Xarray Basics	Helpful	Basic Dataset/DataArray
Xarray Plotting	Helpful	Basic Plotting/Faceting

Time to learn: 7.5 minutes

Imports

import glob
import os

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from tqdm import tqdm_notebook as tqdm

import wradlib as wrl

Import Australian Radar Data

It is assumed, that data from IDR71 (Terrey Hills, Sidney) from 20th of December 2018 is used in this notebook.

fglob = "data/hdf5/terrey_*.h5"
idr71 = glob.glob(fglob)
idr71.sort()
print("Files available: {}".format(len(idr71)))

Files available: 40

TimeSeries QVP

All wradlib xarray backends have the capability to read multiple sweeps/volumes in one go. We have to prepare the list of files a bit, though.

We have a time dimension for every ray (the coordinate time depends on azimuth), since xradar tries to adhere as much as possible to the cfradial standard. Then, we need to concatenate into a new dimension other than “time”, otherwise it will fail. Another option would be to reduce the time dimension to only one value, we can use a preprocess function for that.

def preprocess(ds):
    ds.coords["time"] = ds.coords["time"].values.min()
    return ds

ts = xr.open_mfdataset(
    idr71,
    engine="odim",
    group="sweep_0",
    combine="nested",
    concat_dim="vtime",
#    preprocess=preprocess
)

display(ts)

<xarray.Dataset>
Dimensions:            (vtime: 40, azimuth: 360, range: 200)
Coordinates:
  * azimuth            (azimuth) float32 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    time               (vtime, azimuth) datetime64[ns] dask.array<chunksize=(1, 360), meta=np.ndarray>
  * range              (range) float32 125.0 375.0 625.0 ... 4.962e+04 4.988e+04
    longitude          float64 151.2
    latitude           float64 -33.7
    altitude           float64 195.0
Dimensions without coordinates: vtime
Data variables: (12/13)
    DBZH               (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    VRADH              (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    WRADH              (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    TH                 (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    ZDR                (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    RHOHV              (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    ...                 ...
    SNRH               (vtime, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    sweep_mode         (vtime) <U20 'azimuth_surveillance' ... 'azimuth_surve...
    sweep_number       (vtime) int64 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
    prt_mode           (vtime) <U7 'not_set' 'not_set' ... 'not_set' 'not_set'
    follow_mode        (vtime) <U7 'not_set' 'not_set' ... 'not_set' 'not_set'
    sweep_fixed_angle  (vtime) float64 32.0 32.0 32.0 32.0 ... 32.0 32.0 32.0

Georeference and add height coordinate

Here we have another issue if we want to use wradlib’s georeference method. The method reads “sweep_mode” to separate volume scans from RHIs, but since the task information comes as variables in cfradial, it gets concatenated over time and the method will fail

ts.pipe(add_height)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[12], line 1
----> 1 ts.pipe(add_height)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/common.py:775, in DataWithCoords.pipe(self, func, *args, **kwargs)
    773     return func(*args, **kwargs)
    774 else:
--> 775     return func(self, *args, **kwargs)

Cell In[5], line 2, in add_height(ds)
      1 def add_height(ds):
----> 2     ds = ds.pipe(wrl.georef.georeference_dataset)
      3     height = ds.z.mean("azimuth")
      4     ds = ds.assign_coords({"height": (["range"], height.data)})

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/common.py:775, in DataWithCoords.pipe(self, func, *args, **kwargs)
    773     return func(*args, **kwargs)
    774 else:
--> 775     return func(self, *args, **kwargs)

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/wradlib/georef/xarray.py:232, in georeference_dataset(obj, **kwargs)
    229 obj.coords["gr"] = (dimlist, gr)
    231 # adding rays, bins coordinates
--> 232 if obj.sweep_mode == "azimuth_surveillance":
    233     bins, rays = np.meshgrid(obj["range"], obj["azimuth"], indexing="xy")
    234 else:

File /home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/xarray/core/common.py:153, in AbstractArray.__bool__(self)
    152 def __bool__(self: Any) -> bool:
--> 153     return bool(self.values)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

We could reduce it’s dimension to only one value, or we can set is as a coordinate in the preprocessing function to avoid it’s concatenation over time (as long as it does not change)

def preprocess(ds):
    ds.coords["time"] = ds.coords["time"].values.min()
    ds = ds.set_coords(["sweep_mode"])
    return ds

ts = xr.open_mfdataset(
    idr71,
    engine="odim",
    group="sweep_0",
    combine="nested",
    concat_dim="time",
    preprocess=preprocess
)

ts

<xarray.Dataset>
Dimensions:            (time: 40, azimuth: 360, range: 200)
Coordinates:
  * azimuth            (azimuth) float32 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
  * time               (time) datetime64[ns] 2018-12-20T05:04:32.015277824 .....
  * range              (range) float32 125.0 375.0 625.0 ... 4.962e+04 4.988e+04
    sweep_mode         <U20 'azimuth_surveillance'
    longitude          float64 151.2
    latitude           float64 -33.7
    altitude           float64 195.0
Data variables:
    DBZH               (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    VRADH              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    WRADH              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    TH                 (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    ZDR                (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    RHOHV              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    PHIDP              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    SNRH               (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    sweep_number       (time) int64 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
    prt_mode           (time) <U7 'not_set' 'not_set' ... 'not_set' 'not_set'
    follow_mode        (time) <U7 'not_set' 'not_set' ... 'not_set' 'not_set'
    sweep_fixed_angle  (time) float64 32.0 32.0 32.0 32.0 ... 32.0 32.0 32.0

ts = ts.pipe(add_height)
display(ts)

<xarray.Dataset>
Dimensions:            (time: 40, azimuth: 360, range: 200)
Coordinates: (12/15)
  * azimuth            (azimuth) float32 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation          (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
  * time               (time) datetime64[ns] 2018-12-20T05:04:32.015277824 .....
  * range              (range) float32 125.0 375.0 625.0 ... 4.962e+04 4.988e+04
    sweep_mode         <U20 'azimuth_surveillance'
    longitude          float64 151.2
    ...                 ...
    y                  (azimuth, range) float32 106.0 318.0 ... 4.216e+04
    z                  (azimuth, range) float32 262.0 395.0 ... 2.673e+04
    gr                 (azimuth, range) float32 106.0 318.0 ... 4.216e+04
    rays               (azimuth, range) float32 0.5 0.5 0.5 ... 359.5 359.5
    bins               (azimuth, range) float32 125.0 375.0 ... 4.988e+04
    height             (range) float32 262.0 395.0 527.0 ... 2.66e+04 2.673e+04
Data variables:
    DBZH               (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    VRADH              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    WRADH              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    TH                 (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    ZDR                (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    RHOHV              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    PHIDP              (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    SNRH               (time, azimuth, range) float32 dask.array<chunksize=(1, 360, 200), meta=np.ndarray>
    sweep_number       (time) int64 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
    prt_mode           (time) <U7 'not_set' 'not_set' ... 'not_set' 'not_set'
    follow_mode        (time) <U7 'not_set' 'not_set' ... 'not_set' 'not_set'
    sweep_fixed_angle  (time) float64 32.0 32.0 32.0 32.0 ... 32.0 32.0 32.0

Calculate Statistics

stats = ["median", "mean", "min", "max"]
stat = [
    getattr(ts.where(ts.RHOHV > 0.8), st)("azimuth", skipna=True, keep_attrs=True)
    for st in stats
]
ts_stats = xr.concat(stat, dim="stats")
ts_stats = ts_stats.assign_coords({"stats": stats})

display(ts_stats)

<xarray.Dataset>
Dimensions:            (stats: 4, time: 40, range: 200)
Coordinates:
  * time               (time) datetime64[ns] 2018-12-20T05:04:32.015277824 .....
  * range              (range) float32 125.0 375.0 625.0 ... 4.962e+04 4.988e+04
    sweep_mode         <U20 'azimuth_surveillance'
    longitude          float64 151.2
    latitude           float64 -33.7
    altitude           float64 195.0
    height             (range) float32 262.0 395.0 527.0 ... 2.66e+04 2.673e+04
  * stats              (stats) <U6 'median' 'mean' 'min' 'max'
Data variables:
    DBZH               (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    VRADH              (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    WRADH              (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    TH                 (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    ZDR                (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    RHOHV              (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    PHIDP              (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    SNRH               (stats, time, range) float32 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    sweep_number       (stats, time, range) float64 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    sweep_fixed_angle  (stats, time, range) float64 dask.array<chunksize=(1, 1, 200), meta=np.ndarray>
    prt_mode           (stats, time, range) object dask.array<chunksize=(3, 1, 200), meta=np.ndarray>
    follow_mode        (stats, time, range) object dask.array<chunksize=(3, 1, 200), meta=np.ndarray>

Plot QVP’s

levels = np.arange(-30, 80, 5)
facet = ts_stats.TH.plot(
    x="time",
    y="height",
    col="stats",
    col_wrap=2,
    cmap="turbo",
    figsize=(12, 10),
    levels=levels,
)

/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,
/home/jgiles/mambaforge/envs/wradlib4/lib/python3.11/site-packages/dask/array/reductions.py:618: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)