我如何使用scipy进行二维插值？

本问答旨在作为关于使用scipy进行二维（和多维）插值的规范（-ish）。 关于各种多维插值方法的基本语法常常有问题，我希望也可以直接设置这些。

我有一组分散的二维数据点，我想绘制它们作为一个很好的表面，最好是在matplotlib.pyplot使用类似contourf或plot_surface东西。 我怎么能插入我的二维或多维数据使用scipy网格？

我find了scipy.interpolate子包，但是在使用interp2d或者bisplrep或者griddata或者bisplrep时候，我总是收到错误。 这些方法的正确语法是什么？

Answer 1

免责声明：我主要是写这篇文章的语法考虑和一般行为。 我对所描述的方法的内存和CPU方面并不熟悉，而且我将这个答案瞄准那些拥有相当小的数据集的人，这样插值的质量就成为需要考虑的主要方面。 我知道，当处理非常大的数据集时，性能较好的方法（即griddata和Rbf ）可能不可行。

我将比较三种多维插值方法（ interp2d / splines， griddata和Rbf ）。 我将对它们进行两种内插任务和两种基本function（要插入的点）。 具体例子将演示二维插值，但可行的方法适用于任意维度。 每种方法提供各种插值; 在所有情况下，我将使用三次插值（或closures¹ ）。 需要注意的是，无论何时使用插值，都会引入与原始数据相比较的偏差，并且所使用的具体方法会影响您最终得到的工件。 一定要意识到这一点，并负责任地插入。

这两个插值任务将是

1. 上采样（input数据位于矩形网格上，输出数据位于更密集的网格上）
2. 将散乱的数据插入一个规则的网格

这两个函数（ [x,y] in [-1,1]x[-1,1]的域[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]

以下是他们的样子：

我将首先演示三种方法在这四种testing下的行为，然后我将详细介绍这三种方法的语法。 如果你知道你应该从某种方法中得到什么，你可能不想浪费时间学习它的语法（看着你， interp2d ）。

testing数据

为了明确起见，这里是我生成input数据的代码。 虽然在这个特定的情况下，我明显知道数据的基础function，我只会用它来为插值方法生成input。 为了方便起见，我使用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) 

平滑function和上采样

让我们从最简单的任务开始。 下面是从形状网格[6,7]到[20,21]之一的上采样对平滑testing函数的作用：

尽pipe这是一个简单的任务，但是输出之间已经有了细微的差别。 乍一看，这三个输出都是合理的。 有两个特点需要注意，根据我们对底层函数的先验知识： griddata的中间情况最griddata扭曲数据。 注意plot的y==-1边界（最接近x标签）：函数应该严格为零（因为y==-1是平滑函数的交点线），但griddata不是这种情况。 还要注意曲线的x==-1边界（在左边后面）：底层函数在[-1, -0.5]处有一个局部极大值（意味着在边界附近有零梯度），但griddata输出显示清楚该区域的非零梯度。 效果是微妙的，但是这是一种偏见。 （ Rbf的保真度与径向函数的默认select相比更好，被称为multiquadratic函数。）

邪恶的function和上采样

有一点困难的任务是对我们的恶function进行升频采样：

这三种方法已经开始显现出明显的差异。 观察曲面图， interp2d的输出中出现了明显的虚假极值（注意绘制曲面右侧的两个interp2d ）。 虽然griddata和Rbf看起来似乎乍一看产生了相似的结果，但后者似乎在[0.4, -0.4]附近产生了一个更深的最小值，而这个最小值不存在于底层函数中。

然而，有一个关键的方面，其中Rbf是远远优越的：它尊重底层函数的对称性（当然也可以通过样本网格的对称性来实现）。 griddata的输出打破了采样点的对称性，在平滑的情况下，这些对称点已经很弱。

平滑的function和分散的数据

大多数人想要对分散的数据进行插值。 为此，我期望这些testing更重要。 如上所示，样本点在感兴趣的区域中被伪均匀地select。 在实际情况下，每次测量可能会有额外的噪音，您应该考虑插入原始数据是否合理。

输出为平滑function：

现在已经有一个恐怖节目正在进行。 我将interp2d的输出interp2d到[-1, 1]之间[-1, 1]专门用于绘图，以便保留至lessless量的信息。 很明显，虽然存在一些潜在的形状，但是存在噪声很大的区域，其中该方法完全失效。 griddata的第二种情况很好地重现了形状，但是注意到等高线图边界处的白色区域。 这是因为griddata只能在input数据点的凸包内运行（换句话说，它不会执行任何外插 ）。 我保留了位于凸包外的输出点的默认NaN值。 考虑到这些特点， Rbf似乎performance最好。

邪恶的function和分散的数据

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

