问题描述：

这篇问答旨在成为使用 scipy 进行二维（和多维）插值的规范问答。关于各种多维插值方法的基本语法，经常会有人提出问题，我希望能够解答这些问题。

我有一组散布的二维数据点，我想将它们绘制成一个漂亮的曲面，最好使用类似contourf或plot_surface之类的函数matplotlib.pyplot。如何使用 scipy 将二维或多维数据插值到网格中？

我找到了子包，但在使用或或或（或旧版）scipy.interpolate时总是出错。这些方法的正确语法是什么？interp2d`bisplrepgriddataRBFInterpolator`Rbf

解决方案 1：

免责声明：我撰写本文主要考虑的是语法方面的考虑和一般行为。我对所述方法的内存和 CPU 方面不太熟悉，因此我的目标是针对那些拥有相对较少数据集的用户，因为插值的质量可能是主要考虑因素。我知道，在处理非常大的数据集时，性能更好的方法（即griddata不RBFInterpolator带neighbors关键字参数的方法）可能并不可行。

请注意，此答案使用了1.7.0RBFInterpolator中引入的新类。有关旧类，请参阅此答案的先前版本。SciPy`Rbf`

我将比较三种多维插值方法（interp2d/splines、griddata和RBFInterpolator）。我将对它们进行两种插值任务和两种底层函数（要从中进行插值的点）。具体示例将演示二维插值，但可行的方法适用于任意维度。每种方法都提供各种插值；在所有情况下，我都会使用三次插值（或接近1 的插值）。需要注意的是，无论何时使用插值，都会引入与原始数据相比的偏差，而所使用的具体方法会影响最终得到的工件。请始终注意这一点，并负责任地进行插值。

两个插值任务将是

1. 上采样（输入数据在矩形网格上，输出数据在更密集的网格上）

2. 将散点数据插值到规则网格上

这两个函数（在域上[x, y] in [-1, 1]x[-1, 1]）将是

1. 流畅友好的功能：cos(pi*x)*sin(pi*y)；范围[-1, 1]

2. 一个邪恶的（特别是非连续的）函数：x*y / (x^2 + y^2)在原点附近值为 0.5；范围在[-0.5, 0.5]

它们的外观如下：

我将首先演示这三种方法在这四个测试下的表现，然后详细介绍这三种方法的语法。如果你已经了解某个方法的用途，那么你可能就不想浪费时间学习它的语法了（说的就是你interp2d）。

测试数据

为了清晰起见，以下是我生成输入数据的代码。虽然在这个特定案例中，我显然知道数据背后的函数，但我只会用它来生成插值方法的输入。我使用 numpy 是为了方便（主要是为了生成数据），但单独使用 scipy 也足够了。

import numpy as np
import scipy.interpolate as interp

# auxiliary function for mesh generation
def gimme_mesh(n):
    minval = -1
    maxval =  1
    # produce an asymmetric shape in order to catch issues with transpositions
    return np.meshgrid(np.linspace(minval, maxval, n),
                       np.linspace(minval, maxval, n + 1))

# set up underlying test functions, vectorized
def fun_smooth(x, y):
    return np.cos(np.pi*x) * np.sin(np.pi*y)

def fun_evil(x, y):
    # watch out for singular origin; function has no unique limit there
    return np.where(x**2 + y**2 > 1e-10, x*y/(x**2+y**2), 0.5)

# sparse input mesh, 6x7 in shape
N_sparse = 6
x_sparse, y_sparse = gimme_mesh(N_sparse)
z_sparse_smooth = fun_smooth(x_sparse, y_sparse)
z_sparse_evil = fun_evil(x_sparse, y_sparse)

# scattered input points, 10^2 altogether (shape (100,))
N_scattered = 10
rng = np.random.default_rng()
x_scattered, y_scattered = rng.random((2, N_scattered**2))*2 - 1
z_scattered_smooth = fun_smooth(x_scattered, y_scattered)
z_scattered_evil = fun_evil(x_scattered, y_scattered)

