ADD: ZDR bias code from RadTraQ

examples/correct/plot_zdr_check.py

@@ -0,0 +1,50 @@
+"""
+ZDR Bias Calculation
+---------------------
+
+This example shows how to calculate the ZDR bias from VPT/Birdbath scans.
+The technique here uses a vertically pointing scan in regions of light rain.
+In these regions, raindrops should be approximately spherical and therefore their
+ZDR near zero. Therefore, we want the average ZDR in these regions.
+This code applies reflectivity and cross correlation ratio-based thresholds to the ZDR
+bias calculation to ensure that we are taking the average ZDR in light rain.
+
+"""
+
+import matplotlib.pyplot as plt
+from open_radar_data import DATASETS
+
+import pyart
+
+# Read in example data
+filename = DATASETS.fetch("sgpxsaprcfrvptI4.a1.20200205.100827.nc")
+ds = pyart.io.read(filename)
+
+# Set up a typical filter for ZDR bias calculation in birdbath scan
+# Light rain and RhoHV near 1 ensures that raindrops are close to spherical
+# Therefore ZDR should be zero in these regions
+gatefilter = pyart.filters.GateFilter(ds)
+gatefilter.exclude_below("cross_correlation_ratio_hv", 0.995)
+gatefilter.exclude_above("cross_correlation_ratio_hv", 1)
+gatefilter.exclude_below("reflectivity", 10)
+gatefilter.exclude_above("reflectivity", 30)
+
+results = pyart.correct.calc_zdr_offset(
+    ds,
+    zdr_var="differential_reflectivity",
+    gatefilter=gatefilter,
+    height_range=(1000, 3000),
+)
+
+print("Zdr Bias: " + "{:.2f}".format(results["bias"]))
+
+fig, ax = plt.subplots(1, 3, figsize=(8, 5))
+ax[0].plot(results["profile_zdr"], results["range"])
+ax[0].set_ylabel("Range (m)")
+ax[0].set_xlabel("Zdr (dB)")
+ax[1].plot(results["profile_reflectivity"], results["range"])
+ax[1].set_xlabel("Zh (dBZ)")
+ax[2].plot(results["profile_cross_correlation_ratio_hv"], results["range"])
+ax[2].set_xlabel("RhoHV ()")
+fig.tight_layout()
+plt.show()

pyart/correct/__init__.py

@@ -8,6 +8,7 @@
 from .attenuation import calculate_attenuation  # noqa
 from .attenuation import calculate_attenuation_philinear  # noqa
 from .attenuation import calculate_attenuation_zphi  # noqa
+from .bias_and_noise import calc_zdr_offset  # noqa
 from .bias_and_noise import calc_cloud_mask, calc_noise_floor, correct_bias  # noqa
 from .bias_and_noise import (
     correct_noise_rhohv,  # noqa

pyart/correct/bias_and_noise.py

@@ -2,6 +2,7 @@
 Corrects polarimetric variables for noise
 
 """
+
 import warnings
 
 import dask
@@ -144,6 +145,62 @@ def correct_bias(radar, bias=0.0, field_name=None):
     return corr_field
 
 
+def calc_zdr_offset(radar, gatefilter=None, height_range=None, zdr_var=None):
+    """
+    Function for calculating the ZDR bias from a VPT scan.
+
+    Parameters
+    ----------
+    radar : PyART radar object
+        Radar object with radar data
+    gatefilter: PyART GateFilter
+        Gatefilter for filtering out data for calculating ZDR bias
+    height_range: tuple
+        The minimum and maximum heights in meters for the scan.
+    zdr_var: string or None
+        The name of the ZDR variable. Set to None to have PyART try to determine this automatically.
+
+    Returns
+    -------
+    profiles : dict
+        The mean vertical profiles of each radar moment are extracted along with the ZDR bias.
+
+    """
+    if height_range is None:
+        height_range = (0, 100000.0)
+
+    if zdr_var is None:
+        zdr_var = get_field_name("differential_reflectivity")
+
+    height_mask = np.logical_and(
+        radar.range["data"] >= height_range[0], radar.range["data"] <= height_range[1]
+    )
+
+    mask = gatefilter.gate_included
+    Zdr = radar.fields[zdr_var]["data"]
+    if isinstance(Zdr, np.ma.MaskedArray):
+        Zdr = Zdr.filled(np.nan)
+    Zdr = np.where(mask, Zdr, np.nan)
+    bias = np.nanmean(Zdr[:, height_mask])
+    height_mask_tiled = np.tile(height_mask, (radar.nrays, 1))
+    Zdr = np.where(height_mask_tiled, Zdr, np.nan)
+    results = {
+        "bias": bias,
+        "profile_zdr": np.nanmean(Zdr[:, :], axis=0),
+        "range": radar.range["data"],
+    }
+    for k in radar.fields.keys():
+        if k != "range":
+            field_data = radar.fields[k]["data"].astype(float)
+            if isinstance(field_data, np.ma.MaskedArray):
+                field_data = field_data.filled(np.nan)
+
+            field_data = np.where(mask, field_data, np.nan)
+            field_data = np.where(height_mask_tiled, field_data, np.nan)
+            results["profile_" + k] = np.nanmean(field_data[:, :], axis=0)
+    return results
+
+
 def calc_cloud_mask(
     radar,
     field,

tests/correct/test_correct_bias.py

@@ -31,6 +31,25 @@ def test_correct_bias():
     pytest.raises(KeyError, pyart.correct.correct_bias, radar, bias, field_name)
 
 
+def test_calc_zdr_offset():
+    xsapr_test_file = DATASETS.fetch("sgpxsaprcfrvptI4.a1.20200205.100827.nc")
+    ds = pyart.io.read(xsapr_test_file)
+    gatefilter = pyart.filters.GateFilter(ds)
+    gatefilter.exclude_below("cross_correlation_ratio_hv", 0.995)
+    gatefilter.exclude_above("cross_correlation_ratio_hv", 1)
+    gatefilter.exclude_below("reflectivity", 10)
+    gatefilter.exclude_above("reflectivity", 30)
+
+    results = pyart.correct.calc_zdr_offset(
+        ds,
+        zdr_var="differential_reflectivity",
+        gatefilter=gatefilter,
+        height_range=(1000, 3000),
+    )
+    assert_almost_equal(results["bias"], 2.69, decimal=2)
+    assert_almost_equal(results["profile_reflectivity"][15], 14.37, decimal=2)
+
+
 def test_calc_noise_floor():
     expected = [-46.25460013, -46.48371626, -46.3314618, -46.82639895, -46.76403711]