氣象水文科研貓公眾號交流郵箱：leolovehydrometeor@hotmail.com歡迎投稿&批評指正如有侵權且本公眾號未能正確引用原文，請聯繫刪除，謝謝理解、謝謝配合。猶記得研三那年，中山大學某博士建議我棄用Matlab，轉戰R語言和Python。故本文將對Matlab的地圖白化絕口不提，畢竟當年我真的沒搞明白如何利用m_map工具箱。1 我的m_map學習成果：
註：我未能學到精髓，不敢班門弄斧耍大刀。
2 MeteoInfo的idw插值：
註：MeteoInfo的Kriging模塊尚未添加。
fn = 'F:/RMeteoInfo/data/AQI3.txt'ncol = numasciicol(fn)nrow = numasciirow(fn)a = asciiread(fn,shape=(nrow,ncol))lon = a[:,0]lat = a[:,1]pm= a[:,2]x = arange(113, 120, 0.05)y = arange(36, 43, 0.05)gtemp,gx,gy = griddata((lon, lat), pm, xi=(x, y), method='idw', radius=1.6)axesm()bou1_layer = shaperead('F:/Rpeng/31/data3/bound/bound.shp')mlayer = shaperead('F:/Rpeng/31/data3/bound/bound.shp')geoshow(bou1_layer, edgecolor='lightgray')geoshow(mlayer, visible=False)levs = [0, 25, 35, 50, 75, 95, 115, 130, 150,200]cols = [(255,255,255),(0,255,0),(127,255,0),(255,255,0),(255,215,0),(255,128,0),(255,97,0),(255,0,0),(176,23,31),(135,38,87),(255,0,255)]layer = contourfm(x, y, gtemp,levs,colors=cols)slayer = scatterm(lon, lat,pm,levs,colors=cols, size=8)masklayer(mlayer, [layer])xlim(113, 120)ylim(36, 43) colorbar(layer, orientation='horizontal',ticklen=0,extendrect=False, shrink=1, aspect=30)yaxis(tickvisible=False,location='right')xaxis(tickvisible=False,location='top') savefig('F:/RMeteoInfo/plot34.png', dpi=600)3 我的ArcGis成果：
註：繪圖很漂亮但操作較為繁瑣。
4 R語言的Kriging、IDW往期推文超連結：
《R語言基礎地圖構建（12）》
5 Python—Kriging、IDW空間插值：
註：本文僅以basemap為例，如同布朗尼和布萊斯的關係，故cartopy請讀者自行修改。
5.1 Python—Kriging空間插值：
import numpy as npimport pandas as pdfrom pykrige.ok import OrdinaryKrigingfrom pykrige.kriging_tools import write_asc_gridimport pykrige.kriging_tools as ktimport matplotlib.pyplot as pltfrom mpl_toolkits.basemap import Basemapfrom matplotlib.colors import LinearSegmentedColormapfrom matplotlib.patches import Path, PathPatchimport maskout2shp_path=r'F:/Rpython/lp3/guangxi/map/guangxi.shp'filename=r'F:/Rpython/lp16/data/ex3.xlsx'df=pd.read_excel(filename)lons=df['lon']lats=df['lat']data=df['temp2']extent=[104,113,21,27]grid_space = 0.01grid_lon = np.arange(104,113,0.01)grid_lat = np.arange(21,27,0.01)OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian', verbose=True, enable_plotting=False,nlags=5)z1, ss1 = OK.execute('grid', grid_lon, grid_lat)xintrp, yintrp = np.meshgrid(grid_lon, grid_lat)fig, ax = plt.subplots(figsize=(10,10))m=Basemap(llcrnrlon=104-0.1,llcrnrlat=21-0.1,urcrnrlon=113+0.1,urcrnrlat=27+0.1,projection = 'cyl')m.readshapefile('F:/Rpython/lp3/guangxi/map/guangxi','whatevername',color='k',linewidth=1)x,y=m(xintrp, yintrp)cs1=ax.contourf(x, y, z1,levels=np.arange(1300,2000,100),extend='both',cmap='gist_rainbow') cs2= ax.contour(x, y, z1,colors='red',linewidths=0.6)b=plt.colorbar(cs1,shrink=0.55,orientation='vertical',extend='both',pad=0.035,aspect=20) parallels = np.arange(21,27.1,1)m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.2,color='red',dashes=[1,4],size=12)meridians = np.arange(104,113.1,1)m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.2,color='red',dashes=[1,4],size=12)plt.yticks(parallels,len(parallels)*[''])ax.tick_params(axis="y",which="both",direction="in")plt.xticks(meridians,len(meridians)*[''])ax.tick_params(axis="x",which="both",direction="in")clip1=maskout2.shp2clip(cs1,ax,r'F:/Rpython/lp3/guangxi/map/guangxi.shp') clip2=maskout2.shp2clip(cs2,ax,r'F:/Rpython/lp3/guangxi/map/guangxi.shp') plt.savefig(r"F:/Rpython/lp27/plot65.1.2.1.png",dpi=600)plt.show()5.2 Python—IDW空間插值：
讀者可嘗試定義idw(本文先不給出idw)。
import numpy as npimport matplotlib.pyplot as pltfrom scipy.interpolate import Rbfdef main(): n = 10 nx, ny = 50, 50 x, y, z = map(np.random.random, [n, n, n]) xi = np.linspace(x.min(), x.max(), nx) yi = np.linspace(y.min(), y.max(), ny) xi, yi = np.meshgrid(xi, yi) xi, yi = xi.flatten(), yi.flatten() grid1 = simple_idw(x,y,z,xi,yi) grid1 = grid1.reshape((ny, nx)) grid2 = scipy_idw(x,y,z,xi,yi) grid2 = grid2.reshape((ny, nx)) grid3 = linear_rbf(x,y,z,xi,yi) grid3 = grid3.reshape((ny, nx)) plot(x,y,z,grid1) plt.title('Homemade IDW') plot(x,y,z,grid2) plt.title("Scipy's Rbf with function=linear") plot(x,y,z,grid3) plt.title('Homemade linear Rbf') plt.show()def simple_idw(x, y, z, xi, yi): dist = distance_matrix(x,y, xi,yi) weights = 1.0 / dist weights /= weights.sum(axis=0) zi = np.dot(weights.T, z) return zidef linear_rbf(x, y, z, xi, yi): dist = distance_matrix(x,y, xi,yi) internal_dist = distance_matrix(x,y, x,y) weights = np.linalg.solve(internal_dist, z) zi = np.dot(dist.T, weights) return zidef scipy_idw(x, y, z, xi, yi): interp = Rbf(x, y, z, function='linear') return interp(xi, yi)def distance_matrix(x0, y0, x1, y1): obs = np.vstack((x0, y0)).T interp = np.vstack((x1, y1)).T d0 = np.subtract.outer(obs[:,0], interp[:,0]) d1 = np.subtract.outer(obs[:,1], interp[:,1]) return np.hypot(d0, d1)def plot(x,y,z,grid): plt.figure() plt.imshow(grid, extent=(x.min(), x.max(), y.max(), y.min())) plt.hold(True) plt.scatter(x,y,c=z) plt.colorbar()if __name__ == '__main__': main()func=Rbf(lon,lat,tem,function='linear')func=Rbf(lon,lat,tem,function='cubic')綜上所述，不同的插值軟體（本文僅選用五種）、不同的插值方法（本文僅選用Kriging、IDW、Rbf三種）所得到的插值結果的差異還是很大的！具體是什麼原因造成的呢？做研究時我們應該選用哪種方法呢？大家一起在留言區討論下吧，因筆者時間有限，本期暫不做深入探討，下期我們再一起探討、交流、學習吧。關注氣象水文科研貓公眾號一起交流學習吧.
只要你點右下方 【贊】 和 【在看】
我們就是好朋友，Thanks♪(･ω･)ﾉ