# dense output mesh, 20x21 in shape
N_dense = 20
x_dense, y_dense = gimme_mesh(N_dense)

平滑函数和上采样

让我们从最简单的任务开始。以下是平滑测试函数中从 形状网格[6, 7]到 形状网格的上采样过程：[20, 21]

尽管这是一项简单的任务，但输出之间已经存在细微的差别。乍一看，所有三个输出都是合理的。根据我们对底层函数的先验知识，有两个特点需要注意：中间的情况对griddata数据的扭曲最大。请注意y == -1图的边界（最靠近x标签）：函数应该严格为零（因为y == -1是平滑函数的节点线），但情况并非如此griddata。还要注意x == -1图的边界（后方，左侧）：底层函数在处有一个局部最大值（意味着边界附近的梯度为零）[-1, -0.5]，但griddata输出在此区域明显显示非零梯度。这种影响很微妙，但它仍然是一种偏差。

邪恶函数和上采样

稍微困难一点的任务是对我们的邪恶函数执行上采样：

这三种方法之间的差异开始显现。观察表面图，在输出中出现了明显的伪极值interp2d（注意绘制表面右侧的两个凸起）。虽然griddata和RBFInterpolator乍一看似乎产生了相似的结果，但[0.4, -0.4]在底层函数中不存在的局部最小值附近产生了局部最小值。

然而，它有一个关键方面RBFInterpolator远胜于其他方法：它尊重底层函数的对称性（当然，这也得益于样本网格的对称性）。 的输出griddata破坏了样本点的对称性，这在平滑情况下已经相当明显。

平滑函数和散乱数据

人们通常希望对散点数据进行插值。因此，我认为这些测试更为重要。如上所示，样本点是在目标域内伪均匀选取的。在实际情况下，每次测量都可能存在额外的噪声，因此您应该首先考虑对原始数据进行插值是否合理。

平滑函数的输出：

现在已经开始有点恐怖了。我专门将输出从interp2d到 之间的部分剪切[-1, 1]用于绘图，以便至少保留少量信息。很明显，虽然一些底层形状仍然存在，但仍存在大量噪声区域，该方法在这些区域中完全失效。第二种情况griddata相当好地再现了形状，但请注意轮廓图边界处的白色区域。这是因为griddata仅在输入数据点的凸包内起作用（换句话说，它不执行任何外推）。我保留了位于凸包外部的输出点的默认 NaN 值。2考虑到这些特性，RBFInterpolator似乎表现最佳。

邪恶函数和散乱数据

我们一直在等待的时刻到了：

放弃并不令人意外interp2d。事实上，在调用过程中，interp2d你应该预料到会有人友好地RuntimeWarning抱怨无法构造样条曲线。至于其他两种方法，RBFInterpolator即使在结果外推域的边界附近，它们似乎也能产生最佳输出。

因此，让我按优先顺序（最差的方法最不可能被任何人阅读）对这三种方法说几句话。

scipy.interpolate.RBFInterpolator

类名中的 RBFRBFInterpolator代表“径向基函数”。说实话，在我开始研究这篇文章之前，我从未考虑过这种方法，但我很确定将来我会用到它。

与基于样条线的方法（见后文）类似，使用方法分为两步：首先，RBFInterpolator根据输入数据创建一个可调用的类实例；然后，针对给定的输出网格调用此对象以获取插值结果。平滑上采样测试的示例：

import scipy.interpolate as interp

sparse_points = np.stack([x_sparse.ravel(), y_sparse.ravel()], -1)  # shape (N, 2) in 2d
dense_points = np.stack([x_dense.ravel(), y_dense.ravel()], -1)  # shape (N, 2) in 2d

zfun_smooth_rbf = interp.RBFInterpolator(sparse_points, z_sparse_smooth.ravel(),
                                         smoothing=0, kernel='cubic')  # explicit default smoothing=0 for interpolation
