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

Question

该Q&A旨在作为关于使用scipy的二维(和多维)插值的规范(-ish).关于各种多维插值方法的基本语法经常存在问题,我希望也能这样做.

我有一组散二维数据点,我想他们绘制作为一个很好的表面,最好使用类似contourf或plot_surface在matplotlib.pyplot.如何使用scipy将我的二维或多维数据插入到网格中？

我已经找到了scipy.interpolate子包,但我一直使用时收到错误interp2d或bisplrep或griddata或rbf.这些方法的正确语法是什么？

Answer 1

免责声明:我主要是在写这篇文章时考虑到语法因素和一般行为.我不熟悉所描述的方法的内存和CPU方面,我将这个答案的目标对准那些拥有相当小的数据集的人,这样插值的质量可能是需要考虑的主要方面.我知道,在处理非常大的数据集时,性能更好的方法(即griddata和Rbf)可能不可行.

我将比较三种多维插值方法(interp2d/ splines,griddata和Rbf).我将使它们受到两种插值任务和两种底层函数(要插入的点).具体的例子将展示二维插值,但可行的方法适用于任意维度.每种方法都提供各种插值; 在所有情况下,我将使用三次插值(或接近¹).重要的是要注意,无论何时使用插值,都会引入与原始数据相比的偏差,并且所使用的特定方法会影响您最终会得到的工件.始终注意这一点,并负责任地插入.

两个插值任务将是

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
x_scattered,y_scattered = np.random.rand(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输出显示该区域内的非零梯度.效果很微妙,但它仍然是一种偏见.(Rbf配音默认选择径向功能,保真度更高multiquadratic.)

邪恶的功能和上采样

一个更难的任务是对我们的邪恶功能进行上采样:

三种方法之间开始显示出明显的差异.观察表面图,输出中出现明显的虚假极值interp2d(注意绘制表面右侧的两个凸起).虽然griddata并且Rbf乍一看似乎产生类似的结果,但后者似乎产生了一个更深的最小值[0.4, -0.4],而不是基本功能.

然而,有一个关键方面Rbf是远远优越的:它尊重底层函数的对称性(当然也可以通过样本网格的对称性实现).输出griddata打破了采样点的对称性,这在光滑的情况下已经很弱.

平滑的功能和分散的数据

通常,人们希望对分散的数据执行插值.出于这个原因,我希望这些测试更重要.如上所示,在感兴趣的域中伪均匀地选择样本点.在实际场景中,每次测量可能会产生额外的噪声,您应该考虑是否有必要插入原始数据.

平滑功能的输出:

现在已经有一些恐怖表演正在进行中.我将输出剪切interp2d到两者之间[-1, 1]专门用于绘图,以便至少保留少量信息.很明显,虽然存在一些潜在的形状,但是存在巨大的嘈杂区域,其中该方法完全破坏.第二种情况griddata很好地再现了形状,但注意了等高线边界处的白色区域.这是因为griddata只在输入数据点的凸包内部工作(换句话说,它不执行任何外推).我保留了位于凸包外面的输出点的默认NaN值.²考虑到这些功能,Rbf似乎表现最佳.

邪恶的功能和分散的数据

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

interp2d放弃这并不奇怪.事实上,在打电话给interp2d你的时候,应该期待一些友好RuntimeWarning的人抱怨样条的不可能性.至于其他两种方法,Rbf似乎产生最佳输出,甚至在结果被外推的域的边界附近.

因此,让我对三种方法说几句话,按优先顺序递减(这样最差的人最不可能被任何人阅读).

scipy.interpolate.Rbf

该Rbf级代表"径向基函数".老实说,在我开始研究这篇文章之前,我从未考虑过这种方法,但我很确定将来会使用这些方法.

就像基于样条的方法(见下文)一样,用法分为两步:第一步Rbf根据输入数据创建一个可调用的类实例,然后为给定的输出网格调用该对象以获得插值结果.平滑上采样测试的示例:

import scipy.interpolate as interp
zfun_smooth_rbf = interp.Rbf(x_sparse, y_sparse, z_sparse_smooth, function='cubic', smooth=0)  # default smooth=0 for interpolation
z_dense_smooth_rbf = zfun_smooth_rbf(x_dense, y_dense)  # not really a function, but a callable class instance

请注意,在这种情况下,输入和输出点都是2d数组,并且输出z_dense_smooth_rbf具有相同的形状,x_dense并且y_dense不需要任何努力.另请注意,Rbf支持插值的任意尺寸.

所以, scipy.interpolate.Rbf

• 即使对于疯狂的输入数据,也能产生良好的输出
• 支持更高维度的插值
• 在输入点的凸包外面推断(当然外推总是赌博,你通常不应该完全依赖它)
• 创建插补器作为第一步,因此在各种输出点进行评估是额外的努力
• 可以有任意形状的输出点(而不是被约束到矩形网格,见后面)
• 倾向于保持输入数据的对称性
• 支持多种径向函数关键字function:multiquadric,inverse,gaussian,linear,cubic,quintic,thin_plate和用户自定义的任意

scipy.interpolate.griddata

我以前的最爱,griddata是任意尺寸插值的一般主力.除了为节点的凸包外部的点设置单个预设值之外,它不执行外推,但由于外推是非常易变和危险的事情,因此这不一定是骗局.用法示例:

z_dense_smooth_griddata = interp.griddata(np.array([x_sparse.ravel(),y_sparse.ravel()]).T,
                                          z_sparse_smooth.ravel(),
                                          (x_dense,y_dense), method='cubic')   # default method is linear

请注意稍微克隆的语法.必须[N, D]在D尺寸的形状数组中指定输入点.为此,我们首先必须展平我们的2d坐标数组(使用ravel),然后连接数组并转置结果.有多种方法可以做到这一点,但所有这些方法看起来都很笨重.输入z数据也必须展平.当涉及输出点时,我们有更多的自由:由于某种原因,这些也可以被指定为多维数组的元组.注意help的griddata是误导,因为它表明,同样是为真正的输入点(至少在0.17.0或更新版本):

griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)
    Interpolate unstructured D-dimensional data.

    Parameters
    ----------
    points : ndarray of floats, shape (n, D)
        Data point coordinates. Can either be an array of
        shape (n, D), or a tuple of `ndim` arrays.
    values : ndarray of float or complex, shape (n,)
        Data values.
    xi : ndarray of float, shape (M, D)
        Points at which to interpolate data.

