grid_from_radars generates misleading gridded arrays when applied on the last sweep of a given file

[isilber]

• Py-ART version: 1.16.1
• Python version: 3.11.5
• Operating System: Red Hat 7.9

Description

When converting radar objects to a grid object (following a single sweep extraction) using the grid_from_radars with the default map_gates_to_grid algorithm, gridded fields have an apparent maximum range much smaller than the true range. This artifact only occurs when the last sweep is processed, so if nsweeps is 4, the artifact will be generated only when gridding is applied to the 4th sweep.
I saw this issue in multiple files in multiple datastreams (e.g., ARM's anxkasacr and houkasacr) and this issue is consistent.
The ROI field appears consistent in all cases, while the map module does not appear to use any sweep information directly. I have no leads on this, but I suspect the solution should be straightforward (an indexing issue?) to someone more familiar with the code base.

Working example shown below.

What I Did

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pyart

# path, file, and fields
file_path = '/path/to/datastream/anxkasacrcfrcorppivM1.c1/'
file_name = 'anxkasacrcfrcorppivM1.c1.20200312.003009.nc'
processed_field = 'reflectivity_at_cor'

# load file
radar_obj = pyart.io.read(f"{file_path}{file_name}")  # load file
print(f"Total sweeps = {radar_obj.nsweeps}")

# Plot polar coordinate (raw) data
print(f"now displaying Ze data for 2nd and 4th (final) sweeps in polar coordinates\n"
      f"(both approaching the 40 km range)")
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.RadarMapDisplay(sweep_obj)
    display.plot_ppi_map(
        processed_field,
        sweep=0,
        ax=ax,
    )
    display.plot_range_rings([10, 20, 30, 40])
plt.show()

print(f"now displaying gridded Ze data for 2nd and 4th (final) sweeps")
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.), [-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))
plt.show()

[mgrover1]

@isilber - sorry for the delay here, do you have a sample file we can test this with?

[isilber]

anxkasacrcfrcorppivM1.c1.20200312.003009.zip
@mgrover1 I'm attaching the file you can use in the example I added above. Thanks!

[mgrover1]

Thank you!!

[mgrover1]

@isilber - with your current settings, you are configuring the elevation to grid to as 0 meters (above ground level). The scan strategy outlined in the data would look something like this:

You are using the default gridding parameters, which is the barnes scheme, as indicated in the documentation
https://arm-doe.github.io/pyart/API/generated/pyart.map.map_to_grid.html#pyart.map.map_to_grid

The default minimum radius for the scheme is 500 meter. Looking at the scan strategy visualization here, it would make sense for the 2 degree elevation scan to be beyond 500 meters from the desired grid level, at that ~15 km point we see in the data.

I attached the original code snippet, using 0 meters as the grid level

If we change the following snippet to use 500 meters as the lower bound instead:

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),
    # Use 500 meters as the lower bound for z
    grid_limits=((500, 10000.), [-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))

we see the following:

Does this make sense? Where in the documentation do you think we could be more clear here?

[mgrover1]

@isilber - is it okay if we transfer this over to a discussion?

[isilber]

@mgrover1, my pull request #1509 resolves this issue. The additions to the code are minimal (no new methods) and are consistent with the general structure of the gridding module. These edits do not require us to set a case-dependent value to the grid_limits input parameter but by simply setting the h_factor parameter (which was already in the code base) and the zdist_factor parameter (which was added now) to 0.0, for example:

    grid = pyart.map.grid_from_radars(
    (sweep_obj,),
    grid_shape=(1, 1601, 1601),
    grid_limits=((0., 10000.), [-40000, 40000], [-40000, 40000]),
    fields=[processed_field],
    h_factor=0.0,
    zdist_factor=0.0,
    )

I think an intuitive way to think about that "2-D" option (or better say, case) is for gridding a single PPI sweep of a single radar to a grid based on the distance of points from the sweep's (elevation-angle-dependent) cone.
This will also resolve other cases in which the elevation angle might be higher, such as in ARM's EPCAPE KaSACR PPI files, in which setting z-dimension in grid_limits to (500, 10000) will no longer work and would therefore require some additional processing step of vertical dimension manipulation such as traditional composites.