z_dense_smooth_rbf = zfun_smooth_rbf(dense_points).reshape(x_dense.shape)  # not really a function, but a callable class instance

zfun_evil_rbf = interp.RBFInterpolator(sparse_points, z_sparse_evil.ravel(),
                                       smoothing=0, kernel='cubic')  # explicit default smoothing=0 for interpolation
z_dense_evil_rbf = zfun_evil_rbf(dense_points).reshape(x_dense.shape)  # not really a function, but a callable class instance

请注意，为了让 API 顺利运行，我们必须进行一些数组构建工作RBFInterpolator。由于我们必须将二维点作为形状为 的数组传递(N, 2)，因此我们必须展平输入网格，并将两个展平后的数组堆叠起来。构造的插值器也需要这种格式的查询点，结果将是一个形状为 的一维数组，(N,)我们必须将其重新整形以匹配我们的二维网格进行绘图。由于RBFInterpolator它对输入点的维数没有任何假设，因此它支持任意维度的插值。

所以，scipy.interpolate.RBFInterpolator

• 即使输入数据疯狂，也能产生良好的输出

• 支持更高维度的插值

• 在输入点的凸包之外进行外推（当然，外推总是一种赌博，通常你不应该完全依赖它）

• 第一步是创建一个插值器，这样在各个输出点对其进行评估可以减少额外的工作量

• 可以具有任意形状的输出点数组（而不是限制为矩形网格，见下文）

• 更有可能保持输入数据的对称性

• 支持关键字的多种径向函数kernel：multiquadric，，，，，，，， （默认） 。从SciPy 1.7.0 开始inverse_multiquadric，由于技术原因，该类不允许传递自定义可调用函数，但这可能会在未来的版本中添加。inverse_quadratic`gaussianlinearcubicquinticthin_plate_spline`

• 可以通过增加smoothing参数给出不精确的插值

RBF 插值的一个缺点是插值N数据点需要对N x N矩阵求逆。这种二次复杂度会迅速增加大量数据点的内存需求。不过，新RBFInterpolator类还支持一个neighbors关键字参数，该参数将每个径向基函数的计算限制在k最近邻域内，从而减少内存需求。

scipy.interpolate.griddata

我以前最喜欢的，griddata是 ，它是用于任意维度插值的通用工具。它不会对节点凸包之外的点进行外推，除非设置一个预设值。但由于外推本身就非常不稳定且危险，所以这不一定是个缺点。使用示例：

sparse_points = np.stack([x_sparse.ravel(), y_sparse.ravel()], -1)  # shape (N, 2) in 2d
z_dense_smooth_griddata = interp.griddata(sparse_points, z_sparse_smooth.ravel(),
                                          (x_dense, y_dense), method='cubic')  # default method is linear

请注意，输入数组需要与 相同的数组变换RBFInterpolator。输入点必须指定为维度[N, D]为形状的数组D，或者指定为一维数组的元组：

z_dense_smooth_griddata = interp.griddata((x_sparse.ravel(), y_sparse.ravel()),
                                          z_sparse_smooth.ravel(), (x_dense, y_dense), method='cubic')

输出点数组可以指定为任意维度数组的元组（如以上两个片段所示），这为我们提供了更多的灵活性。

简而言之，scipy.interpolate.griddata

• 即使输入数据疯狂，也能产生良好的输出

• 支持更高维度的插值

• 不执行外推，可以为输入点凸包之外的输出设置单个值（参见fill_value）

• 在一次调用中计算插值，因此探测多组输出点从头开始

• 可以有任意形状的输出点

• 支持任意维度的最近邻插值和线性插值，以及一维和二维的三次插值。最近邻插值和线性插值分别使用NearestNDInterpolator和LinearNDInterpolator。一维三次插值使用样条函数，二维三次插值用于CloughTocher2DInterpolator构建连续可微的分段三次插值器。

