大矩阵的 SciPy 插值答案

【问题标题】：SciPy interpolation of large matrix大矩阵的 SciPy 插值
【发布时间】：2011-07-16 17:37:41
【问题描述】：

我有一个 ndarray (Z)，在矩形网格 (X, Y) 上有大约 500000 个元素。

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

我有一些代码在 Matlab 中工作：

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

但是，当我尝试对 scipy.interpolate 使用相同的方法时，我会收到各种错误，具体取决于方法。例如，如果我指定 kind = 'linear'，interp2d 会因 MemoryError 而失败，如果我指定 kind='cubic'，则出现“OverflowError: Too many data points to interpolate”。我也尝试过Rbf 和bisplev，但它们也失败了。

所以问题是：是否有允许对大型矩阵进行插值的插值函数？是否有其他解决方案？（或者我是否必须编写一个函数来选择点周围的合适区域进行插值，然后调用 interp2d？）

另外：如何处理复数？

【问题讨论】：

想展示您的代码吗？ 500000并不是那么大。谢谢

标签： python matlab scipy interpolation

【解决方案1】：

由于您的数据在网格上，您可以使用RectBivariateSpline。

要处理复数，您可以分别插入 data.real 和 data.imag（FITPACK 例程 IIRC 不处理复数数据）。

【讨论】：

【解决方案2】：

编辑： 哎呀。刚刚意识到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])

【讨论】：

【解决方案3】：

当通过传递一组网格坐标进行初始化时，在相对较大的数据集上构建 scipy interp2d 插值器可能需要很长时间。如果数据位于矩形网格上，您可以考虑另一种初始化 interp2d 的方法：

一）

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 也很好用。

【讨论】：