From 92aa9b870ab3850f9890707c4ab6f5603a2031d1 Mon Sep 17 00:00:00 2001 From: Bobby Jackson Date: Fri, 30 Aug 2024 09:34:55 -0500 Subject: [PATCH] ADD: ZDR bias code from RadTraQ (#1630) * ADD: ZDR bias calculations from Radtraq * ADD: Example for gallery * ADD: ZDR bias calculations * FIX: Revert back for spectra_calculations.py * FIX: Restore phase proc from main branch * STY: Black and ruff changes. * FIX: Duplicate entries in __init__ and obj to radar * FIX: Reinsert calc_bias * FIX: Never mind, there were two calc_biases --------- Co-authored-by: Zach Sherman <19153455+zssherman@users.noreply.github.com> --- examples/correct/plot_zdr_check.py | 50 ++++++++++++++++++++++++++ pyart/correct/__init__.py | 1 + pyart/correct/bias_and_noise.py | 57 ++++++++++++++++++++++++++++++ tests/correct/test_correct_bias.py | 19 ++++++++++ 4 files changed, 127 insertions(+) create mode 100644 examples/correct/plot_zdr_check.py diff --git a/examples/correct/plot_zdr_check.py b/examples/correct/plot_zdr_check.py new file mode 100644 index 0000000000..df98780e0c --- /dev/null +++ b/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() diff --git a/pyart/correct/__init__.py b/pyart/correct/__init__.py index 6a9812ec77..7b050057ad 100644 --- a/pyart/correct/__init__.py +++ b/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 diff --git a/pyart/correct/bias_and_noise.py b/pyart/correct/bias_and_noise.py index 295a6689fa..866a290cc9 100755 --- a/pyart/correct/bias_and_noise.py +++ b/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, diff --git a/tests/correct/test_correct_bias.py b/tests/correct/test_correct_bias.py index 3be079411f..aec10fa576 100644 --- a/tests/correct/test_correct_bias.py +++ b/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]