2011-03-16 85 views
2

我有一个在矩形网格(X,Y)上有大约500000个元素的ndarray(Z)。大矩阵的SciPy插值

现在我想在x,y中的大约100个位置插值,这些位置不一定在网格上。

我有一些代码在Matlab工作:

data = interp2(X,Y,Z, x,y); 

然而,当我尝试使用与scipy.interpolate同样的方法,我得到取决于方法的各种错误。例如,如果我指定kind = 'linear'和“溢出错误:太多的数据点进行插值”,如果我指定了kind='cubic',则interp2d将失败并显示MemoryError。我也试过Rbfbisplev,但他们也失败了。

所以问题是:是否有插值函数允许插值大矩阵?有没有解决问题的另一种方法? (或者我必须编码一个函数,选择点周围的适当区域来插入并调用interp2d?)

另外:如何用复数来做到这一点?

+0

小心显示您的代码? 500000不是那么大。谢谢 – eat 2011-03-16 16:18:32

回答

2

由于您的数据位于网格上,因此您可以使用RectBivariateSpline

要处理复杂数字,可以分别插入data.realdata.imag(FITPACK例程IIRC不处理复杂数据)。

1

编辑:哎呦。刚刚意识到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]) 
+0

感谢您的示例代码! – Joma 2011-04-05 15:18:01