From 70c314c4640a999357d8cc0c7cc804af8c9337db Mon Sep 17 00:00:00 2001 From: Israel Silber <55999320+isilber@users.noreply.github.com> Date: Thu, 17 Aug 2023 07:17:23 -0700 Subject: [PATCH] ENH: optional rectangular window dims for velocity texture analysis and bug fixes to the load_kazr method (#1446) * 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 * FIX: Add in linting fixes --------- Co-authored-by: mgrover1 --- pyart/aux_io/arm_vpt.py | 2 +- pyart/retrieve/simple_moment_calculations.py | 18 +++++++++++------- pyart/util/sigmath.py | 16 +++++++++++----- .../test_simple_moment_calculations.py | 18 ++++++++++++++++++ 4 files changed, 41 insertions(+), 13 deletions(-) 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) diff --git a/pyart/retrieve/simple_moment_calculations.py b/pyart/retrieve/simple_moment_calculations.py index 9b7c8f5f7e..61beb2db04 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. @@ -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) @@ -300,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, wind_size) - ) + vel_texture_field["data"] = ndimage.median_filter(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 diff --git a/tests/retrieve/test_simple_moment_calculations.py b/tests/retrieve/test_simple_moment_calculations.py index 40708f7aee..b222a317cf 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,17 @@ 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.0 + ) + texture_field3 = pyart.retrieve.calculate_velocity_texture( + radar2, vel_field="velocity", wind_size=(3, 3), nyq=10.0 + ) + assert_allclose(texture_field2["data"], texture_field3["data"], atol=1e-10)