简而言之, scipy.interpolate.griddata

• 即使对于疯狂的输入数据,也能产生良好的输出
• 支持更高维度的插值
• 不执行外推,可以为输入点的凸包外部的输出设置单个值(请参阅参考资料fill_value)
• 计算单个调用中的插值,因此探测多组输出点从头开始
• 可以有任意形状的输出点
• 支持任意维度的最近邻和线性插值,1d和2d中的立方.最近邻和线性插值分别使用NearestNDInterpolator和LinearNDInterpolator引擎盖下.1d三次插值使用样条,2d三次插值用于CloughTocher2DInterpolator构造连续可微分段 - 三次插值器.
• 可能会违反输入数据的对称性

scipy.interpolate.interp2d/scipy.interpolate.bisplrep

我正在讨论的唯一原因interp2d及其亲属是它有一个欺骗性的名字,人们可能会尝试使用它.扰流警报:不要使用它(从scipy版本0.17.0开始).它已经比以前的主题更特殊,因为它专门用于二维插值,但我怀疑这是多变量插值的最常见情况.

就语法而言,interp2d类似于Rbf它首先需要构造一个插值实例,可以调用该插值实例来提供实际的插值.然而,有一个问题:输出点必须位于矩形网格上,因此进入内插器调用的输入必须是跨越输出网格的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包装器,它们又是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的版本0.17.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网格上定义的输入点有限)
• 执行外推
• 创建插补器作为第一步,因此在各种输出点进行评估是额外的努力
• 只能在矩形网格上产生输出,对于分散输出,您必须在循环中调用插值器
• 支持线性,立方和五次插值
• 可能会违反输入数据的对称性

¹我很确定基本函数cubic和linear基本函数Rbf与其他同名内插器并不完全对应.
²这些NaN也是表面图看起来如此奇怪的原因:matplotlib历史上很难用复杂的深度信息绘制复杂的3d对象.数据中的NaN值会混淆渲染器,因此应该在后面的曲面部分被绘制在前面.这是可视化问题,而不是插值问题.