Hydrometeor partitioning ratio retrievals for GPM#

In this notebook, GPM Dual Frequency Radar (DPR) measurements are used to derive Hydrometeor Partitioning Ratios (HPR) according to [Pejcic et al., 2026]. This requires the measured Ku-band reflectivity, the dual-frequency ratios (Ku-band - Ka-band) and the DPR temperature and rain type information. The HPRs for the different hydrometeor classes are then presented.

import warnings

import matplotlib.pyplot as plt
import numpy as np
import wradlib as wrl
import wradlib_data
import xarray as xr
from dask.diagnostics import ProgressBar
from IPython.display import display

warnings.filterwarnings("ignore")

/home/docs/checkouts/readthedocs.org/user_builds/wradlib/conda/latest/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Read dual-frequency satellite observations (GPM)#

path_gpm = wradlib_data.DATASETS.fetch(
    "gpm/2A-CS-VP-24.GPM.DPR.V9-20211125.20180625-S050710-E051028.024557.V07A.HDF5"
)
# Read GPM data
sr_data = wrl.io.open_gpm_dataset(path_gpm, group="FS").chunk(nray=1)
sr_data = sr_data.set_coords(["Longitude", "Latitude"])
sr_data = xr.decode_cf(sr_data)

Downloading file 'gpm/2A-CS-VP-24.GPM.DPR.V9-20211125.20180625-S050710-E051028.024557.V07A.HDF5' from 'https://github.com/wradlib/wradlib-data/raw/main/data/gpm/2A-CS-VP-24.GPM.DPR.V9-20211125.20180625-S050710-E051028.024557.V07A.HDF5' to '/home/docs/.cache/wradlib-data'.

Plot GPM overpass#

plt.figure(figsize=(5, 4))
sr_data.zFactorFinalNearSurface.isel(nfreq=0).plot(
    x="Longitude",
    y="Latitude",
    vmin=0,
    vmax=40,
    cmap="turbo",
)

<matplotlib.collections.QuadMesh at 0x75f4316b4c20>

../../../_images/0f35af052ad2a563e54c622cce459a938edcc6e4f627861bd264301869444503.png

Assign coordinates#

sr_data = sr_data.set_coords("height")
sr_data = sr_data.assign_coords(nbin=sr_data.nbin.data)
sr_data = sr_data.assign_coords(nscan=sr_data.nscan.data)
sr_data = sr_data.assign_coords(nray=sr_data.nray.data)

Plot overview along track#

zlvl = np.arange(10, 57.5, 2.5)
zlvl2 = np.arange(10, 57.5, 5)
dpr_lvl = np.array([-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30])

ff = 10
lw = 2.5
x1, x2 = -101, -98.5
y1, y2 = 0, 15000

fig, axs = plt.subplots(1, 3, figsize=(20, 5))  # , sharey='row', sharex='col'

# Ku_measured
KU = sr_data.zFactorMeasured.sel(nfreq=0, nray=19)
plot = KU.plot(
    ax=axs[0],
    x="Longitude",
    y="height",
    cmap="HomeyerRainbow",
    levels=zlvl,
    cbar_kwargs={"extend": "neither", "label": "", "pad": 0.01, "ticks": zlvl2},
    xlim=(x1, x2),
    ylim=(y1, y2),
)

colorbar = plot.colorbar
colorbar.ax.tick_params(labelsize=ff)

# Ka_measured
KA = sr_data.zFactorMeasured.sel(nfreq=1, nray=19)
plot = KA.plot(
    ax=axs[1],
    x="Longitude",
    y="height",
    cmap="HomeyerRainbow",
    levels=zlvl,
    cbar_kwargs={"extend": "neither", "label": "", "pad": 0.01, "ticks": zlvl2},
    xlim=(x1, x2),
    ylim=(y1, y2),
)

colorbar = plot.colorbar
colorbar.ax.tick_params(labelsize=ff)


# DFR_measured
DFR = sr_data.zFactorMeasured.sel(nfreq=0, nray=19) - sr_data.zFactorMeasured.sel(
    nfreq=1, nray=19
)

plot = DFR.plot(
    ax=axs[2],
    x="Longitude",
    y="height",
    cmap="HomeyerRainbow",
    levels=dpr_lvl,
    cbar_kwargs={"extend": "neither", "label": "", "pad": 0.01, "ticks": dpr_lvl},
    xlim=(x1, x2),
    ylim=(y1, y2),
)

colorbar = plot.colorbar
colorbar.ax.tick_params(labelsize=ff)

T = [r"$Z_m^{K_u}$ in dBZ", r"$Z_m^{K_a}$ in dBZ", r"$DFR_m^{K_u-K_a}$ in dB"]
for i in range(len(T)):
    axs[i].set_title("", fontsize=ff)
    axs[i].set_title(T[i], fontsize=ff, loc="right")
    axs[i].set_ylabel("Height in m", fontsize=ff)
    axs[i].set_xlabel("Longitude in deg", fontsize=ff)
    axs[i].grid(ls=":", zorder=-100)
    axs[i].tick_params(axis="both", labelsize=ff)

