Hydrometeor partitioning ratio retrievals for Ground Radar#

In this notebook, measurements from NEXRAD’s KDDC ground radar are used to derive Hydrometeor Partitioning Ratios (HPR) following Pejcic et al. 2025 (in review). This requires the horizontal reflectivity, differential reflectivity, specific differential phase, cross correlation coefficient, temperature information and rain type. The temperature information is derived from sounding and a rain type classification is applied following Park et al. 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 datetime as dt
import xradar as xd
import scipy

Plot Temperature Profile#

[7]:

fig = plt.figure(figsize=(5, 10))
ax = fig.add_subplot(111)
rs_ds.TEMP.plot(y="HGHT", ax=ax, zorder=0, c="r")
ax.grid(True)

get freezing level height#

We need to obtain the freezing level height, which is needed for an ad-hoc retrieval of raintype.

[8]:

fl = np.abs(rs_ds).argmin("HGHT").TEMP
display(fl)

<xarray.DataArray 'TEMP' ()> Size: 8B
array(4914)

Interpolate Temperature onto sweeps#

The following function interpolates the vertical temperature profile onto the radar sweeps.

[11]:

def merge_radar_profile(rds, cds):
    if "z" in rds.coords:
        cds = cds.interp({"HGHT": rds.z}, method="linear")
        rds = rds.assign({"TEMP": cds})
    return rds


gr_data3 = gr_data2.map_over_datasets(merge_radar_profile, rs_ds.TEMP)

[12]:

gr_data3["sweep_1"].TEMP.plot(x="x", y="y")

[12]:

<matplotlib.collections.QuadMesh at 0x7f25dbc0f4d0>

../../_images/notebooks_classify_hmcp_nexrad_20_1.png

Ad-hoc retrieval of raintype#

The following algorithm of raintype estimation is derived after (Park et al.).

keep all radar bins >= 45 dBZ
keep all radar bins > 30 dBZ and height > fl + 1600m
combine 1 and 2
iterate over x,y pairs and fetch from whole tree to set as convective.

[13]:

def mask_data(rds, fl):
    if "z" in rds.coords:
        # Thresholding and smoothing (Park et al.)
        # -----------------------------------------
        xwin_zh = 5
        rds = rds.where(rds.RH > 0.8)
        rds["CZ"] = rds.CZ.rolling(
            range=xwin_zh, min_periods=xwin_zh // 2, center=True
        ).mean(skipna=True)
        mask = (rds.CZ >= 45) | ((rds.CZ > 30) & (rds.z > (fl + 1600)))
        rds = rds.assign(mask=mask)
    return rds


gr_data4 = gr_data3.map_over_datasets(mask_data, fl)

Extract xyz bin coordinates#

This iterates over the whole DataTree and extracts the RainType-mask as 1-dimensional array. This keeps only valid values.

[14]:

def get_xyz(tree):
    swp_list = []
    for key in list(tree.children):
        if "sweep" in key:
            ds = tree[key].ds.stack(npoints=("azimuth", "range"))
            ds = ds.reset_coords().where(ds.mask, drop=True)
            swp_list.append(ds.mask)
    return xr.concat(swp_list, "npoints")

Calculate hydrometeor partitioning ratios (HPR)#

This uses the loaded weights and centroids to retrieve the hydrometeor partitioning ratio from the observations.

[21]:

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

Plotting all Hydrometeor-Classes#

For better plotting we transfrom to 100% and drop NaN data.

[22]:

hmpr = hmpr.dropna("range", how="all") * 100
hpr_bins = [0, 1, 2.5, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]

hmpr.plot(
    col="hmc",
    col_wrap=3,
    x="x",
    y="y",
    cmap="HomeyerRainbow",
    levels=hpr_bins,
    cbar_kwargs={"ticks": hpr_bins},
)

[22]:

<xarray.plot.facetgrid.FacetGrid at 0x7f25c0232120>

../../_images/notebooks_classify_hmcp_nexrad_39_1.png

Hydrometeor partitioning ratio retrievals for Ground Radar

Contents

Hydrometeor partitioning ratio retrievals for Ground Radar#

Read centroids, covariances and weights#

Read polarimetric radar observations#

Get Temperature Profile#

Plot Temperature Profile#

get freezing level height#

georeference DataTree#

Interpolate Temperature onto sweeps#

Ad-hoc retrieval of raintype#

Extract xyz bin coordinates#

Interpolation of RainType mask#

ZDR Offset retrieval#

Extract sweep 2 for further processing#

Combine observations into xr.DataArray#

Calculate hydrometeor partitioning ratios (HPR)#

Plotting all Hydrometeor-Classes#