Very Slow Interpolation Using `scipy.interpolate.griddata`

I am experiencing excruciatingly slow performance of scipy.interpolate.griddata when trying to interpolate 'almost' regularly gridded data to map coordinates so that both map and d

Solution 1:

After a long time of putting up with excruciatingly slow performance of scipy.interpolate.griddata I've decided to ditch griddata in favor of image transformation with OpenCV. Specifically, Perspective Transformation.

So for the above example, the one in the question above, for which you can get the input files here, this is a piece of code which takes 1.1 ms as opposed to the 692 ms which takes the regridding part in the example above.

import cv2
new_data = data.T[::-1]

# calculate the pixel coordinates of the
# computational domain corners in the data array
w,e,s,n = map_extent
dx = float(e-w)/new_data.shape[1]
dy = float(n-s)/new_data.shape[0]
x = (lon.ravel()-w)/dx
y = (n-lat.ravel())/dy

computational_domain_corners = np.float32(zip(x,y))

data_array_corners = np.float32([[0,new_data.shape[0]],
                                 [0,0],
                                 [new_data.shape[1],new_data.shape[0]],
                                 [new_data.shape[1],0]])

# Compute the transformation matrix which places
# the corners of the data array at the corners of
# the computational domain in data array pixel coordinates
tranformation_matrix = cv2.getPerspectiveTransform(data_array_corners,
                                                   computational_domain_corners)

# Make the transformation making the final array the same shape
# as the data array, cubic interpolate the data placing NaN's
# outside the new array geometry
mapped_data = cv2.warpPerspective(new_data,tranformation_matrix,
                                  (new_data.shape[1],new_data.shape[0]),
                                  flags=2,
                                  borderMode=0,
                                  borderValue=np.nan)

The only drawback I see to this solution is a slight offset in the data as illustrated by the non-overlapping contours in the attached image. regridded data contours (probably more accurate) in black and warpPerspective data contours in 'jet' colorscale.

At the moment, I live just fine with the discrepancy at the advantage of performance and I hope this solution helps others as-well.

Someone (not me...) should find a way to improve the performance of griddata :) Enjoy!

Solution 2:

I used numpy ndimage.map_coordinates. It worked great!

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html

copied from the above link:

scipy.ndimage.interpolation.map_coordinates(input, coordinates, output=None, order=3, mode='constant', cval=0.0, prefilter=True)

Map the input array to new coordinates by interpolation.

The array of coordinates is used to find, for each point in the output, the corresponding coordinates in the input. The value of the input at those coordinates is determined by spline interpolation of the requested order.

The shape of the output is derived from that of the coordinate array by dropping the first axis. The values of the array along the first axis are the coordinates in the input array at which the output value is found.

    from scipy import ndimage
    a = np.arange(12.).reshape((4, 3))
    a
    array([[  0.,   1.,   2.],
            [  3.,   4.,   5.],
            [  6.,   7.,   8.],
            [  9.,  10.,  11.]])
    ndimage.map_coordinates(a, [[0.5, 2], [0.5, 1]], order=1)
    [ 2.  7.]