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

/srv/conda/envs/notebook/lib/python3.9/site-packages/requests/__init__.py:102: RequestsDependencyWarning: urllib3 (1.26.8) or chardet (5.2.0)/charset_normalizer (2.0.10) doesn't match a supported version!
  warnings.warn("urllib3 ({}) or chardet ({})/charset_normalizer ({}) doesn't match a supported "

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.

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

display(ts)

<xarray.Dataset>
Dimensions:     (time: 40, azimuth: 360, range: 200)
Coordinates:
  * azimuth     (azimuth) float32 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5
    elevation   (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    rtime       (time, 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
  * time        (time) datetime64[ns] 2018-12-20T05:04:32 ... 2018-12-20T08:5...
    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>
Attributes:
    fixed_angle:  32.0

Georeference and add height coordinate

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

<xarray.Dataset>
Dimensions:     (time: 40, azimuth: 360, range: 200)
Coordinates: (12/16)
  * azimuth     (azimuth) float32 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5
    elevation   (azimuth) float32 dask.array<chunksize=(360,), meta=np.ndarray>
    rtime       (time, 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
  * time        (time) datetime64[ns] 2018-12-20T05:04:32 ... 2018-12-20T08:5...
    sweep_mode  <U20 'azimuth_surveillance'
    ...          ...
    y           (azimuth, range) float32 106.0 318.0 ... 4.195e+04 4.216e+04
    z           (azimuth, range) float32 262.0 395.0 ... 2.66e+04 2.673e+04
    gr          (azimuth, range) float32 106.0 318.0 ... 4.195e+04 4.216e+04
    rays        (azimuth, range) float32 0.5 0.5 0.5 0.5 ... 359.5 359.5 359.5
    bins        (azimuth, range) float32 125.0 375.0 ... 4.962e+04 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>
Attributes:
    fixed_angle:  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:
  * range       (range) float32 125.0 375.0 625.0 ... 4.962e+04 4.988e+04
  * time        (time) datetime64[ns] 2018-12-20T05:04:32 ... 2018-12-20T08:5...
    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>
Attributes:
    fixed_angle:  32.0

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,
)

/srv/conda/envs/notebook/lib/python3.9/site-packages/dask/array/reductions.py:640: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(x_chunk, axis=axis, keepdims=keepdims)
/srv/conda/envs/notebook/lib/python3.9/site-packages/dask/array/reductions.py:611: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)
/srv/conda/envs/notebook/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out,