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 <[email protected]>

* Update pyart/util/radar_utils.py

Co-authored-by: Max Grover <[email protected]>

---------

Co-authored-by: mgrover1 <[email protected]>

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,
     )
 

pyart/util/__init__.py

@@ -29,6 +29,7 @@
     image_mute_radar,
     is_vpt,
     join_radar,
+    determine_sweeps,
     subset_radar,
     to_vpt,
 )

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,

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

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))}