Move functions and change math to numpy

pyart/retrieve/__init__.py

@@ -11,14 +11,14 @@
 from .echo_class import steiner_conv_strat  # noqa
 from .gate_id import fetch_radar_time_profile, map_profile_to_gates  # noqa
 from .kdp_proc import kdp_maesaka, kdp_schneebeli, kdp_vulpiani  # noqa
-from .precip_rate import ZtoR  # noqa
 from .qpe import est_rain_rate_a  # noqa
 from .qpe import est_rain_rate_hydro  # noqa
 from .qpe import est_rain_rate_kdp  # noqa
 from .qpe import est_rain_rate_z  # noqa
 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

pyart/retrieve/qpe.py

@@ -684,3 +684,57 @@ 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

tests/retrieve/test_qpe.py

@@ -1,6 +1,7 @@
 """ Unit tests for rainfall rate estimation module. """
 
 from numpy.testing import assert_allclose
+import numpy as np
 
 import pyart
 
@@ -153,3 +154,27 @@ 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()
+
+    # check calculations are within 10^-4 orders of magnitude
+    assert (
+        np.floor(
+            np.log10(
+                grid.fields["NWS_primary_prate"]["data"][0, 10, 10]
+                - (
+                    (
+                        (10 ** (grid.fields["reflectivity"]["data"][0, 10, 10] / 10))
+                        / 300
+                    )
+                    ** (1 / 1.4)
+                )
+            )
+        )
+        < -4
+    )