../../../_images/88cc848a3d78297a3ecd61becb98c31ebd73f1df747f5a90b1db70904abc9d3a.png

Fetch Observations#

obs = sr_data.pipe(wrl.classify.create_gpm_observations)
display(obs)

<xarray.DataArray 'data' (obs: 4, nscan: 284, nray: 49, nbin: 176)> Size: 78MB
dask.array<where, shape=(4, 284, 49, 176), dtype=float64, chunksize=(1, 284, 1, 176), chunktype=numpy.ndarray>
Coordinates:
  * obs        (obs) <U4 64B 'ZKUM' 'DFRM' 'RT' 'TEMP'
  * nscan      (nscan) int64 2kB 0 1 2 3 4 5 6 7 ... 277 278 279 280 281 282 283
    date       (nscan) datetime64[us] 2kB dask.array<chunksize=(284,), meta=np.ndarray>
  * nray       (nray) int64 392B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48
    Latitude   (nscan, nray) float32 56kB dask.array<chunksize=(284, 1), meta=np.ndarray>
    Longitude  (nscan, nray) float32 56kB dask.array<chunksize=(284, 1), meta=np.ndarray>
  * nbin       (nbin) int64 1kB 0 1 2 3 4 5 6 7 ... 169 170 171 172 173 174 175
    height     (nscan, nray, nbin) float32 10MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
Attributes:
    DimensionNames:    nscan,nray,nbin,nfreq
    Units:             dBZ
    units:             dBZ
    CodeMissingValue:  -9999.9

Apply classifier#

hmpr = obs.wrl.classify.calculate_hmpr(cw.weights, cdp)
display(hmpr)

<xarray.DataArray 'HPR' (hmc: 9, nscan: 284, nray: 49, nbin: 176)> Size: 176MB
dask.array<truediv, shape=(9, 284, 49, 176), dtype=float64, chunksize=(9, 284, 49, 176), chunktype=numpy.ndarray>
Coordinates:
  * hmc        (hmc) <U2 72B 'LR' 'MR' 'HR' 'BD' 'RH' 'GR' 'IC' 'WS' 'SN'
  * nscan      (nscan) int64 2kB 0 1 2 3 4 5 6 7 ... 277 278 279 280 281 282 283
    date       (nscan) datetime64[us] 2kB dask.array<chunksize=(284,), meta=np.ndarray>
  * nray       (nray) int64 392B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48
    Latitude   (nscan, nray) float32 56kB dask.array<chunksize=(284, 1), meta=np.ndarray>
    Longitude  (nscan, nray) float32 56kB dask.array<chunksize=(284, 1), meta=np.ndarray>
  * nbin       (nbin) int64 1kB 0 1 2 3 4 5 6 7 ... 169 170 171 172 173 174 175
    temp       (nscan, nray, nbin) float64 20MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
    height     (nscan, nray, nbin) float32 10MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
Attributes:
    standard_name:  hydrometeor_partitioning_ratio
    units:          %

hmpr = hmpr.chunk(hmc=1, nray=1)
display(hmpr)

<xarray.DataArray 'HPR' (hmc: 9, nscan: 284, nray: 49, nbin: 176)> Size: 176MB
dask.array<rechunk-merge, shape=(9, 284, 49, 176), dtype=float64, chunksize=(1, 284, 1, 176), chunktype=numpy.ndarray>
Coordinates:
  * hmc        (hmc) <U2 72B 'LR' 'MR' 'HR' 'BD' 'RH' 'GR' 'IC' 'WS' 'SN'
  * nscan      (nscan) int64 2kB 0 1 2 3 4 5 6 7 ... 277 278 279 280 281 282 283
    date       (nscan) datetime64[us] 2kB dask.array<chunksize=(284,), meta=np.ndarray>
  * nray       (nray) int64 392B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48
    Latitude   (nscan, nray) float32 56kB dask.array<chunksize=(284, 1), meta=np.ndarray>
    Longitude  (nscan, nray) float32 56kB dask.array<chunksize=(284, 1), meta=np.ndarray>
  * nbin       (nbin) int64 1kB 0 1 2 3 4 5 6 7 ... 169 170 171 172 173 174 175
    temp       (nscan, nray, nbin) float64 20MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
    height     (nscan, nray, nbin) float32 10MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
Attributes:
    standard_name:  hydrometeor_partitioning_ratio
    units:          %

Plot results#

hpr_bins = [0, 1, 2.5, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]
x1, x2 = -101, -98.5
y1, y2 = 0, 15000
with ProgressBar():
    hmpr_sel.plot(
        col="hmc",
        col_wrap=3,
        x="Longitude",
        y="height",
        cmap="HomeyerRainbow",
        levels=hpr_bins,
        xlim=(x1, x2),
        ylim=(y1, y2),
        cbar_kwargs={"ticks": hpr_bins},
    )

../../../_images/fc0128197bf83b233b1bb95e59ee1909c6840861c928fdda59d285a84cac366e.png

Hydrometeor partitioning ratio retrievals for GPM

Contents