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 2025 (in review). 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.

[1]:

import wradlib as wrl
import wradlib_data
import numpy as np

import warnings
warnings.filterwarnings("ignore")

import matplotlib.pyplot as plt
try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()


import glob
import xarray as xr
import scipy
import pandas as dp
from dask.diagnostics import ProgressBar

Read dual-frequency satellite observations (GPM)#

[2]:

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/runner/work/wradlib/wradlib/wradlib-data'.

Plot GPM overpass#

[3]:

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

[3]:

<matplotlib.collections.QuadMesh at 0x7f590c655cd0>

../../_images/notebooks_classify_hmcp_gpm_6_1.png

Assign coordinates#

[4]:

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#

[5]:

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/notebooks_classify_hmcp_gpm_10_0.png

[6]:

# centroids and covariances
cdp_file = wradlib_data.DATASETS.fetch("misc/hmcp_centroids_df.nc")
with xr.open_dataset(cdp_file) as cdp:
    cdp
cdp

Downloading file 'misc/hmcp_centroids_df.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/misc/hmcp_centroids_df.nc' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

[7]:

# weights
weights_file = wradlib_data.DATASETS.fetch("misc/hmcp_weights.nc")
with xr.open_dataset(weights_file) as cw:
    display(cw)
    pass

Downloading file 'misc/hmcp_weights.nc' from 'https://github.com/wradlib/wradlib-data/raw/main/data/misc/hmcp_weights.nc' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

[8]:

with ProgressBar():
    obs = wrl.classify.create_gpm_observations(sr_data)
obs

[8]:

<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:
    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>
    date       (nscan) datetime64[ns] 2kB dask.array<chunksize=(284,), meta=np.ndarray>
    height     (nscan, nray, nbin) float32 10MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
  * nbin       (nbin) int64 1kB 0 1 2 3 4 5 6 7 ... 169 170 171 172 173 174 175
  * nscan      (nscan) int64 2kB 0 1 2 3 4 5 6 7 ... 277 278 279 280 281 282 283
  * nray       (nray) int64 392B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48
  * obs        (obs) <U4 64B 'ZKUM' 'DFRM' 'RT' 'TEMP'
Attributes:
    DimensionNames:    nscan,nray,nbin,nfreq
    Units:             dBZ
    units:             dBZ
    CodeMissingValue:  -9999.9

[9]:

%%time
with ProgressBar():
    hmpr = wrl.classify.calculate_hmpr(obs, cw.weights, cdp)  # .compute()
hmpr

CPU times: user 206 ms, sys: 0 ns, total: 206 ms
Wall time: 206 ms

[9]:

<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'
    temp       (nscan, nray, nbin) float64 20MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
    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>
    date       (nscan) datetime64[ns] 2kB dask.array<chunksize=(284,), meta=np.ndarray>
    height     (nscan, nray, nbin) float32 10MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
  * nbin       (nbin) int64 1kB 0 1 2 3 4 5 6 7 ... 169 170 171 172 173 174 175
  * nscan      (nscan) int64 2kB 0 1 2 3 4 5 6 7 ... 277 278 279 280 281 282 283
  * nray       (nray) int64 392B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48
Attributes:
    standard_name:  hydrometeor_partitioning_ratio
    units:          %

[10]:

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

[10]:

<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'
    temp       (nscan, nray, nbin) float64 20MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
    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>
    date       (nscan) datetime64[ns] 2kB dask.array<chunksize=(284,), meta=np.ndarray>
    height     (nscan, nray, nbin) float32 10MB dask.array<chunksize=(284, 1, 176), meta=np.ndarray>
  * nbin       (nbin) int64 1kB 0 1 2 3 4 5 6 7 ... 169 170 171 172 173 174 175
  * nscan      (nscan) int64 2kB 0 1 2 3 4 5 6 7 ... 277 278 279 280 281 282 283
  * nray       (nray) int64 392B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48
Attributes:
    standard_name:  hydrometeor_partitioning_ratio
    units:          %

[11]:

with ProgressBar():
    hmpr_sel = hmpr.sel(nray=19) * 100
    hmpr_sel = hmpr_sel.compute()
hmpr_sel

[########################################] | 100% Completed | 14.08 s

[12]:

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/notebooks_classify_hmcp_gpm_17_0.png

Hydrometeor partitioning ratio retrievals for GPM

Contents

Hydrometeor partitioning ratio retrievals for GPM#

Read dual-frequency satellite observations (GPM)#

Plot GPM overpass#

Assign coordinates#

Plot overview along track#