python - 如何使用 scipy 执行二维插值？

Question

This Q&A is intended as a canonical(-ish) concerning two-dimensional (and multi-dimensional) interpolation using scipy. There are often questions concerning the basic syntax of various multidimensional interpolation methods, I hope to set these straight too.

我有一组分散的二维数据点，我想将它们绘制为一个漂亮的表面，最好使用 contourf 或 plot_surface 之类的东西 matplotlib.pyplot 。如何使用 scipy 将二维或多维数据插入网格？
我找到了 scipy.interpolate 子包，但在使用 interp2d 或 bisplrep 或 griddata 或 1x21812231343141 或 0x21812231343141 或 0x21812318 或 0x21813141 或 0x2181318 时不断出错。这些方法的正确语法是什么？

Answer 1

免责声明:我写这篇文章主要是考虑到语法和一般行为。我不熟悉所描述方法的内存和 CPU 方面，我将这个答案针对那些拥有相当小的数据集的人，因此插值的质量可以成为要考虑的主要方面。我知道在处理非常大的数据集时，性能更好的方法(即 griddata 和 RBFInterpolator 没有 neighbors 关键字参数)可能不可行。
请注意，此答案使用 RBFInterpolator 中引入的新 SciPy 1.7.0 类。对于旧的 Rbf 类，请参阅 the previous version of this answer 。
我将比较三种多维插值方法( interp2d /splines, griddata 和 RBFInterpolator )。我将让它们接受两种插值任务和两种底层函数(要插值的点)。具体例子将演示二维插值，但可行的方法适用于任意维度。每种方法都提供了各种插值；在所有情况下，我都会使用三次插值(或类似的方法)。重要的是要注意，无论何时使用插值，与原始数据相比都会引入偏差，并且使用的特定方法会影响最终得到的工件。始终注意这一点，并负责任地进行插值。
这两个插值任务将是

上采样(输入数据在矩形网格上，输出数据在更密集的网格上)

将分散数据插值到规则网格上

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

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

一个邪恶的(特别是非连续的)函数: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 是平滑函数的节点线)，但 14232318143232323232323232323232323323323232323232323232323232323232323232323232323232323323232332332323都都是还要注意图的 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 RBFInterpolator 类名称中的 RBF 代表“径向基函数”。老实说，在我开始研究这篇文章之前，我从未考虑过这种方法，但我很确定我将来会使用这些方法。
就像基于样条的方法(见下文)，使用分为两步:第一步，根据输入数据创建一个可调用的 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

请注意，我们必须进行一些数组构建操作才能使 RBFInterpolator 的 API 满意。由于我们必须将 2d 点作为形状为 (N, 2) 的数组传递，因此我们必须展平输入网格并堆叠两个展平的阵列。构造的内插器也需要这种格式的查询点，结果将是一个形状为 (N,) 的一维数组，我们必须重新整形以匹配我们的二维网格以进行绘图。由于 RBFInterpolator 不对输入点的维数做任何假设，因此它支持插值的任意维数。
所以，scipy.interpolate.RBFInterpolator

即使对于疯狂的输入数据

也会产生表现良好的输出

支持更高维度的插值

在输入点的凸包外外推(当然外推总是一场赌博，你一般根本不应该依赖它)

作为第一步创建一个插值器，因此在不同的输出点评估它是更少的额外工作

可以有任意形状的输出点阵列(与被限制为矩形网格相反，见下文)

更可能保留输入数据的对称性

支持多种径向函数关键字kernel:multiquadric，inverse_multiquadric，inverse_quadratic，gaussian，linear，cubic，quintic，thin_plate_spline(缺省值)。从 SciPy 1.7.0 开始，由于技术原因，该类不允许传递自定义可调用项，但这可能会在 future 版本中添加。

可以通过增加 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 )

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

可以有任意形状的输出点

支持任意维度的最近邻和线性插值，1d 和 2d 的立方。最近邻和线性插值分别使用 NearestNDInterpolator 和 LinearNDInterpolator。 1d 三次插值使用样条，2d 三次插值使用 CloughTocher2DInterpolator 来构造一个连续可微的分段三次插值器。

可能违反输入数据的对称性

scipy.interpolate.interp2d / scipy.interpolate.bisplrep
我讨论 interp2d 及其亲属的唯一原因是它的名称具有欺骗性，人们可能会尝试使用它。剧透警告:不要使用它(从 scipy 1.7.0 版开始)。它已经比之前的主题更加特殊，因为它专门用于二维插值，但我怀疑这是迄今为止多元插值最常见的情况。
就语法而言，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

现在，这是关于 interp2d 的事情:(在 scipy 版本 1.7.0 中)有一个很好的 comment in 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 历来难以用适当的深度信息绘制复杂的 3d 对象。数据中的 NaN 值会混淆渲染器，因此应该在后面的表面部分被绘制在前面。这是可视化的问题，而不是插值。