Nearest neighbor interpolation issue when gridding

[lauratomkins]

I'm having an issue with the weighting = 'Nearest' option in the pyart.map.grid_from_radars function with NEXRAD Level 2 data. Anytime I use this weighting option, every field except the reflectivity looks wrong. I can use the Cressman weighting with no issues. I've provided some example images and code to reproduce. Any insight would be very helpful! Thanks.

Here is an example set of images using weighting = 'Nearest'

And here are the same example set using weighting = 'Cressman'

Here is the code to reproduce:

import pyart
import matplotlib.pyplot as plt

# import NEXRAD level2 file
radar = pyart.io.read_nexrad_archive('KOKX20200213_060331_V06')
radar = radar.extract_sweeps([0,1]) # extract first 2 0.5deg sweeps

# grid
center_lat, center_lon, _ = pyart.io.nexrad_common.get_nexrad_location('KOKX')

grid_nearestn = pyart.map.grid_from_radars(
        (radar,), grid_shape=(1,241,241),
        grid_limits=((0, 10000), (-200000, 200000), (-200000, 200000)),
        grid_projection={'proj':'lcca', 'lat_0':center_lat, 'lon_0':center_lon},
        weighting_function='Nearest') # nearest neightbor


grid_cressman = pyart.map.grid_from_radars(
        (radar,), grid_shape=(1,241,241),
        grid_limits=((0, 10000), (-200000, 200000), (-200000, 200000)),
        grid_projection={'proj':'lcca', 'lat_0':center_lat, 'lon_0':center_lon},
        weighting_function='Cressman') # cressman weighting

# Plot

# Nearest Neighbor grid
x = grid_nearestn.point_x['data'][0]/1000
y = grid_nearestn.point_y['data'][0]/1000

figs, axs = plt.subplots(2,2, figsize=(10,10)) 

ref = axs[0,0].pcolormesh(x, y, grid_nearestn.fields['reflectivity']['data'][0], 
         cmap='magma_r', vmin=0, vmax=40, shading='nearest')
ref_cbar = plt.colorbar(ref, ax=axs[0,0], fraction=0.040, pad=0.04)
axs[0,0].set_title(radar.metadata['instrument_name'] + ' Reflectivity',fontsize=14)
axs[0,0].set(xlim=[-125,125], ylim=[-125,125])
axs[0,0].set_aspect('equal','box')

vel = axs[0,1].pcolormesh(x, y, grid_nearestn.fields['velocity']['data'][0], 
         cmap='RdBu_r', vmin=-40, vmax=40, shading='nearest')
vel_cbar = plt.colorbar(vel, ax=axs[0,1], fraction=0.040, pad=0.04)
axs[0,1].set_title(radar.metadata['instrument_name'] + ' Velocity',fontsize=14)
axs[0,1].set(xlim=[-125,125], ylim=[-125,125])
axs[0,1].set_aspect('equal','box')

rho = axs[1,0].pcolormesh(x, y, grid_nearestn.fields['cross_correlation_ratio']['data'][0], 
         cmap='cividis', vmin=0.8, vmax=1, shading='nearest')
rho_cbar = plt.colorbar(rho, ax=axs[1,0], fraction=0.040, pad=0.04)
axs[1,0].set_title(radar.metadata['instrument_name'] + ' RHOHV',fontsize=14)
axs[1,0].set(xlim=[-125,125], ylim=[-125,125])
axs[1,0].set_aspect('equal','box')

sw = axs[1,1].pcolormesh(x, y, grid_nearestn.fields['spectrum_width']['data'][0], 
         cmap='plasma', vmin=0, vmax=5, shading='nearest')
sw_cbar = plt.colorbar(sw, ax=axs[1,1], fraction=0.040, pad=0.04)
axs[1,1].set_title(radar.metadata['instrument_name'] + ' SW',fontsize=14)
axs[1,1].set(xlim=[-125,125], ylim=[-125,125])
axs[1,1].set_aspect('equal','box')

figs.savefig('nearest.png', dpi=300, bbox_inches='tight')


# Cressman  grid
x = grid_cressman.point_x['data'][0]/1000
y = grid_cressman.point_y['data'][0]/1000

figs, axs = plt.subplots(2,2, figsize=(10,10)) 

