【问题标题】:How can I perform two-dimensional interpolation using scipy?如何使用 scipy 执行二维插值?
【发布时间】:2016-10-18 17:41:20
【问题描述】:

本问答旨在作为关于使用 scipy 进行二维(和多维)插值的规范(-ish)。各种多维插值方法的基本语法经常有问题,我也希望把这些弄清楚。

我有一组分散的二维数据点,我想将它们绘制成一个漂亮的表面,最好在matplotlib.pyplot 中使用contourfplot_surface 之类的东西。如何使用 scipy 将我的二维或多维数据插入到网格中?

我找到了scipy.interpolate 子包,但在使用interp2dbisplrepgriddataRBFInterpolator(或更旧的Rbf)时,我不断收到错误消息。这些方法的正确语法是什么?

【问题讨论】:

  • 我还建议添加其他插值任务:输入数据是分散的或矩形网格,输出数据是分散的(即任意坐标)。
  • @BrianD 除非我误解了“输入数据是分散的或矩形网格”,这正是我在回答中考虑的两个任务。至于“输出数据是分散的”,我不确定这是否有趣:一旦你有一个插值函数,你就可以在任何地方查询它。结果不应该对输出点的位置在质量上敏感。你有什么想法?
  • 这很公平。我想我在考虑“任意采样”(而不是在网格上进行上采样或插值)并且只是在与 interp2d 苦苦挣扎,现在使用与 Rbf 相同的代码“正常工作”所以......没关系?
  • @BrianD 哈,是的,这几乎是我整理此问答时的经验。还有一个 Rbf 的新实现正在工作中,它非常接近完成。我会在答案发布后更新我的答案,因为它看起来不错(而且更具前瞻性)。
  • @BrianD 我不确定你的意思是什么功能(外推,即域中的缺失值,或者只是钳位插值函数值)。听起来您的意思是后者,但这可以通过在插值后调用np.clip 来完成。但要回答您的问题,不,我认为新功能不会支持这一点。

标签: python scipy interpolation


【解决方案1】:

免责声明:我在写这篇文章时主要考虑了句法考虑和一般行为。我不熟悉所描述方法的内存和 CPU 方面,我将这个答案针对那些拥有相当少量数据集的人,因此插值的质量可能是需要考虑的主要方面。我知道,在处理非常大的数据集时,性能更好的方法(即 griddataRBFInterpolator 没有 neighbors 关键字参数)可能不可行。

请注意,此答案使用 SciPy 1.7.0 中引入的新 RBFInterpolator 类。对于旧版 Rbf 类,请参阅 the previous version of this answer

