From bb1fe24da1112bf7820a3b88febc824fb343ff07 Mon Sep 17 00:00:00 2001
From: lauratomkins <lmtomkin@ncsu.edu>
Date: Thu, 9 Nov 2023 16:19:23 -0500
Subject: [PATCH] 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 <mgroverwx@gmail.com>
---
 pyart/retrieve/__init__.py |  1 +
 pyart/retrieve/qpe.py      | 55 ++++++++++++++++++++++++++++++++++++++
 tests/retrieve/test_qpe.py | 16 +++++++++++
 3 files changed, 72 insertions(+)

diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py
index 4666f6913d..1276e35038 100644
--- 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
index 6a4405af64..c657228d82 100644
--- 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
index e72ebb053c..98b39503e9 100644
--- 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)