ref = axs[0,0].pcolormesh(x, y, grid_cressman.fields['reflectivity']['data'][0], 
         cmap='magma_r', vmin=0, vmax=40, shading='nearest')
ref_cbar = plt.colorbar(ref, ax=axs[0,0], fraction=0.040, pad=0.04)
axs[0,0].set_title(radar.metadata['instrument_name'] + ' Reflectivity',fontsize=14)
axs[0,0].set(xlim=[-125,125], ylim=[-125,125])
axs[0,0].set_aspect('equal','box')

vel = axs[0,1].pcolormesh(x, y, grid_cressman.fields['velocity']['data'][0], 
         cmap='RdBu_r', vmin=-40, vmax=40, shading='nearest')
vel_cbar = plt.colorbar(vel, ax=axs[0,1], fraction=0.040, pad=0.04)
axs[0,1].set_title(radar.metadata['instrument_name'] + ' Velocity',fontsize=14)
axs[0,1].set(xlim=[-125,125], ylim=[-125,125])
axs[0,1].set_aspect('equal','box')

rho = axs[1,0].pcolormesh(x, y, grid_cressman.fields['cross_correlation_ratio']['data'][0], 
         cmap='cividis', vmin=0.8, vmax=1, shading='nearest')
rho_cbar = plt.colorbar(rho, ax=axs[1,0], fraction=0.040, pad=0.04)
axs[1,0].set_title(radar.metadata['instrument_name'] + ' RHOHV',fontsize=14)
axs[1,0].set(xlim=[-125,125], ylim=[-125,125])
axs[1,0].set_aspect('equal','box')

sw = axs[1,1].pcolormesh(x, y, grid_cressman.fields['spectrum_width']['data'][0], 
         cmap='plasma', vmin=0, vmax=5, shading='nearest')
sw_cbar = plt.colorbar(sw, ax=axs[1,1], fraction=0.040, pad=0.04)
axs[1,1].set_title(radar.metadata['instrument_name'] + ' SW',fontsize=14)
axs[1,1].set(xlim=[-125,125], ylim=[-125,125])
axs[1,1].set_aspect('equal','box')

figs.savefig('cressman.png', dpi=300, bbox_inches='tight')

[zssherman]

Hi @lauratomkins, do you have the most recent github clone of pyart or is it a conda forge install? We recently had a fixed for nearest neighbor, but I haven't done a new pyart conda-forge release yet.

[lauratomkins]

I have a conda forge install

[zssherman]

Gotcha, can you try a install of the github clone? The fix isn't in the conda release. I'll need to get going on a new release.

[lauratomkins]

@zssherman I was able to install with github and here's what the nearest interpolation looks like now. Seems to be better (in that the data is better represented) but still has issues. Very psychedelic though!

[zssherman]

Hmm, yeah @rcjackson I wonder if the nearest neighbor is still having issues. We can try to reproduce on our end. @lauratomkins Are you only running on this one case? If so, you could use the parameter in grid_from_radars, gridding_algo='map_to_grid'. This is an unoptimized version of the grid algos so its kinda slow, but it works. If your doing multiple cases, that might be too slow for that, but if its just the one, it should help for now, until we come up with a solution

[lauratomkins]

Hi @zssherman, I updated the gridding_algo parameter, but no difference in the end result. I also plan to process hundreds of cases so not really an option for me right now.

[zssherman]

Gotcha, we'll be taken a look at the code to see if we can figure out whats going on.

[rcjackson]

So, looking at your data, from what I recall when NEXRAD is in dual PRF mode it will store the same tilt in 2 different sweeps. One tilt will not have dual PRF on, so we only get reflectivity and rhohv, and the other tilt will have dual PRF on is meant for velocity and spectral width. Therefore, you need to be using the first tilt for reflectivity and rhohv, and the second for velocity and spectral width. When I just use tilts separately for gridding, I get fields that make much more sense.

[lauratomkins]

@rcjackson Thanks for your comment, I'm getting similar results as well with this workaround. Not super ideal but hopefully should work for what I'm doing.

[lauratomkins]

