From 860612892897370d327049d4b7e2a1cdffc4656e Mon Sep 17 00:00:00 2001 From: isilber Date: Mon, 14 Aug 2023 16:36:38 +0000 Subject: [PATCH 1/8] 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 2/8] 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 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 57535ed819e5e2395f368cf91607459c33396359 Mon Sep 17 00:00:00 2001 From: mgrover1 Date: Wed, 16 Aug 2023 13:16:31 -0500 Subject: [PATCH 8/8] 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)