插值方法很多,本文只是總結下 scipy 庫中插值用法
一維插值 - 拉格朗日插值
import numpy as np from scipy import interpolate import matplotlib.pylab as plt import pylab as mpl mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定默認字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負號'-'顯示為方塊的問題 ##################################### 一維數據插值 拉格朗日插值 x = np.linspace(0, 2*np.pi+np.pi/4, 10) y = np.sin(x) f_linear = interpolate.lagrange(x, y) # 拉格朗日插值 x_new = np.linspace(0, 2*np.pi+np.pi/4, 100) plt.plot(x, y, "o", label=u"原始數據") plt.plot(x_new, f_linear(x_new), label=u"拉格朗日插值") plt.xlabel(u'安培/A') plt.ylabel(u'伏特/V') plt.legend() plt.show()
輸出
一維插值 - B樣條插值
並且與 線性插值 做了對比
x = np.linspace(0, 2*np.pi+np.pi/4, 10) y = np.sin(x) f_linear = interpolate.interp1d(x, y) # 線性插值 tck = interpolate.splrep(x, y) # B-spline插值 x_new = np.linspace(0, 2*np.pi+np.pi/4, 100) y_bspline = interpolate.splev(x_new, tck) plt.plot(x, y, "o", label=u"原始數據") plt.plot(x_new, f_linear(x_new), label=u"線性插值") plt.plot(x_new, y_bspline, label=u"B-spline插值") plt.xlabel(u'安培/A') plt.ylabel(u'伏特/V') plt.legend() plt.show()
輸出
一維插值 - 其他插值匯總
由於用法一樣,這里統一記錄
x = np.linspace(0, 10, 11) y = np.sin(x) plt.figure(figsize=(12,9)) plt.plot(x, y, 'ro') #建立插值數據點 xnew = np.linspace(0, 10, 101) for kind in ["nearest","zero","slinear","quadratic","cubic","linear"]:#插值方式 #"nearest","zero"為階梯插值 #slinear 線性插值 #"quadratic","cubic" 為2階、3階B樣條曲線插值 f = interpolate.interp1d(x, y, kind = kind) ynew = f(xnew)#計算插值結果 plt.plot(xnew, ynew, label = str(kind)) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.legend(loc = 'lower right') plt.show()
輸出
二維插值 - 樣條插值
這里只是拿 樣條 做個 demo,其他類型的插值 類似
def func(x, y): return (x+y)*np.exp(-5.0*(x**2 + y**2)) # X-Y軸分為15*15的網格 y, x = np.mgrid[-1:1:15j, -1:1:15j] fvals = func(x,y) # 計算每個網格點上的函數值 15*15的值 print(len(fvals[0])) # 三次樣條二維插值 newfunc = interpolate.interp2d(x, y, fvals, kind='cubic') # x y z # 計算100*100的網格上的插值 xnew = np.linspace(-1,1,100)#x ynew = np.linspace(-1,1,100)#y fnew = newfunc(xnew, ynew)#僅僅是y值 100*100的值 ##### 繪圖 # 為了更明顯地比較插值前后的區別,使用關鍵字參數interpolation='nearest' # 關閉 imshow()內置的插值運算 plt.subplot(121) im1 = plt.imshow(fvals, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower") # pl.cm.jet # extent=[-1,1,-1,1]為x,y范圍 plt.colorbar(im1) plt.subplot(122) im2 = plt.imshow(fnew, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower") plt.colorbar(im2) plt.show()
輸出
二維插值 - 三維展示
from mpl_toolkits.mplot3d import Axes3D import matplotlib.cm as cm def func(x, y): return (x+y)*np.exp(-5.0*(x**2 + y**2)) # X-Y軸分為20*20的網格 x = np.linspace(-1, 1, 20) y = np.linspace(-1,1,20) x, y = np.meshgrid(x, y) # 20*20的網格數據 fvals = func(x,y) # 計算每個網格點上的函數值 15*15的值 fig = plt.figure(figsize=(9, 6)) ax = plt.subplot(1, 2, 1,projection = '3d') surf = ax.plot_surface(x, y, fvals, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('f(x, y)') plt.colorbar(surf, shrink=0.5, aspect=5)# 標注 ##### 二維插值 newfunc = interpolate.interp2d(x, y, fvals, kind='cubic') # newfunc為一個函數 # 計算100*100的網格上的插值 xnew = np.linspace(-1,1,100)#x ynew = np.linspace(-1,1,100)#y fnew = newfunc(xnew, ynew)#僅僅是y值 100*100的值 np.shape(fnew) is 100*100 xnew, ynew = np.meshgrid(xnew, ynew) ax2 = plt.subplot(1, 2, 2,projection = '3d') surf2 = ax2.plot_surface(xnew, ynew, fnew, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True) ax2.set_xlabel('xnew') ax2.set_ylabel('ynew') ax2.set_zlabel('fnew(x, y)') plt.colorbar(surf2, shrink=0.5, aspect=5)#標注 plt.show()
輸出
參考資料: