From 90a9b93ed63dbe544a2b5b4c7e71d39d30231f6b Mon Sep 17 00:00:00 2001 From: Israel Silber <55999320+isilber@users.noreply.github.com> Date: Thu, 7 Dec 2023 12:33:52 -0800 Subject: [PATCH] ENH: Add method to automatically determine radar sweep indices (#1493) * FIX: change deprecated np.str in arm_vpt to str * Update default wind_size in calculate_velocity_texture from 4 to 3 Using an odd number prevents single index offsets in placement in case of an even size * ADD: optional rectangular window dims for velocity texture analysis. Also, modify the default window size to 3 instead of the previous value of 4 (even number resulting in a potential offset). * FIX: nyquist velocity dtype in radar obj loaded with read_kazr * pep8 * line reformatting per black * TST: unit tests for rectangular texture window * STY: linting (black) * FIX: Add in linting fixes * FIX: metadata in core.radar.get_elevation * FIX: outdated text in masking_data_with_gatefilters.ipynb * ADD: method to automatically determine sweep number and sweep start and end indices * add comment to 'determine_sweeps' method * FIX: indexing bug in 'determine_sweeps' method * ADD: unit tests for 'determine_sweeps' * linting * FIX: nsweeps wasn't updated when calling 'determine_sweeps'; tests were updated accordingly * line reformatting per black * FIX: update multiple radar object attributes to ensure its full utilization with Py-ART after 'determine_sweeps' is called * line reformatting per black * linting * ENH: add antenna_transition to the testing module's sample objects * ENH: update unit tests to reflect antenna_transition allocation in sample_objects * Update pyart/util/radar_utils.py Co-authored-by: Max Grover * Update pyart/util/radar_utils.py Co-authored-by: Max Grover --------- Co-authored-by: mgrover1 --- pyart/testing/sample_objects.py | 3 + pyart/util/__init__.py | 1 + pyart/util/radar_utils.py | 99 +++++++++++++++++++++++++++++++++ tests/core/test_radar.py | 6 +- tests/util/test_radar_utils.py | 24 ++++++++ 5 files changed, 130 insertions(+), 3 deletions(-) diff --git a/pyart/testing/sample_objects.py b/pyart/testing/sample_objects.py index 86e6e0bbed..fb4f1f0212 100644 --- a/pyart/testing/sample_objects.py +++ b/pyart/testing/sample_objects.py @@ -56,6 +56,7 @@ def make_empty_ppi_radar(ngates, rays_per_sweep, nsweeps): sweep_end_ray_index = get_metadata("sweep_end_ray_index") azimuth = get_metadata("azimuth") elevation = get_metadata("elevation") + antenna_transition = get_metadata("antenna_transition") fields = {} scan_type = "ppi" @@ -76,6 +77,7 @@ def make_empty_ppi_radar(ngates, rays_per_sweep, nsweeps): sweep_end_ray_index["data"] = np.arange( rays_per_sweep - 1, nrays, rays_per_sweep, dtype="int32" ) + antenna_transition["data"] = np.zeros(nrays) azimuth["data"] = np.arange(nrays, dtype="float32") elevation["data"] = np.array([0.75] * nrays, dtype="float32") @@ -96,6 +98,7 @@ def make_empty_ppi_radar(ngates, rays_per_sweep, nsweeps): sweep_end_ray_index, azimuth, elevation, + antenna_transition=antenna_transition, instrument_parameters=None, ) diff --git a/pyart/util/__init__.py b/pyart/util/__init__.py index 1f060fc8ad..ee8e1ae47c 100644 --- a/pyart/util/__init__.py +++ b/pyart/util/__init__.py @@ -29,6 +29,7 @@ image_mute_radar, is_vpt, join_radar, + determine_sweeps, subset_radar, to_vpt, ) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 1a0246f1d7..c2e7c405b9 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -6,6 +6,7 @@ import copy import numpy as np +import numpy.ma as ma from netCDF4 import date2num from ..config import get_fillvalue @@ -104,6 +105,104 @@ def to_vpt(radar, single_scan=True): return +def determine_sweeps(radar, max_offset=0.1, running_win_dt=5.0, deg_rng=(-5.0, 360.0)): + """ + Determine the number of sweeps using elevation data (PPI scans) or azimuth + data (RHI scans) and update the input radar object + + Parameters + ---------- + radar : Radar object + The radar object containing the data. + max_offset : float + Maximum elevation offset (if is_ppi is True) or azimuth offset (if + is_ppi is False) allowed to determine sweeps. + running_win_dt: float + running window period (in seconds) used to determine elevation or + azimuth shifts. + Note: set wisely: the method assumes that a single sweep is longer than this + parameter. + deg_rng: float + angle range (azimuth or elevation) to consider for calculations. + Assuming azimuths between 0 to 360, this should be equal to (0., 360.), but + given that there could be ppi scan strategies at negative elevations, + one might consider a negative values (current default), or , for example, + -180 to 180 if the azimuth range goes from -180 to 180. + + """ + # set fixed and variable coordinates depending on scan type + # ====================== + if "rhi" in radar.scan_type.lower(): + var_array = radar.elevation["data"] + fix_array = radar.azimuth["data"] + else: # ppi or vpt + var_array = radar.azimuth["data"] + fix_array = radar.elevation["data"] + + # set bins and parameters and allocate lists + # ====================== + angle_bins = np.arange( + deg_rng[0] - max_offset, deg_rng[1] + max_offset + 1e-10, max_offset * 2.0 + ) + sample_dt = np.nanmean(np.diff(radar.time["data"])) + win_size = int(np.ceil(running_win_dt / sample_dt)) + if win_size < 2: + raise ValueError( + "Window size <= 1; consider decreasing the value of running_win_dt" + ) + sweep_start_index, sweep_end_index = [], [] + in_sweep = False # determine if sweep is underway in current index + + # Loop through coordinate data and detect sweep edges + # ====================== + t = 0 + while t < radar.time["data"].size - win_size + 1: + var_win = var_array[t : t + win_size] + fix_win = fix_array[t : t + win_size] + idle_sweep = np.diff(var_win) == 0 + if idle_sweep[0]: # sweep did not start + t += 1 + continue + bincounts, _ = np.histogram(fix_win, bins=angle_bins) + moving_radar = np.sum(bincounts > 0) > 1 # radar transition to a new sweep + if in_sweep: + if t == radar.time["data"].size - win_size: + sweep_end_index.append(radar.time["data"].size - 1) + elif moving_radar: + in_sweep = False + sweep_end_index.append(t + win_size - 2) + t += win_size - 2 + elif np.all(~idle_sweep) & ~moving_radar: + in_sweep = True + sweep_start_index.append(t) + t += 1 + sweep_number = np.arange(len(sweep_start_index)) + + # Update radar object + # ====================== + radar.sweep_start_ray_index["data"] = ma.array(sweep_start_index, dtype="int32") + radar.sweep_end_ray_index["data"] = ma.array(sweep_end_index, dtype="int32") + radar.sweep_number["data"] = ma.array(sweep_number, dtype="int32") + fixed_angle = [ + np.mean(fix_array[si : ei + 1]) + for si, ei in zip( + radar.sweep_start_ray_index["data"], radar.sweep_end_ray_index["data"] + ) + ] + radar.fixed_angle["data"] = ma.array(fixed_angle, dtype="float32") + radar.nsweeps = len(sweep_number) + transition = np.zeros(radar.nrays) + for i in range(radar.nsweeps): + if i == 0: + transition[: sweep_start_index[i]] = 1 + else: + transition[sweep_end_index[i - 1] : sweep_start_index[i]] = 1 + radar.antenna_transition["data"] = ma.array(transition, dtype="int32") + bstr_entry = np.array([x for x in f"{radar.scan_type:<22}"], dtype="|S1") + radar.sweep_mode = ma.array(np.tile(bstr_entry[np.newaxis, :], (radar.nsweeps, 1))) + return + + def subset_radar( radar, field_names, diff --git a/tests/core/test_radar.py b/tests/core/test_radar.py index 00e3125296..0ed2d1f875 100644 --- a/tests/core/test_radar.py +++ b/tests/core/test_radar.py @@ -156,7 +156,7 @@ def test_get_gate_x_y_z_transitions(): radar.azimuth["data"][:] = [0, 90, 180, 270, 0, 90, 180, 270] radar.elevation["data"][:] = [0, 0, 0, 0, 10, 10, 10, 10] radar.range["data"][:] = [5, 15, 25, 35, 45] - radar.antenna_transition = {"data": np.array([0, 0, 1, 0, 0, 0, 0, 0])} + radar.antenna_transition["data"][:] = [0, 0, 1, 0, 0, 0, 0, 0] gate_x, gate_y, gate_z = radar.get_gate_x_y_z(0, filter_transitions=True) assert gate_x.shape == (3, 5) @@ -213,7 +213,7 @@ def test_get_gate_lat_lon_alt(): def test_get_gate_lat_lon_alt_transitions(): radar = pyart.testing.make_empty_ppi_radar(5, 4, 2) - radar.antenna_transition = {"data": np.array([0, 0, 1, 0, 0, 0, 0, 0])} + radar.antenna_transition["data"][:] = [0, 0, 1, 0, 0, 0, 0, 0] lat, lon, alt = radar.get_gate_lat_lon_alt(0, filter_transitions=True) assert lat.shape == (3, 5) assert_allclose(lat[0], [36.5, 36.502243, 36.50449, 36.506744, 36.50899], atol=1e-3) @@ -401,7 +401,7 @@ def test_extract_sweeps(): assert eradar.azimuth["data"].shape == (720,) assert eradar.elevation["data"].shape == (720,) assert eradar.scan_rate is None - assert eradar.antenna_transition is None + assert eradar.antenna_transition["data"].shape == (720,) assert eradar.instrument_parameters is None assert eradar.radar_calibration is None diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index 442b2cc0ab..f6dd40fdb1 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -50,6 +50,30 @@ def test_to_vpt(): assert len(radar.instrument_parameters["prt_mode"]["data"]) == 108 +def test_determine_sweeps(): + # ppi + radar = pyart.testing.make_empty_ppi_radar(10, 36, 3) + radar.elevation["data"] = radar.elevation["data"] * np.ceil( + np.arange(1, 36 * 3 + 1) / 36 + ) + pyart.util.determine_sweeps(radar) + assert np.all(radar.sweep_end_ray_index["data"] == [35, 71, 107]) + assert np.all(radar.sweep_start_ray_index["data"] == [0, 36, 72]) + assert len(radar.sweep_number["data"]) == 3 + assert radar.nsweeps == 3 + + # rhi + radar = pyart.testing.make_empty_rhi_radar(10, 25, 5) + radar.azimuth["data"] = ( + radar.azimuth["data"] * np.ceil(np.arange(1, 25 * 5 + 1) / 25) * 25 + ) + pyart.util.determine_sweeps(radar) + assert np.all(radar.sweep_end_ray_index["data"] == [24, 49, 74, 99, 124]) + assert np.all(radar.sweep_start_ray_index["data"] == [0, 25, 50, 75, 100]) + assert len(radar.sweep_number["data"]) == 5 + assert radar.nsweeps == 5 + + def test_subset_radar(): radar = pyart.testing.make_empty_ppi_radar(10, 36, 3) field = {"data": np.ones((36 * 3, 10))}