"""
Functions for rainfall rate estimation.
"""
from warnings import warn
import numpy as np
from ..config import get_field_name, get_fillvalue, get_metadata
from .echo_class import get_freq_band
[docs]def est_rain_rate_zpoly(radar, refl_field=None, rr_field=None):
"""
Estimates rainfall rate from reflectivity using a polynomial Z-R relation
developed at McGill University.
Parameters
----------
radar : Radar
Radar object.
refl_field : str, optional
Name of the reflectivity field to use.
rr_field : str, optional
Name of the rainfall rate field.
Returns
-------
rain : dict
Field dictionary containing the rainfall rate.
References
----------
Doelling et al. Systematic variations of Z–R-relationships from drop size
distributions measured in northern Germany during seven years. 1998. Atmos.
Ocean. Technol, 21, 1545-1556.
Joss et al. Operational Use of Radar for Precipitation Measurements
in Switzerland. 1998. Vdf Hochschulverlag AG ETH Zurich: 134.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
radar.check_field_exists(refl_field)
refl = radar.fields[refl_field]["data"]
refl2 = refl * refl
refl3 = refl * refl2
refl4 = refl * refl3
rr_data = np.ma.power(
10.0, -2.3 + 0.17 * refl - 5.1e-3 * refl2 + 9.8e-5 * refl3 - 6e-7 * refl4
)
rain = get_metadata(rr_field)
rain["data"] = rr_data
return rain
[docs]def est_rain_rate_z(radar, alpha=0.0376, beta=0.6112, refl_field=None, rr_field=None):
"""
Estimates rainfall rate from reflectivity using a power law.
Parameters
----------
radar : Radar
Radar object.
alpha, beta : floats, optional
Factor (alpha) and exponent (beta) of the power law.
refl_field : str, optional
Name of the reflectivity field to use.
rr_field : str, optional
Name of the rainfall rate field.
Returns
-------
rain : dict
Field dictionary containing the rainfall rate.
Reference
---------
Fabry, Frédéric. Radar Meterology. 2015. Ch 9. pg 148-165.
https://doi.org/10.1017/CBO9781107707405
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
radar.check_field_exists(refl_field)
refl = radar.fields[refl_field]["data"]
rr_data = alpha * np.ma.power(np.ma.power(10.0, 0.1 * refl), beta)
rain = get_metadata(rr_field)
rain["data"] = rr_data
return rain
[docs]def est_rain_rate_kdp(radar, alpha=None, beta=None, kdp_field=None, rr_field=None):
"""
Estimates rainfall rate from kdp using alpha power law.
Parameters
----------
radar : Radar
Radar object.
alpha, beta : floats, optional
Factor (alpha) and exponent (beta) of the power law. If not set the
factors are going to be determined according to the radar frequency.
kdp_field : str, optional
Name of the specific differential phase field to use.
rr_field : str, optional
Name of the rainfall rate field.
Returns
-------
rain : dict
Field dictionary containing the rainfall rate.
Reference
---------
Figueras et al. Long-term monitoring of French polarimetric radar data
quality and evaluation of several polarimetric quantitative precipitation
estimators in ideal conditions for operational implementation at C-band.
Quarterly Journal of the Royal Meteorological Society. 2012.
https://doi.org/10.1002/qj.1934
"""
# select the coefficients as alpha function of frequency band
if alpha is None or beta is None:
# assign coefficients according to radar frequency
if "frequency" in radar.instrument_parameters:
alpha, beta = _get_coeff_rkdp(
radar.instrument_parameters["frequency"]["data"][0]
)
else:
alpha, beta = _coeff_rkdp_table()["C"]
warn(
"Radar frequency unknown. "
"Default coefficients for C band will be applied."
)
# parse the field parameters
if kdp_field is None:
kdp_field = get_field_name("specific_differential_phase")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
radar.check_field_exists(kdp_field)
kdp = radar.fields[kdp_field]["data"]
kdp[kdp < 0] = 0.0
rr_data = alpha * np.ma.power(kdp, beta)
rain = get_metadata(rr_field)
rain["data"] = rr_data
return rain
[docs]def est_rain_rate_a(radar, alpha=None, beta=None, a_field=None, rr_field=None):
"""
Estimates rainfall rate from specific attenuation using alpha power law.
Parameters
----------
radar : Radar
Radar object.
alpha, beta : floats, optional
Factor (alpha) and exponent (beta) of the power law. If not set the
factors are going to be determined according to the radar frequency.
a_field : str, optional
Name of the specific attenuation field to use.
rr_field : str, optional
Name of the rainfall rate field.
Returns
-------
rain : dict
Field dictionary containing the rainfall rate.
References
----------
Diederich M., Ryzhkov A., Simmer C., Zhang P. and Tromel S., 2015: Use of
Specific Attenuation for Rainfall Measurement at X-Band Radar Wavelenghts.
Part I: Radar Calibration and Partial Beam Blockage Estimation. Journal of
Hydrometeorology, 16, 487-502.
Ryzhkov A., Diederich M., Zhang P. and Simmer C., 2014: Potential
Utilization of Specific Attenuation for Rainfall Estimation, Mitigation of
Partial Beam Blockage, and Radar Networking. Journal of Atmospheric and
Oceanic Technology, 31, 599-619.
"""
# select the coefficients as alpha function of frequency band
if alpha is None or beta is None:
# assign coefficients according to radar frequency
if "frequency" in radar.instrument_parameters:
alpha, beta = _get_coeff_ra(
radar.instrument_parameters["frequency"]["data"][0]
)
else:
alpha, beta = _coeff_ra_table()["C"]
warn(
"Radar frequency unknown. "
"Default coefficients for C band will be applied."
)
# parse the field parameters
if a_field is None:
a_field = get_field_name("specific_attenuation")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
radar.check_field_exists(a_field)
att = radar.fields[a_field]["data"]
rr_data = alpha * np.ma.power(att, beta)
rain = get_metadata(rr_field)
rain["data"] = rr_data
return rain
[docs]def est_rain_rate_zkdp(
radar,
alphaz=0.0376,
betaz=0.6112,
alphakdp=None,
betakdp=None,
refl_field=None,
kdp_field=None,
rr_field=None,
main_field=None,
thresh=None,
thresh_max=True,
):
"""
Estimates rainfall rate from a blending of power law r-kdp and r-z
relations.
Parameters
----------
radar : Radar
Radar object.
alphaz, betaz : floats, optional
Factor (alpha) and exponent (beta) of the z-r power law.
alphakdp, betakdp : floats, optional
Factor (alpha) and exponent (beta) of the kdp-r power law.
If not set the factors are going to be determined according
to the radar frequency.
refl_field : str, optional
Name of the reflectivity field to use.
kdp_field : str, optional
Name of the specific differential phase field to use.
rr_field : str, optional
Name of the rainfall rate field.
main_field : str, optional
Name of the field that is going to act as main. Has to be
either refl_field or kdp_field. Default is refl_field.
thresh : float, optional
Value of the threshold that determines when to use the secondary
field.
thresh_max : Bool, optional
If true the main field is used up to the thresh value maximum.
Otherwise the main field is not used below thresh value.
Returns
-------
rain_main : dict
Field dictionary containing the rainfall rate.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if kdp_field is None:
kdp_field = get_field_name("specific_differential_phase")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
rain_z = est_rain_rate_z(
radar, alpha=alphaz, beta=betaz, refl_field=refl_field, rr_field=rr_field
)
rain_kdp = est_rain_rate_kdp(
radar, alpha=alphakdp, beta=betakdp, kdp_field=kdp_field, rr_field=rr_field
)
if main_field == refl_field:
rain_main = rain_z
rain_secondary = rain_kdp
elif main_field == kdp_field:
rain_main = rain_kdp
rain_secondary = rain_z
elif main_field is None:
main_field = refl_field
rain_main = rain_z
rain_secondary = rain_kdp
else:
main_field = refl_field
rain_main = rain_z
rain_secondary = rain_kdp
thresh = 40.0
thresh_max = True
warn(
"Unknown main field. Using " + refl_field + " with threshold " + str(thresh)
)
if thresh_max:
is_secondary = rain_main["data"] > thresh
else:
is_secondary = rain_main["data"] < thresh
rain_main["data"][is_secondary] = rain_secondary["data"][is_secondary]
return rain_main
[docs]def est_rain_rate_za(
radar,
alphaz=0.0376,
betaz=0.6112,
alphaa=None,
betaa=None,
refl_field=None,
a_field=None,
rr_field=None,
main_field=None,
thresh=None,
thresh_max=False,
):
"""
Estimates rainfall rate from a blending of power law r-alpha and r-z
relations.
Parameters
----------
radar : Radar
Radar object
alphaz, betaz : floats, optional
Factor (alpha) and exponent (beta) of the z-r power law.
alphaa,betaa : floats, optional
Factor (alpha) and exponent (beta) of the a-r power law. If not set
the factors are going to be determined according to the radar frequency.
refl_field : str, optional
Name of the reflectivity field to use.
a_field : str, optional
Name of the specific attenuation field to use.
rr_field : str, optional
Name of the rainfall rate field.
main_field : str, optional
Name of the field that is going to act as main. Has to be
either refl_field or kdp_field. Default is refl_field.
thresh : float, optional
Value of the threshold that determines when to use the secondary
field.
thresh_max : Bool, optional
If true the main field is used up to the thresh value maximum.
Otherwise the main field is not used below thresh value.
Returns
-------
rain_main : dict
Field dictionary containing the rainfall rate.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if a_field is None:
a_field = get_field_name("specific_attenuation")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
rain_z = est_rain_rate_z(
radar, alpha=alphaz, beta=betaz, refl_field=refl_field, rr_field=rr_field
)
rain_a = est_rain_rate_a(
radar, alpha=alphaa, beta=betaa, a_field=a_field, rr_field=rr_field
)
if main_field == refl_field:
rain_main = rain_z
rain_secondary = rain_a
elif main_field == a_field:
rain_main = rain_a
rain_secondary = rain_z
elif main_field is None:
main_field = a_field
rain_main = rain_a
rain_secondary = rain_z
else:
main_field = a_field
rain_main = rain_a
rain_secondary = rain_z
thresh = 0.04
thresh_max = False
warn("Unknown main field. Using " + a_field + " with threshold " + str(thresh))
if thresh_max:
is_secondary = rain_main["data"] > thresh
else:
is_secondary = rain_main["data"] < thresh
rain_main["data"][is_secondary] = rain_secondary["data"][is_secondary]
return rain_main
[docs]def est_rain_rate_hydro(
radar,
alphazr=0.0376,
betazr=0.6112,
alphazs=0.1,
betazs=0.5,
alphaa=None,
betaa=None,
mp_factor=0.6,
refl_field=None,
a_field=None,
hydro_field=None,
rr_field=None,
main_field=None,
thresh=None,
thresh_max=False,
):
"""
Estimates rainfall rate using different relations between R and the
polarimetric variables depending on the hydrometeor type.
Parameters
----------
radar : Radar
Radar object.
alphazr, betazr : floats, optional
Factor (alpha) and exponent (beta) of the z-r power law for rain.
alphazs, betazs : floats, optional
Factor (alpha) and exponent (beta) of the z-s power law for snow.
alphaa, betaa : floats, optional
Factor (alpha) and exponent (beta) of the a-r power law.
If not set the factors are going to be determined according
to the radar frequency.
mp_factor : float, optional
Factor applied to z-r relation in the melting layer.
refl_field : str, optional
Name of the reflectivity field to use.
a_field : str, optional
Name of the specific attenuation field to use.
hydro_field : str, optional
Name of the hydrometeor classification field to use.
rr_field : str, optional
Name of the rainfall rate field.
main_field : str, optional
Name of the field that is going to act as main. Has to be
either refl_field or kdp_field. Default is refl_field.
thresh : float, optional
Value of the threshold that determines when to use the secondary
field.
thresh_max : Bool, optional
If true the main field is used up to the thresh value maximum.
Otherwise the main field is not used below thresh value.
Returns
-------
rain : dict
Field dictionary containing the rainfall rate.
References
----------
Besic, N., Figueras i Ventura, J., Grazioli, J., Gabella, M., Germann, U., and Berne, A.: Hydrometeor classification
through statistical clustering of polarimetric radar measurements: a semi-supervised approach,
Atmos. Meas. Tech., 9, 4425–4445, https://doi.org/10.5194/amt-9-4425-2016, 2016.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if a_field is None:
a_field = get_field_name("specific_attenuation")
if hydro_field is None:
hydro_field = get_field_name("radar_echo_classification")
if rr_field is None:
rr_field = get_field_name("radar_estimated_rain_rate")
# extract fields and parameters from radar
if hydro_field in radar.fields:
hydroclass = radar.fields[hydro_field]["data"]
else:
raise KeyError("Field not available: " + hydro_field)
# get the location of each hydrometeor class
is_ds = hydroclass == 1
is_cr = hydroclass == 2
is_lr = hydroclass == 3
is_gr = hydroclass == 4
is_rn = hydroclass == 5
is_vi = hydroclass == 6
is_ws = hydroclass == 7
is_mh = hydroclass == 8
is_ih = hydroclass == 9
# compute z-r (in rain) z-r in snow and z-a relations
rain_z = est_rain_rate_z(
radar, alpha=alphazr, beta=betazr, refl_field=refl_field, rr_field=rr_field
)
snow_z = est_rain_rate_z(
radar, alpha=alphazs, beta=betazs, refl_field=refl_field, rr_field=rr_field
)
rain_a = est_rain_rate_a(
radar, alpha=alphaa, beta=betaa, a_field=a_field, rr_field=rr_field
)
# initialize rainfall rate field
rr_data = np.ma.zeros(hydroclass.shape, dtype="float32")
rr_data[:] = np.ma.masked
rr_data.set_fill_value(get_fillvalue())
# apply the relations for each type
# solid phase
rr_data[is_ds] = snow_z["data"][is_ds]
rr_data[is_cr] = snow_z["data"][is_cr]
rr_data[is_vi] = snow_z["data"][is_vi]
rr_data[is_gr] = snow_z["data"][is_gr]
rr_data[is_ih] = snow_z["data"][is_ih]
# rain
if main_field == refl_field:
rain_main = rain_z
rain_secondary = rain_a
elif main_field == a_field:
rain_main = rain_a
rain_secondary = rain_z
elif main_field is None:
main_field = a_field
rain_main = rain_a
rain_secondary = rain_z
else:
main_field = a_field
rain_main = rain_a
rain_secondary = rain_z
thresh = 0.04
thresh_max = False
warn("Unknown main field. Using " + a_field + " with threshold " + str(thresh))
if thresh_max:
is_secondary = rain_main["data"] > thresh
else:
is_secondary = rain_main["data"] < thresh
rain_main["data"][is_secondary] = rain_secondary["data"][is_secondary]
rr_data[is_lr] = rain_main["data"][is_lr]
rr_data[is_rn] = rain_main["data"][is_rn]
# mixed phase
rr_data[is_ws] = mp_factor * rain_z["data"][is_ws]
rr_data[is_mh] = mp_factor * rain_z["data"][is_mh]
rain = get_metadata(rr_field)
rain["data"] = rr_data
return rain
def _get_coeff_rkdp(freq):
"""
Get the R(kdp) power law coefficients for a particular frequency.
Parameters
----------
freq : float
Radar frequency [Hz].
Returns
-------
alpha, beta : floats
The coefficient and exponent of the power law.
"""
coeff_rkdp_dict = _coeff_rkdp_table()
freq_band = get_freq_band(freq)
if (freq_band is not None) and (freq_band in coeff_rkdp_dict):
return coeff_rkdp_dict[freq_band]
if freq < 2e9:
freq_band_aux = "S"
elif freq > 12e9:
freq_band_aux = "X"
warn(
"Radar frequency out of range. "
+ "Coefficients only applied to S, C or X band. "
+ freq_band_aux
+ " band coefficients will be used."
)
return coeff_rkdp_dict[freq_band_aux]
def _coeff_rkdp_table():
"""
Defines the R(kdp) power law coefficients for each frequency band.
Returns
-------
coeff_rkdp_dict : dict
A dictionary with the coefficients at each band.
"""
coeff_rkdp_dict = dict()
# S band: Beard and Chuang coefficients
coeff_rkdp_dict.update({"S": (50.70, 0.8500)})
# C band: Beard and Chuang coefficients
coeff_rkdp_dict.update({"C": (29.70, 0.8500)})
# X band: Brandes coefficients
coeff_rkdp_dict.update({"X": (15.81, 0.7992)})
return coeff_rkdp_dict
def _get_coeff_ra(freq):
"""
Get the R(A) power law coefficients for a particular frequency.
Parameters
----------
freq : float
Radar frequency [Hz].
Returns
-------
alpha, beta : floats
The coefficient and exponent of the power law.
"""
coeff_ra_dict = _coeff_ra_table()
freq_band = get_freq_band(freq)
if (freq_band is not None) and (freq_band in coeff_ra_dict):
return coeff_ra_dict[freq_band]
if freq < 2e9:
freq_band_aux = "S"
elif freq > 12e9:
freq_band_aux = "X"
warn(
"Radar frequency out of range. "
+ "Coefficients only applied to S, C or X band. "
+ freq_band_aux
+ " band coefficients will be used."
)
return coeff_ra_dict[freq_band_aux]
def _coeff_ra_table():
"""
Defines the R(A) power law coefficients for each frequency band.
Returns
-------
coeff_ra_dict : dict
A dictionary with the coefficients at each band.
"""
coeff_ra_dict = dict()
# S band: at 10 deg C according to tables from Ryzhkov et al. 2014
coeff_ra_dict.update({"S": (3100.0, 1.03)})
# C band: at 10 deg C according to tables from Diederich et al. 2015
coeff_ra_dict.update({"C": (250.0, 0.91)})
# X band: at 10 deg C according to tables from Diederich et al. 2015
coeff_ra_dict.update({"X": (45.5, 0.83)})
return coeff_ra_dict
[docs]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