ENH: Add 'zdist_factor' input parameter to 'map_gates_to_grid'. (#1509)

* 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]>

* 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.

* DEL: remove a redundant (liekly residual) input parameter (toa) from map_gates_to_grid

* 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.

* DOC: Add example for gridding data

* Update pyart/map/gates_to_grid.py

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

---------

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

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

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):
+            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,18 +293,19 @@ 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)
     @cython.boundscheck(False)
     @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
+        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
 
@@ -340,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 + (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
@@ -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

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
 ):
     """
@@ -172,9 +173,9 @@ def map_gates_to_grid(
             field_data,
             field_mask,
             excluded_gates,
-            toa,
             roi_func,
             cy_weighting_function,
+            zdist_factor,
         )
 
     # create and return the grid dictionary