Scipy-interpolate

Scipy中有几个通用的插值工具，用于1D、2D以及高维插值：

• interp1d：一维插值函数的类，提供几种插值方法。
• griddata函数提供了一个简单的接口用于N维插值（N=1, 2, 3, 4, …）。也可以使用面向对象的接口来使用底层例程（underlying routines）。
• 1D和2D（光滑的）cubic-spline插值函数是基于FORTRAN库FITPACK的。它们都是FITPACK库的面向过程和面向对象接口。
• 使用径向基（radial basis functions）函数插值。

1D插值（interp1d）

scipy.interpolate中的interp1d类是一个方便的方法用于创建一个基于固定数据点的函数，
能使用线性插值来计算所给定数据范围内任意位置处的值。

此类的实例是通过传递包含数据的一维向量来创建的。类的实例定义了一个__call__方法，因此可以像函数一样处理，用于在已知数据值之间进行插值以获得未知位置的值（可查看docstring获取帮助）。可以在实例化的时候指定边界处的行为。下面的实例演示了它在线性插值（linear interpolation）和三次样条插值（cubic spline interpolation）中的用法。

import numpy as np
from scipy.interpolate import interp1d
# 定义插值函数
x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f = interp1d(x, y)
f2 = interp1d(x, y, kind='cubic')
# 插值
xnew = np.linspace(0, 10, num=41, endpoint=True)
import matplotlib.pyplot as plt
plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()

其他1D插值方法还有nearest、previous和next，它们分别返回沿着x轴最近的、前一点以及后一点的值。
下面的例子展现了它们的用法，使用与前一个例子相同的数据：

import numpy as np
from scipy.interpolate import interp1d
# 定义插值函数
x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f1 = interp1d(x, y, kind='nearest')
f2 = interp1d(x, y, kind='previous')
f3 = interp1d(x, y, kind='next')
# 插值
xnew = np.linspace(0, 10, num=1001, endpoint=True)
import matplotlib.pyplot as plt
plt.plot(x, y, 'o')
plt.plot(xnew, f1(xnew), '-', xnew, f2(xnew), '--', xnew, f3(xnew), ':')
plt.legend(['data','nearest','previous','next'],loc='best')
plt.show()

多维数据插值（griddata）

样条（Spline）插值

Spline interpolation in 1-D: Procedural (interpolate.splXXX)

样条插值需要两个基本步骤：
(1) 计算曲线的样条表示；
(2) 样条曲线在所需点处求值。
为了找到样条表示，有两种不同的方法来表示曲线以及获得（平滑）样条系数：直接法或者参数法。

直接法是在一个2D平面中使用函数splrep来寻找样条曲线的表达式。只需要给定前面两个参数，来给定曲线的$x$和$y$分量。通常的输出是一个三元组$(t,c,k)$，包含钮点（knot-points）、系数$c$和样条曲线的阶数$k$（order）。默认的样条阶数是cubic（三阶），但是可以通过输入参数$k$来修改。

对于N维空间中的曲线，函数splperp允许使用参数化的方法来定义。对于这个函数，只需要1个输入参数。该输入是一个表示N维空间中的曲线的N维数组。这个参变量通过关键字参数$u$来给定，默认为0~1之间的灯具单调序列。默认输出有两个对象构成：一个三元组$(t,c,k)$，它包含了样条曲线的表达参数；参变量$u$。

关键字参数$s$被用来指定样条拟合时的平滑量。默认取值为$s=m-\sqrt{2m}$，其中$m$为插值数据点的个数。因此，如果不需要光滑，则$s=0$。

一旦确定了数据的样条表示，就可以使用函数来评估任何点的样条(splev)及其导数(splev,spalde) 以及任意两点之间的样条积分 (splint)。 此外，对于具有 8 个或更多节的三次样条 ($k=3$)，可以估计样条的根 ( sproot)。 下面的示例演示了这些功能

示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
# cubic-spline 三阶样条插值
x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8) # 插值点：间距np.pi/4，很大
y = np.sin(x)
tck = interpolate.splrep(x,y,s=0)          # 构造插值函数
xnew = np.arange(0, 2*np.pi, np.pi/50)
ynew = interpolate.splev(xnew, tck, der=0) # 插值
# 可视化
plt.figure()
plt.plot(x,y,'xb-', xnew, ynew, xnew,np.sin(xnew))
plt.legend(['Linear','Cubic Spline', 'True'])
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('Cubic-spline interpolation')
plt.show()

可以看到，三阶样条插值的结果与真值几乎重合。

三阶样条的一阶导数：

yder = interpolate.splev(xnew, tck, der=1)
plt.figure()
plt.plot(xnew, yder, xnew, np.cos(xnew),'--')
plt.legend(['Cubic Spline','True'])
plt.title('Derivative estimation from spline')
plt.show()

三阶样条的所有导数：

yders = interpolate.spalde(xnew, tck)
plt.figure()
for i in range(len(yders[0])):
    plt.plot(xnew, [d[i] for d in yders], '--', label=f'{i} derivative')
plt.legend()
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('All derivatives of a B-spline')
plt.show()

样条的积分：

def integ(x, tck, constant=-1):
    x = np.atleast_1d(x)
    out = np.zeros(x.shape, dtype=x.dtype)
    for n in range(len(out)):
        out[n] = interpolate.splint(0, x[n], tck)
    out += constant 
    return out

yint = integ(xnew, tck)
plt.figure()
plt.plot(xnew, yint, xnew, -np.cos(xnew), '--')
plt.legend(['Cubic Spline', 'True'])
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('Integral estimation from spline')
plt.show()

样条的根：

>>> interpolate.sproot(tck)
array([3.14159265])

注意sproot在区间端点$x=0$的附近找根会失败。如果把区间放大一点，就能得到$x=0$和$x=2\pi$这两个根了：

x = np.linspace(-np.pi/4, 2*np.pi + np.pi/4, 21)
y = np.sin(x)
tck = interpolate.splrep(x, y, s=0)

>>> interpolate.sproot(tck)
array([-2.22044605e-16,  3.14159265e+00,  6.28318531e+00])

参数样条曲线：

t = np.arange(0, 1.1, .1)
x = np.sin(2*np.pi*t)
y = np.cos(2*np.pi*t)
tck, u = interpolate.splprep([x, y], s=0)
unew = np.arange(0, 1.01, 0.01)
out = interpolate.splev(unew, tck)
plt.figure()
plt.plot(x, y 'x', out[0], out[1], np.sin(2*np.pi*unew), np.cos(2*np.pi*unew), x,y,'b')
plt.legend(['Linear','Cubic Spline', 'True'])
plt.axis([-1.05, 1.05,-1.05,1.05])
plt.title('Spline of parametrically-defined curve')
plt.show()