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() diff --git a/pyart/map/_gate_to_grid_map.pyx b/pyart/map/_gate_to_grid_map.pyx index 2bb7531937..b7eb8292b3 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): + 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,10 +302,10 @@ 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 + 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 diff --git a/pyart/map/gates_to_grid.py b/pyart/map/gates_to_grid.py index fd2e0ddabe..3031dc87be 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 ): """ @@ -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