SciPy大矩阵插值

Question

我有一个包含大约500000个元素的ndarray（Z），这些元素在一个矩形网格（X, Y）上。

现在我想在一些不一定在这个网格上的位置（大约100个位置）进行插值，也就是想要在这些点上估算出值。

我在Matlab里有一些可以用的代码：

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

但是，当我尝试用scipy.interpolate来做同样的事情时，会遇到各种错误，具体取决于我使用的方法。例如，如果我指定kind = 'linear'，就会出现内存错误（MemoryError）；如果我指定kind='cubic'，则会出现“数据点太多，无法插值”的溢出错误（OverflowError）。我还尝试了Rbf和bisplev，但它们也都失败了。

所以我的问题是：有没有一种插值函数可以处理大矩阵的插值？有没有其他解决这个问题的方法？（或者我是否需要编写一个函数，选择合适的区域来进行插值，然后再调用interp2d？）

另外：如何处理复数的插值？

Answer 1

在使用scipy的interp2d进行插值时，如果数据集比较大，通过传递一组网格坐标来初始化可能会花费很长时间。

如果你的数据是在一个矩形网格上，可以考虑另一种初始化interp2d的方法：

a)

from scipy.interpolate import interp2d    
x = [0,1,2]
y = [0,3]
z = [[1,2,3], [4,5,6]]
i = interp2d(x, y, z)
i(0, 0)[0]

而不是

b)

from scipy.interpolate import interp2d 
x = [0, 1, 2, 0, 1, 2]
y = [0, 0, 0, 3, 3, 3]
z = [1, 2, 3, 4, 5, 6]
i = interp2d(x, y, z)
i(0, 0)[0]

这个在interp2d的实现中是有考虑到的。方法a)启动得明显更快，但只适合矩形网格。我在一个有227000个点的网格上应用这个技巧后，性能从6分钟提升到了3秒。

RectBivariateSpline也表现得很好。

Answer 2

编辑：哎呀，刚发现提问者在问题中已经提到这个解决方案了！

我不知道为什么插值的过程需要花费这么多时间和内存来找到结构化数据的节点，但因为你只是在使用整个网格的一小部分，你可以把插值过程分成几个小块，这样会更高效。

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])

Answer 3

因为你的数据是在一个网格上，所以你可以使用 RectBivariateSpline 这个工具。

如果你要处理复杂的数字，可以分别对 data.real 和 data.imag 进行插值（我记得FITPACK的程序不支持复杂数据）。