Hi @zssherman and @rcjackson. I'm wanted to document some additional issues that I'm having with the nearest neighbor interpolation in the grid_from_radars function. Since I initially posted my question, I separated the sweeps so I'm now only interpolating one sweep (which fixed my initial issue, Thanks!). But now I'm having some slightly different issues.

For reference, here is the dealiased velocity in polar coordinates (zoomed in on the right)

Here is what the gridded dealiased velocity looks like using the default gridding algorithm. The main issues are the really low values in the velocity (close to -5000) and the missing data (which on closer look seems slightly reminiscent of the pattern in the previous plots in this thread).

Then, here is what the field looks like when I change the gridding algorithm to'map_to_grid. It fixes both issues but is too slow for how much data I need to process.

Any thoughts on this? I tried with the updated pyart package from the source instead of conda-forge, but no luck either. Any insight would be greatly appreciate.

[lauratomkins]

The -5000 m/s values show up *after* gridding. Before gridding the lowest value is ~-60 m/s so this seems to be an error that shows up during the interpolation process.

[zssherman]

@lauratomkins Yeah after discussion with @rcjackson seems there might be an issue on how the mask is being applied. We will be taking a look at it.

[millercommamatt]

Unless I'm mistaken, it looks like the gates values are being mapped to a grid instead of a grid being populated with the nearest gate values. This would explain the empty grid boxes in the data shown by @lauratomkins. Basically, no data are being mapped to those empty grid boxes since there are other grid boxes closer to the surrounding gates. I'm not sure what's going on where the values are really low.

[rcjackson]

Do you have the data file for this scan? I was unable to reproduce your behavior even after dealiasing the file you gave me and setting the grid spacing to 2 km. Rather, I get a much more reasonable CAPPI: Going by your plot, I think your grid resolution looks to be >> 2 km.

[lauratomkins]

Here is the file. The second example I've shown with the dealiased velocity is a different file/radar/time compared to the first example. Sorry for not clarifying that.

And I'm actually using a grid spacing of 500 m, that's a typo on my figures above.

[rcjackson]

I still am unable to reproduce your error with map_gates_to_grid at a 500 m scale. When I use this grid specification, I get a very reasonable looking plot.

grid_nearestn = pyart.map.grid_from_radars( (radar,), grid_shape=(1,801,801), griding_algo='map_gates_to_grid', grid_limits=((0, 9000), (-200000, 200000), (-200000, 200000)), grid_projection={'proj':'lcca', 'lat_0':center_lat, 'lon_0':center_lon}, roi_func='constant_roi', constant_roi=2000.0, weighting_function='Nearest') # nearest neightbor

What is the exact code you are using to make the grid? You may need to change the ROI since the nearest neighbor algorithm looks for the closest gate within a certain ROI.

[lauratomkins]

This is the code I'm using. I tried adding in the roi_func='constant_roi' and constant_roi=2000.0 but no luck there either. Are you using the most up to date pyart version? I'm still using the latest conda-forge install.

center_lat, center_lon, _ = pyart.io.nexrad_common.get_nexrad_location('KDIX') # gets central radar lat, lon, altitude

# define grid for interpolating
dx = 0.5; dy = 0.5
xpts = (dx*1000)*np.arange(-(150/dx),(150/dx)+1,1); ypts = (dy*1000)*np.arange(-(150/dy),(150/dy)+1,1); zpts = np.zeros([1,len(xpts)]); # array for interpolation

grid = pyart.map.grid_from_radars(
        (radar,), grid_origin = (center_lat, center_lon), grid_shape=(1, len(ypts), len(xpts)),
        grid_limits=((0, 10000), (np.min(ypts), np.max(ypts)), (np.min(xpts), np.max(xpts))), 
        grid_projection = {'proj':'lcca', 'lat_0': center_lat,'lon_0': center_lon},
        weighting_function = 'Nearest')

[rcjackson]

You may want to install pyart from source. We have made a few fixes to the nearest neighbor code since the conda forge release.
See Issue #942.

To install pyart from the latest build, you'll want to uninstall your conda pyart and then do:

pip install git+https://github.com/ARM-DOE/pyart

This should include this latest fix.

[zssherman]

And I'll get a new release out shortly.

[zssherman]

A new release is now out, I would give it 20 minutes or so for conda-forge to update and have the new downloads.