interp2d放弃并不令人惊讶。 实际上，在调用interp2d期间，您应该期望一些友好的RuntimeWarning抱怨样条不可能被构build。 至于其他两种方法， Rbf似乎产生最好的输出，即使是在外推结果的边界附近。

---

所以，让我对这三种方法说几句话，按照偏好的顺序排列（所以最差的是最不可能被人读的）。

scipy.interpolate.Rbf

Rbf类代表“径向基函数”。 说实话，我从来没有考虑过这个方法，直到我开始研究这个职位，但我很确定我将在未来使用这些。

就像基于样条函数的方法一样（参见后面的内容），用法分两步进行：首先根据input数据创build一个可调用的Rbf类实例，然后为给定的输出网格调用该对象以获得插值结果。 平滑上采样testing的例子：

 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 

请注意，在这种情况下，input点和输出点都是二维数组，输出z_dense_smooth_rbf与x_dense和y_dense具有相同的形状。 另请注意， Rbf支持插值的任意尺寸。

所以， scipy.interpolate.Rbf

• 即使对于疯狂的input数据也能产生良好的输出
• 支持更高维度的插值
• 外推投影点的凸包外（当然，外推总是一场赌博，你一般不应该依赖它）
• 创build一个内插器作为第一步，所以评估它在各个输出点是更less的额外努力
• 可以具有任意形状的输出点（而不是被限制为矩形网格，稍后参见）
• 容易保持input数据的对称性
• 支持关键词function多种径向函数： multiquadric ， inverse ， gaussian ， linear ， cubic ， thin_plate ， thin_plate – thin_plate和自定义任意

scipy.interpolate.griddata

我以前最喜欢的griddata是任意维度插值的一般主力。 除了为节点的凸包外部的点设置单个预设值之外，它不执行外推，但是由于外推是一个非常易变和危险的事情，所以这不一定是con。 用法示例：

 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 

注意稍微有点语法。 input点必须在D维中的形状为[N, D]的数组中指定。 为此，我们首先必须拼合我们的2D坐标数组（使用ravel ），然后连接数组并转置结果。 有多种方法可以做到这一点，但他们都似乎体积庞大。 input的z数据也必须变平。 当涉及到输出点时，我们有更多的自由：由于某些原因，这些也可以被指定为multidimensional array的元组。 请注意， griddata的help是误导性的，因为它表明input点也是如此（至less在版本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

• 即使对于疯狂的input数据也能产生良好的输出
• 支持更高维度的插值
• 不执行外推，可以为input点的凸包外的输出设置单个值（请参阅fill_value ）
• 在一个调用中计算内插值，因此从头开始探测多组输出点
• 可以有任意形状的输出点
• 支持最近邻和任意维的线性插值，立方体在1d和2d。 最近邻和线性插值分别使用NearestNDInterpolator和LinearNDInterpolator 。 1d三次插值采用样条，2d三次插值采用CloughTocher2DInterpolator构造一个连续可微分CloughTocher2DInterpolator三次插值器。
• 可能会违反input数据的对称性

scipy.interpolate.interp2d / scipy.interpolate.bisplrep

我讨论interp2d及其亲属的唯一原因是它有一个欺骗性的名字，人们可能会尝试使用它。 扰stream警报：不要使用它（从scipy版本0.17.0）。 它比以前的主题更专门，它专门用于二维插值，但是我怀疑这是迄今为止多variables插值最常见的情况。

就语法而言， interp2d类似于Rbf ，因为它首先需要构造一个插值实例，可以调用它来提供实际的插值。 但是有一个问题：输出点必须位于一个矩形网格中，所以input到插补器的input必须是跨越输出网格的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最常见的错误interp2d是把你的完整的2D网格插入插值调用，这导致爆炸性的内存消耗，并希望仓促MemoryError 。

现在， interp2d最大的问题是它往往不起作用。 为了理解这一点，我们必须仔细研究。 事实certificate， 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）interp2d的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.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

• 甚至可能导致文物暴躁的数据
• 是专门针对双variables问题的（尽pipe在网格上定义的input点的interpn有限）
• 执行外推
• 创build一个内插器作为第一步，所以评估它在各个输出点是更less的额外努力
• 只能通过矩形网格产生输出，对于分散输出，您必须在循环中调用插补器
• 支持线性，三次和五次插值
• 可能会违反input数据的对称性

---

¹我相当肯定Rbf的cubic和linear基函数并不完全对应于其他相同名称的插值器。
²这些NaN也是表面情节显得如此古怪的原因：matplotlib在绘制具有适当深度信息的复杂3d对象方面历来有困难。 数据中的NaN值混淆了渲染器，所以应该在后面的部分曲面被绘制在前面。 这是可视化问题，而不是插值。