使用`scipy.interpolate.griddata'进行插值的速度非常慢

Question 1

当我试图将 "几乎 "有规律的网格化数据插值到地图坐标时，我遇到了 scipy.interpolate.griddata 的极度缓慢的性能，以便用 matplotlib.pyplot.imshow 绘制地图和数据。因为 matplotlib.pyplot.pcolormesh 花费的时间太长，而且与 alpha 的表现也不太一样。

最好展示一个例子（输入文件可以下载 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这不能通过简单地用仿生变换旋转数据来实现。
在我的真实数据中，griddata每次通话需要80多秒，pcolormesh需要的时间更长（超过2分钟！）。我已经看了Jaimi的答案here和乔-金顿的回答here但我想不出让它为我工作的方法。
我的所有数据集都有完全相同的lons、lats，所以基本上我需要把这些东西映射到地图的坐标上，然后对数据本身应用同样的转换。问题是我如何做到这一点？

Question 2


          
           
            在长时间忍受了
            
             scipy.interpolate.griddata
            
            的极度缓慢的性能之后我决定放弃
            
             griddata
            
            ，改用
            
             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)
            
            进行图像转换。
            
             OpenCV
            
            .具体而言。
            
             观点的转变
            
            .
           
           
            因此，对于上面的例子，也就是上面问题中的例子，你可以得到输入文件
            
             here
            
            ，这是一段需要1.1毫秒的代码，而上面的例子中重新网格化部分需要692毫秒。
           
           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)
我认为这个解决方案的唯一缺点是数据的轻微偏移，如附图中不重叠的等值线所示。黑色的重测数据等值线（可能更准确）和'喷射'色标的warpPerspective数据等值线。
目前，我在性能方面的优势是可以接受这种差异的，我希望这个解决方案也能帮助其他人。
有人（不是我......）应该找到一种方法来提高griddata的性能 :)

Question 3


          
           
            
             我使用了numpy
             
              ndimage.map_coordinates
             
             。它工作得很好!
            
            
             
              http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html
             
            
            
             
              从上述链接中复制的。
             
            
            scipy.ndimage.interpolation.map_coordinates(input, coordinates, output=None, order=3, mode='constant', cval=0.0, prefilter=True)
            
             通过内插法将输入数组映射到新的坐标。
            
            
             坐标数组被用来为输出中的每个点找到输入中的相应坐标。在这些坐标上的输入值是由所要求的花样插值决定的。
            
            
             输出的形状是通过放弃第一轴而从坐标阵列的形状得到的。数组中沿第一轴的值是输入数组中的坐标，在这些坐标上可以找到输出值。
            
                from scipy import ndimage
    a = np.arange(12.).reshape((4, 3))
    array([[  0.,   1.,   2.],
            [  3.,   4.,   5.],