免责声明: 我写这篇文章主要是考虑到句法和一般行为。我不熟悉所描述的方法的内存和 CPU 方面，我的目标是那些拥有相当小的数据集的人，因此插值的质量可以成为主要考虑的方面。我知道，当处理非常大的数据集时，性能更好的方法(即没有 neighbors关键字参数的 griddata和 RBFInterpolator)可能是不可行的。

注意，这个答案使用了 SciPy1.7.0中引入的新的 RBFInterpolator类。

I'm going to compare three kinds of multi-dimensional interpolation methods (interp2d/splines, griddata and RBFInterpolator). I will subject them to two kinds of interpolation tasks and two kinds of underlying functions (points from which are to be interpolated). The specific examples will demonstrate two-dimensional interpolation, but the viable methods are applicable in arbitrary dimensions. Each method provides various kinds of interpolation; in all cases I will use cubic interpolation (or something close¹). It's important to note that whenever you use interpolation you introduce bias compared to your raw data, and the specific methods used affect the artifacts that you will end up with. Always be aware of this, and interpolate responsibly.

两个插值任务将是

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 是为了方便(主要是为了生成数据) ，但是只使用 scypy 也足够了。

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输出在这个区域显示出明显的非零梯度。这种影响是微妙的，但它仍然是一种偏见。

邪恶的功能和升级

A bit harder task is to perform upsampling on our evil function:

这三种方法之间开始显示出明显的差异。观察表面图，在 interp2d的输出中出现了明显的伪极值(注意在绘制表面的右侧有两个驼峰)。虽然 griddata和 RBFInterpolator乍一看似乎产生了相似的结果，在 [0.4, -0.4]附近产生了局部极小值，而这在基础功能中是不存在的。

然而，有一个至关重要的方面，RBFInterpolator远远优越: 它尊重底层函数的对称性(这当然也是由样本网格的对称性使之成为可能)。griddata的输出打破了采样点的对称性，在光滑的情况下，这种对称性已经很弱了。

平滑函数和散乱数据

大多数情况下，人们希望对散乱的数据进行插值。因此，我认为这些测试更加重要。如上所示，样本点是在感兴趣的领域选择伪一致。在现实的场景中，每次测量都可能产生额外的干扰，您应该考虑从一开始就插入原始数据是否有意义。

Output for the smooth function:

现在已经有点恐怖了。我将 interp2d的输出剪切到 [-1, 1]之间，专门用于绘图，以便至少保留最少量的信息。很明显，虽然一些潜在的形状是存在的，有巨大的噪声区域的方法完全崩溃。第二种情况下的 griddata复制的形状相当不错，但注意在轮廓曲线边界的白色区域。这是因为 griddata只在输入数据点的凸包内工作(换句话说，它不执行任何 extrapolation)。我保留了凸壳外输出点的默认 NaN 值。考虑到这些特性，RBFInterpolator似乎性能最好。

邪恶的功能和分散的数据

And the moment we've all been waiting for:

interp2d放弃也就不足为奇了。事实上，在对 interp2d的调用期间，你应该期待一些友好的 RuntimeWarning抱怨不可能构造样条。至于其他两种方法，RBFInterpolator似乎产生了最好的输出，甚至在结果被外推的域的边界附近。

因此，让我说几句关于这三种方法的话，按照偏好的减少顺序(因此，最糟糕的是最不可能被任何人阅读)。

scipy.interpolate.RBFInterpolator

The RBF in the name of the RBFInterpolator class stands for "radial basis functions". To be honest I've never considered this approach until I started researching for this post, but I'm pretty sure I'll be using these in the future.

Just like the spline-based methods (see later), usage comes in two steps: first one creates a callable RBFInterpolator class instance based on the input data, and then calls this object for a given output mesh to obtain the interpolated result. Example from the smooth upsampling test:

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,)的1d 数组，我们必须将其重新形状以匹配绘图的2d 网格。由于 RBFInterpolator对输入点的维数不做任何假设，它支持任意维数的插值。

那么，scipy.interpolate.RBFInterpolator

• produces well-behaved output even for crazy input data
• 支持高维插值
• 外推输入点的凸外壳(当然外推总是一个赌博，你一般不应该依赖它)
• 创建一个插值器作为第一步，因此在不同的输出点评估它是较少的额外工作
• can have output point arrays of arbitrary shape (as opposed to being constrained to rectangular meshes, see later)
• 更有可能保持输入数据的对称性
• supports multiple kinds of radial functions for keyword kernel: multiquadric, inverse_multiquadric, inverse_quadratic, gaussian, linear, cubic, quintic, thin_plate_spline (the default). As of SciPy 1.7.0 the class doesn't allow passing a custom callable due to technical reasons, but this is likely to be added in a future version.
• 通过增加 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维数组中指定，或者作为1d 数组的元组指定:

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

输出点数组可以指定为任意维度的数组元组(如上面两个代码片段所示) ，这使我们具有更大的灵活性。

In a nutshell, scipy.interpolate.griddata

• 即使对于疯狂的输入数据也能产生良好的输出
• 支持高维插值
• 不执行外推，一个单一的值可以设置为输出外凸壳的输入点(见 fill_value)
• 计算单个调用中的插值值，因此从头开始探测多组输出点
• 可以有任意形状的输出点
• 支持任意尺寸的最近邻和线性插值，1 d 和2 d 的立方体。最近的邻居和线性插值分别在引擎盖下使用 NearestNDInterpolator和 LinearNDInterpolator。一维三次插值采用样条函数，二维三次插值采用 CloughTocher2DInterpolator构造连续可微的分段三次插值器。
• 可能违反输入数据的对称性

scipy.interpolate.interp2d /< a href = “ https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.interpolate.bisplrep.html # scipy.interpolate.bisplrep”rel = “ nofollow noReferrer”> scipy.interpolate.bisplrep

我讨论 interp2d和它的亲戚的唯一原因是它有一个欺骗性的名字，人们可能会尝试使用它。剧透一下，别用它。interp2d在 SciPy 版本1.10中已被弃用，并将在 SciPy 版本1.12中删除.见 这个邮件列表讨论的细节。它也比以前的主题更特别，因为它是专门用于二维插值，但我怀疑这是迄今为止最常见的情况下，多元插值。

就语法而言，interp2d类似于 RBFInterpolator，因为它首先需要构造一个插值实例，可以调用该实例来提供实际的插值值。然而，这里有一个问题: 输出点必须位于一个矩形网格上，因此输入到内插器调用中的输入必须是跨越输出网格的1d 向量，就像从 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的包装器，而 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版本中)有一个很好的 在 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有限)
• performs extrapolation
• 创建一个插值器作为第一步，因此在不同的输出点评估它是较少的额外工作
• 只能在矩形网格上产生输出，对于散乱的输出，你必须在一个循环中调用插值器
• supports linear, cubic and quintic interpolation
• 可能违反输入数据的对称性

^{我相当肯定 RBFInterpolator的 cubic和 linear类基函数并不完全对应于同名的其他插值器。
这些 NaNs 也是为什么表面图看起来如此奇怪的原因: matplotlib 在绘制具有适当深度信息的复杂3D 物体方面存在历史上的困难。数据中的 NaN 值混淆了渲染器，因此应该在后面的表面部分被绘制为在前面。这是可视化的问题，而不是插值的问题。}