python⼆维插值_SciPy⼆元样条插值
16. SciPy⼆元样条插值
之前的章节的插值基本上都是⼀个变量x，本章就⼆元(⼆维)的插值问题进⾏展开。
在SciPy模块库的scipy.interpolate⼦模块⾥提供了interp2d⽅法函数可以实现⼆维插值。
16.1 绘制⼆元可视化数据
Python⾥可视化输出的模块很多，但⽐较常⽤的是matplotlib模块来实现数据的可视化的输出。在matplotlib模块⾥的⼦模块pyplot⾥提供很多的类似matlib的⼀些绘制函数，本章、本⽹并没有意愿详述matplotlib模块库，有意学习pyplot的可以访问其官⽅⽹站系统学习。
1). ⽤matplotlib绘制⼆维可视化数据⾸先要引⼊相应的模块库。
import numpy as np, matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
这⾥还引⼊了mpl_toolkits.mplot3d.axes3d模块，⽬前是⽀持三维显⽰。
2). 接下来可以定义⼀个⼆元的$f(x, y)$函数。
def f(x,y):
return np.sin(x) + np.sin(y)
3). 定义⽹格数据：
t = np.linspace(-3, 3, 100)
domain = np.meshgrid(t, t)
有关meshgrid函数可以参考本站的NumPy部分的meshgird函数内容。domain变量是⼀个列表，⾥⾯有两个numpy的array。再将⽹格数据提取出来作为X和Y。在将X和Y代⼊$f(x, y)$函数得到Z值。
X, Y = domain
Z = f(*domain)
4). 绘制$f(x, y)$，这⾥使⽤了subplot2grid函数。
fig = plt.figure()
ax1 = plt.subplot2grid((2,2), (0,0), aspect='equal')
p = ax1.pcolor(X, Y, Z)
CP = ur(X, Y, Z, colors='k')
ax1.clabel(CP)
ax1.set_title('Contour plot')
语句subplot2grid((2,2), (0,0), aspect='equal')的含义是整个绘制图是2⾏2列，此次绘制的ax1图像位于0⾏0列，⽽aspect参数传给内部调⽤的add_subplot函数使⽤,'equal'参数的含义可以查阅该函数。
语句 ax1.pcolor(X, Y, Z)，这是在给此图配置颜⾊。
语句lorbar(p)的colorbar函数的作⽤则是产⽣⾊度条。
语句CP = ur(X, Y, Z, colors='k')的contour函数绘制包络曲线。
语句ax1.clabel(CP)⾥的clalel函数的作⽤则是给绘制的包络曲线加注⽂字信息。
语句ax1.set_title('Contour plot')的set_title函数则是给⼦图添加⼀个标签名。
到此的代码的绘制的图像为下图⾥位置为左上的那个图。
16.2 绘制需要插值的点
1). 产⽣插值点
nodes = 6 * np.random.rand(100, 2) - 3
#print np.random.rand(100, 2)
xi = nodes[:, 0]
yi = nodes[:, 1]
zi = f(xi, yi)
NumPy的random模块下的rand函数可以产⽣$[0, 1)$区间的均匀分布，参数(100, 2)给出函数结构是⼀个100⾏2列的数组(ndarray)。第⼀列作为插值点的$x$，第⼆列作为$y$坐标。将$x$和$y$代⼊$f(x, y)$函数得到$z$值。
2). 绘制插值点的可视化图ax2
ax2 = plt.subplot2grid((2,2), (0,1), aspect='equal')
p2 = ax2.pcolor(X, Y, Z)
ax2.scatter(xi, yi, 25, zi)
ax2.set_xlim(-3, 3)
ax2.set_ylim(-3, 3)
ax2.set_title('Node selection')
这⾥⽤到了scatter函数，他的作⽤是绘制插值点nodes。set_xlim函数是设置显⽰图的范围。⽽其他函数在上⼀节⾥介绍过了。这部分代码被执⾏时会得到下图右上⼦图。
右上⼦图ax2⼩圆圈就是nodes的可视化结果。
16.3 ⼆元3D可视化
这部分的内容是得到上图的下⾯的⼦图，即数据的3D可视化。
ax3 = plt.subplot2grid((2,2), (1,0), projection='3d', colspan=2, rowspan=2)
ax3.plot_surface(X, Y, Z, alpha=0.25)
ax3.scatter(xi, yi, zi, s=25)
cset = ur(X, Y, Z, zdir='z', offset=-4)
cset = ur(X, Y, Z, zdir='x', offset=-5)
ax3.set_xlim3d(-5, 3)
ax3.set_ylim3d(-3, 5)
ax3.set_zlim3d(-4, 2)
ax3.set_title('Surface plot')
fig.tight_layout()
plt.show()
subplot2grid⾥的projection='3d'形参的作⽤是说此图ax3是3D的。位置为(1, 0)由于后续没有(1, 1)位置的⼦图，那么此图ax3占整个第2⾏输出，绘制的X、Y、Z。
语句ax3.plot_surface(X, Y, Z, alpha=0.25)绘制的曲⾯。
⽽语句ax3.scatter(xi, yi, zi, s=25)则是3D绘制的nodes各个插值点。
cset = ur(X, Y, Z, zdir='z', offset=-4)
cset = ur(X, Y, Z, zdir='x', offset=-5)
这两条语句则是绘制的X、Y、Z在xy平⾯、yz平⾯的投影曲线。
16.4 节点插值可视化
节点nodes准备好以后可以使⽤SciPy模块库的scipy.interpolate⼦模块⾥提供了interp2d⽅法函数对nodes进⾏线性插值。
from scipy.interpolate import interp2d
interpolant = interp2d(xi, yi, zi, kind='linear')
然后便可可视化输出。
plt.figure()
plt.axes().set_aspect('equal')
plt.pcolor(X, Y, interpolant(t, t))
plt.scatter(xi, yi, 25, zi)
CP = ur(X, Y, interpolant(t, t), colors='k')
plt.clabel(CP)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.title('Piecewise linear interpolation')
plt.show()
得到结果：
需要注意的是nodes是随机产⽣的，每次运⾏程序的结果可能不⼀样。
16.5 完整的程序代码
import numpy as np, matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
def f(x,y):
return np.sin(x) + np.sin(y)
用subplot函数
t = np.linspace(-3, 3, 100)
domain = np.meshgrid(t, t)
X, Y = domain
Z = f(*domain)
fig = plt.figure()
ax1 = plt.subplot2grid((2,2), (0,0), aspect='equal')
p = ax1.pcolor(X, Y, Z)
CP = ur(X, Y, Z, colors='k')
ax1.clabel(CP)
ax1.set_title('Contour plot')
nodes = 6 * np.random.rand(100, 2) - 3
print np.random.rand(100, 2)
xi = nodes[:, 0]
yi = nodes[:, 1]
zi = f(xi, yi)
ax2 = plt.subplot2grid((2,2), (0,1), aspect='equal')
p2 = ax2.pcolor(X, Y, Z)
ax2.scatter(xi, yi, 25, zi)
ax2.set_xlim(-3, 3)
ax2.set_ylim(-3, 3)
ax2.set_title('Node selection')
ax3 = plt.subplot2grid((2,2), (1,0), projection='3d', colspan=2, rowspan=2) #ax3.plot_surface(X, Y, Z, alpha=0.25)
#ax3.scatter(xi, yi, zi, s=25)
cset = ur(X, Y, Z, zdir='z', offset=-4)
cset = ur(X, Y, Z, zdir='x', offset=-5)
ax3.set_xlim3d(-5, 3)
ax3.set_ylim3d(-3, 5)
ax3.set_zlim3d(-4, 2)
ax3.set_title('Surface plot')
fig.tight_layout()
plt.show()
from scipy.interpolate import interp2d
interpolant = interp2d(xi, yi, zi, kind='linear')
plt.figure()
plt.axes().set_aspect('equal')
plt.pcolor(X, Y, interpolant(t, t))
plt.scatter(xi, yi, 25, zi)
CP = ur(X, Y, interpolant(t, t), colors='k') plt.clabel(CP)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.title('Piecewise linear interpolation')
plt.show()

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。