From 860612892897370d327049d4b7e2a1cdffc4656e Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 16:36:38 +0000 Subject: [PATCH 01/30] FIX: change deprecated np.str in arm_vpt to str --- pyart/aux_io/arm_vpt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyart/aux_io/arm_vpt.py b/pyart/aux_io/arm_vpt.py index 0df3a38f86..94ee51c77b 100644 --- a/pyart/aux_io/arm_vpt.py +++ b/pyart/aux_io/arm_vpt.py @@ -118,7 +118,7 @@ def read_kazr( sweep_number["data"] = np.array([0], dtype=np.int32) sweep_mode = filemetadata("sweep_mode") - sweep_mode["data"] = np.array(["vertical_pointing"], dtype=np.str) + sweep_mode["data"] = np.array(["vertical_pointing"], dtype=str) fixed_angle = filemetadata("fixed_angle") fixed_angle["data"] = np.array([90.0], dtype=np.float32) @@ -168,7 +168,7 @@ def read_kazr( frequency["data"] = np.array([omega / 1e9], dtype=np.float32) prt_mode = filemetadata("prt_mode") - prt_mode["data"] = np.array(["fixed"], dtype=np.str) + prt_mode["data"] = np.array(["fixed"], dtype=str) prf = float(ncobj.pulse_repetition_frequency.split()[0]) prt = filemetadata("prt") From 1fef627a18b10fef4a02eb0e1acbc3b9737a7c54 Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 18:44:50 +0000 Subject: [PATCH 02/30] 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 --- pyart/retrieve/simple_moment_calculations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/retrieve/simple_moment_calculations.py b/pyart/retrieve/simple_moment_calculations.py index 9b7c8f5f7e..43c053ed94 100644 --- a/pyart/retrieve/simple_moment_calculations.py +++ b/pyart/retrieve/simple_moment_calculations.py @@ -240,7 +240,7 @@ def compute_cdr(radar, rhohv_field=None, zdr_field=None, cdr_field=None): def calculate_velocity_texture( - radar, vel_field=None, wind_size=4, nyq=None, check_nyq_uniform=True + radar, vel_field=None, wind_size=3, nyq=None, check_nyq_uniform=True ): """ Derive the texture of the velocity field. From 001596b80e77583781c68951187efa719a73ba9b Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 19:37:41 +0000 Subject: [PATCH 03/30] 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). --- pyart/retrieve/simple_moment_calculations.py | 14 ++++++++++---- pyart/util/sigmath.py | 16 +++++++++++----- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/pyart/retrieve/simple_moment_calculations.py b/pyart/retrieve/simple_moment_calculations.py index 43c053ed94..e9d8bdca53 100644 --- a/pyart/retrieve/simple_moment_calculations.py +++ b/pyart/retrieve/simple_moment_calculations.py @@ -252,9 +252,11 @@ def calculate_velocity_texture( vel_field : str, optional Name of the velocity field. A value of None will force Py-ART to automatically determine the name of the velocity field. - wind_size : int, optional - The size of the window to calculate texture from. The window is - defined to be a square of size wind_size by wind_size. + wind_size : int or 2-element tuple, optional + The size of the window to calculate texture from. + If an integer, the window is defined to be a square of size wind_size + by wind_size. If tuple, defines the m x n dimensions of the filter + window. nyq : float, optional The nyquist velocity of the radar. A value of None will force Py-ART to try and determine this automatically. @@ -274,6 +276,10 @@ def calculate_velocity_texture( if vel_field is None: vel_field = get_field_name("velocity") + # Set wind_size as a tuple if input is int + if isinstance(wind_size, int): + wind_size = (wind_size, wind_size) + # Allocate memory for texture field vel_texture = np.zeros(radar.fields[vel_field]["data"].shape) @@ -301,6 +307,6 @@ def calculate_velocity_texture( "texture_of_radial_velocity" + "_of_scatters_away_from_instrument" ) vel_texture_field["data"] = ndimage.median_filter( - vel_texture, size=(wind_size, wind_size) + vel_texture, size=wind_size ) return vel_texture_field diff --git a/pyart/util/sigmath.py b/pyart/util/sigmath.py index 10205f8ffb..232cf2a25b 100644 --- a/pyart/util/sigmath.py +++ b/pyart/util/sigmath.py @@ -18,9 +18,11 @@ def angular_texture_2d(image, N, interval): image : 2D array of floats The array containing the velocities in which to calculate texture from. - N : int - This is the window size for calculating texture. The texture will be - calculated from an N by N window centered around the gate. + N : int or 2-element tuple + If int, this is the window size for calculating texture. The + texture will be calculated from an N by N window centered + around the gate. If tuple N defines the m x n dimensions of + the window centered around the gate. interval : float The absolute value of the maximum velocity. In conversion to radial coordinates, pi will be defined to be interval @@ -33,6 +35,10 @@ def angular_texture_2d(image, N, interval): Texture of the radial velocity field. """ + # Set N as a tuple if input is int + if isinstance(N, int): + N = (N, N) + # transform distribution from original interval to [-pi, pi] interval_max = interval interval_min = -interval @@ -45,10 +51,10 @@ def angular_texture_2d(image, N, interval): y = np.sin(im) # Calculate convolution - kernel = np.ones((N, N)) + kernel = np.ones(N) xs = signal.convolve2d(x, kernel, mode="same", boundary="symm") ys = signal.convolve2d(y, kernel, mode="same", boundary="symm") - ns = N**2 + ns = np.prod(N) # Calculate norm over specified window xmean = xs / ns From 3ea3a31e05e4addcd64f0eb047dde4ac9d047e51 Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 20:44:24 +0000 Subject: [PATCH 04/30] FIX: nyquist velocity dtype in radar obj loaded with read_kazr --- pyart/aux_io/arm_vpt.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/aux_io/arm_vpt.py b/pyart/aux_io/arm_vpt.py index 94ee51c77b..341dce7bf3 100644 --- a/pyart/aux_io/arm_vpt.py +++ b/pyart/aux_io/arm_vpt.py @@ -176,7 +176,7 @@ def read_kazr( v_nq = float(ncobj.nyquist_velocity.split()[0]) nyquist_velocity = filemetadata("nyquist_velocity") - nyquist_velocity["data"] = (v_nq * np.ones(ncvars["time"].size, dtype=np.float32),) + nyquist_velocity["data"] = v_nq * np.ones(ncvars["time"].size, dtype=np.float32) samples = int(ncobj.num_spectral_averages) n_samples = filemetadata("n_samples") n_samples["data"] = samples * np.ones(ncvars["time"].size, dtype=np.int32) From 64a6b3cc15f4370290a1fd388a921d09105816c4 Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 21:17:30 +0000 Subject: [PATCH 05/30] pep8 --- pyart/retrieve/simple_moment_calculations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/retrieve/simple_moment_calculations.py b/pyart/retrieve/simple_moment_calculations.py index e9d8bdca53..fe1d1e0cd3 100644 --- a/pyart/retrieve/simple_moment_calculations.py +++ b/pyart/retrieve/simple_moment_calculations.py @@ -253,7 +253,7 @@ def calculate_velocity_texture( Name of the velocity field. A value of None will force Py-ART to automatically determine the name of the velocity field. wind_size : int or 2-element tuple, optional - The size of the window to calculate texture from. + The size of the window to calculate texture from. If an integer, the window is defined to be a square of size wind_size by wind_size. If tuple, defines the m x n dimensions of the filter window. From 6571b7876ad5d904e66d0e959b6d8907eb6a440f Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 21:32:28 +0000 Subject: [PATCH 06/30] line reformatting per black --- pyart/retrieve/simple_moment_calculations.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pyart/retrieve/simple_moment_calculations.py b/pyart/retrieve/simple_moment_calculations.py index fe1d1e0cd3..61beb2db04 100644 --- a/pyart/retrieve/simple_moment_calculations.py +++ b/pyart/retrieve/simple_moment_calculations.py @@ -306,7 +306,5 @@ def calculate_velocity_texture( vel_texture_field["standard_name"] = ( "texture_of_radial_velocity" + "_of_scatters_away_from_instrument" ) - vel_texture_field["data"] = ndimage.median_filter( - vel_texture, size=wind_size - ) + vel_texture_field["data"] = ndimage.median_filter(vel_texture, size=wind_size) return vel_texture_field From 4227bf7db58938044504f7d950282616365317a1 Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 16 Aug 2023 18:07:10 +0000 Subject: [PATCH 07/30] TST: unit tests for rectangular texture window --- .../retrieve/test_simple_moment_calculations.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/retrieve/test_simple_moment_calculations.py b/tests/retrieve/test_simple_moment_calculations.py index 40708f7aee..b80f97949d 100644 --- a/tests/retrieve/test_simple_moment_calculations.py +++ b/tests/retrieve/test_simple_moment_calculations.py @@ -75,6 +75,10 @@ def test_calculate_velocity_texture(): radar, vel_field, wind_size=4, nyq=10 ) assert np.all(texture_field["data"] == 0) + texture_field = pyart.retrieve.calculate_velocity_texture( + radar, vel_field, wind_size=(5, 7), nyq=10 + ) + assert np.all(texture_field["data"] == 0) # Test none parameters radar2 = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) @@ -86,3 +90,16 @@ def test_calculate_velocity_texture(): [0.000363, 0.000363, 0.000363, 0.000363, 0.000363], atol=1e-03, ) + + # Test rectangular filter option consistency with default scalar + radar2 = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) + radar2.fields['velocity']['data'] = \ + np.random.normal(scale=0.01, size=radar2.fields['velocity']['data'].shape) * \ + radar2.fields['velocity']['data'] # Generate value variability + texture_field2 = pyart.retrieve.calculate_velocity_texture( + radar2, vel_field='velocity', wind_size=3, nyq=10. + ) + texture_field3 = pyart.retrieve.calculate_velocity_texture( + radar2, vel_field='velocity', wind_size=(3, 3), nyq=10. + ) + assert_allclose(texture_field2["data"], texture_field3["data"], atol=1e-10) From 30ffeccc73903ad837c00d2a8bf9746603a631d8 Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 16 Aug 2023 18:16:12 +0000 Subject: [PATCH 08/30] STY: linting (black) --- tests/retrieve/test_simple_moment_calculations.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/retrieve/test_simple_moment_calculations.py b/tests/retrieve/test_simple_moment_calculations.py index b80f97949d..d7bb2a3a88 100644 --- a/tests/retrieve/test_simple_moment_calculations.py +++ b/tests/retrieve/test_simple_moment_calculations.py @@ -93,13 +93,14 @@ def test_calculate_velocity_texture(): # Test rectangular filter option consistency with default scalar radar2 = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) - radar2.fields['velocity']['data'] = \ - np.random.normal(scale=0.01, size=radar2.fields['velocity']['data'].shape) * \ - radar2.fields['velocity']['data'] # Generate value variability + radar2.fields['velocity']['data'] = ( + np.random.normal(scale=0.01, size=radar2.fields['velocity']['data'].shape) * + radar2.fields['velocity']['data'] + ) # Generate velocity field variability texture_field2 = pyart.retrieve.calculate_velocity_texture( - radar2, vel_field='velocity', wind_size=3, nyq=10. + radar2, vel_field='velocity', wind_size=3, nyq=10.0 ) texture_field3 = pyart.retrieve.calculate_velocity_texture( - radar2, vel_field='velocity', wind_size=(3, 3), nyq=10. + radar2, vel_field='velocity', wind_size=(3, 3), nyq=10.0 ) assert_allclose(texture_field2["data"], texture_field3["data"], atol=1e-10) From 57535ed819e5e2395f368cf91607459c33396359 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Wed, 16 Aug 2023 13:16:31 -0500 Subject: [PATCH 09/30] FIX: Add in linting fixes --- tests/retrieve/test_simple_moment_calculations.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/retrieve/test_simple_moment_calculations.py b/tests/retrieve/test_simple_moment_calculations.py index b80f97949d..b222a317cf 100644 --- a/tests/retrieve/test_simple_moment_calculations.py +++ b/tests/retrieve/test_simple_moment_calculations.py @@ -93,13 +93,14 @@ def test_calculate_velocity_texture(): # Test rectangular filter option consistency with default scalar radar2 = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) - radar2.fields['velocity']['data'] = \ - np.random.normal(scale=0.01, size=radar2.fields['velocity']['data'].shape) * \ - radar2.fields['velocity']['data'] # Generate value variability + radar2.fields["velocity"]["data"] = ( + np.random.normal(scale=0.01, size=radar2.fields["velocity"]["data"].shape) + * radar2.fields["velocity"]["data"] + ) # Generate value variability texture_field2 = pyart.retrieve.calculate_velocity_texture( - radar2, vel_field='velocity', wind_size=3, nyq=10. + radar2, vel_field="velocity", wind_size=3, nyq=10.0 ) texture_field3 = pyart.retrieve.calculate_velocity_texture( - radar2, vel_field='velocity', wind_size=(3, 3), nyq=10. + radar2, vel_field="velocity", wind_size=(3, 3), nyq=10.0 ) assert_allclose(texture_field2["data"], texture_field3["data"], atol=1e-10) From c999ccd3cdd6f4c7a7b93bd131355c188723e18e Mon Sep 17 00:00:00 2001 From: isilber Date: Fri, 1 Dec 2023 05:28:34 +0000 Subject: [PATCH 10/30] FIX: metadata in core.radar.get_elevation --- pyart/core/radar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/core/radar.py b/pyart/core/radar.py index ef96e6791c..85a0398e9a 100644 --- a/pyart/core/radar.py +++ b/pyart/core/radar.py @@ -449,7 +449,7 @@ def get_elevation(self, sweep, copy=False): Returns ------- - azimuths : array + elevation : array Array containing the elevation angles for a given sweep. """ From 48f74aa2e0890f07a645732c21d37330cc00e1cb Mon Sep 17 00:00:00 2001 From: isilber Date: Fri, 1 Dec 2023 22:44:45 +0000 Subject: [PATCH 11/30] FIX: outdated text in masking_data_with_gatefilters.ipynb --- doc/source/notebooks/masking_data_with_gatefilters.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/notebooks/masking_data_with_gatefilters.ipynb b/doc/source/notebooks/masking_data_with_gatefilters.ipynb index 79a86d41b0..c9a1335ea5 100644 --- a/doc/source/notebooks/masking_data_with_gatefilters.ipynb +++ b/doc/source/notebooks/masking_data_with_gatefilters.ipynb @@ -135,7 +135,7 @@ "# The below example shows how you can use a gatefilter to exclude values where gate_id is\n", "# equal to a certain value. (You can exclude values based on conditions from multiple\n", "# fields. In this example, we are excluding regions identified as second trip echoes\n", - "# (gate_id = 0), nonhydrometeors (gate_id = 1), and clutter (gate_id = 2).\n", + "# (gate_id = 0), nonhydrometeors (gate_id = 3), and clutter (gate_id = 5).\n", "gatefilter.exclude_equal(\"gate_id\", 0)\n", "gatefilter.exclude_equal(\"gate_id\", 3)\n", "gatefilter.exclude_equal(\"gate_id\", 5)" From cb6324e25136ceef7b023c86b15c2600c5f8a87d Mon Sep 17 00:00:00 2001 From: isilber Date: Tue, 5 Dec 2023 22:53:20 +0000 Subject: [PATCH 12/30] ADD: method to automatically determine sweep number and sweep start and end indices --- pyart/util/__init__.py | 1 + pyart/util/radar_utils.py | 74 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+) 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..d0e4d114b4 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,79 @@ def to_vpt(radar, single_scan=True): return +def determine_sweeps(radar, + max_offset=0.1, + running_win_dt=5., + deg_rng=(-5., 360.)): + """ + 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.) + sample_dt = np.nanmean(np.diff(radar.time['data'])) + win_size = int(np.ceil(running_win_dt / sample_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] == True: # sweep did not start + t += 1 + continue + bincounts, _ = np.histogram(fix_win, bins=angle_bins) + moving_radar = np.sum(bincounts > 0) > 1 # radar is likely moving to a new sweep position + if in_sweep == True: + if moving_radar | (t == radar.time['data'].size - win_size): + in_sweep = False + sweep_end_index.append(t + win_size - 1) + t += win_size - 1 + elif np.all(idle_sweep == False) & (moving_radar == False): + in_sweep = True + sweep_start_index.append(t) + t += 1 + sweep_number = np.arange(len(sweep_start_index)) + 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') + return + + def subset_radar( radar, field_names, From 26039ce175b00f5713b23bafecb952982172522b Mon Sep 17 00:00:00 2001 From: isilber Date: Tue, 5 Dec 2023 22:58:05 +0000 Subject: [PATCH 13/30] add comment to 'determine_sweeps' method --- pyart/util/radar_utils.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index d0e4d114b4..96350f2e3f 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -172,6 +172,9 @@ def determine_sweeps(radar, 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') From 8ae5a44f3983a6853eb9596e137ff858a1e84e48 Mon Sep 17 00:00:00 2001 From: isilber Date: Tue, 5 Dec 2023 23:52:13 +0000 Subject: [PATCH 14/30] FIX: indexing bug in 'determine_sweeps' method --- pyart/util/radar_utils.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 96350f2e3f..9d6e8991f6 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -147,6 +147,8 @@ def determine_sweeps(radar, angle_bins = np.arange(deg_rng[0] - max_offset, deg_rng[1] + max_offset + 1e-10, max_offset * 2.) 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 @@ -163,10 +165,12 @@ def determine_sweeps(radar, bincounts, _ = np.histogram(fix_win, bins=angle_bins) moving_radar = np.sum(bincounts > 0) > 1 # radar is likely moving to a new sweep position if in_sweep == True: - if moving_radar | (t == radar.time['data'].size - win_size): + 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 - 1) - t += win_size - 1 + sweep_end_index.append(t + win_size - 2) + t += win_size - 2 elif np.all(idle_sweep == False) & (moving_radar == False): in_sweep = True sweep_start_index.append(t) From 47966fed4c51deb0e83ef2fd7dcc2aca801df018 Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 00:09:34 +0000 Subject: [PATCH 15/30] ADD: unit tests for 'determine_sweeps' --- tests/util/test_radar_utils.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index 442b2cc0ab..eefff844e0 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -50,6 +50,24 @@ 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 + + # 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 + + def test_subset_radar(): radar = pyart.testing.make_empty_ppi_radar(10, 36, 3) field = {"data": np.ones((36 * 3, 10))} From ea1a94b1a0461072ba7e4161efba94a50246e8ab Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 00:27:02 +0000 Subject: [PATCH 16/30] linting --- pyart/util/radar_utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 9d6e8991f6..389333971e 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -159,19 +159,19 @@ def determine_sweeps(radar, 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] == True: # sweep did not start + 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 is likely moving to a new sweep position - if in_sweep == True: + 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 == False) & (moving_radar == False): + elif np.all(~idle_sweep) & ~moving_radar: in_sweep = True sweep_start_index.append(t) t += 1 From 1d338c10c86fbf391a7f63b8d5f4c60963fab99f Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 00:47:26 +0000 Subject: [PATCH 17/30] FIX: nsweeps wasn't updated when calling 'determine_sweeps'; tests were updated accordingly --- pyart/util/radar_utils.py | 1 + tests/util/test_radar_utils.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 389333971e..a027bedb2b 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -182,6 +182,7 @@ def determine_sweeps(radar, 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') + radar.nsweeps = len(sweep_number) return diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index eefff844e0..be7050c47a 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -58,6 +58,7 @@ def test_determine_sweeps(): 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) @@ -66,6 +67,7 @@ def test_determine_sweeps(): 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(): From 2324233f4fabefadaa9f6d6137d7ccd4c35a877f Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 04:41:17 +0000 Subject: [PATCH 18/30] line reformatting per black --- pyart/util/radar_utils.py | 43 +++++++++++++++++----------------- tests/util/test_radar_utils.py | 20 +++++++++------- 2 files changed, 34 insertions(+), 29 deletions(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index a027bedb2b..f5ba20b975 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -105,10 +105,7 @@ def to_vpt(radar, single_scan=True): return -def determine_sweeps(radar, - max_offset=0.1, - running_win_dt=5., - deg_rng=(-5., 360.)): +def determine_sweeps(radar, max_offset=0.1, running_win_dt=5., deg_rng=(-5., 360.)): """ determine the number of sweeps using elevation data (PPI scans) or azimuth data (RHI scans) and update the input radar object @@ -135,38 +132,42 @@ def determine_sweeps(radar, """ # 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'] + 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'] + 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.) - sample_dt = np.nanmean(np.diff(radar.time['data'])) + angle_bins = np.arange( + deg_rng[0] - max_offset, deg_rng[1] + max_offset + 1e-10, max_offset * 2. + ) + 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') + 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] + 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 is likely moving to a new sweep position + 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) + 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) @@ -179,9 +180,9 @@ def determine_sweeps(radar, # 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') + 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") radar.nsweeps = len(sweep_number) return diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index be7050c47a..1d2eabf422 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -53,20 +53,24 @@ def test_to_vpt(): 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) + 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 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 + 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 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 From 7c6d52205a5aa24c8f81cef28d8edbfa2290b677 Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 05:55:24 +0000 Subject: [PATCH 19/30] FIX: update multiple radar object attributes to ensure its full utilization with Py-ART after 'determine_sweeps' is called --- pyart/util/radar_utils.py | 12 ++++++++++++ tests/util/test_radar_utils.py | 4 ++-- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index f5ba20b975..c9a0dfb4ba 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -183,7 +183,19 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5., deg_rng=(-5., 360 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.antenna_transition["data"].size) + 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 "{:<22}".format(radar.scan_type)], dtype='|S1') + radar.sweep_mode = ma.array(np.tile(bstr_entry[np.newaxis, :], (radar.nsweeps, 1))) return diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index 1d2eabf422..f6dd40fdb1 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -54,7 +54,7 @@ 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 + np.arange(1, 36 * 3 + 1) / 36 ) pyart.util.determine_sweeps(radar) assert np.all(radar.sweep_end_ray_index["data"] == [35, 71, 107]) @@ -65,7 +65,7 @@ def test_determine_sweeps(): # 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 + 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]) From 11999c71f36c6b32ac66635ebdb6d0b49fe1777e Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 06:08:26 +0000 Subject: [PATCH 20/30] line reformatting per black --- pyart/util/radar_utils.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index c9a0dfb4ba..54db73836b 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -105,7 +105,7 @@ def to_vpt(radar, single_scan=True): return -def determine_sweeps(radar, max_offset=0.1, running_win_dt=5., deg_rng=(-5., 360.)): +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 @@ -142,13 +142,13 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5., deg_rng=(-5., 360 # 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. + 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' + "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 @@ -183,9 +183,13 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5., deg_rng=(-5., 360 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") + 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.antenna_transition["data"].size) for i in range(radar.nsweeps): @@ -194,7 +198,7 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5., deg_rng=(-5., 360 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 "{:<22}".format(radar.scan_type)], dtype='|S1') + 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 From e197ab221419c4343d37f37820f84d5ecdd59bfe Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 06:11:29 +0000 Subject: [PATCH 21/30] linting --- pyart/util/radar_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 54db73836b..a8da47cf83 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -198,7 +198,7 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5.0, deg_rng=(-5.0, 3 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') + 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 From 83e2e166cdbd2e67ca2bec31fe1cc8b802722d9e Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 06:30:27 +0000 Subject: [PATCH 22/30] ENH: add antenna_transition to the testing module's sample objects --- pyart/testing/sample_objects.py | 3 +++ pyart/util/radar_utils.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) 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/radar_utils.py b/pyart/util/radar_utils.py index a8da47cf83..64735a8a50 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -191,7 +191,7 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5.0, deg_rng=(-5.0, 3 ] radar.fixed_angle["data"] = ma.array(fixed_angle, dtype="float32") radar.nsweeps = len(sweep_number) - transition = np.zeros(radar.antenna_transition["data"].size) + transition = np.zeros(radar.nrays) for i in range(radar.nsweeps): if i == 0: transition[: sweep_start_index[i]] = 1 From f4592b9545d41a65559fa26d46c189f174257997 Mon Sep 17 00:00:00 2001 From: isilber Date: Wed, 6 Dec 2023 06:47:30 +0000 Subject: [PATCH 23/30] ENH: update unit tests to reflect antenna_transition allocation in sample_objects --- tests/core/test_radar.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 From 312c50845eb92d1c917e6acae77b100073a61c1b Mon Sep 17 00:00:00 2001 From: Israel Silber <55999320+isilber@users.noreply.github.com> Date: Wed, 6 Dec 2023 08:03:13 -0800 Subject: [PATCH 24/30] Update pyart/util/radar_utils.py Co-authored-by: Max Grover --- pyart/util/radar_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 64735a8a50..b4097b62fd 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -112,7 +112,7 @@ def determine_sweeps(radar, max_offset=0.1, running_win_dt=5.0, deg_rng=(-5.0, 3 Parameters ---------- - radar : radar object + radar : Radar object The radar object containing the data. max_offset : float Maximum elevation offset (if is_ppi is True) or azimuth offset (if From 17412796b9c6a88de205eab8cd8d8b98353c13fa Mon Sep 17 00:00:00 2001 From: Israel Silber <55999320+isilber@users.noreply.github.com> Date: Wed, 6 Dec 2023 08:03:24 -0800 Subject: [PATCH 25/30] Update pyart/util/radar_utils.py Co-authored-by: Max Grover --- pyart/util/radar_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index b4097b62fd..c2e7c405b9 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -107,7 +107,7 @@ def to_vpt(radar, single_scan=True): 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 + Determine the number of sweeps using elevation data (PPI scans) or azimuth data (RHI scans) and update the input radar object Parameters From 9f72439063a82c85a2d1a73e0833ad89b1c303ca Mon Sep 17 00:00:00 2001 From: isilber Date: Thu, 18 Jan 2024 22:11:43 +0000 Subject: [PATCH 26/30] ENH: add 'zdist_factor' input parameter to 'map_gates_to_grid'. This parameter scales the distance in the z-dimension when calculating weights upon gridding. It is a semi-equivalent to the scaling factors in the RoI class methods. Example for using this parameter would be setting a 'zdist_factor' value of 0.0 combined with an h_factor=0.0 (if calling DistBeamRoI) or z_factor=0.0 (if calling DistRoI), resulting in the exclusion of the z dimension in gridding weighting, which could serve as a potential solution for gridding a single PPI sweep from a single radar. In that case, the distance in the weighting function effectively serves as the distance of a given point from the PPI sweep's cone. --- pyart/map/_gate_to_grid_map.pyx | 18 +++++++++++++----- pyart/map/gates_to_grid.py | 2 ++ 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/pyart/map/_gate_to_grid_map.pyx b/pyart/map/_gate_to_grid_map.pyx index 2bb7531937..512af344a2 100644 --- a/pyart/map/_gate_to_grid_map.pyx +++ b/pyart/map/_gate_to_grid_map.pyx @@ -232,7 +232,8 @@ cdef class GateToGridMapper: float[:, ::1] gate_z, float[:, ::1] gate_y, float[:, ::1] gate_x, float[:, :, ::1] field_data, char[:, :, ::1] field_mask, char[:, ::1] excluded_gates, - float toa, RoIFunction roi_func, int weighting_function): + float toa, RoIFunction roi_func, int weighting_function, + float zdist_factor): """ Map radar gates unto the regular grid. @@ -264,6 +265,12 @@ cdef class GateToGridMapper: Function to use for weighting gates based upon distance. 0 for Barnes, 1 for Cressman, 2 for Nearest and 3 for Barnes 2 neighbor weighting. + zdist_factor: float + Scaling factor for squared z difference in distance calculation. + A value of 0.0 combined with an h_factor=0.0 (if calling + DistBeamRoI) or z_factor=0.0 (if calling DistRoI) results in + the exclusion of the z dimension in gridding weighting and could + serve as a potential solution for gridding a single PPI sweep. """ @@ -286,7 +293,8 @@ cdef class GateToGridMapper: values = field_data[nray, ngate] masks = field_mask[nray, ngate] - self.map_gate(x, y, z, roi, values, masks, weighting_function) + self.map_gate(x, y, z, roi, values, masks, weighting_function, + zdist_factor) @cython.initializedcheck(False) @cython.cdivision(True) @@ -294,7 +302,7 @@ cdef class GateToGridMapper: @cython.wraparound(False) cdef int map_gate(self, float x, float y, float z, float roi, float[:] values, char[:] masks, - int weighting_function): + int weighting_function, float zdist_factor): """ Map a single gate to the grid. """ cdef float xg, yg, zg, dist, weight, roi2, dist2, min_dist2 @@ -340,7 +348,7 @@ cdef class GateToGridMapper: xg = self.x_step * xi yg = self.y_step * yi zg = self.z_step * zi - dist = ((xg - x)**2 + (yg - y)**2 + (zg - z)**2) + dist = ((xg - x)**2 + (yg - y)**2 + zdist_factor * (zg - z)**2) if dist >= roi2: continue for i in range(self.nfields): @@ -362,7 +370,7 @@ cdef class GateToGridMapper: xg = self.x_step * xi yg = self.y_step * yi zg = self.z_step * zi - dist2 = (xg-x)*(xg-x) + (yg-y)*(yg-y) + (zg-z)*(zg-z) + dist2 = (xg-x)*(xg-x) + (yg-y)*(yg-y) + zdist_factor * (zg-z)*(zg-z) if dist2 > roi2: continue diff --git a/pyart/map/gates_to_grid.py b/pyart/map/gates_to_grid.py index fd2e0ddabe..49771e71fa 100644 --- a/pyart/map/gates_to_grid.py +++ b/pyart/map/gates_to_grid.py @@ -39,6 +39,7 @@ def map_gates_to_grid( h_factor=1.0, nb=1.5, bsp=1.0, + zdist_factor=1.0, **kwargs ): """ @@ -175,6 +176,7 @@ def map_gates_to_grid( toa, roi_func, cy_weighting_function, + zdist_factor, ) # create and return the grid dictionary From 9c470876e5bc8f60755e7edc0464628ece46811a Mon Sep 17 00:00:00 2001 From: isilber Date: Thu, 18 Jan 2024 22:28:33 +0000 Subject: [PATCH 27/30] DEL: remove a redundant (liekly residual) input parameter (toa) from map_gates_to_grid --- pyart/map/_gate_to_grid_map.pyx | 2 +- pyart/map/gates_to_grid.py | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/pyart/map/_gate_to_grid_map.pyx b/pyart/map/_gate_to_grid_map.pyx index 512af344a2..6dd22f4521 100644 --- a/pyart/map/_gate_to_grid_map.pyx +++ b/pyart/map/_gate_to_grid_map.pyx @@ -232,7 +232,7 @@ cdef class GateToGridMapper: float[:, ::1] gate_z, float[:, ::1] gate_y, float[:, ::1] gate_x, float[:, :, ::1] field_data, char[:, :, ::1] field_mask, char[:, ::1] excluded_gates, - float toa, RoIFunction roi_func, int weighting_function, + RoIFunction roi_func, int weighting_function, float zdist_factor): """ Map radar gates unto the regular grid. diff --git a/pyart/map/gates_to_grid.py b/pyart/map/gates_to_grid.py index 49771e71fa..882f4c6785 100644 --- a/pyart/map/gates_to_grid.py +++ b/pyart/map/gates_to_grid.py @@ -30,7 +30,6 @@ def map_gates_to_grid( gatefilters=False, map_roi=True, weighting_function="Barnes2", - toa=17000.0, roi_func="dist_beam", constant_roi=None, z_factor=0.05, @@ -173,7 +172,6 @@ def map_gates_to_grid( field_data, field_mask, excluded_gates, - toa, roi_func, cy_weighting_function, zdist_factor, From d8899cea6bc0ce7c8e8a08156ab1ae61ec78ae95 Mon Sep 17 00:00:00 2001 From: isilber Date: Thu, 18 Jan 2024 22:54:02 +0000 Subject: [PATCH 28/30] MNT: set a consistent nomenclature for the distance squared parameter. In parts of the 'map_grid' method it was referred to as 'dist' whereas in another as 'dist2' (parameter calculation lines differ in style but identical in value), resulting in an unnecessary scalar allocation, and a confusing nomenclature. --- pyart/map/_gate_to_grid_map.pyx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyart/map/_gate_to_grid_map.pyx b/pyart/map/_gate_to_grid_map.pyx index 6dd22f4521..b7eb8292b3 100644 --- a/pyart/map/_gate_to_grid_map.pyx +++ b/pyart/map/_gate_to_grid_map.pyx @@ -305,7 +305,7 @@ cdef class GateToGridMapper: int weighting_function, float zdist_factor): """ Map a single gate to the grid. """ - cdef float xg, yg, zg, dist, weight, roi2, dist2, min_dist2 + cdef float xg, yg, zg, weight, roi2, dist2, min_dist2 cdef int x_min, x_max, y_min, y_max, z_min, z_max cdef int xi, yi, zi, x_argmin, y_argmin, z_argmin @@ -348,12 +348,12 @@ cdef class GateToGridMapper: xg = self.x_step * xi yg = self.y_step * yi zg = self.z_step * zi - dist = ((xg - x)**2 + (yg - y)**2 + zdist_factor * (zg - z)**2) - if dist >= roi2: + dist2 = ((xg - x)**2 + (yg - y)**2 + zdist_factor * (zg - z)**2) + if dist2 >= roi2: continue for i in range(self.nfields): - if dist < self.min_dist2[zi, yi, xi, i]: - self.min_dist2[zi, yi, xi, i] = dist + if dist2 < self.min_dist2[zi, yi, xi, i]: + self.min_dist2[zi, yi, xi, i] = dist2 x_argmin = xi y_argmin = yi z_argmin = zi From 137098e25d21b60299304ba034ed5e75b048a933 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Tue, 23 Jan 2024 14:08:48 -0500 Subject: [PATCH 29/30] DOC: Add example for gridding data --- .../mapping/plot_grid_single_sweep_ppi.py | 105 ++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 examples/mapping/plot_grid_single_sweep_ppi.py diff --git a/examples/mapping/plot_grid_single_sweep_ppi.py b/examples/mapping/plot_grid_single_sweep_ppi.py new file mode 100644 index 0000000000..bb22864c03 --- /dev/null +++ b/examples/mapping/plot_grid_single_sweep_ppi.py @@ -0,0 +1,105 @@ +""" +====================================== +Map a single PPI sweep to a Cartesian grid using 2D weighting +====================================== + +Map the reflectivity field of a single PPI sweep from Antenna (polar) coordinates +to a Cartesian grid, while using a 2D weighting. +This solution is valid only for this case (single PPI sweep), yet it is a useful one +(an arguably global solution since it overlooks the z-dimension grid limits). +The exclusion of the z-dimension from the RoI and weighting calculations results in +minor errors, especially considering the high co-variance of neighboring radar +volumes and the typically small Cartesian grid separation. Errors are effectively +negligible at low elevation angles common to PPI sweeps. + +""" +print(__doc__) + +# ===================== +# Author: Israel Silber +# ===================== + +import cartopy.crs as ccrs +import matplotlib.pyplot as plt +from open_radar_data import DATASETS + +import pyart + +# file, and fields +# ======================= +file_name = DATASETS.fetch("example_plot_ppi_single_sweep.nc") +processed_field = "reflectivity_at_cor" + +# load file +# ======================= +radar_obj = pyart.io.read(file_name) +print(f"Total sweeps = {radar_obj.nsweeps}") + +# Plot polar coordinate (raw) data +# ======================= +print( + "===\n\nnow displaying Ze data for 2nd and 4th sweeps in polar coordinates\n" + "(both approaching the 40 km range denoted by the purple ring)\n\n===" +) +fig = plt.figure(figsize=(12, 6), tight_layout=True) +fig.suptitle("polar") +for sweep, ax_ind in zip([1, 3], range(2)): + ax = fig.add_subplot(1, 2, ax_ind + 1, projection=ccrs.Mercator()) + sweep_obj = radar_obj.extract_sweeps((sweep,)) + display = pyart.graph.RadarDisplay(sweep_obj) + display.plot( + processed_field, + sweep=0, + ax=ax, + ) + display.plot_range_rings([10, 20, 30]) + display.plot_range_rings([40], col="purple") +plt.show() + +print( + "===\n\nnow displaying gridded Ze data for 2nd and 4th (final) sweeps; note the " + "truncated max range in the case of the 4th sweep\n\n===" +) +fig2 = plt.figure(figsize=(12, 6), tight_layout=True) +fig2.suptitle("Cartesian gridded") +for sweep, ax_ind in zip([1, 3], range(2)): + sweep_obj = radar_obj.extract_sweeps((sweep,)) + grid = pyart.map.grid_from_radars( + (sweep_obj,), + grid_shape=(1, 1601, 1601), + grid_limits=((0, 10000.0), [-40000, 40000], [-40000, 40000]), + fields=[processed_field], + ) + ax = fig2.add_subplot(1, 2, ax_ind + 1) + ax.imshow( + grid.fields[processed_field]["data"][0], + origin="lower", + extent=(-40, 40, -40, 40), + ) + ax.set_title(f"sweep #{sweep + 1}") +plt.show() + +print( + "===\n\nUsing 2D weighting by " + "setting h_factor and zdist_factor to 0.0, the max range looks OK now\n\n===" +) +fig2 = plt.figure(figsize=(12, 6), tight_layout=True) +fig2.suptitle("Cartesian gridded") +for sweep, ax_ind in zip([1, 3], range(2)): + sweep_obj = radar_obj.extract_sweeps((sweep,)) + grid = pyart.map.grid_from_radars( + (sweep_obj,), + grid_shape=(1, 1601, 1601), + grid_limits=((0, 10000.0), [-40000, 40000], [-40000, 40000]), + fields=[processed_field], + h_factor=0.0, + zdist_factor=0.0, + ) + ax = fig2.add_subplot(1, 2, ax_ind + 1) + ax.imshow( + grid.fields[processed_field]["data"][0], + origin="lower", + extent=(-40, 40, -40, 40), + ) + ax.set_title(f"sweep #{sweep + 1}") +plt.show() From 1d960f51f127360847e0059e5808e4f122ff0fd6 Mon Sep 17 00:00:00 2001 From: Israel Silber <55999320+isilber@users.noreply.github.com> Date: Tue, 23 Jan 2024 13:00:38 -0800 Subject: [PATCH 30/30] Update pyart/map/gates_to_grid.py Co-authored-by: Max Grover --- pyart/map/gates_to_grid.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyart/map/gates_to_grid.py b/pyart/map/gates_to_grid.py index 882f4c6785..3031dc87be 100644 --- a/pyart/map/gates_to_grid.py +++ b/pyart/map/gates_to_grid.py @@ -30,6 +30,7 @@ def map_gates_to_grid( gatefilters=False, map_roi=True, weighting_function="Barnes2", + toa=17000.0, roi_func="dist_beam", constant_roi=None, z_factor=0.05,