ADD: function to rescale reflectivity to a precipitation rate (#1486)

* Create precip_rate.py * Update __init__.py * Create test_precip_rate.py * Update precip_rate.py * Linting formatting * Update test_precip_rate.py * Update test_precip_rate.py * Move functions and change math to numpy * Update formatting * Update formatting * FIX: Fix the testing suite to avoid nans --------- Co-authored-by: mgrover1 <[email protected]>
ARM-DOE · Nov 9, 2023 · bb1fe24 · bb1fe24
1 parent d58b687
commit bb1fe24
Show file tree

Hide file tree

Showing 3 changed files with 72 additions and 0 deletions.
diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py
@@ -18,6 +18,7 @@
 from .qpe import est_rain_rate_za  # noqa
 from .qpe import est_rain_rate_zkdp  # noqa
 from .qpe import est_rain_rate_zpoly  # noqa
+from .qpe import ZtoR  # noqa
 from .qvp import quasi_vertical_profile  # noqa
 from .simple_moment_calculations import calculate_snr_from_reflectivity  # noqa
 from .simple_moment_calculations import calculate_velocity_texture  # noqa

diff --git a/pyart/retrieve/qpe.py b/pyart/retrieve/qpe.py
@@ -684,3 +684,58 @@ def _coeff_ra_table():
     coeff_ra_dict.update({"X": (45.5, 0.83)})
 
     return coeff_ra_dict
+
+
+def ZtoR(radar, ref_field="reflectivity", a=300, b=1.4, save_name="NWS_primary_prate"):
+    """
+    Convert reflectivity (dBZ) to precipitation rate (mm/hr)
+
+    Author: Laura Tomkins
+
+    Parameters
+    ----------
+    radar : Radar
+        Radar object used.
+    ref_field : str
+        Reflectivity field name to use to look up reflectivity data. In the
+        radar object. Default field name is 'reflectivity'. Units are expected
+        to be dBZ.
+    a : float
+        a value (coefficient) in the Z-R relationship
+    b: float
+        b value (exponent) in the Z-R relationship
+
+    Returns
+    -------
+    radar : Radar
+        The radar object containing the precipitation rate field
+
+    References
+    ----------
+    American Meteorological Society, 2022: "Z-R relation". Glossary of Meteorology,
+    https://glossary.ametsoc.org/wiki/Z-r_relation
+
+    """
+
+    # get reflectivity data
+    ref_data = radar.fields[ref_field]["data"]
+    ref_data = np.ma.masked_invalid(ref_data)
+
+    # convert to linear reflectivity
+    ref_linear = 10 ** (ref_data / 10)
+    precip_rate = (ref_linear / a) ** (1 / b)
+
+    # create dictionary
+    prate_dict = {
+        "data": precip_rate,
+        "standard_name": save_name,
+        "long_name": f"{save_name} rescaled from linear reflectivity",
+        "units": "mm/hr",
+        "valid_min": 0,
+        "valid_max": 10000,
+    }
+
+    # add field to radar object
+    radar.add_field(save_name, prate_dict, replace_existing=True)
+
+    return radar
diff --git a/tests/retrieve/test_qpe.py b/tests/retrieve/test_qpe.py
@@ -153,3 +153,19 @@ def test_get_coeff_rkdp():
 
     coeff_rkdp_use_x = pyart.retrieve.qpe._get_coeff_rkdp(13e9)
     assert coeff_rkdp_use_x == (15.81, 0.7992)
+
+
+def test_precip_rate():
+    grid = pyart.testing.make_storm_grid()
+    grid = pyart.retrieve.ZtoR(grid)
+
+    # check that field is in grid object
+    assert "NWS_primary_prate" in grid.fields.keys()
+
+    # Calculate the estimated value
+    correct = (
+        (10 ** (grid.fields["reflectivity"]["data"][0, 10, 10] / 10.0)) / 300.0
+    ) ** (1 / 1.4)
+
+    # Check for correctness
+    assert_allclose(grid.fields["NWS_primary_prate"]["data"][0, 10, 10], correct)