我将比较三种多维插值方法(interp2d/splines、griddataRBFInterpolator)。我将让他们接受两种插值任务和两种底层函数(要插值的点)。具体示例将演示二维插值,但可行的方法适用于任意维度。每种方法都提供各种插值;在所有情况下,我都会使用三次插值(或接近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 输出清楚地显示非零梯度在这个地区。效果是微妙的,但它仍然是一种偏见。

邪恶函数和上采样

一个更难的任务是对我们的 evil 函数执行上采样:

三种方法之间开始出现明显的差异。查看表面图,interp2d 的输出中出现了明显的虚假极值(注意绘制表面右侧的两个驼峰)。虽然griddataRBFInterpolator 乍一看似乎产生了相似的结果,但在[0.4, -0.4] 附近产生了底层函数中不存在的局部最小值。

但是,RBFInterpolator 在一个关键方面要优越得多:它尊重底层函数的对称性(这当然也可以通过示例网格的对称性实现)。 griddata 的输出打破了样本点的对称性,这在平滑情况下已经很微弱了。

平滑函数和分散数据

大多数情况下,人们希望对分散的数据执行插值。出于这个原因,我希望这些测试更加重要。如上所示,样本点是在感兴趣的域中伪均匀地选择的。在实际场景中,每次测量可能会产生额外的噪音,您应该考虑从一开始就插入原始数据是否有意义。

平滑函数的输出:

现在已经有一些恐怖节目正在上演。我将输出从interp2d 剪裁到[-1, 1] 之间,专门用于绘图,以便至少保留最少量的信息。很明显,虽然存在一些基本形状,但存在巨大的噪声区域,该方法完全失效。 griddata 的第二种情况很好地再现了形状,但请注意等高线图边界处的白色区域。这是因为griddata 仅在输入数据点的凸包内起作用(换句话说,它不执行任何外推)。我为凸包外的输出点保留了默认的 NaN 值。2考虑到这些特性,RBFInterpolator 似乎表现最好。

邪恶的功能和分散的数据

我们一直在等待的那一刻:

interp2d 放弃并不奇怪。事实上,在调用interp2d 期间,您应该期待一些友好的RuntimeWarnings 抱怨无法构造样条线。至于其他两种方法,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 满意。由于我们必须将二维点作为形状(N, 2) 的数组传递,我们必须展平输入网格并将两个展平的数组堆叠起来。构造的插值器也需要这种格式的查询点,结果将是一个形状为(N,) 的一维数组,我们必须对其进行整形以匹配我们的二维网格以进行绘图。由于RBFInterpolator 不对输入点的维数做任何假设,因此它支持任意维数进行插值。

所以,scipy.interpolate.RBFInterpolator

  • 即使对于疯狂的输入数据也能产生良好的输出
  • 支持更高维度的插值
  • 在输入点的凸包之外进行外推(当然外推始终是一场赌博,您通常根本不应该依赖它)
  • 第一步创建一个插值器,因此在不同的输出点上评估它的额外工作量更少
  • 可以有任意形状的输出点数组(而不是被限制为矩形网格,见下文)
  • 更有可能保持输入数据的对称性
  • 支持关键字kernel的多种径向函数:multiquadricinverse_multiquadricinverse_quadraticgaussianlinearcubicquinticthin_plate_spline(默认) .从 SciPy 1.7.0 开始,由于技术原因,该类不允许传递自定义可调用对象,但这很可能会在未来的版本中添加。
  • 可以通过增加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 相同的数组转换。输入点必须在 D 维度的形状为 [N, D] 的数组中指定,或者作为一维数组的元组:

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

输出点数组可以指定为任意维度数组的元组(如上面的 sn-ps 中),这给了我们更多的灵活性。

简而言之,scipy.interpolate.griddata

  • 即使对于疯狂的输入数据也能产生良好的输出
  • 支持更高维度的插值
  • 不进行外插,可以为输入点凸包外的输出设置单个值(见fill_value
  • 在一次调用中计算内插值,因此从头开始探测多组输出点
  • 可以有任意形状的输出点
  • 支持任意维度的最近邻和线性插值,1d 和 2d 立方。最近邻和线性插值分别在底层使用NearestNDInterpolatorLinearNDInterpolator。一维三次插值使用样条曲线,二维三次插值使用CloughTocher2DInterpolator 构造连续可微分的分段三次插值器。
  • 可能会违反输入数据的对称性

scipy.interpolate.interp2d/scipy.interpolate.bisplrep

我讨论interp2d 及其亲属的唯一原因是它有一个欺骗性的名称,人们可能会尝试使用它。剧透警告:不要使用它(从 scipy 版本 1.7.0 开始)。它已经比之前的主题更特别了,因为它专门用于二维插值,但我怀疑这是迄今为止多元插值最常见的情况。

就语法而言,interp2dRBFInterpolator 相似之处在于它首先需要构造一个插值实例,可以调用该实例来提供实际的插值。但是有一个问题:输出点必须位于矩形网格上,因此进入插值器调用的输入必须是跨越输出网格的一维向量,就像来自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 中)interp2d 有一个不错的 comment in interpolate/interpolate.py

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.pybisplrep 中有一些设置,最终

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我相当确定cubiclinear 的基函数类型与RBFInterpolator 不完全对应于其他同名插值器。
2这些 NaN 也是表面图看起来如此奇怪的原因:matplotlib 历来难以绘制具有适当深度信息的复杂 3d 对象。数据中的 NaN 值混淆了渲染器,因此应该在后面的表面部分被绘制在前面。这是可视化问题,而不是插值问题。

【讨论】:

  • Rbf 可以比 griddata 消耗更多的内存,具体取决于数据点的数量和维度的数量。 griddata 也有底层对象 LinearNDInterpolator,可以像 Rbf 一样分两步使用。
  • Griddata 的三次插值仅限于 2 (?) 维。对于更高维度,基于 chebfun 的 smolyak 稀疏网格值得考虑。
  • 让我用这个链接结束我的 cmets,在那里我研究了所有分散的插值选项:scicomp.stackexchange.com/questions/19137/…
  • griddata 线性插值是局部的,griddata 三次插值是全局的。不支持外推,因为我没有时间弄清楚如何保持连续性/可微性。 Rbf 适用于小型数据集,但要插入 n 个数据点,它需要反转 n x n 矩阵,这最终在 n>5000 后变得不可能。 Rbf 也可能对数据的分布很敏感,您可能需要手动微调其参数。可以为大型数据集做 Rbf,但这在 scipy 中没有实现。
  • 这里是大型数据集的 rbf:github.com/scipy/scipy/issues/5180
猜你喜欢
  • 2018-09-03
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2019-06-01
  • 1970-01-01
  • 1970-01-01
  • 2021-10-01
相关资源
最近更新 更多