From 3653c2eacf9cc207286a30f7999a1efdf3ae26b0 Mon Sep 17 00:00:00 2001 From: Bhupendra Raut Date: Tue, 27 Aug 2024 09:44:16 -0500 Subject: [PATCH 1/7] Updated Documentation for retrieve/hydroclass_semisupervised function (#1616) * class names added to the documentation * example added * description * description refined * Corrected spelling/punctuations --- pyart/retrieve/echo_class.py | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index ac19ce5e826..8057704ac22 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -612,8 +612,10 @@ def hydroclass_semisupervised( radar_freq=None, ): """ - Classifies precipitation echoes following the approach by Besic et al - (2016). + Classifies precipitation echoes into hydrometeor types. + + The `hydroclass_semisupervised` function classifies precipitation echoes in the polarimetric radar data + into 9 hydrometeor types using a semi-supervised approach (Besic et al., 2016). Parameters ---------- @@ -642,7 +644,17 @@ def hydroclass_semisupervised( Returns ------- hydro : dict - Hydrometeor classification field. + Hydrometeor classification. + - 0: Not classified + - 1: Aggregates + - 2: Ice crystals + - 3: Light rain + - 4: Rimed particles + - 5: Rain + - 6: Vertically oriented ice + - 7: Wet snow + - 8: Melting hail + - 9: Dry hail or high-density graupel References ---------- @@ -651,6 +663,18 @@ def hydroclass_semisupervised( of polarimetric radar measurements: a semi-supervised approach, Atmos. Meas. Tech., 9, 4425-4445, doi:10.5194/amt-9-4425-2016, 2016 + Usage + ----- + .. code-block:: python + hydro_class = pyart.retrieve.hydroclass_semisupervised( + radar, + refl_field="corrected_reflectivity", + zdr_field="corrected_differential_reflectivity", + kdp_field="specific_differential_phase", + rhv_field="uncorrected_cross_correlation_ratio", + temp_field="temperature", + ) + Notes ----- The default hydrometeor classification is valid for C-band radars. For X-band radars, @@ -660,8 +684,8 @@ def hydroclass_semisupervised( If the radar frequency information is missing from the radar object, you can add it in `radar.instrument_parameters`, as follows: - .. code-block:: python + .. code-block:: python radar.instrument_parameters["frequency"] = { "long_name": "Radar frequency", "units": "Hz", From 664865b552092c944a877db8e3644a0feaa9adbf Mon Sep 17 00:00:00 2001 From: Max Grover Date: Tue, 27 Aug 2024 14:52:22 -0500 Subject: [PATCH 2/7] Add more references (#1620) * ADD: Add reference for compute_cdf * ADD: Add reference for computations * ADD: Add qpe references --- pyart/retrieve/qpe.py | 6 ++++++ pyart/retrieve/simple_moment_calculations.py | 14 ++++++++++++++ 2 files changed, 20 insertions(+) diff --git a/pyart/retrieve/qpe.py b/pyart/retrieve/qpe.py index c657228d826..ad8b28856ba 100644 --- a/pyart/retrieve/qpe.py +++ b/pyart/retrieve/qpe.py @@ -475,6 +475,12 @@ def est_rain_rate_hydro( 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: diff --git a/pyart/retrieve/simple_moment_calculations.py b/pyart/retrieve/simple_moment_calculations.py index 61beb2db04e..0dfd32c6291 100644 --- a/pyart/retrieve/simple_moment_calculations.py +++ b/pyart/retrieve/simple_moment_calculations.py @@ -165,6 +165,11 @@ def compute_l(radar, rhohv_field=None, l_field=None): l_field_out : dict L field. + References + ---------- + Ryzhkov, A. V., 2001: Interpretation of Polarimetric Radar Covariance Matrix for Meteorological Scatterers: Theoretical Analysis. + J. Atmos. Oceanic Technol., 18, 315–328, https://doi.org/10.1175/1520-0426(2001)018<0315:IOPRCM>2.0.CO;2. + """ # parse the field parameters if rhohv_field is None: @@ -191,6 +196,8 @@ def compute_cdr(radar, rhohv_field=None, zdr_field=None, cdr_field=None): """ Computes the Circular Depolarization Ratio. + + Parameters ---------- radar : Radar @@ -207,6 +214,13 @@ def compute_cdr(radar, rhohv_field=None, zdr_field=None, cdr_field=None): cdr : dict CDR field. + References + ---------- + Matrosov, S. Y., 2004: Depolarization Estimates from Linear H and V Measurements with Weather + Radars Operating in Simultaneous Transmission–Simultaneous Receiving Mode. + J. Atmos. Oceanic Technol., 21, 574–583, + https://doi.org/10.1175/1520-0426(2004)021<0574:DEFLHA>2.0.CO;2. + """ # parse the field parameters if rhohv_field is None: From 48892899b2cdce1f9274224250e3d93d3de20157 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Thu, 29 Aug 2024 14:37:50 +0200 Subject: [PATCH 3/7] Addition of qvp functions from MeteoSwiss Py-ART + support of demixing in hydroclass_semisupervised (#1618) * ENH: addition of qvp functions from MeteoSwiss Py-ART + entropy computation in hydroclass_semisupervised * FIX: fixing est_rain_rate_hydro to match ARM original * ENH: black formatting * FIX: modified pyart examples to handle new hydroclass_semisupervised output * FIX: fix in plot_hydrometeor_class_x_band.py example which was using wrong hydroclass labels * ENH: updated defaut_config.py to account for new outputs of hydroclass_semisupervised * FIX: adapted hydroclass_semisupervised to account for specified radar_freq * FIX: adapted hydroclass_semisupervised to account for specified radar_freq * ADD: new tests for missing fcts * ENH: add support for more scan modes in qvp fcts * FIX: docstring correction in cartesian_to_antenna * FIX: missing hydrometeor variables in default_config.py * fixes to make tests pass * fixes to make tests pass * FIX: changed value in ZDR standardization from -5 to 5 to -1.5 to 5 according to Besic (2016) --- examples/correct/plot_attenuation.py | 1 + examples/correct/plot_dealias.py | 1 + examples/io/plot_nexrad_data_aws.py | 1 + examples/io/plot_older_nexrad_data_aws.py | 1 + .../plot_compare_two_radars_gatemapper.py | 1 + .../mapping/plot_grid_single_sweep_ppi.py | 1 + .../mapping/plot_map_one_radar_to_grid.py | 1 + .../mapping/plot_map_two_radars_to_grid.py | 1 + examples/plotting/plot_choose_a_colormap.py | 1 + examples/plotting/plot_cross_section.py | 1 + examples/plotting/plot_modify_colorbar.py | 1 + .../plotting/plot_nexrad_multiple_moments.py | 1 + examples/plotting/plot_nexrad_reflectivity.py | 1 + examples/plotting/plot_ppi_cfradial.py | 1 + examples/plotting/plot_ppi_mdv.py | 1 + examples/plotting/plot_rhi_cfradial.py | 1 + .../plotting/plot_rhi_cfradial_singlescan.py | 1 + ..._rhi_contours_differential_reflectivity.py | 1 + examples/plotting/plot_rhi_data_overlay.py | 1 + examples/plotting/plot_rhi_mdv.py | 1 + examples/plotting/plot_rhi_two_panel.py | 1 + examples/plotting/plot_xsect.py | 1 + examples/plotting/radar-cross-section.ipynb | 5 +- examples/retrieve/column-example.ipynb | 7 +- examples/retrieve/plot_hydrometeor.py | 2 +- .../retrieve/plot_hydrometeor_class_x_band.py | 22 +- .../retrieve/wavelet_echo_class_example.ipynb | 10 +- pyart/__check_build/__init__.py | 1 + pyart/core/__init__.py | 1 + pyart/core/transforms.py | 27 + pyart/core/wind_profile.py | 1 - pyart/default_config.py | 100 + pyart/graph/_cm_colorblind.py | 1 + pyart/graph/gridmapdisplay_basemap.py | 1 - pyart/io/mdv_grid.py | 1 - pyart/io/nexrad_common.py | 1 + pyart/io/nexrad_level2.py | 1 - pyart/io/nexrad_level3.py | 2 +- pyart/map/gates_to_grid.py | 1 + pyart/retrieve/__init__.py | 3 +- pyart/retrieve/_echo_class_wt.py | 1 - pyart/retrieve/echo_class.py | 462 ++++- pyart/retrieve/qvp.py | 1744 ++++++++++++++++- pyart/util/__init__.py | 2 + pyart/util/circular_stats.py | 35 + pyart/util/radar_utils.py | 30 + pyart/xradar/accessor.py | 1 - tests/bridge/test_wradlib_bridge.py | 1 - tests/core/test_transforms.py | 12 + tests/graph/test_cm.py | 1 - tests/graph/test_cm_colorblind.py | 1 - tests/graph/test_gridmapdisplay.py | 1 + tests/graph/test_gridmapdisplay_basemap.py | 1 + tests/graph/test_radarmapdisplay.py | 1 + tests/graph/test_radarmapdisplay_basemap.py | 1 + tests/retrieve/test_echo_class.py | 57 +- tests/retrieve/test_qvp.py | 383 ++++ tests/util/test_circular_stats.py | 30 + tests/util/test_radar_utils.py | 10 + 59 files changed, 2829 insertions(+), 154 deletions(-) diff --git a/examples/correct/plot_attenuation.py b/examples/correct/plot_attenuation.py index 469b2d049fa..6bdbb3ee745 100644 --- a/examples/correct/plot_attenuation.py +++ b/examples/correct/plot_attenuation.py @@ -7,6 +7,7 @@ for a polarimetric radar using a Z-PHI method implemented in Py-ART. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/correct/plot_dealias.py b/examples/correct/plot_dealias.py index 8687b22c81b..e33a3acb1f8 100644 --- a/examples/correct/plot_dealias.py +++ b/examples/correct/plot_dealias.py @@ -6,6 +6,7 @@ In this example doppler velocities are dealiased using the ial condition of the dealiasing, using the region-based dealiasing algorithm in Py-ART. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/io/plot_nexrad_data_aws.py b/examples/io/plot_nexrad_data_aws.py index 7a0a055191c..a9eabc55412 100644 --- a/examples/io/plot_nexrad_data_aws.py +++ b/examples/io/plot_nexrad_data_aws.py @@ -7,6 +7,7 @@ and plot quick looks of the datasets. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/io/plot_older_nexrad_data_aws.py b/examples/io/plot_older_nexrad_data_aws.py index b50f8d81c64..e650da6fb3f 100644 --- a/examples/io/plot_older_nexrad_data_aws.py +++ b/examples/io/plot_older_nexrad_data_aws.py @@ -7,6 +7,7 @@ to 2008 that are missing some coordinate metadata. """ + print(__doc__) diff --git a/examples/mapping/plot_compare_two_radars_gatemapper.py b/examples/mapping/plot_compare_two_radars_gatemapper.py index fe57ea4557f..8b8274218ee 100644 --- a/examples/mapping/plot_compare_two_radars_gatemapper.py +++ b/examples/mapping/plot_compare_two_radars_gatemapper.py @@ -7,6 +7,7 @@ another radar in Antenna coordinates and compare the fields. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) and Bobby Jackson (rjackson@anl.gov) diff --git a/examples/mapping/plot_grid_single_sweep_ppi.py b/examples/mapping/plot_grid_single_sweep_ppi.py index f7c5eca8ac6..beac728837d 100644 --- a/examples/mapping/plot_grid_single_sweep_ppi.py +++ b/examples/mapping/plot_grid_single_sweep_ppi.py @@ -13,6 +13,7 @@ negligible at low elevation angles common to PPI sweeps. """ + print(__doc__) # ===================== diff --git a/examples/mapping/plot_map_one_radar_to_grid.py b/examples/mapping/plot_map_one_radar_to_grid.py index fc2b9dfd561..9b5280ee7dd 100644 --- a/examples/mapping/plot_map_one_radar_to_grid.py +++ b/examples/mapping/plot_map_one_radar_to_grid.py @@ -7,6 +7,7 @@ Cartesian grid. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/mapping/plot_map_two_radars_to_grid.py b/examples/mapping/plot_map_two_radars_to_grid.py index 7298ea43568..b4b0c1cbfbb 100644 --- a/examples/mapping/plot_map_two_radars_to_grid.py +++ b/examples/mapping/plot_map_two_radars_to_grid.py @@ -7,6 +7,7 @@ coordinates to a Cartesian grid. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_choose_a_colormap.py b/examples/plotting/plot_choose_a_colormap.py index acc1e9880df..c322a93442f 100644 --- a/examples/plotting/plot_choose_a_colormap.py +++ b/examples/plotting/plot_choose_a_colormap.py @@ -7,6 +7,7 @@ and how to add them to your own plots. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_cross_section.py b/examples/plotting/plot_cross_section.py index aafa9e3b790..420d4b40153 100644 --- a/examples/plotting/plot_cross_section.py +++ b/examples/plotting/plot_cross_section.py @@ -7,6 +7,7 @@ of your radar grid using the GridMapDisplay """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_modify_colorbar.py b/examples/plotting/plot_modify_colorbar.py index 9451f897ca7..44495da9503 100644 --- a/examples/plotting/plot_modify_colorbar.py +++ b/examples/plotting/plot_modify_colorbar.py @@ -7,6 +7,7 @@ within a Py-ART display object. """ + print(__doc__) # Author: Joe O'Brien (obrienj@anl.gov) diff --git a/examples/plotting/plot_nexrad_multiple_moments.py b/examples/plotting/plot_nexrad_multiple_moments.py index b6c03fa9a26..7e6cc6041d7 100644 --- a/examples/plotting/plot_nexrad_multiple_moments.py +++ b/examples/plotting/plot_nexrad_multiple_moments.py @@ -7,6 +7,7 @@ NEXRAD Archive file. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_nexrad_reflectivity.py b/examples/plotting/plot_nexrad_reflectivity.py index a3e73f0edca..49f12bccf39 100644 --- a/examples/plotting/plot_nexrad_reflectivity.py +++ b/examples/plotting/plot_nexrad_reflectivity.py @@ -7,6 +7,7 @@ NEXRAD file. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_ppi_cfradial.py b/examples/plotting/plot_ppi_cfradial.py index 982e6602261..33f18ea8542 100644 --- a/examples/plotting/plot_ppi_cfradial.py +++ b/examples/plotting/plot_ppi_cfradial.py @@ -6,6 +6,7 @@ An example which creates a PPI plot of a Cfradial file. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_ppi_mdv.py b/examples/plotting/plot_ppi_mdv.py index 706495bf3a6..2b1036e41f6 100644 --- a/examples/plotting/plot_ppi_mdv.py +++ b/examples/plotting/plot_ppi_mdv.py @@ -6,6 +6,7 @@ An example which creates a PPI plot of a MDV file using a RadarDisplay object. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_rhi_cfradial.py b/examples/plotting/plot_rhi_cfradial.py index 5878a287064..fa85dd5fe09 100755 --- a/examples/plotting/plot_rhi_cfradial.py +++ b/examples/plotting/plot_rhi_cfradial.py @@ -7,6 +7,7 @@ a RadarDisplay object. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_rhi_cfradial_singlescan.py b/examples/plotting/plot_rhi_cfradial_singlescan.py index 088b74d2c37..9f32eeedcf5 100755 --- a/examples/plotting/plot_rhi_cfradial_singlescan.py +++ b/examples/plotting/plot_rhi_cfradial_singlescan.py @@ -7,6 +7,7 @@ a RadarDisplay object. """ + print(__doc__) import matplotlib.pyplot as plt diff --git a/examples/plotting/plot_rhi_contours_differential_reflectivity.py b/examples/plotting/plot_rhi_contours_differential_reflectivity.py index bde4afbb58a..f64a9d0b8ca 100644 --- a/examples/plotting/plot_rhi_contours_differential_reflectivity.py +++ b/examples/plotting/plot_rhi_contours_differential_reflectivity.py @@ -7,6 +7,7 @@ and adding differnential Reflectivity contours from the same MDV file. """ + print(__doc__) # Author: Cory Weber (cweber@anl.gov) diff --git a/examples/plotting/plot_rhi_data_overlay.py b/examples/plotting/plot_rhi_data_overlay.py index 2e6a901e29d..c65b95f9843 100644 --- a/examples/plotting/plot_rhi_data_overlay.py +++ b/examples/plotting/plot_rhi_data_overlay.py @@ -7,6 +7,7 @@ and adding Reflectivity contours from the same MDV file. """ + print(__doc__) # Author: Cory Weber (cweber@anl.gov) diff --git a/examples/plotting/plot_rhi_mdv.py b/examples/plotting/plot_rhi_mdv.py index c5cfcfea67a..e6e5ca9cf4b 100644 --- a/examples/plotting/plot_rhi_mdv.py +++ b/examples/plotting/plot_rhi_mdv.py @@ -6,6 +6,7 @@ An example which creates a RHI plot of a MDV file using a RadarDisplay object. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_rhi_two_panel.py b/examples/plotting/plot_rhi_two_panel.py index 9d6b84ae890..9efb26a60b8 100755 --- a/examples/plotting/plot_rhi_two_panel.py +++ b/examples/plotting/plot_rhi_two_panel.py @@ -7,6 +7,7 @@ included in the two panels are reflectivity and doppler velocity. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_xsect.py b/examples/plotting/plot_xsect.py index a64fcca70c4..5bc92beed34 100644 --- a/examples/plotting/plot_xsect.py +++ b/examples/plotting/plot_xsect.py @@ -7,6 +7,7 @@ of PPI scans and plots both cross sections. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/radar-cross-section.ipynb b/examples/plotting/radar-cross-section.ipynb index e5389e7d99b..4fe48a2ca0a 100644 --- a/examples/plotting/radar-cross-section.ipynb +++ b/examples/plotting/radar-cross-section.ipynb @@ -8,6 +8,7 @@ "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", + "\n", "import pyart\n", "from pyart.testing import get_test_data\n", "\n", @@ -344,7 +345,7 @@ "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 8))\n", - "title = \"($^\\circ$N)\"\n", + "title = r\"($^\\circ$N)\"\n", "cross.corrected_reflectivity_horizontal.plot(\n", " cmap=\"pyart_HomeyerRainbow\",\n", " x=None,\n", @@ -362,7 +363,7 @@ "metadata": {}, "outputs": [], "source": [ - "from pyart.graph.common import generate_grid_time_begin, generate_field_name\n", + "from pyart.graph.common import generate_field_name, generate_grid_time_begin\n", "\n", "\n", "def generate_cross_section_title(grid, field, start, end):\n", diff --git a/examples/retrieve/column-example.ipynb b/examples/retrieve/column-example.ipynb index eeeca093759..5f6beb51fe1 100644 --- a/examples/retrieve/column-example.ipynb +++ b/examples/retrieve/column-example.ipynb @@ -16,10 +16,11 @@ "metadata": {}, "outputs": [], "source": [ - "import pyart\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import pyart\n", "from pyart.testing import get_test_data" ] }, diff --git a/examples/retrieve/plot_hydrometeor.py b/examples/retrieve/plot_hydrometeor.py index 43702ca962b..f3f062c0391 100644 --- a/examples/retrieve/plot_hydrometeor.py +++ b/examples/retrieve/plot_hydrometeor.py @@ -53,7 +53,7 @@ kdp_field="specific_differential_phase", rhv_field="uncorrected_cross_correlation_ratio", temp_field="temperature", -) +)["hydro"] radar.add_field("radar_echo_classification", hydro) diff --git a/examples/retrieve/plot_hydrometeor_class_x_band.py b/examples/retrieve/plot_hydrometeor_class_x_band.py index 1841562fd89..2e8d40b7692 100644 --- a/examples/retrieve/plot_hydrometeor_class_x_band.py +++ b/examples/retrieve/plot_hydrometeor_class_x_band.py @@ -51,7 +51,7 @@ kdp_field="filtered_corrected_specific_diff_phase", rhv_field="RHOHV", temp_field="sounding_temperature", -) +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -69,7 +69,7 @@ rhv_field="RHOHV", temp_field="sounding_temperature", radar_freq=9.2e9, -) +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -98,7 +98,7 @@ rhv_field="RHOHV", temp_field="sounding_temperature", radar_freq=9.2e9, -) +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -121,7 +121,6 @@ "Lime", "Yellow", "Red", - "Fuchsia", ] cmaphid = colors.ListedColormap(hid_colors) cmapmeth = colors.ListedColormap(hid_colors[0:6]) @@ -129,19 +128,18 @@ def adjust_fhc_colorbar_for_pyart(cb): - cb.set_ticks(np.arange(1.4, 10, 0.9)) + cb.set_ticks(np.arange(1.4, 9, 0.9)) cb.ax.set_yticklabels( [ - "Drizzle", - "Rain", - "Ice Crystals", "Aggregates", - "Wet Snow", + "Ice Crystals", + "Light rain", + "Rain", "Vertical Ice", - "LD Graupel", + "Wet snow", "HD Graupel", - "Hail", - "Big Drops", + "Melting hail", + "Dry hail or high-density graupel", ] ) cb.ax.set_ylabel("") diff --git a/examples/retrieve/wavelet_echo_class_example.ipynb b/examples/retrieve/wavelet_echo_class_example.ipynb index 1b92486415c..64c1870c0d8 100644 --- a/examples/retrieve/wavelet_echo_class_example.ipynb +++ b/examples/retrieve/wavelet_echo_class_example.ipynb @@ -74,8 +74,9 @@ } ], "source": [ - "import pyart\n", - "import numpy as np" + "import numpy as np\n", + "\n", + "import pyart" ] }, { @@ -281,10 +282,9 @@ "source": [ "# Required imports\n", "\n", - "import matplotlib.pyplot as plt\n", - "import matplotlib.colors as mcolors\n", "import cartopy.crs as ccrs\n", - "\n", + "import matplotlib.colors as mcolors\n", + "import matplotlib.pyplot as plt\n", "\n", "display = pyart.graph.GridMapDisplay(grid)\n", "\n", diff --git a/pyart/__check_build/__init__.py b/pyart/__check_build/__init__.py index 7c050370481..610b381934c 100644 --- a/pyart/__check_build/__init__.py +++ b/pyart/__check_build/__init__.py @@ -1,6 +1,7 @@ """ Module to give helpful messages to the user that did not compile Py-ART properly. """ + import os INPLACE_MSG = """ diff --git a/pyart/core/__init__.py b/pyart/core/__init__.py index 15788cac001..abcf1e0c930 100644 --- a/pyart/core/__init__.py +++ b/pyart/core/__init__.py @@ -13,6 +13,7 @@ from .transforms import cartesian_vectors_to_geographic # noqa from .transforms import geographic_to_cartesian # noqa from .transforms import geographic_to_cartesian_aeqd # noqa +from .transforms import cartesian_to_antenna # noqa from .wind_profile import HorizontalWindProfile # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/pyart/core/transforms.py b/pyart/core/transforms.py index e9581f5c0eb..4dff4b16e08 100644 --- a/pyart/core/transforms.py +++ b/pyart/core/transforms.py @@ -78,6 +78,33 @@ def antenna_to_cartesian(ranges, azimuths, elevations): return x, y, z +def cartesian_to_antenna(x, y, z): + """ + Returns antenna coordinates from Cartesian coordinates. + + Parameters + ---------- + x, y, z : array + Cartesian coordinates in meters from the radar. + + Returns + ------- + ranges : array + Distances to the center of the radar gates (bins) in m. + azimuths : array + Azimuth angle of the radar in degrees. [-180., 180] + elevations : array + Elevation angle of the radar in degrees. + + """ + ranges = np.sqrt(x**2.0 + y**2.0 + z**2.0) + elevations = np.rad2deg(np.arctan(z / np.sqrt(x**2.0 + y**2.0))) + azimuths = np.rad2deg(np.arctan2(x, y)) # [-180, 180] + azimuths[azimuths < 0.0] += 360.0 # [0, 360] + + return ranges, azimuths, elevations + + def antenna_vectors_to_cartesian(ranges, azimuths, elevations, edges=False): """ Calculate Cartesian coordinate for gates from antenna coordinate vectors. diff --git a/pyart/core/wind_profile.py b/pyart/core/wind_profile.py index 2a8520a79a0..2ebfad767cb 100644 --- a/pyart/core/wind_profile.py +++ b/pyart/core/wind_profile.py @@ -3,7 +3,6 @@ """ - import numpy as np diff --git a/pyart/default_config.py b/pyart/default_config.py index 876d654f475..cba8690fc82 100644 --- a/pyart/default_config.py +++ b/pyart/default_config.py @@ -91,6 +91,16 @@ rain_rate = "rain_rate" radar_estimated_rain_rate = "radar_estimated_rain_rate" radar_echo_classification = "radar_echo_classification" +hydroclass_entropy = "hydroclass_entropy" +proportion_AG = "proportion_AG" +proportion_CR = "proportion_CR" +proportion_LR = "proportion_LR" +proportion_RP = "proportion_RP" +proportion_RN = "proportion_RN" +proportion_VI = "proportion_VI" +proportion_WS = "proportion_WS" +proportion_MH = "proportion_MH" +proportion_IH = "proportion_IH" specific_attenuation = "specific_attenuation" specific_differential_attenuation = "specific_differential_attenuation" clutter_filter_power_removed = "clutter_filter_power_removed" @@ -186,6 +196,16 @@ "rain_rate": rain_rate, "radar_estimated_rain_rate": radar_estimated_rain_rate, "radar_echo_classification": radar_echo_classification, + "hydroclass_entropy": hydroclass_entropy, + "proportion_AG": proportion_AG, + "proportion_CR": proportion_CR, + "proportion_LR": proportion_LR, + "proportion_RP": proportion_RP, + "proportion_RN": proportion_RN, + "proportion_VI": proportion_VI, + "proportion_WS": proportion_WS, + "proportion_MH": proportion_MH, + "proportion_IH": proportion_IH, "specific_attenuation": specific_attenuation, "differential_phase_texture": differential_phase_texture, "eastward_wind_component": eastward_wind_component, @@ -640,6 +660,66 @@ "long_name": "Radar Echo classification", "coordinates": "elevation azimuth range", }, + hydroclass_entropy: { + "units": "-", + "standard_name": "hydroclass_entropy", + "long_name": "Semi-supervised hydrometeor classification entropy", + "coordinates": "elevation azimuth range", + }, + proportion_AG: { + "units": "percent", + "standard_name": "proportion_AG", + "long_name": "Aggregates proportion", + "coordinates": "elevation azimuth range", + }, + proportion_CR: { + "units": "percent", + "standard_name": "proportion_CR", + "long_name": "Crystals proportion", + "coordinates": "elevation azimuth range", + }, + proportion_LR: { + "units": "percent", + "standard_name": "proportion_LR", + "long_name": "Light Rain proportion", + "coordinates": "elevation azimuth range", + }, + proportion_RP: { + "units": "percent", + "standard_name": "proportion_RP", + "long_name": "Rimed particles proportion", + "coordinates": "elevation azimuth range", + }, + proportion_RN: { + "units": "percent", + "standard_name": "proportion_RN", + "long_name": "Rain proportion", + "coordinates": "elevation azimuth range", + }, + proportion_VI: { + "units": "percent", + "standard_name": "proportion_VI", + "long_name": "Vertical Ice Crystals proportion", + "coordinates": "elevation azimuth range", + }, + proportion_WS: { + "units": "percent", + "standard_name": "proportion_WS", + "long_name": "Wet Snow proportion", + "coordinates": "elevation azimuth range", + }, + proportion_MH: { + "units": "percent", + "standard_name": "proportion_MH", + "long_name": "Melting Hail proportion", + "coordinates": "elevation azimuth range", + }, + proportion_IH: { + "units": "percent", + "standard_name": "proportion_IH", + "long_name": "Ice Hail proportion", + "coordinates": "elevation azimuth range", + }, specific_attenuation: { "units": "dB/km", "standard_name": "specific_attenuation", @@ -1442,6 +1522,16 @@ def spectrum_width_limit(container=None, selection=0): rain_rate: "pyart_RRate11", radar_estimated_rain_rate: "pyart_RRate11", radar_echo_classification: "pyart_LangRainbow12", + hydroclass_entropy: "pyart_LangRainbow12", + proportion_AG: "pyart_LangRainbow12", + proportion_CR: "pyart_LangRainbow12", + proportion_LR: "pyart_LangRainbow12", + proportion_RP: "pyart_LangRainbow12", + proportion_RN: "pyart_LangRainbow12", + proportion_VI: "pyart_LangRainbow12", + proportion_WS: "pyart_LangRainbow12", + proportion_MH: "pyart_LangRainbow12", + proportion_IH: "pyart_LangRainbow12", specific_attenuation: "pyart_Carbone17", differential_phase_texture: "pyart_BlueBrown11", height: "pyart_SCook18", @@ -1492,6 +1582,16 @@ def spectrum_width_limit(container=None, selection=0): rain_rate: (0.0, 50.0), radar_estimated_rain_rate: (0.0, 50.0), radar_echo_classification: (0, 11), + hydroclass_entropy: (0.0, 1.0), + proportion_AG: (0.0, 100.0), + proportion_CR: (0.0, 100.0), + proportion_LR: (0.0, 100.0), + proportion_RP: (0.0, 100.0), + proportion_RN: (0.0, 100.0), + proportion_VI: (0.0, 100.0), + proportion_WS: (0.0, 100.0), + proportion_MH: (0.0, 100.0), + proportion_IH: (0.0, 100.0), specific_attenuation: (0.0, 10.0), differential_phase_texture: (0, 180.0), height: (0, 20000), diff --git a/pyart/graph/_cm_colorblind.py b/pyart/graph/_cm_colorblind.py index 8b078ea58e9..164fcd02dec 100644 --- a/pyart/graph/_cm_colorblind.py +++ b/pyart/graph/_cm_colorblind.py @@ -2,6 +2,7 @@ Data for colorblind friendly radar colormaps """ + import os import numpy as np diff --git a/pyart/graph/gridmapdisplay_basemap.py b/pyart/graph/gridmapdisplay_basemap.py index 1b6f914aaee..997de186ea3 100644 --- a/pyart/graph/gridmapdisplay_basemap.py +++ b/pyart/graph/gridmapdisplay_basemap.py @@ -3,7 +3,6 @@ """ - import matplotlib.pyplot as plt import numpy as np diff --git a/pyart/io/mdv_grid.py b/pyart/io/mdv_grid.py index ebcedf3bb87..e8619665001 100644 --- a/pyart/io/mdv_grid.py +++ b/pyart/io/mdv_grid.py @@ -3,7 +3,6 @@ """ - import datetime import warnings diff --git a/pyart/io/nexrad_common.py b/pyart/io/nexrad_common.py index e34c27751e2..13f3c6acd8a 100644 --- a/pyart/io/nexrad_common.py +++ b/pyart/io/nexrad_common.py @@ -2,6 +2,7 @@ Data and functions common to all types of NEXRAD files. """ + # The functions in this module are intended to be used in other # nexrad related modules. The functions are not and should not be # exported into the pyart.io namespace. diff --git a/pyart/io/nexrad_level2.py b/pyart/io/nexrad_level2.py index 38a3e88530b..0facc446dbd 100644 --- a/pyart/io/nexrad_level2.py +++ b/pyart/io/nexrad_level2.py @@ -2,7 +2,6 @@ Functions for reading NEXRAD level 2 files. """ - # This file is part of the Py-ART, the Python ARM Radar Toolkit # https://github.com/ARM-DOE/pyart diff --git a/pyart/io/nexrad_level3.py b/pyart/io/nexrad_level3.py index 0b2ed137f47..72605fb3a6e 100644 --- a/pyart/io/nexrad_level3.py +++ b/pyart/io/nexrad_level3.py @@ -761,7 +761,7 @@ def _int16_to_float16(val): ("block_length", INT4), # Length of block in bytes ("layers", INT2), # Number of data layers ("layer_divider", INT2), # Delineate data layers, -1 - ("layer_length", INT4) # Length of data layer in bytes + ("layer_length", INT4), # Length of data layer in bytes # Display data packets ) diff --git a/pyart/map/gates_to_grid.py b/pyart/map/gates_to_grid.py index bcc470d55a9..95a6e888248 100644 --- a/pyart/map/gates_to_grid.py +++ b/pyart/map/gates_to_grid.py @@ -2,6 +2,7 @@ Generate a Cartesian grid by mapping from radar gates onto the grid. """ + import gc import warnings diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py index 9be0b022e98..563ced9a2cb 100644 --- a/pyart/retrieve/__init__.py +++ b/pyart/retrieve/__init__.py @@ -21,7 +21,8 @@ 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 .qvp import quasi_vertical_profile, compute_qvp, compute_rqvp # noqa +from .qvp import compute_evp, compute_svp, compute_vp, compute_ts_along_coord # noqa from .simple_moment_calculations import calculate_snr_from_reflectivity # noqa from .simple_moment_calculations import calculate_velocity_texture # noqa from .simple_moment_calculations import compute_cdr # noqa diff --git a/pyart/retrieve/_echo_class_wt.py b/pyart/retrieve/_echo_class_wt.py index e9d38365e34..1e1c30615f7 100644 --- a/pyart/retrieve/_echo_class_wt.py +++ b/pyart/retrieve/_echo_class_wt.py @@ -13,7 +13,6 @@ atwt2d """ - import numpy as np diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index 8057704ac22..327acb9c00b 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -3,6 +3,7 @@ """ +from copy import deepcopy from warnings import warn import numpy as np @@ -10,6 +11,7 @@ # Local imports from ..config import get_field_name, get_fillvalue, get_metadata from ..core import Grid +from ..util.radar_utils import ma_broadcast_to from ._echo_class import _feature_detection, steiner_class_buff from ._echo_class_wt import calc_scale_break, wavelet_reclass @@ -601,15 +603,25 @@ def feature_detection( def hydroclass_semisupervised( radar, + hydro_names=("AG", "CR", "LR", "RP", "RN", "VI", "WS", "MH", "IH/HDG"), + var_names=("Zh", "ZDR", "KDP", "RhoHV", "relH"), mass_centers=None, weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), + value=50.0, + lapse_rate=-6.5, refl_field=None, zdr_field=None, rhv_field=None, kdp_field=None, temp_field=None, + iso0_field=None, hydro_field=None, + entropy_field=None, radar_freq=None, + temp_ref="temperature", + compute_entropy=False, + output_distances=False, + vectorize=False, ): """ Classifies precipitation echoes into hydrometeor types. @@ -621,28 +633,57 @@ def hydroclass_semisupervised( ---------- radar : radar Radar object. + hydro_names : array of str + name of the types of hydrometeors + var_names : array of str + name of the variables mass_centers : ndarray 2D, optional The centroids for each variable and hydrometeor class in (nclasses, nvariables). weights : ndarray 1D, optional - The weight given to each variable. - refl_field, zdr_field, rhv_field, kdp_field, temp_field : str, optional + The weight given to each variable. Ordered by [dBZ, ZDR, KDP, RhoHV, + H_ISO0] + value : float + The value controlling the rate of decay in the distance transformation + lapse_rate : float + The decrease in temperature for each vertical km [deg/km] + refl_field, zdr_field, rhv_field, kdp_field, temp_field, iso0_field : str Inputs. Field names within the radar object which represent the horizonal reflectivity, the differential reflectivity, the copolar - correlation coefficient, the specific differential phase and the - temperature field. A value of None for any of these parameters will - use the default field name as defined in the Py-ART configuration - file. - hydro_field : str, optional + correlation coefficient, the specific differential phase, the + temperature (in deg celsius) and the height respect to the iso0 fields. + A value of None for any of these parameters will use the default field + name as defined in the Py-ART configuration file. + hydro_field : str Output. Field name which represents the hydrometeor class field. A value of None will use the default field name as defined in the Py-ART configuration file. + entropy_field : str + Output. Field name which represents the entropy class field. + A value of None will use the default field name as defined in the + Py-ART configuration file. radar_freq : str, optional Radar frequency in Hertz (Hz) used for classification. This parameter will be ignored, if the radar object has frequency information. + temp_ref : str + the field use as reference for temperature. Can be either temperature + or height_over_iso0 + compute_entropy : bool + If true, the entropy is computed + output_distances : bool + If true, the normalized distances to the centroids for each + hydrometeor are provided as output + vectorize : bool + If true, a vectorized version of the class assignation is going to be + used Returns ------- + fields_dict : dict + Dictionary containing the retrieved fields + + The output directionary field_dict has the following keys: + hydro : dict Hydrometeor classification. - 0: Not classified @@ -656,6 +697,15 @@ def hydroclass_semisupervised( - 8: Melting hail - 9: Dry hail or high-density graupel + if compute_entropy is True: + entropy : dict + Shannon entropy of the hydrometeor demixing + + if output_distances is True: + propX : dict + Proportion of a given hydrometeor class in the polarimetric + decomposition of a radar volume + References ---------- Besic, N., Figueras i Ventura, J., Grazioli, J., Gabella, M., Germann, U., @@ -663,6 +713,11 @@ def hydroclass_semisupervised( of polarimetric radar measurements: a semi-supervised approach, Atmos. Meas. Tech., 9, 4425-4445, doi:10.5194/amt-9-4425-2016, 2016 + Besic, N., Gehring, J., Praz, C., Figueras i Ventura, J., Grazioli, J., Gabella, + M., Germann, U., and Berne, A.: Unraveling hydrometeor mixtures in polarimetric + radar measurements, Atmos. Meas. Tech., 11, 4847–4866, + https://doi.org/10.5194/amt-11-4847-2018, 2018. + Usage ----- .. code-block:: python @@ -692,7 +747,6 @@ def hydroclass_semisupervised( "data": [9.2e9] } """ - lapse_rate = -6.5 # select the centroids as a function of frequency band if mass_centers is None: @@ -713,61 +767,108 @@ def hydroclass_semisupervised( "So frequency is unknown. Default coefficients for C band will be applied." ) - # parse the field parameters - if refl_field is None: - refl_field = get_field_name("reflectivity") - if zdr_field is None: - zdr_field = get_field_name("differential_reflectivity") - if rhv_field is None: - rhv_field = get_field_name("cross_correlation_ratio") - if kdp_field is None: - kdp_field = get_field_name("specific_differential_phase") - if temp_field is None: - temp_field = get_field_name("temperature") if hydro_field is None: hydro_field = get_field_name("radar_echo_classification") + if compute_entropy: + if entropy_field is None: + entropy_field = get_field_name("hydroclass_entropy") + + # Get the data fields + fields_dict = {} + for var_name in var_names: + if var_name == "Zh": + if refl_field is None: + refl_field = get_field_name("reflectivity") + radar.check_field_exists(refl_field) + fields_dict.update({var_name: deepcopy(radar.fields[refl_field]["data"])}) + elif var_name == "ZDR": + if zdr_field is None: + zdr_field = get_field_name("differential_reflectivity") + radar.check_field_exists(zdr_field) + fields_dict.update({var_name: deepcopy(radar.fields[zdr_field]["data"])}) + elif var_name == "KDP": + if kdp_field is None: + kdp_field = get_field_name("specific_differential_phase") + radar.check_field_exists(kdp_field) + fields_dict.update({var_name: deepcopy(radar.fields[kdp_field]["data"])}) + elif var_name == "RhoHV": + if rhv_field is None: + rhv_field = get_field_name("cross_correlation_ratio") + radar.check_field_exists(rhv_field) + fields_dict.update({var_name: deepcopy(radar.fields[rhv_field]["data"])}) + elif var_name == "relH": + if temp_ref == "temperature": + if temp_field is None: + temp_field = get_field_name("temperature") + radar.check_field_exists(temp_field) + # convert temp in relative height respect to iso0 + temp = deepcopy(radar.fields[temp_field]["data"]) + fields_dict.update({var_name: temp * (1000.0 / lapse_rate)}) + else: + if iso0_field is None: + iso0_field = get_field_name("height_over_iso0") + radar.check_field_exists(iso0_field) + fields_dict.update( + {var_name: deepcopy(radar.fields[iso0_field]["data"])} + ) + else: + raise ValueError( + "Variable " + + var_name + + " unknown. " + + "Valid variable names for hydrometeor classification are: " + + "relH, Zh, ZDR, KDP and RhoHV" + ) + + # standardize data and centroids + mc_std = np.empty(np.shape(mass_centers), dtype=fields_dict[var_names[0]].dtype) + for i, var_name in enumerate(var_names): + mc_std[:, i] = _standardize(mass_centers[:, i], var_name) + fields_dict[var_name] = _standardize(fields_dict[var_name], var_name) - # extract fields and parameters from radar - radar.check_field_exists(refl_field) - radar.check_field_exists(zdr_field) - radar.check_field_exists(rhv_field) - radar.check_field_exists(kdp_field) - radar.check_field_exists(temp_field) - - refl = radar.fields[refl_field]["data"] - zdr = radar.fields[zdr_field]["data"] - rhohv = radar.fields[rhv_field]["data"] - kdp = radar.fields[kdp_field]["data"] - temp = radar.fields[temp_field]["data"] - - # convert temp in relative height respect to iso0 - relh = temp * (1000.0 / lapse_rate) - - # standardize data - refl_std = _standardize(refl, "Zh") - zdr_std = _standardize(zdr, "ZDR") - kdp_std = _standardize(kdp, "KDP") - rhohv_std = _standardize(rhohv, "RhoHV") - relh_std = _standardize(relh, "relH") - - # standardize centroids - mc_std = np.zeros(np.shape(mass_centers)) - mc_std[:, 0] = _standardize(mass_centers[:, 0], "Zh") - mc_std[:, 1] = _standardize(mass_centers[:, 1], "ZDR") - mc_std[:, 2] = _standardize(mass_centers[:, 2], "KDP") - mc_std[:, 3] = _standardize(mass_centers[:, 3], "RhoHV") - mc_std[:, 4] = _standardize(mass_centers[:, 4], "relH") + # if entropy has to be computed get transformation parameters + t_vals = None + if compute_entropy: + t_vals = _compute_coeff_transform(mc_std, weights=weights, value=value) # assign to class - hydroclass_data, min_dist = _assign_to_class( - refl_std, zdr_std, kdp_std, rhohv_std, relh_std, mc_std, weights=weights - ) + if vectorize: + hydroclass_data, entropy_data, prop_data = _assign_to_class_scan( + fields_dict, mc_std, var_names=var_names, weights=weights, t_vals=t_vals + ) + else: + hydroclass_data, entropy_data, prop_data = _assign_to_class( + fields_dict, mc_std, weights=weights, t_vals=t_vals + ) # prepare output fields + fields_dict = dict() hydro = get_metadata(hydro_field) hydro["data"] = hydroclass_data - - return hydro + hydro.update({"_FillValue": 0}) + labels = ["NC"] + ticks = [1] + boundaries = [-0.5, 1.5] + for i, hydro_name in enumerate(hydro_names): + labels.append(hydro_name) + ticks.append(i + 1) + boundaries.append(i + 1.5) + hydro.update({"labels": labels, "ticks": ticks, "boundaries": boundaries}) + fields_dict.update({"hydro": hydro}) + + if compute_entropy: + entropy = get_metadata(entropy_field) + entropy["data"] = entropy_data + fields_dict.update({"entropy": entropy}) + + if output_distances: + for field_name in hydro_names: + field_name = "proportion_" + field_name + prop = get_metadata(field_name) + prop["data"] = prop_data[:, :, i] + fields_dict.update({field_name: prop}) + + return fields_dict def _standardize(data, field_name, mx=None, mn=None): @@ -809,7 +910,9 @@ def _standardize(data, field_name, mx=None, mn=None): data[data < -0.5] = -0.5 data = 10.0 * np.ma.log10(data + 0.6) elif field_name == "RhoHV": - data = 10.0 * np.ma.log10(1.0 - data) + # avoid infinite result + data[data > 1.0] = 1.0 + data = 10.0 * np.ma.log10(1.0000000000001 - data) mask = np.ma.getmaskarray(data) field_std = 2.0 * (data - mn) / (mx - mn) - 1.0 @@ -821,13 +924,11 @@ def _standardize(data, field_name, mx=None, mn=None): def _assign_to_class( - zh, - zdr, - kdp, - rhohv, - relh, + fields_dict, mass_centers, + var_names=("Zh", "ZDR", "KDP", "RhoHV", "relH"), weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), + t_vals=None, ): """ Assigns an hydrometeor class to a radar range bin computing @@ -835,57 +936,201 @@ def _assign_to_class( Parameters ---------- - zh, zdr, kdp, rhohv, relh : radar fields - Variables used for assignment normalized to [-1, 1] values. + fields_dict : dict + Dictionary containg the variables used for assigment normalized to + [-1, 1] values + mass_centers : matrix + centroids normalized to [-1, 1] values (nclasses, nvariables) + var_names : array of str + Name of the variables + weights : array + optional. The weight given to each variable (nvariables) + t_vals : array + transformation values for the distance to centroids (nclasses) + + Returns + ------- + hydroclass : int array + the index corresponding to the assigned class + entropy : float array + the entropy + t_dist : float matrix + if entropy is computed, the transformed distances of each class + (proxy for proportions of each hydrometeor) (nrays, nbins, nclasses) + + """ + # prepare data + nrays = fields_dict[var_names[0]].shape[0] + nbins = fields_dict[var_names[0]].shape[1] + nclasses = mass_centers.shape[0] + nvariables = mass_centers.shape[1] + dtype = fields_dict[var_names[0]].dtype + + hydroclass = np.ma.empty((nrays, nbins), dtype=np.uint8) + entropy = None + t_dist = None + if t_vals is not None: + entropy = np.ma.empty((nrays, nbins), dtype=dtype) + t_dist = np.ma.masked_all((nrays, nbins, nclasses), dtype=dtype) + + for ray in range(nrays): + data = [] + for var_name in var_names: + data.append(fields_dict[var_name][ray, :]) + data = np.ma.array(data, dtype=dtype) + weights_mat = np.broadcast_to( + weights.reshape(nvariables, 1), (nvariables, nbins) + ) + dist = np.ma.zeros((nclasses, nbins), dtype=dtype) + + # compute distance: masked entries will not contribute to the distance + mask = np.ma.getmaskarray(fields_dict[var_names[0]][ray, :]) + for i in range(nclasses): + centroids_class = mass_centers[i, :] + centroids_class = np.broadcast_to( + centroids_class.reshape(nvariables, 1), (nvariables, nbins) + ) + dist_ray = np.ma.sqrt( + np.ma.sum(((centroids_class - data) ** 2.0) * weights_mat, axis=0) + ) + dist_ray[mask] = np.ma.masked + dist[i, :] = dist_ray + + # Get hydrometeor class + class_vec = dist.argsort(axis=0, fill_value=10e40) + hydroclass_ray = (class_vec[0, :] + 1).astype(np.uint8) + hydroclass_ray[mask] = 0 + hydroclass[ray, :] = hydroclass_ray + + if t_vals is None: + continue + + # Transform the distance using the coefficient of the dominant class + t_vals_ray = np.ma.masked_where(mask, t_vals[class_vec[0, :]]) + t_vals_ray = ma_broadcast_to(t_vals_ray.reshape(1, nbins), (nclasses, nbins)) + t_dist_ray = np.ma.exp(-t_vals_ray * dist) + + # set transformed distances to a value between 0 and 1 + dist_total = np.ma.sum(t_dist_ray, axis=0) + dist_total = ma_broadcast_to(dist_total.reshape(1, nbins), (nclasses, nbins)) + t_dist_ray /= dist_total + + # Compute entropy + entropy_ray = -np.ma.sum( + t_dist_ray * np.ma.log(t_dist_ray) / np.ma.log(nclasses), axis=0 + ) + entropy_ray[mask] = np.ma.masked + entropy[ray, :] = entropy_ray + + t_dist[ray, :, :] = np.ma.transpose(t_dist_ray) + + if t_vals is not None: + t_dist *= 100.0 + + return hydroclass, entropy, t_dist + + +def _assign_to_class_scan( + fields_dict, + mass_centers, + var_names=("Zh", "ZDR", "KDP", "RhoHV", "relH"), + weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), + t_vals=None, +): + """ + assigns an hydrometeor class to a radar range bin computing + the distance between the radar variables an a centroid. + Computes the entire radar volume at once + + Parameters + ---------- + fields_dict : dict + Dictionary containg the variables used for assigment normalized to + [-1, 1] values mass_centers : matrix - Centroids normalized to [-1, 1] values. - weights : array, optional - The weight given to each variable. + centroids normalized to [-1, 1] values + var_names : array of str + Name of the variables + weights : array + optional. The weight given to each variable + t_vals : matrix + transformation values for the distance to centroids + (nclasses, nvariables) Returns ------- hydroclass : int array - The index corresponding to the assigned class. - mind_dist : float array - The minimum distance to the centroids. + the index corresponding to the assigned class + entropy : float array + the entropy + t_dist : float matrix + if entropy is computed, the transformed distances of each class + (proxy for proportions of each hydrometeor) (nrays, nbins, nclasses) """ # prepare data - nrays = zh.shape[0] - nbins = zdr.shape[1] + nrays = fields_dict[var_names[0]].shape[0] + nbins = fields_dict[var_names[0]].shape[1] nclasses = mass_centers.shape[0] nvariables = mass_centers.shape[1] + dtype = fields_dict[var_names[0]].dtype - data = np.ma.array([zh, zdr, kdp, rhohv, relh]) + data = [] + for var_name in var_names: + data.append(fields_dict[var_name]) + data = np.ma.array(data, dtype=dtype) weights_mat = np.broadcast_to( weights.reshape(nvariables, 1, 1), (nvariables, nrays, nbins) ) - dist = np.ma.zeros((nclasses, nrays, nbins), dtype="float64") # compute distance: masked entries will not contribute to the distance + mask = np.ma.getmaskarray(fields_dict[var_names[0]]) + dist = np.ma.zeros((nrays, nbins, nclasses), dtype=dtype) + t_dist = None + entropy = None for i in range(nclasses): centroids_class = mass_centers[i, :] centroids_class = np.broadcast_to( centroids_class.reshape(nvariables, 1, 1), (nvariables, nrays, nbins) ) - dist[i, :, :] = np.ma.sqrt( + dist_aux = np.ma.sqrt( np.ma.sum(((centroids_class - data) ** 2.0) * weights_mat, axis=0) ) + dist_aux[mask] = np.ma.masked + dist[:, :, i] = dist_aux - # use very large fill_value so that masked entries will be sorted at the - # end. There should not be any masked entry anyway - class_vec = dist.argsort(axis=0, fill_value=10e40) - - # get minimum distance. Acts as a confidence value - dist.sort(axis=0, fill_value=10e40) - min_dist = dist[0, :, :] + del data + del weights_mat - # Entries with non-valid reflectivity values are set to 0 (No class) - mask = np.ma.getmaskarray(zh) - hydroclass = class_vec[0, :, :] + 1 + # Get hydrometeor class + class_vec = dist.argsort(axis=-1, fill_value=10e40) + hydroclass = np.ma.asarray(class_vec[:, :, 0] + 1, dtype=np.uint8) hydroclass[mask] = 0 - return hydroclass, min_dist + if t_vals is not None: + # Transform the distance using the coefficient of the dominant class + t_vals_aux = np.ma.masked_where(mask, t_vals[class_vec[:, :, 0]]) + t_vals_aux = ma_broadcast_to( + t_vals_aux.reshape(nrays, nbins, 1), (nrays, nbins, nclasses) + ) + t_dist = np.ma.exp(-t_vals_aux * dist) + del t_vals_aux + + # set distance to a value between 0 and 1 + dist_total = np.ma.sum(t_dist, axis=-1) + dist_total = ma_broadcast_to( + dist_total.reshape(nrays, nbins, 1), (nrays, nbins, nclasses) + ) + t_dist /= dist_total + del dist_total + + # compute entroy + entropy = -np.ma.sum(t_dist * np.ma.log(t_dist) / np.ma.log(nclasses), axis=-1) + entropy[mask] = np.ma.masked + + t_dist *= 100.0 + + return hydroclass, entropy, t_dist def _get_mass_centers(freq): @@ -985,9 +1230,10 @@ def _data_limits_table(): """ dlimits_dict = dict() dlimits_dict.update({"Zh": (60.0, -10.0)}) - dlimits_dict.update({"ZDR": (5.0, -5.0)}) + dlimits_dict.update({"ZDR": (5.0, -1.5)}) dlimits_dict.update({"KDP": (7.0, -10.0)}) dlimits_dict.update({"RhoHV": (-5.23, -50.0)}) + dlimits_dict.update({"RelH": (5000.0, -5000.0)}) return dlimits_dict @@ -1189,3 +1435,49 @@ def conv_strat_raut( } return reclass_dict + + +def _compute_coeff_transform( + mass_centers, weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), value=50.0 +): + """ + get the transformation coefficients + + Parameters + ---------- + mass_centers : ndarray 2D + The centroids for each class and variable (nclasses, nvariables) + weights : array + optional. The weight given to each variable (nvariables) + value : float + parameter controlling the rate of decay of the distance transformation + + Returns + ------- + t_vals : ndarray 1D + The coefficients used to transform the distances to each centroid for + each class (nclasses) + + """ + nclasses, nvariables = np.shape(mass_centers) + t_vals = np.empty((nclasses, nclasses), dtype=mass_centers.dtype) + for i in range(nclasses): + weights_mat = np.broadcast_to( + weights.reshape(1, nvariables), (nclasses, nvariables) + ) + centroids_class = mass_centers[i, :] + centroids_class = np.broadcast_to( + centroids_class.reshape(1, nvariables), (nclasses, nvariables) + ) + t_vals[i, :] = np.sqrt( + np.sum( + weights_mat * np.power(np.abs(centroids_class - mass_centers), 2.0), + axis=1, + ) + ) + + # pick the second lowest value (the first is 0) + t_vals = np.sort(t_vals, axis=-1)[:, 1] + t_vals = np.log(value) / t_vals + + return t_vals diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 991b43d5930..b4a6b4a63ee 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -1,16 +1,47 @@ """ +pyart.retrieve.quasi_vertical_profile +===================================== + Retrieval of QVPs from a radar object -""" +.. autosummary:: + :toctree: generated/ + + quasi_vertical_profile + compute_qvp + compute_rqvp + compute_evp + compute_svp + compute_vp + compute_ts_along_coord + project_to_vertical + find_rng_index + get_target_elevations + find_nearest_gate + find_neighbour_gates + _create_qvp_object + _create_along_coord_object + _update_qvp_metadata + _update_along_coord_metadata + +""" "" "" "" "" "" "" "" "" + +from copy import deepcopy +from warnings import warn import numpy as np +from netCDF4 import num2date +from scipy.interpolate import interp1d from ..core.transforms import antenna_to_cartesian +from ..io.common import make_time_unit_str +from ..util.circular_stats import compute_directional_stats +from ..util.datetime_utils import datetime_from_radar +from ..util.radar_utils import ma_broadcast_to +from ..util.xsect import cross_section_ppi, cross_section_rhi -def quasi_vertical_profile( - radar, desired_angle=None, fields=None, gatefilter=None, verbose=False -): +def quasi_vertical_profile(radar, desired_angle=None, fields=None, gatefilter=None): """ Quasi Vertical Profile. @@ -26,8 +57,6 @@ def quasi_vertical_profile( desired_angle : float Radar tilt angle to use for indexing radar field data. None will result in wanted_angle = 20.0 - verbose : bool - boolean flag to turn the printing of radar tilts on or off. Other Parameters ---------------- @@ -38,38 +67,38 @@ def quasi_vertical_profile( Returns ------- qvp : Dictonary - A quasai vertical profile object containing fields + A quasi vertical profile object containing fields from a radar object References ---------- Troemel, S., M. Kumjian, A. Ryzhkov, and C. Simmer, 2013: Backscatter - differential phase - estimation and variability. J Appl. Meteor. Clim.. + differential phase - estimation and variability. J Appl. Meteor. Clim. 52, 2529 - 2548. - Troemel, S., A. Ryzhkov, P. Zhang, and C. Simmer, 2014: Investigations - of backscatter differential phase in the melting layer. J. Appl. Meteorol. - Clim. 54, 2344 - 2359. + Troemel, S., A. Ryzhkov, P. Zhang, and C. Simmer, 2014: + Investigations of backscatter differential phase in the melting layer. J. + Appl. Meteorol. Clim. 54, 2344 - 2359. Ryzhkov, A., P. Zhang, H. Reeves, M. Kumjian, T. Tschallener, S. Tromel, - C. Simmer, 2015: Quasi-vertical profiles - a new way to look at polarimetric - radar data. Submitted to J. Atmos. Oceanic Technol. + C. Simmer, 2015: Quasi-vertical profiles - a new way to look at + polarimetric radar data. Submitted to J. Atmos. Oceanic Technol. """ # Creating an empty dictonary qvp = {} - # Setting the desired radar angle and getting index value for desired radar angle + # Setting the desired radar angle and getting index value for desired + # radar angle if desired_angle is None: desired_angle = 20.0 index = abs(radar.fixed_angle["data"] - desired_angle).argmin() radar_slice = radar.get_slice(index) - if verbose: - # Printing radar tilt angles and radar elevation - print(radar.fixed_angle["data"]) - print(radar.elevation["data"][-1]) + # Printing radar tilt angles and radar elevation + print(radar.fixed_angle["data"]) + print(radar.elevation["data"][-1]) # Setting field parameters # If fields is None then all radar fields pulled else defined field is used @@ -91,7 +120,7 @@ def quasi_vertical_profile( qvp.update({field: radar_fields}) else: - # Filtereing data based on defined gatefilter + # Filtering data based on defined gatefilter # If none is defined goes to else statement if gatefilter is not None: get_field = radar.get_field(index, fields) @@ -111,3 +140,1680 @@ def quasi_vertical_profile( ) qvp.update({"height": z}) return qvp + + +def compute_qvp( + radar, + field_names, + ref_time=None, + angle=0.0, + ang_tol=1.0, + hmax=10000.0, + hres=50.0, + avg_type="mean", + nvalid_min=30, + interp_kind="none", + qvp=None, +): + """ + Computes quasi vertical profiles. + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of str + list of field names to add to the QVP + ref_time : datetime object + reference time for current radar volume + angle : int or float + If the radar object contains a PPI volume, the sweep number to + use, if it contains an RHI volume the elevation angle. + ang_tol : float + If the radar object contains an RHI volume, the tolerance in the + elevation angle for the conversion into PPI + hmax : float + The maximum height to plot [m]. + hres : float + The height resolution [m]. + avg_type : str + The type of averaging to perform. Can be either "mean" or "median" + nvalid_min : int + Minimum number of valid points to accept average. + interp_kind : str + type of interpolation when projecting to vertical grid: 'none', + or 'nearest', etc. + 'none' will select from all data points within the regular grid + height bin the closest to the center of the bin. + 'nearest' will select the closest data point to the center of the + height bin regardless if it is within the height bin or not. + Data points can be masked values + If another type of interpolation is selected masked values will be + eliminated from the data points before the interpolation + qvp : QVP object or None + If it is not None this is the QVP object where to store the data from + the current time step. Otherwise a new QVP object will be created + + + Returns + ------- + qvp : qvp object + The computed QVP object + + Reference + --------- + Ryzhkov A., Zhang P., Reeves H., Kumjian M., Tschallener T., Trömel S., + Simmer C. 2016: Quasi-Vertical Profiles: A New Way to Look at Polarimetric + Radar Data. JTECH vol. 33 pp 551-562 + + """ + if avg_type not in ("mean", "median"): + warn("Unsuported statistics " + avg_type) + return None + + radar_aux = deepcopy(radar) + # transform radar into ppi over the required elevation + if radar_aux.scan_type == "rhi": + radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol) + elif radar_aux.scan_type == "ppi": + radar_aux = radar_aux.extract_sweeps([int(angle)]) + else: + warn("Error: unsupported scan type.") + return None + + if qvp is None: + qvp = _create_qvp_object( + radar_aux, + field_names, + qvp_type="qvp", + start_time=ref_time, + hmax=hmax, + hres=hres, + ) + + # modify metadata + if ref_time is None: + ref_time = datetime_from_radar(radar_aux) + qvp = _update_qvp_metadata( + qvp, + ref_time, + qvp.longitude["data"][0], + qvp.latitude["data"][0], + elev=qvp.fixed_angle["data"][0], + ) + + for field_name in field_names: + # compute QVP data + if field_name not in radar_aux.fields: + warn("Field " + field_name + " not in radar object") + qvp_data = np.ma.masked_all(qvp.ngates) + else: + values, _ = compute_directional_stats( + radar_aux.fields[field_name]["data"], + avg_type=avg_type, + nvalid_min=nvalid_min, + axis=0, + ) + + # Project to vertical grid: + qvp_data = project_to_vertical( + values, + radar_aux.gate_altitude["data"][0, :], + qvp.range["data"], + interp_kind=interp_kind, + ) + + # Put data in radar object + if np.size(qvp.fields[field_name]["data"]) == 0: + qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates) + else: + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) + + return qvp + + +def compute_rqvp( + radar, + field_names, + ref_time=None, + hmax=10000.0, + hres=2.0, + avg_type="mean", + nvalid_min=30, + interp_kind="nearest", + rmax=50000.0, + weight_power=2.0, + qvp=None, +): + """ + Computes range-defined quasi vertical profiles. + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of str + list of field names to add to the QVP + ref_time : datetime object + reference time for current radar volume + hmax : float + The maximum height to plot [m]. + hres : float + The height resolution [m]. + avg_type : str + The type of averaging to perform. Can be either "mean" or "median" + nvalid_min : int + Minimum number of valid points to accept average. + interp_kind : str + type of interpolation when projecting to vertical grid: 'none', + or 'nearest', etc. + 'none' will select from all data points within the regular grid + height bin the closest to the center of the bin. + 'nearest' will select the closest data point to the center of the + height bin regardless if it is within the height bin or not. + Data points can be masked values + If another type of interpolation is selected masked values will be + eliminated from the data points before the interpolation + rmax : float + ground range up to which the data is intended for use [m]. + weight_power : float + Power p of the weighting function 1/abs(grng-(rmax-1))**p given to + the data outside the desired range. -1 will set the weight to 0. + qvp : QVP object or None + If it is not None this is the QVP object where to store the data from + the current time step. Otherwise a new QVP object will be created + + + Returns + ------- + qvp : qvp object + The computed range defined quasi vertical profile + + Reference + --------- + Tobin D.M., Kumjian M.R. 2017: Polarimetric Radar and Surface-Based + Precipitation-Type Observations of ice Pellet to Freezing Rain + Transitions. Weather and Forecasting vol. 32 pp 2065-2082 + + """ + if avg_type not in ("mean", "median"): + warn("Unsuported statistics " + avg_type) + return None + + radar_aux = deepcopy(radar) + # transform radar into ppi over the required elevation + if radar_aux.scan_type == "rhi": + target_elevations, el_tol = get_target_elevations(radar_aux) + radar_ppi = cross_section_rhi(radar_aux, target_elevations, el_tol=el_tol) + elif radar_aux.scan_type == "ppi": + radar_ppi = radar_aux + else: + warn("Error: unsupported scan type.") + return None + + if qvp is None: + qvp = _create_qvp_object( + radar_aux, + field_names, + qvp_type="rqvp", + start_time=ref_time, + hmax=hmax, + hres=hres, + ) + + # modify metadata + if ref_time is None: + ref_time = datetime_from_radar(radar_ppi) + qvp = _update_qvp_metadata( + qvp, ref_time, qvp.longitude["data"][0], qvp.latitude["data"][0], elev=90.0 + ) + + rmax_km = rmax / 1000.0 + grng_interp = np.ma.masked_all((radar_ppi.nsweeps, qvp.ngates)) + val_interp = dict() + for field_name in field_names: + val_interp.update( + {field_name: np.ma.masked_all((radar_ppi.nsweeps, qvp.ngates))} + ) + + for sweep in range(radar_ppi.nsweeps): + radar_aux = deepcopy(radar_ppi) + radar_aux = radar_aux.extract_sweeps([sweep]) + height = radar_aux.gate_altitude["data"][0, :] + + # compute ground range [Km] + grng = ( + np.sqrt( + np.power(radar_aux.gate_x["data"][0, :], 2.0) + + np.power(radar_aux.gate_y["data"][0, :], 2.0) + ) + / 1000.0 + ) + + # Project ground range to grid + f = interp1d( + height, grng, kind=interp_kind, bounds_error=False, fill_value="extrapolate" + ) + grng_interp[sweep, :] = f(qvp.range["data"]) + + for field_name in field_names: + if field_name not in radar_aux.fields: + warn("Field " + field_name + " not in radar object") + continue + + # Compute QVP for this sweep + values, _ = compute_directional_stats( + radar_aux.fields[field_name]["data"], + avg_type=avg_type, + nvalid_min=nvalid_min, + axis=0, + ) + + # Project to grid + val_interp[field_name][sweep, :] = project_to_vertical( + values, height, qvp.range["data"], interp_kind=interp_kind + ) + + # Compute weight + weight = np.ma.abs(grng_interp - (rmax_km - 1.0)) + weight[grng_interp <= rmax_km - 1.0] = 1.0 / np.power( + weight[grng_interp <= rmax_km - 1.0], 0.0 + ) + + if weight_power == -1: + weight[grng_interp > rmax_km - 1.0] = 0.0 + else: + weight[grng_interp > rmax_km - 1.0] = 1.0 / np.power( + weight[grng_interp > rmax_km - 1.0], weight_power + ) + + for field_name in field_names: + # mask weights where there is no data + mask = np.ma.getmaskarray(val_interp[field_name]) + weight_aux = np.ma.masked_where(mask, weight) + + # Weighted average + qvp_data = np.ma.sum(val_interp[field_name] * weight_aux, axis=0) / np.ma.sum( + weight_aux, axis=0 + ) + + # Put data in radar object + if np.size(qvp.fields[field_name]["data"]) == 0: + qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates) + else: + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) + + return qvp + + +def compute_evp( + radar, + field_names, + lon, + lat, + ref_time=None, + latlon_tol=0.0005, + delta_rng=15000.0, + delta_azi=10, + hmax=10000.0, + hres=250.0, + avg_type="mean", + nvalid_min=1, + interp_kind="none", + qvp=None, +): + """ + Computes enhanced vertical profiles. + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of str + list of field names to add to the QVP + lat, lon : float + latitude and longitude of the point of interest [deg] + ref_time : datetime object + reference time for current radar volume + latlon_tol : float + tolerance in latitude and longitude in deg. + delta_rng, delta_azi : float + maximum range distance [m] and azimuth distance [degree] from the + central point of the evp containing data to average. + hmax : float + The maximum height to plot [m]. + hres : float + The height resolution [m]. + avg_type : str + The type of averaging to perform. Can be either "mean" or "median" + nvalid_min : int + Minimum number of valid points to accept average. + interp_kind : str + type of interpolation when projecting to vertical grid: 'none', + or 'nearest', etc. + 'none' will select from all data points within the regular grid + height bin the closest to the center of the bin. + 'nearest' will select the closest data point to the center of the + height bin regardless if it is within the height bin or not. + Data points can be masked values + If another type of interpolation is selected masked values will be + eliminated from the data points before the interpolation + qvp : QVP object or None + If it is not None this is the QVP object where to store the data from + the current time step. Otherwise a new QVP object will be created + + + Returns + ------- + qvp : qvp object + The computed enhanced vertical profile + + Reference + --------- + Kaltenboeck R., Ryzhkov A. 2016: A freezing rain storm explored with a + C-band polarimetric weather radar using the QVP methodology. + Meteorologische Zeitschrift vol. 26 pp 207-222 + + """ + if avg_type not in ("mean", "median"): + warn("Unsuported statistics " + avg_type) + return None + + radar_aux = deepcopy(radar) + # transform radar into ppi over the required elevation + if radar_aux.scan_type == "rhi": + target_elevations, el_tol = get_target_elevations(radar_aux) + radar_ppi = cross_section_rhi(radar_aux, target_elevations, el_tol=el_tol) + elif radar_aux.scan_type == "ppi": + radar_ppi = radar_aux + else: + warn("Error: unsupported scan type.") + return None + + radar_aux = radar_ppi.extract_sweeps([0]) + + if qvp is None: + qvp = _create_qvp_object( + radar_aux, + field_names, + qvp_type="evp", + start_time=ref_time, + hmax=hmax, + hres=hres, + ) + + # modify metadata + if ref_time is None: + ref_time = datetime_from_radar(radar_aux) + qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0) + + height = dict() + values = dict() + for field_name in field_names: + height.update({field_name: np.array([], dtype=float)}) + values.update({field_name: np.ma.array([], dtype=float)}) + + for sweep in range(radar_ppi.nsweeps): + radar_aux = deepcopy(radar_ppi) + radar_aux = radar_aux.extract_sweeps([sweep]) + + # find nearest gate to lat lon point + ind_ray, _, azi, rng = find_nearest_gate( + radar_aux, lat, lon, latlon_tol=latlon_tol + ) + + if ind_ray is None: + continue + + # find neighbouring gates to be selected + inds_ray, inds_rng = find_neighbour_gates( + radar_aux, azi, rng, delta_azi=delta_azi, delta_rng=delta_rng + ) + + for field_name in field_names: + if field_name not in radar_aux.fields: + warn("Field " + field_name + " not in radar object") + continue + + height[field_name] = np.append( + height[field_name], radar_aux.gate_altitude["data"][ind_ray, inds_rng] + ) + + # keep only data we are interested in + field = radar_aux.fields[field_name]["data"][:, inds_rng] + field = field[inds_ray, :] + + vals, _ = compute_directional_stats( + field, avg_type=avg_type, nvalid_min=nvalid_min, axis=0 + ) + values[field_name] = np.ma.append(values[field_name], vals) + + for field_name in field_names: + # Project to vertical grid: + qvp_data = project_to_vertical( + values[field_name], + height[field_name], + qvp.range["data"], + interp_kind=interp_kind, + ) + + # Put data in radar object + if np.size(qvp.fields[field_name]["data"]) == 0: + qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates) + else: + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) + + return qvp + + +def compute_svp( + radar, + field_names, + lon, + lat, + angle, + ref_time=None, + ang_tol=1.0, + latlon_tol=0.0005, + delta_rng=15000.0, + delta_azi=10, + hmax=10000.0, + hres=250.0, + avg_type="mean", + nvalid_min=1, + interp_kind="none", + qvp=None, +): + """ + Computes slanted vertical profiles. + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of str + list of field names to add to the QVP + lat, lon : float + latitude and longitude of the point of interest [deg] + angle : int or float + If the radar object contains a PPI volume, the sweep number to + use, if it contains an RHI volume the elevation angle. + ref_time : datetime object + reference time for current radar volume + ang_tol : float + If the radar object contains an RHI volume, the tolerance in the + elevation angle for the conversion into PPI + latlon_tol : float + tolerance in latitude and longitude in deg. + delta_rng, delta_azi : float + maximum range distance [m] and azimuth distance [degree] from the + central point of the evp containing data to average. + hmax : float + The maximum height to plot [m]. + hres : float + The height resolution [m]. + avg_type : str + The type of averaging to perform. Can be either "mean" or "median" + nvalid_min : int + Minimum number of valid points to accept average. + interp_kind : str + type of interpolation when projecting to vertical grid: 'none', + or 'nearest', etc. + 'none' will select from all data points within the regular grid + height bin the closest to the center of the bin. + 'nearest' will select the closest data point to the center of the + height bin regardless if it is within the height bin or not. + Data points can be masked values + If another type of interpolation is selected masked values will be + eliminated from the data points before the interpolation + qvp : QVP object or None + If it is not None this is the QVP object where to store the data from + the current time step. Otherwise a new QVP object will be created + + + Returns + ------- + qvp : qvp object + The computed slanted vertical profile + + Reference + --------- + Bukovcic P., Zrnic D., Zhang G. 2017: Winter Precipitation Liquid-Ice + Phase Transitions Revealed with Polarimetric Radar and 2DVD Observations + in Central Oklahoma. JTECH vol. 56 pp 1345-1363 + + """ + if avg_type not in ("mean", "median"): + warn("Unsuported statistics " + avg_type) + return None + + radar_aux = deepcopy(radar) + # transform radar into ppi over the required elevation + if radar_aux.scan_type == "rhi": + radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol) + elif radar_aux.scan_type == "ppi": + radar_aux = radar_aux.extract_sweeps([int(angle)]) + else: + warn("Error: unsupported scan type.") + return None + + if qvp is None: + qvp = _create_qvp_object( + radar_aux, + field_names, + qvp_type="svp", + start_time=ref_time, + hmax=hmax, + hres=hres, + ) + + # modify metadata + if ref_time is None: + ref_time = datetime_from_radar(radar_aux) + qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=qvp.fixed_angle["data"][0]) + + # find nearest gate to lat lon point + ind_ray, _, azi, rng = find_nearest_gate(radar_aux, lat, lon, latlon_tol=latlon_tol) + + if ind_ray is None: + warn("No data in selected point") + qvp_data = np.ma.masked_all(qvp.ngates) + + for field_name in field_names: + # Put data in radar object + if np.size(qvp.fields[field_name]["data"]) == 0: + qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates) + else: + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) + return qvp + + # find neighbouring gates to be selected + inds_ray, inds_rng = find_neighbour_gates( + radar_aux, azi, rng, delta_azi=delta_azi, delta_rng=delta_rng + ) + height = radar_aux.gate_altitude["data"][ind_ray, inds_rng] + + for field_name in field_names: + if field_name not in radar_aux.fields: + warn("Field " + field_name + " not in radar object") + qvp_data = np.ma.masked_all(qvp.ngates) + else: + # keep only data we are interested in + field = radar_aux.fields[field_name]["data"][:, inds_rng] + field = field[inds_ray, :] + + # compute values + values, _ = compute_directional_stats( + field, avg_type=avg_type, nvalid_min=nvalid_min, axis=0 + ) + + # Project to vertical grid: + qvp_data = project_to_vertical( + values, height, qvp.range["data"], interp_kind=interp_kind + ) + + # Put data in radar object + if np.size(qvp.fields[field_name]["data"]) == 0: + qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates) + else: + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) + + return qvp + + +def compute_vp( + radar, + field_names, + lon, + lat, + ref_time=None, + latlon_tol=0.0005, + hmax=10000.0, + hres=50.0, + interp_kind="none", + qvp=None, +): + """ + Computes vertical profiles. + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of str + list of field names to add to the QVP + lat, lon : float + latitude and longitude of the point of interest [deg] + ref_time : datetime object + reference time for current radar volume + latlon_tol : float + tolerance in latitude and longitude in deg. + hmax : float + The maximum height to plot [m]. + hres : float + The height resolution [m]. + interp_kind : str + type of interpolation when projecting to vertical grid: 'none', + or 'nearest', etc. + 'none' will select from all data points within the regular grid + height bin the closest to the center of the bin. + 'nearest' will select the closest data point to the center of the + height bin regardless if it is within the height bin or not. + Data points can be masked values + If another type of interpolation is selected masked values will be + eliminated from the data points before the interpolation + qvp : QVP object or None + If it is not None this is the QVP object where to store the data from + the current time step. Otherwise a new QVP object will be created + + + Returns + ------- + qvp : qvp object + The computed vertical profile + + """ + radar_aux = deepcopy(radar) + # transform radar into ppi over the required elevation + if radar_aux.scan_type == "rhi": + target_elevations, el_tol = get_target_elevations(radar_aux) + radar_ppi = cross_section_rhi(radar_aux, target_elevations, el_tol=el_tol) + elif radar_aux.scan_type == "ppi": + radar_ppi = radar_aux + else: + warn("Error: unsupported scan type.") + return None + + if qvp is None: + qvp = _create_qvp_object( + radar_ppi, + field_names, + qvp_type="vp", + start_time=ref_time, + hmax=hmax, + hres=hres, + ) + + # modify metadata + if ref_time is None: + ref_time = datetime_from_radar(radar_ppi) + qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0) + + height = dict() + values = dict() + for field_name in field_names: + height.update({field_name: np.array([], dtype=float)}) + values.update({field_name: np.ma.array([], dtype=float)}) + + for sweep in range(radar_ppi.nsweeps): + radar_aux = deepcopy(radar_ppi.extract_sweeps([sweep])) + + # find nearest gate to lat lon point + ind_ray, ind_rng, _, _ = find_nearest_gate( + radar_aux, lat, lon, latlon_tol=latlon_tol + ) + + if ind_ray is None: + continue + + for field_name in field_names: + if field_name not in radar_aux.fields: + warn("Field " + field_name + " not in radar object") + continue + + height[field_name] = np.append( + height[field_name], radar_aux.gate_altitude["data"][ind_ray, ind_rng] + ) + values[field_name] = np.ma.append( + values[field_name], + radar_aux.fields[field_name]["data"][ind_ray, ind_rng], + ) + + for field_name in field_names: + # Project to vertical grid: + qvp_data = project_to_vertical( + values[field_name], + height[field_name], + qvp.range["data"], + interp_kind=interp_kind, + ) + + # Put data in radar object + if np.size(qvp.fields[field_name]["data"]) == 0: + qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates) + else: + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) + + return qvp + + +def compute_ts_along_coord( + radar, + field_names, + mode="ALONG_AZI", + fixed_range=None, + fixed_azimuth=None, + fixed_elevation=None, + ang_tol=1.0, + rng_tol=50.0, + value_start=None, + value_stop=None, + ref_time=None, + acoord=None, +): + """ + Computes time series along a particular antenna coordinate, i.e. along + azimuth, elevation or range + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : str + list of field names + mode : str + coordinate to extract data along. Can be ALONG_AZI, ALONG_ELE or + ALONG_RNG + fixed_range, fixed_azimuth, fixed_elevation : float + The fixed range [m], azimuth [deg] or elevation [deg] to extract. + In each mode two of these parameters have to be defined. If they are + not defined they default to 0. + ang_tol, rng_tol : float + The angle tolerance [deg] and range tolerance [m] around the fixed + range or azimuth/elevation + value_start, value_stop : float + The minimum and maximum value at which the data along a coordinate + start and stop + ref_time : datetime object + reference time for current radar volume + acoord : acoord object or None + If it is not None this is the object where to store the data from the + current time step. Otherwise a new acoord object will be created + + Returns + ------- + acoord : acoord object + The computed data along a coordinate + + """ + if mode == "ALONG_RNG": + if value_start is None: + value_start = 0.0 + if value_stop is None: + value_stop = radar.range["data"][-1] + if fixed_azimuth is None: + fixed_azimuth = 0.0 + if fixed_elevation is None: + fixed_elevation = 0.0 + elif mode == "ALONG_AZI": + if value_start is None: + value_start = np.min(radar.azimuth["data"]) + if value_stop is None: + value_stop = np.max(radar.azimuth["data"]) + if fixed_range is None: + fixed_range = 0.0 + if fixed_elevation is None: + fixed_elevation = 0.0 + elif mode == "ALONG_ELE": + if value_start is None: + value_start = np.min(radar.elevation["data"]) + if value_stop is None: + value_stop = np.max(radar.elevation["data"]) + if fixed_range is None: + fixed_range = 0.0 + if fixed_azimuth is None: + fixed_azimuth = 0.0 + else: + warn("Unknown time series of coordinate mode " + mode) + return None + + if mode == "ALONG_RNG": + # rng_values : range + # fixed_angle: elevation + # elevation: elevation + # azimuth: azimuth + vals_dict = {} + for field_name in field_names: + rng_values, vals, _, _ = get_data_along_rng( + radar, + field_name, + [fixed_elevation], + [fixed_azimuth], + ang_tol=ang_tol, + rmin=value_start, + rmax=value_stop, + ) + + if vals.size == 0: + warn("No data found") + return None + vals_dict.update({field_name: vals}) + fixed_angle = fixed_elevation + elevation = fixed_elevation + azimuth = fixed_azimuth + atol = rng_tol + elif mode == "ALONG_AZI": + # rng_values : azimuth + # fixed_angle : elevation + # elevation : elevation + # azimuth : range + vals_dict = {} + for field_name in field_names: + rng_values, vals, _, _ = get_data_along_azi( + radar, + field_name, + [fixed_range], + [fixed_elevation], + rng_tol=rng_tol, + ang_tol=ang_tol, + azi_start=value_start, + azi_stop=value_stop, + ) + + if vals.size == 0: + warn("No data found") + return None + vals_dict.update({field_name: vals}) + fixed_angle = fixed_elevation + elevation = fixed_elevation + azimuth = fixed_range + atol = ang_tol + else: + # rng_values : elevation + # fixed_angle : azimuth + # elevation : range + # azimuth : azimuth + vals_dict = {} + for field_name in field_names: + rng_values, vals, _, _ = get_data_along_ele( + radar, + field_name, + [fixed_range], + [fixed_azimuth], + rng_tol=rng_tol, + ang_tol=ang_tol, + ele_min=value_start, + ele_max=value_stop, + ) + + if vals.size == 0: + warn("No data found") + return None + vals_dict.update({field_name: vals}) + fixed_angle = fixed_azimuth + elevation = fixed_range + azimuth = fixed_azimuth + atol = ang_tol + + if acoord is None: + acoord = _create_along_coord_object( + radar, field_names, rng_values, fixed_angle, mode, start_time=ref_time + ) + + if not np.allclose(rng_values, acoord.range["data"], rtol=1e5, atol=atol): + warn("Unable to add data. xvalues different from previous ones") + return None + + # modify metadata + if ref_time is None: + ref_time = datetime_from_radar(radar) + acoord = _update_along_coord_metadata(acoord, ref_time, elevation, azimuth) + + # Put data in radar object + for field_name in field_names: + if np.size(acoord.fields[field_name]["data"]) == 0: + acoord.fields[field_name]["data"] = vals_dict[field_name].reshape( + 1, acoord.ngates + ) + else: + acoord.fields[field_name]["data"] = np.ma.concatenate( + ( + acoord.fields[field_name]["data"], + vals_dict[field_name].reshape(1, acoord.ngates), + ) + ) + + return acoord + + +def project_to_vertical( + data_in, data_height, grid_height, interp_kind="none", fill_value=-9999.0 +): + """ + Projects radar data to a regular vertical grid + + Parameters + ---------- + data_in : ndarray 1D + the radar data to project + data_height : ndarray 1D + the height of each radar point + grid_height : ndarray 1D + the regular vertical grid to project to + interp_kind : str + The type of interpolation to use: 'none' or 'nearest' + fill_value : float + The fill value used for interpolation + + Returns + ------- + data_out : ndarray 1D + The projected data + + """ + if data_in.size == 0: + data_out = np.ma.masked_all(grid_height.size) + return data_out + + if interp_kind == "none": + hres = grid_height[1] - grid_height[0] + data_out = np.ma.masked_all(grid_height.size) + for ind_r, h in enumerate(grid_height): + ind_h = find_rng_index(data_height, h, rng_tol=hres / 2.0) + if ind_h is None: + continue + data_out[ind_r] = data_in[ind_h] + elif interp_kind == "nearest": + data_filled = data_in.filled(fill_value=fill_value) + f = interp1d( + data_height, + data_filled, + kind=interp_kind, + bounds_error=False, + fill_value=fill_value, + ) + data_out = np.ma.masked_values(f(grid_height), fill_value) + else: + valid = np.logical_not(np.ma.getmaskarray(data_in)) + height_valid = data_height[valid] + data_valid = data_in[valid] + f = interp1d( + height_valid, + data_valid, + kind=interp_kind, + bounds_error=False, + fill_value=fill_value, + ) + data_out = np.ma.masked_values(f(grid_height), fill_value) + + return data_out + + +def find_rng_index(rng_vec, rng, rng_tol=0.0): + """ + Find the range index corresponding to a particular range + + Parameters + ---------- + rng_vec : float array + The range data array where to look for + rng : float + The range to search + rng_tol : float + Tolerance [m] + + Returns + ------- + ind_rng : int + The range index + + """ + dist = np.abs(rng_vec - rng) + ind_rng = np.argmin(dist) + if dist[ind_rng] > rng_tol: + return None + + return ind_rng + + +def get_target_elevations(radar_in): + """ + Gets RHI target elevations + + Parameters + ---------- + radar_in : Radar object + current radar object + + Returns + ------- + target_elevations : 1D-array + Azimuth angles + el_tol : float + azimuth tolerance + """ + sweep_start = radar_in.sweep_start_ray_index["data"][0] + sweep_end = radar_in.sweep_end_ray_index["data"][0] + target_elevations = np.sort(radar_in.elevation["data"][sweep_start : sweep_end + 1]) + el_tol = np.median(target_elevations[1:] - target_elevations[:-1]) + + return target_elevations, el_tol + + +def find_nearest_gate(radar, lat, lon, latlon_tol=0.0005): + """ + Find the radar gate closest to a lat,lon point + + Parameters + ---------- + radar : radar object + the radar object + lat, lon : float + The position of the point + latlon_tol : float + The tolerance around this point + + Returns + ------- + ind_ray, ind_rng : int + The ray and range index + azi, rng : float + the range and azimuth position of the gate + + """ + # find gates close to lat lon point + inds_ray_aux, inds_rng_aux = np.where( + np.logical_and( + np.logical_and( + radar.gate_latitude["data"] < lat + latlon_tol, + radar.gate_latitude["data"] > lat - latlon_tol, + ), + np.logical_and( + radar.gate_longitude["data"] < lon + latlon_tol, + radar.gate_longitude["data"] > lon - latlon_tol, + ), + ) + ) + + if inds_ray_aux.size == 0: + warn( + "No data found at point lat " + + str(lat) + + " +- " + + str(latlon_tol) + + " lon " + + str(lon) + + " +- " + + str(latlon_tol) + + " deg" + ) + + return None, None, None, None + + # find closest latitude + ind_min = np.argmin( + np.abs(radar.gate_latitude["data"][inds_ray_aux, inds_rng_aux] - lat) + ) + ind_ray = inds_ray_aux[ind_min] + ind_rng = inds_rng_aux[ind_min] + + azi = radar.azimuth["data"][ind_ray] + rng = radar.range["data"][ind_rng] + + return ind_ray, ind_rng, azi, rng + + +def find_neighbour_gates(radar, azi, rng, delta_azi=None, delta_rng=None): + """ + Find the neighbouring gates within +-delta_azi and +-delta_rng + + Parameters + ---------- + radar : radar object + the radar object + azi, rng : float + The azimuth [deg] and range [m] of the central gate + delta_azi, delta_rng : float + The extend where to look for + + Returns + ------- + inds_ray_aux, ind_rng_aux : int + The indices (ray, rng) of the neighbouring gates + + """ + # find gates close to lat lon point + if delta_azi is None: + inds_ray = np.ma.arange(radar.azimuth["data"].size) + else: + azi_max = azi + delta_azi + azi_min = azi - delta_azi + if azi_max > 360.0: + azi_max -= 360.0 + if azi_min < 0.0: + azi_min += 360.0 + if azi_max > azi_min: + inds_ray = np.where( + np.logical_and( + radar.azimuth["data"] < azi_max, radar.azimuth["data"] > azi_min + ) + )[0] + else: + inds_ray = np.where( + np.logical_or( + radar.azimuth["data"] > azi_min, radar.azimuth["data"] < azi_max + ) + )[0] + if delta_rng is None: + inds_rng = np.ma.arange(radar.range["data"].size) + else: + inds_rng = np.where( + np.logical_and( + radar.range["data"] < rng + delta_rng, + radar.range["data"] > rng - delta_rng, + ) + )[0] + + return inds_ray, inds_rng + + +def get_data_along_rng( + radar, field_name, fix_elevations, fix_azimuths, ang_tol=1.0, rmin=None, rmax=None +): + """ + Get data at particular (azimuths, elevations) + + Parameters + ---------- + radar : radar object + the radar object where the data is + field_name : str + name of the field to filter + fix_elevations, fix_azimuths: list of floats + List of elevations, azimuths couples [deg] + ang_tol : float + Tolerance between the nominal angle and the radar angle [deg] + rmin, rmax: float + Min and Max range of the obtained data [m] + + Returns + ------- + xvals : 1D float array + The ranges of each azi, ele pair + yvals : 1D float array + The values + valid_azi, valid_ele : float arrays + The azi, ele pairs + + """ + if rmin is None: + rmin = 0.0 + if rmax is None: + rmax = np.max(radar.range["data"]) + + rng_mask = np.logical_and(radar.range["data"] >= rmin, radar.range["data"] <= rmax) + + x = radar.range["data"][rng_mask] + + xvals = [] + yvals = [] + valid_azi = [] + valid_ele = [] + if radar.scan_type in ("ppi", "vertical_pointing"): + for ele, azi in zip(fix_elevations, fix_azimuths): + if radar.scan_type == "vertical_pointing": + dataset_line = deepcopy(radar) + xvals.extend(x) + else: + ind_sweep = find_ang_index( + radar.fixed_angle["data"], ele, ang_tol=ang_tol + ) + if ind_sweep is None: + warn("No elevation angle found for fix_elevation " + str(ele)) + continue + + new_dataset = radar.extract_sweeps([ind_sweep]) + + try: + dataset_line = cross_section_ppi(new_dataset, [azi], az_tol=ang_tol) + except OSError: + warn( + " No data found at azimuth " + + str(azi) + + " and elevation " + + str(ele) + ) + continue + xvals.append(x) + yvals.append(dataset_line.fields[field_name]["data"][0, rng_mask]) + valid_azi.append(dataset_line.azimuth["data"][0]) + valid_ele.append(dataset_line.elevation["data"][0]) + else: + for ele, azi in zip(fix_elevations, fix_azimuths): + ind_sweep = find_ang_index(radar.fixed_angle["data"], azi, ang_tol=ang_tol) + if ind_sweep is None: + warn("No azimuth angle found for fix_azimuth " + str(azi)) + continue + new_dataset = radar.extract_sweeps([ind_sweep]) + + try: + dataset_line = cross_section_rhi(new_dataset, [ele], el_tol=ang_tol) + except OSError: + warn( + " No data found at azimuth " + + str(azi) + + " and elevation " + + str(ele) + ) + continue + yvals.extend(dataset_line.fields[field_name]["data"][0, rng_mask]) + xvals.extend(x) + valid_azi.append(dataset_line.azimuth["data"][0]) + valid_ele.append(dataset_line.elevation["data"][0]) + + return np.array(xvals), np.ma.array(yvals), valid_azi, valid_ele + + +def get_data_along_azi( + radar, + field_name, + fix_ranges, + fix_elevations, + rng_tol=50.0, + ang_tol=1.0, + azi_start=None, + azi_stop=None, +): + """ + Get data at particular (ranges, elevations) + + Parameters + ---------- + radar : radar object + the radar object where the data is + field_name : str + name of the field to filter + fix_ranges, fix_elevations: list of floats + List of ranges [m], elevations [deg] couples + rng_tol : float + Tolerance between the nominal range and the radar range [m] + ang_tol : float + Tolerance between the nominal angle and the radar angle [deg] + azi_start, azi_stop: float + Start and stop azimuth angle of the data [deg] + + Returns + ------- + xvals : 1D float array + The ranges of each rng, ele pair + yvals : 1D float array + The values + valid_rng, valid_ele : float arrays + The rng, ele pairs + + """ + if azi_start is None: + azi_start = np.min(radar.azimuth["data"]) + if azi_stop is None: + azi_stop = np.max(radar.azimuth["data"]) + + yvals = [] + xvals = [] + valid_rng = [] + valid_ele = [] + for rng, ele in zip(fix_ranges, fix_elevations): + ind_rng = find_rng_index(radar.range["data"], rng, rng_tol=rng_tol) + if ind_rng is None: + warn("No range gate found for fix_range " + str(rng)) + continue + + if radar.scan_type == "ppi": + ind_sweep = find_ang_index(radar.fixed_angle["data"], ele, ang_tol=ang_tol) + if ind_sweep is None: + warn("No elevation angle found for fix_elevation " + str(ele)) + continue + new_dataset = radar.extract_sweeps([ind_sweep]) + else: + try: + new_dataset = cross_section_rhi(radar, [ele], el_tol=ang_tol) + except OSError: + warn( + " No data found at range " + str(rng) + " and elevation " + str(ele) + ) + continue + if azi_start < azi_stop: + azi_mask = np.logical_and( + new_dataset.azimuth["data"] >= azi_start, + new_dataset.azimuth["data"] <= azi_stop, + ) + else: + azi_mask = np.logical_or( + new_dataset.azimuth["data"] >= azi_start, + new_dataset.azimuth["data"] <= azi_stop, + ) + yvals.extend(new_dataset.fields[field_name]["data"][azi_mask, ind_rng]) + xvals.extend(new_dataset.azimuth["data"][azi_mask]) + valid_rng.append(new_dataset.range["data"][ind_rng]) + valid_ele.append(new_dataset.elevation["data"][0]) + + return np.array(xvals), np.ma.array(yvals), valid_rng, valid_ele + + +def get_data_along_ele( + radar, + field_name, + fix_ranges, + fix_azimuths, + rng_tol=50.0, + ang_tol=1.0, + ele_min=None, + ele_max=None, +): + """ + Get data at particular (ranges, azimuths) + + Parameters + ---------- + radar : radar object + the radar object where the data is + field_name : str + name of the field to filter + fix_ranges, fix_azimuths: list of floats + List of ranges [m], azimuths [deg] couples + rng_tol : float + Tolerance between the nominal range and the radar range [m] + ang_tol : float + Tolerance between the nominal angle and the radar angle [deg] + ele_min, ele_max: float + Min and max elevation angle [deg] + + Returns + ------- + xvals : 1D float array + The ranges of each rng, ele pair + yvals : 1D float array + The values + valid_rng, valid_ele : float arrays + The rng, ele pairs + + """ + if ele_min is None: + ele_min = np.min(radar.elevation["data"]) + if ele_max is None: + ele_max = np.max(radar.elevation["data"]) + + yvals = [] + xvals = [] + valid_rng = [] + valid_azi = [] + for rng, azi in zip(fix_ranges, fix_azimuths): + ind_rng = find_rng_index(radar.range["data"], rng, rng_tol=rng_tol) + if ind_rng is None: + warn("No range gate found for fix_range " + str(rng)) + continue + + if radar.scan_type == "ppi": + try: + new_dataset = cross_section_ppi(radar, [azi], az_tol=ang_tol) + except OSError: + warn( + " No data found at range " + str(rng) + " and elevation " + str(azi) + ) + continue + else: + ind_sweep = find_ang_index(radar.fixed_angle["data"], azi, ang_tol=ang_tol) + if ind_sweep is None: + warn("No azimuth angle found for fix_azimuth " + str(azi)) + continue + new_dataset = radar.extract_sweeps([ind_sweep]) + + ele_mask = np.logical_and( + new_dataset.elevation["data"] >= ele_min, + new_dataset.elevation["data"] <= ele_max, + ) + yvals.extend(new_dataset.fields[field_name]["data"][ele_mask, ind_rng]) + xvals.extend(new_dataset.elevation["data"][ele_mask]) + valid_rng.append(new_dataset.range["data"][ind_rng]) + valid_azi.append(new_dataset.elevation["data"][0]) + + return np.array(xvals), np.ma.array(yvals), valid_rng, valid_azi + + +def find_ang_index(ang_vec, ang, ang_tol=0.0): + """ + Find the angle index corresponding to a particular fixed angle + + Parameters + ---------- + ang_vec : float array + The angle data array where to look for + ang : float + The angle to search + ang_tol : float + Tolerance [deg] + + Returns + ------- + ind_ang : int + The angle index + + """ + dist = np.abs(ang_vec - ang) + ind_ang = np.argmin(dist) + if dist[ind_ang] > ang_tol: + return None + + return ind_ang + + +def _create_qvp_object( + radar, field_names, qvp_type="qvp", start_time=None, hmax=10000.0, hres=200.0 +): + """ + Creates a QVP object containing fields from a radar object that can + be used to plot and produce the quasi vertical profile + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of strings + Radar fields to use for QVP calculation. + qvp_type : str + Type of QVP. Can be qvp, rqvp, evp + start_time : datetime object + the QVP start time + hmax : float + The maximum height of the QVP [m]. Default 10000. + hres : float + The QVP range resolution [m]. Default 50 + + Returns + ------- + qvp : Radar-like object + A quasi vertical profile object containing fields + from a radar object + + """ + qvp = deepcopy(radar) + + # prepare space for field + qvp.fields = dict() + for field_name in field_names: + qvp.add_field(field_name, deepcopy(radar.fields[field_name])) + qvp.fields[field_name]["data"] = np.array([], dtype="float64") + + # fixed radar objects parameters + qvp.range["data"] = np.arange(hmax / hres) * hres + hres / 2.0 + qvp.ngates = len(qvp.range["data"]) + + if start_time is None: + qvp.time["units"] = radar.time["units"] + else: + qvp.time["units"] = make_time_unit_str(start_time) + + qvp.time["data"] = np.array([], dtype="float64") + qvp.scan_type = qvp_type + qvp.sweep_number["data"] = np.array([0], dtype="int32") + qvp.nsweeps = 1 + qvp.sweep_mode["data"] = np.array([qvp_type]) + qvp.sweep_start_ray_index["data"] = np.array([0], dtype="int32") + + if qvp.rays_are_indexed is not None: + qvp.rays_are_indexed["data"] = np.array([qvp.rays_are_indexed["data"][0]]) + if qvp.ray_angle_res is not None: + qvp.ray_angle_res["data"] = np.array([qvp.ray_angle_res["data"][0]]) + + if qvp_type in ("rqvp", "evp", "vp"): + qvp.fixed_angle["data"] = np.array([90.0], dtype="float64") + + # ray dependent radar objects parameters + qvp.sweep_end_ray_index["data"] = np.array([-1], dtype="int32") + qvp.rays_per_sweep["data"] = np.array([0], dtype="int32") + qvp.azimuth["data"] = np.array([], dtype="float64") + qvp.elevation["data"] = np.array([], dtype="float64") + qvp.nrays = 0 + + return qvp + + +def _create_along_coord_object( + radar, field_names, rng_values, fixed_angle, mode, start_time=None +): + """ + Creates an along coord object containing fields from a radar object that + can be used to plot and produce the time series along a coordinate + + Parameters + ---------- + radar : Radar + Radar object used. + field_names : list of strings + Radar fields to use for QVP calculation. + rng_values : 1D-array + The values to put in the range field + fixed_angle : float + the fixed angle + mode : str + The along coord mode, can be ALONG_AZI, ALONG_ELE, ALONG_RNG + start_time : datetime object + the acoord start time + + Returns + ------- + acoord : Radar-like object + An along coordinate object containing fields from a radar object + + """ + acoord = deepcopy(radar) + + # prepare space for field + acoord.fields = dict() + for field_name in field_names: + acoord.add_field(field_name, deepcopy(radar.fields[field_name])) + acoord.fields[field_name]["data"] = np.array([], dtype="float64") + + # fixed radar objects parameters + acoord.range["data"] = rng_values + acoord.ngates = len(acoord.range["data"]) + + if start_time is None: + acoord.time["units"] = radar.time["units"] + else: + acoord.time["units"] = make_time_unit_str(start_time) + + acoord.time["data"] = np.array([], dtype="float64") + acoord.scan_type = mode + acoord.sweep_number["data"] = np.array([0], dtype="int32") + acoord.nsweeps = 1 + acoord.sweep_mode["data"] = np.array([mode]) + acoord.sweep_start_ray_index["data"] = np.array([0], dtype="int32") + + if acoord.rays_are_indexed is not None: + acoord.rays_are_indexed["data"] = np.array([acoord.rays_are_indexed["data"][0]]) + if acoord.ray_angle_res is not None: + acoord.ray_angle_res["data"] = np.array([acoord.ray_angle_res["data"][0]]) + + acoord.fixed_angle["data"] = np.array([fixed_angle], dtype="float64") + + # ray dependent radar objects parameters + acoord.sweep_end_ray_index["data"] = np.array([-1], dtype="int32") + acoord.rays_per_sweep["data"] = np.array([0], dtype="int32") + acoord.azimuth["data"] = np.array([], dtype="float64") + acoord.elevation["data"] = np.array([], dtype="float64") + acoord.nrays = 0 + + return acoord + + +def _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0): + """ + updates a QVP object metadata with data from the current radar volume + + Parameters + ---------- + qvp : QVP object + QVP object + ref_time : datetime object + the current radar volume reference time + + Returns + ------- + qvp : QVP object + The updated QVP object + + """ + start_time = num2date(0, qvp.time["units"], qvp.time["calendar"]) + qvp.time["data"] = np.append( + qvp.time["data"], (ref_time - start_time).total_seconds() + ) + qvp.sweep_end_ray_index["data"][0] += 1 + qvp.rays_per_sweep["data"][0] += 1 + qvp.nrays += 1 + + qvp.azimuth["data"] = np.ones((qvp.nrays,), dtype="float64") * 0.0 + qvp.elevation["data"] = np.ones((qvp.nrays,), dtype="float64") * elev + + qvp.gate_longitude["data"] = np.ones((qvp.nrays, qvp.ngates), dtype="float64") * lon + qvp.gate_latitude["data"] = np.ones((qvp.nrays, qvp.ngates), dtype="float64") * lat + qvp.gate_altitude["data"] = ma_broadcast_to( + qvp.range["data"], (qvp.nrays, qvp.ngates) + ) + + return qvp + + +def _update_along_coord_metadata(acoord, ref_time, elevation, azimuth): + """ + updates an along coordinate object metadata with data from the current + radar volume + + Parameters + ---------- + acoord : along coordinate object + along coordinate object + ref_time : datetime object + the current radar volume reference time + elevation, azimuth : 1D-array + the elevation and azimuth value of the data selected + + Returns + ------- + acoord : along coordinate object + The updated along coordinate object + + """ + start_time = num2date(0, acoord.time["units"], acoord.time["calendar"]) + acoord.time["data"] = np.append( + acoord.time["data"], (ref_time - start_time).total_seconds() + ) + acoord.sweep_end_ray_index["data"][0] += 1 + acoord.rays_per_sweep["data"][0] += 1 + acoord.nrays += 1 + + acoord.azimuth["data"] = np.append(acoord.azimuth["data"], azimuth) + acoord.elevation["data"] = np.append(acoord.elevation["data"], elevation) + + return acoord diff --git a/pyart/util/__init__.py b/pyart/util/__init__.py index ee8e1ae47cd..b3591df6892 100644 --- a/pyart/util/__init__.py +++ b/pyart/util/__init__.py @@ -14,6 +14,7 @@ from .circular_stats import interval_std # noqa from .circular_stats import mean_of_two_angles # noqa from .circular_stats import mean_of_two_angles_deg # noqa +from .circular_stats import compute_directional_stats # noqa from .columnsect import for_azimuth # noqa from .columnsect import get_column_rays # noqa from .columnsect import get_field_location # noqa @@ -32,6 +33,7 @@ determine_sweeps, subset_radar, to_vpt, + ma_broadcast_to, ) from .sigmath import ( # noqa angular_texture_2d, diff --git a/pyart/util/circular_stats.py b/pyart/util/circular_stats.py index 639cbb83b26..c69359c978b 100644 --- a/pyart/util/circular_stats.py +++ b/pyart/util/circular_stats.py @@ -10,6 +10,41 @@ # https://en.wikipedia.org/wiki/Mean_of_circular_quantities +def compute_directional_stats(field, avg_type="mean", nvalid_min=1, axis=0): + """ + Computes the mean or the median along one of the axis (ray or range) + + Parameters + ---------- + field : ndarray + the radar field data + avg_type :str + the type of average: 'mean' or 'median' + nvalid_min : int + the minimum number of points to consider the stats valid. Default 1 + axis : int + the axis along which to compute (0=ray, 1=range) + + Returns + ------- + values : ndarray 1D + The resultant statistics + nvalid : ndarray 1D + The number of valid points used in the computation + """ + if avg_type == "mean": + values = np.ma.mean(field, axis=axis) + else: + values = np.ma.median(field, axis=axis) + + # Set to non-valid if there is not a minimum number of valid gates + valid = np.logical_not(np.ma.getmaskarray(field)) + nvalid = np.sum(valid, axis=axis, dtype=int) + values[nvalid < nvalid_min] = np.ma.masked + + return values, nvalid + + def mean_of_two_angles(angles1, angles2): """ Compute the element by element mean of two sets of angles. diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index a8f056728c6..ba7063a98cb 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -689,3 +689,33 @@ def image_mute_radar(radar, field, mute_field, mute_threshold, field_threshold=N muted_dict["long_name"] = "Muted " + field radar.add_field("muted_" + field, muted_dict) return radar + + +def ma_broadcast_to(array, tup): + """ + Is used to guarantee that a masked array can be broadcasted without + loosing the mask + + Parameters + ---------- + array : Numpy masked array or normal array + tup : shape as tuple + + Returns + ------- + broadcasted_array + The broadcasted numpy array including its mask if available + otherwise only the broadcasted array is returned + + """ + broadcasted_array = np.broadcast_to(array, tup) + + if np.ma.is_masked(array): + initial_mask = np.ma.getmask(array) + initial_fill_value = array.fill_value + broadcasted_mask = np.broadcast_to(initial_mask, tup) + return np.ma.array( + broadcasted_array, mask=broadcasted_mask, fill_value=initial_fill_value + ) + + return broadcasted_array diff --git a/pyart/xradar/accessor.py b/pyart/xradar/accessor.py index 8a374c88ac4..e79a1bcb150 100644 --- a/pyart/xradar/accessor.py +++ b/pyart/xradar/accessor.py @@ -3,7 +3,6 @@ """ - import copy import numpy as np diff --git a/tests/bridge/test_wradlib_bridge.py b/tests/bridge/test_wradlib_bridge.py index e777bdb2d74..d2d6a9cdda5 100644 --- a/tests/bridge/test_wradlib_bridge.py +++ b/tests/bridge/test_wradlib_bridge.py @@ -1,6 +1,5 @@ """ Unit Tests for Py-ART's io/mdv.py module. """ - import numpy as np import pytest diff --git a/tests/core/test_transforms.py b/tests/core/test_transforms.py index ab7527bbf39..d3a5ac20886 100644 --- a/tests/core/test_transforms.py +++ b/tests/core/test_transforms.py @@ -188,3 +188,15 @@ def test_cartesian_to_geographic_aeqd(): lon, lat = transforms.cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R) assert_almost_equal(lon, -100.0, 3) assert_almost_equal(lat, 40.0, 3) + + +def test_cartesian_to_antenna(): + r, az, el = transforms.cartesian_to_antenna( + np.array([1000, 3000, 2000]), + np.array([1000, 2000, 3000]), + np.array([500, 1000, 1500]), + ) + + assert np.allclose(r, [1500, 3741.6573, 3905.1248]) + assert np.allclose(az, [45.0, 56.30993247, 33.69006753]) + assert np.allclose(el, [19.47122063, 15.50135957, 22.5885387]) diff --git a/tests/graph/test_cm.py b/tests/graph/test_cm.py index 562ef1a976b..4ad9841d5bd 100644 --- a/tests/graph/test_cm.py +++ b/tests/graph/test_cm.py @@ -1,6 +1,5 @@ """ Unit Tests for Py-ART's graph/cm.py module. """ - import matplotlib from pyart.graph import cm diff --git a/tests/graph/test_cm_colorblind.py b/tests/graph/test_cm_colorblind.py index 5060ff01433..e71ff0ccee5 100644 --- a/tests/graph/test_cm_colorblind.py +++ b/tests/graph/test_cm_colorblind.py @@ -1,6 +1,5 @@ """ Unit Tests for Py-ART's graph/cm_colorblind.py module. """ - import matplotlib from pyart.graph import cm_colorblind diff --git a/tests/graph/test_gridmapdisplay.py b/tests/graph/test_gridmapdisplay.py index b0ce20c2542..7068bb142cd 100644 --- a/tests/graph/test_gridmapdisplay.py +++ b/tests/graph/test_gridmapdisplay.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/gridmapdisplay.py module. """ + # execute this script to create figure_gridmapdisplay_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/graph/test_gridmapdisplay_basemap.py b/tests/graph/test_gridmapdisplay_basemap.py index 4e31d2095ab..34efc2e34de 100644 --- a/tests/graph/test_gridmapdisplay_basemap.py +++ b/tests/graph/test_gridmapdisplay_basemap.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/gridmapdisplay_basemap.py module. """ + # execute this script to create figure_gridmapdisplay_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/graph/test_radarmapdisplay.py b/tests/graph/test_radarmapdisplay.py index 9a5680b77e1..51d3976d78f 100644 --- a/tests/graph/test_radarmapdisplay.py +++ b/tests/graph/test_radarmapdisplay.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/radarmapdisplay.py module. """ + # execute this script to create figure_plot_radar_display_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/graph/test_radarmapdisplay_basemap.py b/tests/graph/test_radarmapdisplay_basemap.py index dcee1412c5b..a811010644d 100644 --- a/tests/graph/test_radarmapdisplay_basemap.py +++ b/tests/graph/test_radarmapdisplay_basemap.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/radarmapdisplay.py module. """ + # execute this script to create figure_plot_radar_display_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index ce9e276250e..47348447477 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -144,37 +144,66 @@ def test_hydroclass_semisupervised(): mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small) + hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] assert "units" in hydro_nofreq.keys() assert "standard_name" in hydro_nofreq.keys() assert "long_name" in hydro_nofreq.keys() assert "coordinates" in hydro_nofreq.keys() - assert hydro_nofreq["data"].dtype == "int64" + assert hydro_nofreq["data"].dtype == "uint8" assert_allclose(hydro_nofreq["data"][0][0:5], [6, 6, 6, 6, 6]) assert_allclose(hydro_nofreq["data"][0][-5:], [2, 2, 2, 2, 2]) radar_small.instrument_parameters["frequency"] = {"data": np.array([5e9])} - hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small) + hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] assert_allclose(hydro_freq["data"][0][0:5], [6, 6, 6, 6, 6]) assert_allclose(hydro_freq["data"][0][-5:], [2, 2, 2, 2, 2]) hydro = pyart.retrieve.hydroclass_semisupervised( - radar_small, mass_centers=mass_centers + radar_small, + mass_centers=mass_centers, + compute_entropy=True, + output_distances=True, + ) + assert_allclose(hydro["hydro"]["data"][0][0:5], [6, 6, 6, 6, 6]) + assert_allclose(hydro["hydro"]["data"][0][-5:], [2, 2, 2, 2, 2]) + assert_allclose( + hydro["entropy"]["data"][0][0:5].data, + [0.35945634, 0.35945634, 0.35945634, 0.35945634, 0.35945634], + rtol=1e-5, + ) + assert_allclose( + hydro["entropy"]["data"][0][-5:].data, + [0.232788, 0.232788, 0.232788, 0.232788, 0.232788], + rtol=1e-5, + ) + assert_allclose( + hydro["proportion_CR"]["data"][0][0:5].data, + [0.03524214, 0.03524214, 0.03524214, 0.03524214, 0.03524214], + rtol=1e-5, + ) + assert_allclose( + hydro["proportion_CR"]["data"][0][-5:].data, + [ + 8.336829723993246e-05, + 8.336829723993246e-05, + 8.336829723993246e-05, + 8.336829723993246e-05, + 8.336829723993246e-05, + ], + rtol=1e-5, ) - assert_allclose(hydro["data"][0][0:5], [6, 6, 6, 6, 6]) - assert_allclose(hydro["data"][0][-5:], [2, 2, 2, 2, 2]) def test_data_limits_table(): dlimits_dict = pyart.retrieve.echo_class._data_limits_table() test_dict = { "Zh": (60.0, -10.0), - "ZDR": (5.0, -5.0), + "ZDR": (5, -1.5), "KDP": (7.0, -10.0), "RhoHV": (-5.23, -50.0), + "RelH": (5000, -5000), } - assert isinstance(dlimits_dict, dict) assert dlimits_dict == test_dict @@ -259,16 +288,12 @@ def test_assign_to_class(): ) mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydroclass, min_dist = pyart.retrieve.echo_class._assign_to_class( - zh, zdr, kdp, rhohv, relh, mass_centers - ) + field_dict = {"Zh": zh, "ZDR": zdr, "KDP": kdp, "RhoHV": rhohv, "relH": relh} - assert_allclose(hydroclass[0], [7, 7, 7, 7, 7], atol=1e-7) - assert_allclose( - min_dist[0], - [227.0343910, 227.0343910, 227.0343910, 227.0343910, 227.0343910], - atol=1e-7, + hydroclass, _, _ = pyart.retrieve.echo_class._assign_to_class( + field_dict, mass_centers ) + assert_allclose(hydroclass[0], [7, 7, 7, 7, 7], atol=1e-7) def test_standardize(): diff --git a/tests/retrieve/test_qvp.py b/tests/retrieve/test_qvp.py index d292dfa4ffd..2293ae56f7e 100644 --- a/tests/retrieve/test_qvp.py +++ b/tests/retrieve/test_qvp.py @@ -1,11 +1,188 @@ """ Unit Tests for Py-ART's retrieve/qvp.py module. """ +import datetime + import numpy as np +import pytest +from netCDF4 import num2date from numpy.testing import assert_almost_equal import pyart +@pytest.fixture +def test_radar(): + test_radar = pyart.testing.make_empty_ppi_radar(1000, 360, 1) + test_radar.range["data"] *= 100 + test_radar.elevation["data"] = np.ones(test_radar.elevation["data"].shape) * 10.0 + test_radar.fixed_angle["data"] = np.array([10.0]) + refl = 0.1 * np.arange(test_radar.ngates) + refl = np.tile(refl, (test_radar.nrays, 1)) + test_radar.add_field("reflectivity", {"data": refl}) + return test_radar + + +def test_compute_qvp(test_radar): + qvp = pyart.retrieve.compute_qvp(test_radar, ["reflectivity"], hmax=10000) + qvp_refl = np.array( + [ + 1.8999999999999875, + 2.200000000000011, + 2.3999999999999853, + 2.7000000000000175, + 3.0, + 3.299999999999978, + 3.599999999999994, + 3.9000000000000132, + 4.200000000000027, + 4.5, + 4.7000000000000295, + 5.0, + 5.299999999999969, + 5.599999999999963, + 5.900000000000039, + 6.2000000000000135, + 6.5, + 6.7, + 7.0, + 7.300000000000015, + ] + ) + + assert np.allclose(qvp.fields["reflectivity"]["data"][0, 10:30], qvp_refl) + assert len(qvp.range["data"]) == 200 + assert np.all(qvp.range["data"] == np.arange(25, 10000, 50)) + assert qvp.azimuth["data"][0] == 0 + + +def test_compute_rqvp(test_radar): + rqvp = pyart.retrieve.compute_rqvp( + test_radar, ["reflectivity"], hres=50, hmax=5000, rmax=5000, weight_power=-1 + ) + rqvp_refl = np.array( + [ + 1.8999999999999875, + 2.200000000000011, + 2.3999999999999853, + 2.7000000000000175, + 3.0, + 3.299999999999978, + 3.599999999999994, + 3.9000000000000132, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + ] + ) + rqvp_refl_mask = np.array( + [ + False, + False, + False, + False, + False, + False, + False, + False, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + ] + ) + + assert np.allclose(rqvp.fields["reflectivity"]["data"][0, 10:30], rqvp_refl) + assert np.allclose( + rqvp.fields["reflectivity"]["data"][0, 10:30].mask, rqvp_refl_mask + ) + + assert len(rqvp.range["data"]) == 100 + assert np.all(rqvp.range["data"] == np.arange(25, 5000, 50)) + assert rqvp.azimuth["data"][0] == 0 + + +def test_compute_evp(test_radar): + evp = pyart.retrieve.compute_evp( + test_radar, ["reflectivity"], -97, 36, delta_rng=40000 + ) + + evp_refl = np.array( + [ + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + 32.7, + 33.7, + 35.10000000000001, + 36.5, + 37.89999999999999, + 39.29999999999999, + 40.70000000000001, + ] + ) + + evp_refl_mask = np.array( + [ + [ + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + False, + False, + False, + False, + False, + False, + False, + ] + ] + ) + + assert np.allclose(evp.fields["reflectivity"]["data"][0, 10:30], evp_refl) + assert np.allclose(evp.fields["reflectivity"]["data"][0, 10:30].mask, evp_refl_mask) + + assert len(evp.range["data"]) == 40 + assert np.all(evp.range["data"] == np.arange(125, 10000, 250)) + assert evp.azimuth["data"][0] == 0 + + def test_quasi_vertical_profile(): test_radar = pyart.testing.make_target_radar() height = np.arange(0, 1000, 200) @@ -179,3 +356,209 @@ def test_quasi_vertical_profile(): assert_almost_equal(qvp["height"], qvp_height, 3) assert_almost_equal(qvp["range"], qvp_range, 3) assert_almost_equal(qvp["reflectivity"], qvp_reflectivity, 3) + + +def test_project_to_vertical(): + z = np.array([100, 200, 500, 1000, 2000]) + data = np.ma.array(np.linspace(0, 50, len(z))) + znew = np.arange(0, 2300, 100) + data_out = pyart.retrieve.qvp.project_to_vertical(data, z, znew) + data_out2 = pyart.retrieve.qvp.project_to_vertical( + data, z, znew, interp_kind="nearest" + ) + ref_out = np.ma.array( + data=[ + np.nan, + 0.0, + 12.5, + np.nan, + np.nan, + 25.0, + np.nan, + np.nan, + np.nan, + np.nan, + 37.5, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + 50.0, + np.nan, + np.nan, + ], + mask=[ + True, + False, + False, + True, + True, + False, + True, + True, + True, + True, + False, + True, + True, + True, + True, + True, + True, + True, + True, + True, + False, + True, + True, + ], + ) + ref_out2 = np.ma.array( + data=[ + np.nan, + 0.0, + 12.5, + 12.5, + 25.0, + 25.0, + 25.0, + 25.0, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + np.nan, + np.nan, + ], + mask=[ + True, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + True, + True, + ], + ) + + assert np.all(data_out == ref_out) + assert np.all(data_out2 == ref_out2) + + +def test_find_rng_index(): + vec = np.arange(0, 1000, 10) + idx = pyart.retrieve.qvp.find_rng_index(vec, 140) + assert idx == 14 + + +def test_get_target_elevations(): + test_radar = pyart.testing.make_target_radar() + test_radar.elevation["data"][0] = 2 + target_elevations, el_tol = pyart.retrieve.qvp.get_target_elevations(test_radar) + assert target_elevations[-1] == 2 + assert np.all(target_elevations[0:-1] == 0.75) + assert el_tol == 0 + + +def test_find_nearest_gate(test_radar): + ind_ray, ind_rng, azi, rng = pyart.retrieve.qvp.find_nearest_gate( + test_radar, 36.4, -97.4 + ) + + assert ind_ray == 141.0 + assert ind_rng == 145.0 + assert azi == 141.0 + assert rng == 14514.514 + + +def test_find_neighbour_gates(test_radar): + inds_ray, inds_rng = pyart.retrieve.qvp.find_neighbour_gates( + test_radar, 141, 14514, delta_azi=3, delta_rng=200 + ) + + assert np.all(inds_ray == [139, 140, 141, 142, 143]) + assert np.all(inds_rng == [143, 144, 145, 146]) + + +def test_find_ang_index(): + ang_vec = np.arange(0, 360) + idx = pyart.retrieve.qvp.find_ang_index(ang_vec, 36) + assert idx == 36 + + +def test__create_qvp_object(test_radar): + qvp = pyart.retrieve.qvp._create_qvp_object( + test_radar, ["reflectivity"], hmax=5000, hres=500 + ) + assert np.all(qvp.range["data"] == np.arange(250, 5000, 500)) + assert "reflectivity" in qvp.fields + assert qvp.sweep_mode["data"] == ["qvp"] + assert qvp.fixed_angle["data"] == 10 + + +def test__create_along_coord_object(test_radar): + acoord = pyart.retrieve.qvp._create_along_coord_object( + test_radar, ["reflectivity"], np.arange(0, 10000, 100), 10, "ALONG_AZI" + ) + + assert np.all(acoord.range["data"] == np.arange(0, 10000, 100)) + assert acoord.sweep_mode["data"] == ["ALONG_AZI"] + + +def test__update_qvp_metadata(test_radar): + qvp = pyart.retrieve.qvp._create_qvp_object( + test_radar, ["reflectivity"], hmax=5000, hres=500 + ) + new_time = datetime.datetime(2024, 6, 10) + qvp = pyart.retrieve.qvp._update_qvp_metadata(qvp, new_time, 10, 45) + start_time = num2date(0, qvp.time["units"], qvp.time["calendar"]) + time_offset = (new_time - start_time).total_seconds() + + assert np.all(qvp.gate_longitude["data"] == 10) + assert np.all(qvp.gate_latitude["data"] == 45) + assert qvp.time["data"] == time_offset + + +def test__update_along_coord_metadata(test_radar): + acoord = pyart.retrieve.qvp._create_along_coord_object( + test_radar, ["reflectivity"], np.arange(0, 10000, 100), 10, "ALONG_AZI" + ) + new_time = datetime.datetime(2024, 6, 10) + acoord = pyart.retrieve.qvp._update_along_coord_metadata(acoord, new_time, 5, 200) + + start_time = num2date(0, acoord.time["units"], acoord.time["calendar"]) + time_offset = (new_time - start_time).total_seconds() + + assert acoord.time["data"] == time_offset + assert acoord.elevation["data"] == [5] + assert acoord.azimuth["data"] == [200] diff --git a/tests/util/test_circular_stats.py b/tests/util/test_circular_stats.py index 79063ec3848..e5574c719f8 100644 --- a/tests/util/test_circular_stats.py +++ b/tests/util/test_circular_stats.py @@ -6,6 +6,36 @@ from pyart.util import circular_stats +def test_compute_directional_stats(): + field_1d = np.ma.array([1, 1, 3, 4, 5]) + field_1d.mask = [True, False, False, False, False] + field = np.tile(field_1d, (10, 1)) + + mean, nvalid = circular_stats.compute_directional_stats(field, axis=1) + median, nvalid = circular_stats.compute_directional_stats( + field, axis=1, avg_type="median" + ) + + assert np.all(mean == 3.25) + assert np.all(median == 3.5) + assert np.all(nvalid == 4) + + # Test with larger nb of nvalid_min + mean, nvalid = circular_stats.compute_directional_stats(field, axis=1, nvalid_min=5) + assert np.all(mean.mask) + + # Test other direction + mean, nvalid = circular_stats.compute_directional_stats(field, axis=0) + median, nvalid = circular_stats.compute_directional_stats( + field, axis=0, avg_type="median" + ) + + ref = np.ma.array([np.nan, 1, 3, 4, 5], mask=[1, 0, 0, 0, 0]) + assert np.all(mean == ref) + assert np.all(median == ref) + assert np.all(nvalid == [0, 10, 10, 10, 10]) + + def test_mean_of_two_angles(): mean = circular_stats.mean_of_two_angles(np.pi / 4 + 0.2, 7 * np.pi / 4) assert_almost_equal(mean, 0.10, 2) diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index f6dd40fdb1d..0a7afbc61b5 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -105,6 +105,16 @@ def test_subset_radar(): assert list(radarcut.fields) == ["f1"] +def test_ma_broadcast_to(): + buf = np.ma.zeros(5) + buf.mask = [1, 1, 0, 0, 0] + buf_broad = pyart.util.radar_utils.ma_broadcast_to(buf, (10, 5)) + assert buf_broad.shape == (10, 5) + assert buf_broad.mask.shape == (10, 5) + expected_mask = np.tile(np.array([True, True, False, False, False]), [10, 1]) + assert np.all(expected_mask == buf_broad.mask) + + # read in example file radar = pyart.io.read_nexrad_archive(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) From 23ce76744796818a0580b051817ef3341f31aea7 Mon Sep 17 00:00:00 2001 From: Max Grover Date: Thu, 29 Aug 2024 07:41:45 -0500 Subject: [PATCH 4/7] ENH: Improve performance for xradar interface (#1621) --- pyart/xradar/accessor.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/pyart/xradar/accessor.py b/pyart/xradar/accessor.py index e79a1bcb150..e64d30b8021 100644 --- a/pyart/xradar/accessor.py +++ b/pyart/xradar/accessor.py @@ -275,17 +275,26 @@ def __init__(self, xradar, default_sweep="sweep_0", scan_type=None): self.nsweeps = len(self.sweep_group_names) self.combined_sweeps = self._combine_sweeps() self.fields = self._find_fields(self.combined_sweeps) + + # Check to see if the time is multidimensional - if it is, collapse it + if len(self.combined_sweeps.time.dims) > 1: + time = self.combined_sweeps.time.isel(range=0) + else: + time = self.combined_sweeps.time self.time = dict( - data=(self.combined_sweeps.time - self.combined_sweeps.time.min()).astype( - "int64" - ) - / 1e9, + data=(time - time.min()).astype("int64").values / 1e9, units=f"seconds since {pd.to_datetime(self.combined_sweeps.time.min().values).strftime('%Y-%m-%d %H:%M:%S.0')}", calendar="gregorian", ) self.range = dict(data=self.combined_sweeps.range.values) self.azimuth = dict(data=self.combined_sweeps.azimuth.values) self.elevation = dict(data=self.combined_sweeps.elevation.values) + # Check to see if the time is multidimensional - if it is, collapse it + if len(self.combined_sweeps.sweep_fixed_angle.dims) > 1: + self.combined_sweeps["sweep_fixed_angle"] = ( + "sweep_number", + np.unique(self.combined_sweeps.sweep_fixed_angle), + ) self.fixed_angle = dict(data=self.combined_sweeps.sweep_fixed_angle.values) self.antenna_transition = None self.latitude = dict( @@ -590,6 +599,10 @@ def _combine_sweeps(self): # Stack the sweep number and azimuth together stacked = merged.stack(gates=["sweep_number", "azimuth"]).transpose() + # Select the valid azimuths + good_azimuths = stacked.time.dropna("gates", how="all").gates + stacked = stacked.sel(gates=good_azimuths) + # Drop the missing gates cleaned = stacked.where(stacked.time == stacked.time.dropna("gates")) From c810816821589bed7956b4e86c50968546072052 Mon Sep 17 00:00:00 2001 From: Max Grover Date: Thu, 29 Aug 2024 12:50:20 -0500 Subject: [PATCH 5/7] ADD: Add to_radar() method to xradar objects (#1622) * ADD: Add to_radar() method * FIX: Fix the examples to use new to_radar() interface --- examples/xradar/plot_dealias_xradar.py | 2 +- examples/xradar/plot_grid_xradar.py | 2 +- examples/xradar/plot_xradar.py | 2 +- pyart/xradar/accessor.py | 28 +++++++++++++++++++++++++- tests/xradar/test_accessor.py | 13 ++++++++++++ 5 files changed, 43 insertions(+), 4 deletions(-) diff --git a/examples/xradar/plot_dealias_xradar.py b/examples/xradar/plot_dealias_xradar.py index 9e78a49b10d..5ef3ce1327f 100644 --- a/examples/xradar/plot_dealias_xradar.py +++ b/examples/xradar/plot_dealias_xradar.py @@ -23,7 +23,7 @@ tree = xd.io.open_cfradial1_datatree(filename) # Give the tree Py-ART radar methods -radar = pyart.xradar.Xradar(tree) +radar = tree.pyart.to_radar() # Determine the nyquist velocity using the maximum radial velocity from the first sweep nyq = radar["sweep_0"]["mean_doppler_velocity"].max().values diff --git a/examples/xradar/plot_grid_xradar.py b/examples/xradar/plot_grid_xradar.py index 8a6e00f0c85..085fc4d3f45 100644 --- a/examples/xradar/plot_grid_xradar.py +++ b/examples/xradar/plot_grid_xradar.py @@ -21,7 +21,7 @@ tree = xd.io.open_cfradial1_datatree(filename) # Give the tree Py-ART radar methods -radar = pyart.xradar.Xradar(tree) +radar = tree.pyart.to_radar() # Grid using 11 vertical levels, and 101 horizontal grid cells at a resolution on 1 km grid = pyart.map.grid_from_radars( diff --git a/examples/xradar/plot_xradar.py b/examples/xradar/plot_xradar.py index 708adffe725..2fc7ba439ab 100644 --- a/examples/xradar/plot_xradar.py +++ b/examples/xradar/plot_xradar.py @@ -21,7 +21,7 @@ tree = xd.io.open_cfradial1_datatree(filename) # Give the tree Py-ART radar methods -radar = pyart.xradar.Xradar(tree) +radar = tree.pyart.to_radar() # Plot the Reflectivity Field (corrected_reflectivity_horizontal) display = pyart.graph.RadarMapDisplay(radar) diff --git a/pyart/xradar/accessor.py b/pyart/xradar/accessor.py index e64d30b8021..f0749687132 100644 --- a/pyart/xradar/accessor.py +++ b/pyart/xradar/accessor.py @@ -5,12 +5,14 @@ import copy +import datatree import numpy as np import pandas as pd from datatree import DataTree, formatting, formatting_html from datatree.treenode import NodePath from xarray import DataArray, Dataset, concat from xarray.core import utils +from xradar.accessors import XradarAccessor from xradar.util import get_sweep_keys from ..config import get_metadata @@ -591,7 +593,11 @@ def _combine_sweeps(self): # Loop through and extract the different datasets ds_list = [] for sweep in self.sweep_group_names: - ds_list.append(self.xradar[sweep].ds.drop_duplicates("azimuth")) + ds_list.append( + self.xradar[sweep] + .ds.drop_duplicates("azimuth") + .set_coords("sweep_number") + ) # Merge based on the sweep number merged = concat(ds_list, dim="sweep_number") @@ -800,3 +806,23 @@ def _point_altitude_data(): return grid.origin_altitude["data"][0] + grid.point_z["data"] return _point_altitude_data + + +@datatree.register_datatree_accessor("pyart") +class XradarDataTreeAccessor(XradarAccessor): + """Adds a number of pyart specific methods to datatree.DataTree objects.""" + + def to_radar(self, scan_type=None) -> DataTree: + """ + Add pyart radar object methods to the xradar datatree object + Parameters + ---------- + scan_type: string + Scan type (ppi, rhi, etc.) + Returns + ------- + dt: datatree.Datatree + Datatree including pyart.Radar methods + """ + dt = self.xarray_obj + return Xradar(dt, scan_type=scan_type) diff --git a/tests/xradar/test_accessor.py b/tests/xradar/test_accessor.py index 044df28f478..1ec9a2af5db 100644 --- a/tests/xradar/test_accessor.py +++ b/tests/xradar/test_accessor.py @@ -154,3 +154,16 @@ def test_grid_write_read(): assert grid1.radar_name["data"] == grid2.radar_name["data"] assert grid1.nradar == grid2.nradar grid2.ds.close() + + +def test_to_xradar_object(): + """Test the to_radar_object""" + dtree = xd.io.open_cfradial1_datatree( + filename, + optional=False, + ) + radar = dtree.pyart.to_radar() + reflectivity = radar.get_field(0, "DBZ") + assert reflectivity.shape == (480, 996) + assert radar.nsweeps == 9 + assert radar.ngates == 996 From 37ae458550ec7da70d063319bb8177048396fe6d Mon Sep 17 00:00:00 2001 From: Max Grover Date: Thu, 29 Aug 2024 16:03:49 -0500 Subject: [PATCH 6/7] DOC: Enable readthedocs (#1629) --- .readthedocs.yaml | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 .readthedocs.yaml diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000000..33eb10b46fc --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,7 @@ +version: 2 +conda: + environment: doc/environment.yml +build: + os: "ubuntu-20.04" + tools: + python: "mambaforge-4.10" From 2410559ecb78d763876b2e86dd7421a684e30576 Mon Sep 17 00:00:00 2001 From: Zach Sherman <19153455+zssherman@users.noreply.github.com> Date: Thu, 29 Aug 2024 16:08:32 -0500 Subject: [PATCH 7/7] ADD: Adding cloud mask code. (#1628) * ADD: Adding cloud mask code. * MNT: Update change suggested by max. Co-authored-by: Max Grover --------- Co-authored-by: Max Grover --- pyart/correct/__init__.py | 7 +- pyart/correct/bias_and_noise.py | 281 +++++++++++++++++++++++++++++ tests/correct/test_correct_bias.py | 111 +++++++++++- 3 files changed, 397 insertions(+), 2 deletions(-) diff --git a/pyart/correct/__init__.py b/pyart/correct/__init__.py index caa3ac75b9e..6a9812ec775 100644 --- a/pyart/correct/__init__.py +++ b/pyart/correct/__init__.py @@ -8,7 +8,12 @@ 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 correct_bias, correct_noise_rhohv # noqa +from .bias_and_noise import calc_cloud_mask, calc_noise_floor, correct_bias # noqa +from .bias_and_noise import ( + correct_noise_rhohv, # noqa + cloud_threshold, # noqa + range_correction, # noqa +) # noqa from .dealias import dealias_fourdd # noqa from .despeckle import despeckle_field, find_objects # noqa from .phase_proc import phase_proc_lp, phase_proc_lp_gf # noqa diff --git a/pyart/correct/bias_and_noise.py b/pyart/correct/bias_and_noise.py index 03bf7a31618..295a6689faa 100755 --- a/pyart/correct/bias_and_noise.py +++ b/pyart/correct/bias_and_noise.py @@ -2,10 +2,15 @@ Corrects polarimetric variables for noise """ +import warnings +import dask import numpy as np +import pint +from scipy import signal from ..config import get_field_name, get_metadata +from ..core.radar import Radar def correct_noise_rhohv( @@ -137,3 +142,279 @@ def correct_bias(radar, bias=0.0, field_name=None): corr_field["data"] = corr_field_data return corr_field + + +def calc_cloud_mask( + radar, + field, + height=None, + noise_threshold=-45.0, + threshold_offset=5.0, + counts_threshold=12, +): + """ + Primary function for calculating the cloud mask. + + Parameters + ---------- + radar : Radar + Py-ART Radar object. + field : string + Reflectivity field name to calculate. + height : string + Height name to use for calculations. + noise_threshold : float + Threshold value used for noise detection. Greater than this value. + threshold_offset : float + Threshold offset value used for noise detection + counts_threshold : int + Threshold of counts used to determine mask. Greater than or equal to this value. + + Returns + ------- + radar : Radar + Returns an updated Radar object with cloud mask fields. + + References + ---------- + Kollias, P., I. Jo, P. Borque, A. Tatarevic, K. Lamer, N. Bharadwaj, K. Widener, + K. Johnson, and E.E. Clothiaux, 2014: Scanning ARM Cloud Radars. Part II: Data + Quality Control and Processing. J. Atmos. Oceanic Technol., 31, 583–598, + https://doi.org/10.1175/JTECH-D-13-00045.1 + + """ + + if not isinstance(radar, Radar): + raise ValueError("Please use a valid Py-ART Radar object") + + if not isinstance(field, str): + raise ValueError("Please specify a valid field name.") + + noise = calc_noise_floor(radar, field, height=height) + + noise_thresh = ( + np.nanmin( + np.vstack( + [ + noise, + np.full(np.shape(radar.fields[field]["data"])[0], noise_threshold), + ] + ), + axis=0, + ) + + threshold_offset + ) + + data = range_correction(radar, field, height=height) + + task = [] + for i in range(np.shape(data)[0]): + task.append(dask.delayed(_first_mask)(data[i, :], noise_thresh[i])) + + result = dask.compute(task) + mask1 = np.array(result[0]) + + counts = signal.convolve2d(mask1, np.ones((4, 4), dtype=int), mode="same") + mask2 = np.zeros_like(data, dtype=np.int16) + mask2[counts >= counts_threshold] = 1 + + cloud_mask_1 = {} + cloud_mask_1["long_name"] = "Cloud mask 1 (linear profile)" + cloud_mask_1["units"] = "1" + cloud_mask_1["comment"] = ( + "The mask is calculated with a " "linear mask along each time profile." + ) + cloud_mask_1["flag_values"] = [0, 1] + cloud_mask_1["flag_meanings"] = ["no_cloud", "cloud"] + cloud_mask_1["data"] = mask1 + + cloud_mask_2 = {} + cloud_mask_2["long_name"] = "Cloud mask 2 (2D box)" + cloud_mask_2["units"] = "1" + cloud_mask_2["comment"] = "The mask uses a 2D box to " "filter out noise." + cloud_mask_2["flag_values"] = [0, 1] + cloud_mask_2["flag_meanings"] = ["no_cloud", "cloud"] + cloud_mask_2["data"] = mask2 + + radar.add_field("cloud_mask_1", cloud_mask_1, replace_existing=True) + radar.add_field("cloud_mask_2", cloud_mask_2, replace_existing=True) + + return radar + + +def calc_noise_floor(radar, field, height): + """ + Calculation for getting the noise floor + + Parameters + ---------- + radar : Radar + Py-ART Radar object containing data. + field : string + Reflectivity field name to correct. + height : string + Height name to use in correction. + + Returns + ------- + noise : array + Returns the noise floor value for each time sample. + + References + ---------- + Kollias, P., I. Jo, P. Borque, A. Tatarevic, K. Lamer, N. Bharadwaj, K. Widener, + K. Johnson, and E.E. Clothiaux, 2014: Scanning ARM Cloud Radars. Part II: Data + Quality Control and Processing. J. Atmos. Oceanic Technol., 31, 583–598, + https://doi.org/10.1175/JTECH-D-13-00045.1 + + """ + + if not isinstance(radar, Radar): + raise ValueError("Please use a valid Py-ART Radar object.") + + # Range correct data and return the array from the Radar object + data = range_correction(radar, field, height=height) + + # Pass each timestep into task list to calculate cloud threshhold + # with a delayed dask process + task = [dask.delayed(cloud_threshold)(row) for row in data] + + # Perform dask computation + noise = dask.compute(*task) + + # Convert returned dask tuple into numpy array + noise = np.array(noise, dtype=float) + + return noise + + +def cloud_threshold(data, n_avg=1.0, nffts=None): + """ + Calculates the noise floor from a cloud threshold. + + Parameters + ---------- + data : array + Numpy array + n_avg : float + Number of points to average over + nffts : int + Number of heights to iterate over. If None will use the size of data. + + Returns + ------- + result : numpy scalar float + Returns the noise floor value for each time sample. + + References + ---------- + Kollias, P., I. Jo, P. Borque, A. Tatarevic, K. Lamer, N. Bharadwaj, K. Widener, + K. Johnson, and E.E. Clothiaux, 2014: Scanning ARM Cloud Radars. Part II: Data + Quality Control and Processing. J. Atmos. Oceanic Technol., 31, 583–598, + https://doi.org/10.1175/JTECH-D-13-00045.1 + + """ + + if nffts is None: + nffts = data.size + + data = 10.0 ** (data / 10.0) + data = np.sort(data) + + nthld = 10.0**-10.0 + dsum = 0.0 + sumSq = 0.0 + n = 0.0 + numNs = [] + sqrt_n_avg = np.sqrt(n_avg) + for i in range(nffts): + if data[i] > nthld: + dsum += data[i] + sumSq += data[i] ** 2.0 + n += 1.0 + a3 = dsum * dsum + a1 = sqrt_n_avg * (n * sumSq - a3) + if n > nffts / 4.0: + if a1 <= a3: + sumNs = dsum + numNs = [n] + else: + sumNs = dsum + numNs = [n] + + if len(numNs) > 0: + n_mean = sumNs / numNs[0] + else: + n_mean = np.nan + + if n_mean == 0.0: + result = np.nan + else: + result = 10.0 * np.log10(n_mean) + + return result + + +def range_correction(radar, field, height): + """ + Corrects reflectivity for range to help get the + correct noise floor values + + Parameters + ---------- + radar : Radar + Py-ART Radar object containing data. + field : string + Reflectivity field name to correct. + height : string + Height name to use in correction. + + Returns + ------- + data : array + Returns a range corrected array matching reflectivity field. + + """ + + if not isinstance(radar, Radar): + raise ValueError("Please use a valid Py-ART Radar object.") + + try: + height_units = getattr(radar, height)["units"] + except KeyError: + warnings.warn( + f"Height '{height} does not have units attribute. " + "Assuming units are meters." + ) + height_units = "m" + + height = getattr(radar, height)["data"] + desired_unit = "m" + if height_units is not desired_unit: + if isinstance(height, np.ma.MaskedArray): + height = height.filled(np.nan) + unit_registry = pint.UnitRegistry() + height = height * unit_registry.parse_expression(height_units) + height = height.to(desired_unit) + height = height.magnitude + + data = radar.fields[field]["data"] + + if isinstance(data, np.ma.MaskedArray) and not data.mask: + data = data.data + elif isinstance(data, np.ma.MaskedArray) and data.mask: + data = data.filled(np.nan) + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", category=RuntimeWarning, message=".*divide by zero encountered.*" + ) + data = data - 20.0 * np.log10(height / 1000.0) + + return data + + +def _first_mask(data, noise_threshold): + mask = np.zeros_like(data, dtype=np.int16) + mask[data > noise_threshold] = 1 + return mask diff --git a/tests/correct/test_correct_bias.py b/tests/correct/test_correct_bias.py index 9ab077749e1..3be079411f2 100644 --- a/tests/correct/test_correct_bias.py +++ b/tests/correct/test_correct_bias.py @@ -1,11 +1,16 @@ """ Unit tests for bias and noise module. """ +import dask +import numpy as np import pytest -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_almost_equal, assert_array_equal +from open_radar_data import DATASETS import pyart radar = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) +kazr_file = DATASETS.fetch("sgpkazrgeC1.a1.20190529.000002.cdf") +radar_kazr = pyart.aux_io.read_kazr(kazr_file) def test_correct_bias(): @@ -24,3 +29,107 @@ def test_correct_bias(): bias = 0 field_name = "foo" pytest.raises(KeyError, pyart.correct.correct_bias, radar, bias, field_name) + + +def test_calc_noise_floor(): + expected = [-46.25460013, -46.48371626, -46.3314618, -46.82639895, -46.76403711] + + result = pyart.correct.calc_noise_floor(radar_kazr, "reflectivity_copol", "range") + + assert_almost_equal(result[0:5], expected, decimal=3) + + bad_radar = "foo" + pytest.raises( + ValueError, + pyart.correct.calc_noise_floor, + bad_radar, + "reflectivity_copol", + "range", + ) + + +def test_range_correction(): + expected_first = [-40.92386, -39.07553, -29.10794, -26.25786, -26.27254] + expected_last = [-43.42107, -23.74393, -22.77576, -16.700004, -19.68343] + + result = pyart.correct.range_correction(radar_kazr, "reflectivity_copol", "range") + assert_almost_equal(result[0][0:5], expected_first, decimal=3) + assert_almost_equal(result[-1][0:5], expected_last, decimal=3) + + bad_radar = "foo" + pytest.raises( + ValueError, + pyart.correct.range_correction, + bad_radar, + "reflectivity_copol", + "range", + ) + + radar_kazr.range.pop("units") + pytest.warns( + UserWarning, + pyart.correct.range_correction, + radar_kazr, + "reflectivity_copol", + "range", + ) + + +def test_cloud_threshold(): + expected = [-46.25460, -46.48371, -46.33146, -46.82639, -46.76403] + + result = pyart.correct.range_correction(radar_kazr, "reflectivity_copol", "range") + + task = [dask.delayed(pyart.correct.cloud_threshold)(row) for row in result] + noise = dask.compute(*task) + noise = np.array(noise, dtype=float) + + assert_almost_equal(noise[0:5], expected, decimal=3) + + +def test_calc_cloud_mask(): + expected_cloud_1_first = [1, 1, 1, 1, 1] + expected_cloud_1_last = [0, 1, 1, 1, 1] + expected_cloud_2_first = [0, 0, 0, 0, 0] + expected_cloud_2_last = [0, 0, 0, 1, 1] + + radar_mask = pyart.correct.calc_cloud_mask( + radar_kazr, "reflectivity_copol", "range" + ) + + assert "cloud_mask_1" in radar_mask.fields.keys() + assert "cloud_mask_2" in radar_mask.fields.keys() + assert_array_equal( + expected_cloud_1_first, radar_mask.fields["cloud_mask_1"]["data"][0][0:5] + ) + assert_array_equal( + expected_cloud_1_last, radar_mask.fields["cloud_mask_1"]["data"][-1][0:5] + ) + assert_array_equal( + expected_cloud_2_first, radar_mask.fields["cloud_mask_2"]["data"][0][0:5] + ) + assert_array_equal( + expected_cloud_2_last, radar_mask.fields["cloud_mask_2"]["data"][-1][0:5] + ) + + assert ( + radar_mask.fields["cloud_mask_1"]["long_name"] + == "Cloud mask 1 (linear profile)" + ) + assert radar_mask.fields["cloud_mask_2"]["long_name"] == "Cloud mask 2 (2D box)" + + assert radar_mask.fields["cloud_mask_2"]["units"] == "1" + + assert radar_mask.fields["cloud_mask_2"]["flag_meanings"] == ["no_cloud", "cloud"] + + bad_radar = "foo" + pytest.raises( + ValueError, + pyart.correct.calc_cloud_mask, + bad_radar, + "reflectivity_copol", + "range", + ) + + bad_field = 42 + pytest.raises(ValueError, pyart.correct.calc_cloud_mask, radar, bad_field, "range")