reduce_data bug with polar grids

Hello,

This is a bug I have been chasing for a while, but I have the impression it somehow is worse now (1.5.0) than before. It happens when using kd_tree.resample_* and reduce_data=True (which is the default) with polar orbiting satellites, on polar grids. At least the polar grids I am using for my work.

I prepared a script to generate test data and compare the remapping with and without reduce_data. It results in the image below.

To me, the reduce_data feature should only speed-up the resampling, not change the result in any way. So this is a serious flaw.

I hope this can be made rather easily into a unit test for reduce_data since it is self contained (no need for external data files to be loaded).

Where should we go from there? The github tool will not let me attach a .py. I can send it by mail. I try to paste it here but it does not look good.

import numpy as np
import pyresample as pr
import matplotlib.pylab as plt
from mpl_toolkits.basemap import Basemap
import pyproj

print "Pyresample {}".format(pr.__version__,)

# routine to create a polar area_def (NSIDC's 25km SSM/I grid)
def get_areadef(hemis):
    nx, ny = {'nh':(448,304), 'sh':(332,316)}[hemis][::-1]
    a = 25000
    if hemis == 'nh':
        xl = -3850000.00
        xL = xl+0.5*a + (nx-1)*a + 0.5*a
        yL = 5850000.00
        yl = yL-0.5*a - (ny-1)*a - 0.5*a
        pdict = {u'a': u'6378273', u'b': u'6356889.44891', u'lat_ts': u'70', u'lon_0': u'-45', u'proj': u'stere', u'lat_0': u'90', u'units':u'm'}
    elif hemis == 'sh':
        xl = -3950000.00
        xL = xl+0.5*a + (nx-1)*a + 0.5*a
        yL = +4350000.00
        yl = yL-0.5*a - (ny-1)*a - 0.5*a
        pdict = {u'a': u'6378273', u'b': u'6356889.44891', u'lat_ts': u'-70', u'lon_0': u'0', u'proj': u'stere', u'lat_0': u'-90', u'units':u'm'}
    else:
        raise ValueError("Not a valid area ('nh' and 'sh' are valid)")

    area_def  = pr.geometry.AreaDefinition('nsidc_ps_{}'.format(hemis,),
                                          'nsidc_ps_{}'.format(hemis,),
                                          'nsidc_ps_{}'.format(hemis,),
                                          pdict,
                                          nx, ny,
                                          [xl, yl, xL, yL])
    if area_def.pixel_size_x != 25000 or area_def.pixel_size_y != 25000:
        raise ValueError("The grid parameters do not match a grid spacing of 25km! (x:{} y:{})".format(area_def.pixel_size_x, area_def.pixel_size_y))
    
    return area_def


area_def = get_areadef('nh')

#simulate swath data:
#  1) have a single line of data going through (North) pole
slon = 125.
swath_lats_1 = np.arange(30.,90.,1.)
swath_lons_1 = np.array([slon,]*swath_lats_1.size)
swath_lats_2 = swath_lats_1[::-1]
swath_lons_2 = np.array([180+slon,]*swath_lats_2.size)

swath_lats = np.concatenate((swath_lats_1,[90.,],swath_lats_2))
swath_lons = np.concatenate((swath_lons_1,[0.,],swath_lons_2))

#  2) use pyproj to obtain x,y coordinates in the map
P = pyproj.Proj(area_def.proj4_string)
x,y = P(swath_lons, swath_lats)

#  3) use pyproj(inverse=True) to get 5 neighbouring lines that are side by side and not going through Pole. This is our simulated polar orbiting swath
nbscanpos = 5
lons = np.empty((swath_lons.size,nbscanpos))
lats = np.empty((swath_lats.size,nbscanpos))
data = np.empty((swath_lats.size,nbscanpos))
off = 100*1000  # be at least 100km from Pole
for i in range(0,nbscanpos):
    xp = x + off
    yp = y - off
    lons[:,i], lats[:,i] = P(xp,yp,inverse=True)
    data[:,i] = i+5
    off += 50*1000 # distance between scan positions is 50km


# Do 2 kdtree.resample() calls: exactly the same except for reduce_data=False|True (default is True)
swath_def = pr.geometry.SwathDefinition(lons, lats)
res_rdF = pr.kd_tree.resample_nearest(swath_def,data,area_def,radius_of_influence=50000,reduce_data=False)
res_rdT = pr.kd_tree.resample_nearest(swath_def,data,area_def,radius_of_influence=50000,reduce_data=True)

#plot
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(121)
ax.set_title('kdtree data_reduce=False')
bmap = pr.plot.area_def2basemap(area_def)
bmap.drawcoastlines()
bmap.imshow(res_rdF,interpolation='none',cmap=plt.cm.plasma, origin='upper')

ax = fig.add_subplot(122)
ax.set_title('kdtree data_reduce=True')
bmap = pr.plot.area_def2basemap(area_def)
bmap.drawcoastlines()
bmap.imshow(res_rdT,interpolation='none',cmap=plt.cm.plasma, origin='upper')


plt.savefig('./pyresample_reduce_data_bug.png',bbox_inches='tight')
print "./pyresample_reduce_data_bug.png ready"

plt.show()

Start by running the self-contained Python reproduction from the issue and compare kd_tree.resample_nearest with reduce_data=False and reduce_data=True on the polar grid. Trace the kd_tree resampling entry point to determine why the results differ. Done means the behavior is covered by a regression test and reduce_data changes performance without changing the result.