ENH: optional rectangular window dims for velocity texture analysis a…

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

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)

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

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

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)