编辑:哎呦。刚刚意识到OP在问题中提出了这个解决方案!
我不知道为什么插值例程花费很多时间和内存来查找结构化数据的节点,但由于您只使用整个网格的一小部分,因此可以将插值分解为多个补丁让事情更有效率。
from scipy import interpolate
import numpy as np
def my_interp(X, Y, Z, x, y, spn=3):
xs,ys = map(np.array,(x,y))
z = np.zeros(xs.shape)
for i,(x,y) in enumerate(zip(xs,ys)):
# get the indices of the nearest x,y
xi = np.argmin(np.abs(X[0,:]-x))
yi = np.argmin(np.abs(Y[:,0]-y))
xlo = max(xi-spn, 0)
ylo = max(yi-spn, 0)
xhi = min(xi+spn, X[0,:].size)
yhi = min(yi+spn, Y[:,0].size)
# make slices of X,Y,Z that are only a few items wide
nX = X[xlo:xhi, ylo:yhi]
nY = Y[xlo:xhi, ylo:yhi]
nZ = Z[xlo:xhi, ylo:yhi]
intp = interpolate.interp2d(nX, nY, nZ)
z[i] = intp(x,y)[0]
return z
N = 1000
X,Y = np.meshgrid(np.arange(N), np.arange(N))
Z = np.random.random((N, N))
print my_interp(X, Y, Z, [13.2, 999.9], [0.01, 45.3])
小心显示您的代码? 500000不是那么大。谢谢 – eat 2011-03-16 16:18:32