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)中的用法。

1
2
3
4
5
6
7
8
9
10
11
12
13
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插值方法还有nearestpreviousnext,它们分别返回沿着x轴最近的、前一点以及后一点的值。
下面的例子展现了它们的用法,使用与前一个例子相同的数据:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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来寻找样条曲线的表达式。只需要给定前面两个参数,来给定曲线的xxyy分量。通常的输出是一个三元组(t,c,k)(t,c,k),包含钮点(knot-points)、系数cc和样条曲线的阶数kk(order)。默认的样条阶数是cubic(三阶),但是可以通过输入参数kk来修改。

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

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

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

示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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()

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

三阶样条的一阶导数:

1
2
3
4
5
6
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()

结果

三阶样条的所有导数:

1
2
3
4
5
6
7
8
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()

结果

样条的积分:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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()

结果

样条的根:

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

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

1
2
3
x = np.linspace(-np.pi/4, 2*np.pi + np.pi/4, 21)
y = np.sin(x)
tck = interpolate.splrep(x, y, s=0)
1
2
>>> interpolate.sproot(tck)
array([-2.22044605e-16, 3.14159265e+00, 6.28318531e+00])

参数样条曲线:

1
2
3
4
5
6
7
8
9
10
11
12
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()