python中griddata的外插值,使用`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 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?

解决方案

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!