• 可能会破坏输入数据的对称性

scipy.interpolate.interp2d/scipy.interpolate.bisplrep

我讨论interp2d及其相关函数的唯一原因是它的名称具有误导性，人们很可能会尝试使用它。剧透警告：请勿使用它。interp2d在 SciPy 1.10 版中已弃用，并将在 SciPy 1.12 版中移除。有关详细信息，请参阅此邮件列表讨论。它比之前的主题更特殊，因为它专门用于二维插值，但我怀疑这是迄今为止多变量插值最常见的情况。

就语法而言，它interp2d与 类似RBFInterpolator，首先需要构造一个插值实例，然后调用该实例来提供实际的插值值。然而，有一个问题：输出点必须位于矩形网格上，因此调用插值器的输入必须是跨越输出网格的一维向量，就像来自numpy.meshgrid：

# reminder: x_sparse and y_sparse are of shape [6, 7] from numpy.meshgrid
zfun_smooth_interp2d = interp.interp2d(x_sparse, y_sparse, z_sparse_smooth, kind='cubic')   # default kind is 'linear'
# reminder: x_dense and y_dense are of shape (20, 21) from numpy.meshgrid
xvec = x_dense[0,:] # 1d array of unique x values, 20 elements
yvec = y_dense[:,0] # 1d array of unique y values, 21 elements
z_dense_smooth_interp2d = zfun_smooth_interp2d(xvec, yvec)   # output is (20, 21)-shaped array

使用时最常见的错误之一interp2d是将完整的 2d 网格放入插值调用中，这会导致内存消耗激增，并有望导致仓促的MemoryError。

现在， 最大的问题是interp2d它经常不起作用。为了理解这一点，我们必须深入了解。原来，interp2d是低级函数bisplrep+的包装器bisplev，而这些低级函数又是 FITPACK 例程（用 Fortran 编写）的包装器。与上例等效的调用是

kind = 'cubic'
if kind == 'linear':
    kx = ky = 1
elif kind == 'cubic':
    kx = ky = 3
elif kind == 'quintic':
    kx = ky = 5
# bisplrep constructs a spline representation, bisplev evaluates the spline at given points
bisp_smooth = interp.bisplrep(x_sparse.ravel(), y_sparse.ravel(),
                              z_sparse_smooth.ravel(), kx=kx, ky=ky, s=0)
z_dense_smooth_bisplrep = interp.bisplev(xvec, yvec, bisp_smooth).T  # note the transpose

现在，这是关于的事情：（在 scipy 版本 1.7.0 中） for中interp2d有一个很好的注释interpolate/interpolate.py`interp2d`：

if not rectangular_grid:
    # TODO: surfit is really not meant for interpolation!
    self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)

确实interpolate/fitpack.py，有bisplrep一些设置，最终

tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
                                task, s, eps, tx, ty, nxest, nyest,
                                wrk, lwrk1, lwrk2)                 

就是这样。底层例程interp2d实际上并非用于执行插值。它们可能足以处理足够良好的数据，但在实际情况下，您可能需要使用其他方法。

总而言之，interpolate.interp2d

• 即使数据经过良好调整，也可能导致伪影

• interpn专门针对双变量问题（尽管在网格上定义的输入点有限）

• 进行外推

• 第一步是创建一个插值器，这样在各个输出点对其进行评估可以减少额外的工作量

• 只能在矩形网格上产生输出，对于分散的输出，您必须在循环中调用插值器

• 支持线性、三次和五次插值

• 可能会破坏输入数据的对称性

1我相当肯定cubic和 的linear基函数类型RBFInterpolator与其他同名插值器并不完全对应。2这些 NaN 也是导致曲面图看起来如此奇怪的原因：matplotlib 历来难以绘制具有适当深度信息的复杂三维对象。数据中的 NaN 值会使渲染器感到困惑，导致曲面中本应位于后方的部分被绘制到了前方。

这是一个可视化问题，而非插值问题。