前面幾篇推文我們分辨介紹了使用Python和R繪製了二維核密度空間插值方法，並使用了Python視覺化庫plotnine、Basemap以及R的ggplot2完成了相關視覺化教程的繪製推文，詳細內容如下：

• Python-plotnine 核密度空間插值視覺化繪製
• Python-Basemap核密度空間插值視覺化繪製
• R-ggplot2+sf 核密度空間插值視覺化繪製

接下來，我們將繼續介紹空間插值的其他方法，本期推文，我們將介紹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插值我們將使用現有的優秀三方包進行計算操作。文中有很多重複的知識點沒有詳細介紹，大家可以檢視之前的推文，或者等這個系列完成後的詳細原始碼、資料、解釋文件的分享哈！(文中出錯的地方小夥伴們可以私聊指出或者加群討論哈)