前面幾篇推文我們分辨介紹了使用Python和R繪制了二維核密度空間插值方法，並使用了Python可視化庫plotnine、Basemap以及R的ggplot2完成了相關可視化教程的繪制推文，詳細內容如下：
接下來，我們將繼續介紹空間插值的其他方法，本期推文，我們將介紹IDW(反距離加權法（Inverse Distance Weighted）) 插值的Python計算方法及插值結果的可視化繪制過程。主要涉及的知識點如下：
- IDW簡介
- 自定義Python代碼計算空間IDW
- 分別使用plotnine、Basemap進行IDW插值結果可視化繪制
IDW簡介
反距離權重 (IDW) 插值假設：彼此距離較近的事物要比彼此距離較遠的事物更相似。當為任何未測量的位置預測值時，反距離權重法會采用預測位置周圍的測量值與距離預測位置較遠的測量值相比，距離預測位置最近的測量值對預測值的影響更大。反距離權重法假定每個測量點都有一種局部影響，而這種影響會隨着距離的增大而減小。由於這種方法為距離預測位置最近的點分配的權重較大，而權重卻作為距離的函數而減小，因此稱之為反距離權重法。(解釋來源於網絡)，繁瑣的公式也沒放，這里我們給出幾張示意圖即可，原理不解的小伙伴可自行百度。
（基於采樣點距離的IDW插值（左）從高程矢量點插值的IDW曲面（右））
自定義Python代碼計算空間IDW
我們免去了了繁瑣的IDW插值原理部分，這節我們直接根據原理自定義IDW函數，根據已有樣例站點位置及對應值，計算IDW結果。在這之前，我們給出所需樣例的預覽及地圖文件的范圍(構建插值網格所需)，結果如下：
樣例點：
地圖文件范圍信息：
js_box = js.geometry.total_bounds js_box #array([116.36196 , 30.757975, 121.975185, 35.122924])
小伙伴們對上述計算結果有疑惑的地方可以詳細閱讀之前的插值文章(文前鏈接)，或者等我將這系列做完會推出詳細的源碼及解釋文檔(目前在整理中)
定義IDW計算函數
這里主要涉及兩個計算函數，計算經緯度點轉實際距離(km)的haversine方法和計算IDW的函數，定義函數如下：
- haversine方法:
1 import math 2 import numpy as np 3 #更換求距離的函數 4 from math import radians, cos, sin, asin, sqrt 5 6 def haversine(lon1, lat1, lon2, lat2): 7 R = 6372.8 8 dLon = radians(lon2 - lon1) 9 dLat = radians(lat2 - lat1) 10 lat1 = radians(lat1) 11 lat2 = radians(lat2) 12 a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2 13 c = 2*asin(sqrt(a)) 14 d = R * c 15 return d
- IDW
1 def IDW(x, y, z, xi, yi): 2 lstxyzi = [] 3 for p in range(len(xi)): 4 lstdist = [] 5 for s in range(len(x)): 6 d = (haversine(x[s], y[s], xi[p], yi[p])) 7 lstdist.append(d) 8 sumsup = list((1 / np.power(lstdist, 2))) 9 suminf = np.sum(sumsup) 10 sumsup = np.sum(np.array(sumsup) * np.array(z)) 11 u = sumsup / suminf 12 xyzi = [xi[p], yi[p], u] 13 lstxyzi.append(xyzi) 14 return(lstxyzi)
計算所需插值的網格
這里直接給出代碼，階段的結果需要更具上面的函數計算對應網格點出的IDW結果，這樣就可以實現插值操作，代碼如下：
1 js_box = js.geometry.total_bounds 2 #還是插入400*400的網格點 3 grid_lon = np.linspace(js_box[0],js_box[2],400) 4 grid_lat = np.linspace(js_box[1],js_box[3],400) 5 xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
較為簡單，這里就不作多解釋。
計算IDW結果
結合上面兩個部分，我們進行了IDW插值結果，具體計算結果如下：
1 #將插值網格數據整理 2 df_grid =pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten())) 3 #這里將數組轉成列表 4 grid_lon_list = df_grid["long"].tolist() 5 grid_lat_list = df_grid["lat"].tolist() 6 7 pm_idw = IDW(know_lon,know_lat,know_z,grid_lon_list,grid_lat_list) 8 IDW_grid_df = pd.DataFrame(pm_idw,columns=["lon","lat","idw_value"]) 9 IDW_grid_df.head()
這樣就可以得到IDW插值后的DF類型數據了，結果如下(部分）：
可視化繪制
有了規整完的插值結果，那么接下來繪制可視化結果也就非常簡單了，方法和之前的幾篇推文類似，具體如下：
plotnine繪制
首先，我們還是給出樣例點及對應值的映射散點圖，繪圖過程如下：
「散點圖繪制」
1 import plotnine 2 from plotnine import * 3 plotnine.options.figure_size = (5, 4.5) 4 idw_scatter = (ggplot() + 5 geom_map(js,fill='none',color='gray',size=0.4) + 6 geom_point(pm,aes(x='經度',y='緯度',fill='PM2.5'),size=5) + 7 scale_fill_cmap(cmap_name='Spectral_r',name='PM2.5', 8 breaks=[30,40,60,80] 9 )+ 10 scale_x_continuous(breaks=[117,118,119,120,121,122])+ 11 labs(title="Map Charts in Python Exercise 02: Map IDM point", 12 )+ 13 #添加文本信息 14 annotate('text',x=116.5,y=35.3,label="processed map charts with plotnine",ha="left", 15 size=10)+ 16 annotate('text',x=120,y=30.6,label="Visualization by DataCharm",ha="left",size=9)+ 17 theme( 18 text=element_text(family="Roboto Condensed"), 19 #修改背景 20 panel_background=element_blank(), 21 axis_ticks_major_x=element_blank(), 22 axis_ticks_major_y=element_blank(), 23 axis_text=element_text(size=12), 24 axis_title = element_text(size=14,weight="bold"), 25 panel_grid_major_x=element_line(color="gray",size=.5), 26 panel_grid_major_y=element_line(color="gray",size=.5), 27 )) 28 idw_scatter
可視化結果如下：
「IDW插值結果繪制」
1 idw_scatter_inter = (ggplot() + 2 geom_tile(IDW_grid_df,aes(x='lon',y='lat',fill='idw_value'),size=0.1) + 3 geom_map(js,fill='none',color='gray',size=0.4) + 4 geom_point(pm,aes(x='經度',y='緯度',fill='PM2.5'),size=4,stroke=.3,show_legend=False) + 5 scale_fill_cmap(cmap_name='Spectral_r',name='idw_value', 6 breaks=[30,40,60,80] 7 )+ 8 scale_x_continuous(breaks=[117,118,119,120,121,122])+ 9 labs(title="Map Charts in Python Exercise 02: Map IDM point", 10 )+ 11 #添加文本信息 12 annotate('text',x=116.5,y=35.3,label="processed map charts with plotnine",ha="left", 13 size=10)+ 14 annotate('text',x=120,y=30.6,label="Visualization by DataCharm",ha="left",size=9)+ 15 theme( 16 text=element_text(family="Roboto Condensed"), 17 #修改背景 18 panel_background=element_blank(), 19 axis_ticks_major_x=element_blank(), 20 axis_ticks_major_y=element_blank(), 21 axis_text=element_text(size=12), 22 plot_title=element_text(size=15,weight="bold"), 23 axis_title = element_text(size=14), 24 panel_grid_major_x=element_line(color="gray",size=.5), 25 panel_grid_major_y=element_line(color="gray",size=.5), 26 )) 27 idw_scatter_inter
可視化結果如下：
這里加上了散點是為了更好的對比插值結果，不加的效果如下：
裁剪操作
對研究區域的結果進行裁剪，在之前的推文中我們介紹了很多次，這里主要使用geopandas的clip() 方法進行操作，具體過程不再贅述(可以看我之前的推文教程)，我們直接給出裁剪結果：
Basemap繪制
上面介紹了plotnine包進行繪制的，這里我們再使用Basemap進行繪制，直接給出繪圖代碼：
1 from mpl_toolkits.basemap import Basemap 2 3 jiangsu_shp = r"F:\DataCharm\shpfile_data\JS\江蘇省_行政邊界" 4 fig,ax = plt.subplots(figsize=(6,4.5),dpi=130,facecolor="white") 5 map_base = Basemap(llcrnrlon=js_box[0], urcrnrlon=js_box[2], llcrnrlat=js_box[1],urcrnrlat=js_box[3], 6 projection="cyl",lon_0 = 119,lat_0 = 33,ax = ax) 7 # lat = np.arange(30,36,1) 8 # lon = np.arange(116,122,1) 9 map_base.drawparallels(np.arange(30,36,1), labels=[1,0,0,0],fontsize=12,zorder=0) #畫緯度線 10 map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=0) #畫經度線 11 map_base.readshapefile(shapefile = jiangsu_shp, name = "Js", default_encoding="ISO-8859-1", 12 drawbounds=True) 13 cp=map_base.pcolormesh(xgrid, ygrid, data=idw_grid,cmap='Spectral_r') 14 #ct=map_base.contour(xgrid, ygrid, data=idw_grid,colors='w',linewidths=.7) 15 #添加散點 16 vmin = pm["PM2.5"].min() 17 vmax = pm["PM2.5"].max() 18 ax.scatter(pm['經度'],pm["緯度"],c=pm["PM2.5"],s=90,ec="k",lw=0.5,cmap="Spectral_r", 19 vmin=vmin,vmax=vmax) 20 21 22 colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="IDW_inter") 23 #設置colorbar 24 colorbar.outline.set_edgecolor('none') 25 26 for spine in ['top','left','right','bottom']: 27 ax.spines[spine].set_visible(None) #隱去軸脊 28 29 ax.text(.5,1.1,"Map Charts in Python Exercise 02:Map IDW Grid Charts",transform = ax.transAxes,ha='center', 30 va='center',fontweight="bold",fontsize=14) 31 ax.text(.5,1.03, "processed map charts with Basemap", 32 transform = ax.transAxes,ha='center', va='center',fontsize = 10,color='black') 33 ax.text(.83,-.06,'\nVisualization by DataCharm',transform = ax.transAxes, 34 ha='center', va='center',fontsize = 8,color='black')
可視化結果如下：
裁剪操作
裁剪的操作也在之前的推文中介紹多次，這里我們直接給出結果哈：
當然，你也可以通過basemap.contour() 添加二維等值線，可視化結果如下：
總結
這是IDW插值的第一篇推文教程，好多原理的部分也沒有介紹，這里是自定義函數進行插值計算的，當然也是有優秀的第三方包可以完成。下次的R-ggplot2版本的IDW插值我們將使用現有的優秀三方包進行計算操作。文中有很多重復的知識點沒有詳細介紹，大家可以查看之前的推文，或者等這個系列完成后的詳細源碼、數據、解釋文檔的分享哈！(文中出錯的地方小伙伴們可以私聊指出或者加群討論哈)