Very slow interpolation using `scipy.interpolate.griddata`

Question

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 data can be plotted with matplotlib.pyplot.imshow because matplotlib.pyplot.pcolormesh is taking too long and does not behave well with alpha among other things.

Best show an example (input files can be downloaded here):

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata

map_extent = (34.4, 36.2, 30.6, 33.4)
# data corners:
lon = np.array([[34.5,        34.83806236],
                [35.74547079, 36.1173923]])
lat = np.array([[30.8,        33.29936152],
                [30.67890411, 33.17826563]])

# load saved files
topo = np.load('topo.npy')
lons = np.load('lons.npy')
lats = np.load('lats.npy')
data = np.load('data.npy')

# get max res of data
dlon = abs(np.array(np.gradient(lons))).max()
dlat = abs(np.array(np.gradient(lats))).max()

# interpolate the data to the extent of the map
loni,lati = np.meshgrid(np.arange(map_extent[0], map_extent[1]+dlon, dlon),
                        np.arange(map_extent[2], map_extent[3]+dlat, dlat))
zi = griddata((lons.flatten(),lats.flatten()),
              data.flatten(), (loni,lati), method='linear')

Plotting:

fig, (ax1,ax2) = plt.subplots(1,2)
ax1.axis(map_extent)
ax1.imshow(topo,extent=extent,cmap='Greys')

ax2.axis(map_extent)
ax2.imshow(topo,extent=extent,cmap='Greys')

ax1.imshow(zi, vmax=0.1, extent=extent, alpha=0.5, origin='lower')
ax1.plot(lon[0],lat[0], '--k', lw=3, zorder=10)
ax1.plot(lon[-1],lat[-1], '--k', lw=3, zorder=10)
ax1.plot(lon.T[0],lat.T[0], '--k', lw=3, zorder=10)
ax1.plot(lon.T[-1],lat.T[-1], '--k', lw=3, zorder=10)


ax2.pcolormesh(lons,lats,data, alpha=0.5)
ax2.plot(lon[0],lat[0], '--k', lw=3, zorder=10)
ax2.plot(lon[-1],lat[-1], '--k', lw=3, zorder=10)
ax2.plot(lon.T[0],lat.T[0], '--k', lw=3, zorder=10)
ax2.plot(lon.T[-1],lat.T[-1], '--k', lw=3, zorder=10)

Result:

Note, this can not be done by simply rotating the data with affine transformations.

The griddata call takes over 80 seconds per call with my real data and pcolormesh takes even longer (over 2 minutes!). I have looked at both Jaimi's answer here and Joe Kington's answer here but I cant figure out a way to get it to work for me.

All my datasets have exactly the same lons, lats so basically I need to map those once to the map's coordinates and apply the same transformation to the data itself. Question is how do I do that?

Answer 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!

Answer 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.]