前面兩篇推文我們分別介紹了使用Python和R進行IDW(反距離加權法) 插值的計算及結果的可視化過程，詳細內容可見如下：

本期推文，我們將介紹如何使用Python進行克里金(Kriging)插值計算及插值結果的可視化繪制。主要涉及的知識點如下：

克里金(Kriging)插值簡介

Python-pykrige庫克里金插值應用

克里金(Kriging)插值結果可視化繪制

克里金(Kriging)插值簡介

克里金法(Kriging) 是依據協方差函數對隨機過程/隨機場進行空間建模和預測(插值)的回歸算法。在特定的隨機過程，例如固有平穩過程中，克里金法能夠給出最優線性無偏估計(Best Linear Unbiased Prediction, BLUP)，因此在地統計學中也被稱為空間最優無偏估計器(spatial BLUP)(以上定義來自于網絡)。還是IDW插值介紹一樣，我們省去繁瑣的公式推導過程，示意圖如下：

? ? ? ? ? ? ? ? ? ? ? ? (Kriging插值示意圖)

而使用Python進行Kriging插值計算無需自定義復雜的函數，這里我們直接調用pykrige包進行Kriging插值計算，而我們所要做的就是將計算出pykrige包插件計算所需要的參數數據即可。

插值網格制作

無論是自定義還是調用包，我們都需要制作出我們插值區域的網格(grid)，方法也十分簡單，首先根據地圖文件(js)獲取其經緯度范圍，這里我們使用geopandas讀取geojson 地圖文件，并獲取total_bounds屬性，具體代碼如下：

js_box?=?js.geometry.total_bounds

grid_lon?=?np.linspace(js_box[0],js_box[2],400)

grid_lat?=?np.linspace(js_box[1],js_box[3],400)

這里我們還是設置400*400的網格，注意np.linspace()方法和上期中 R的seq() 的使用不同。除此之外，我們還需要獲取已知站點的經緯度信息(lons、lats)和對應值(data)，這里給出點數據預覽，如下：

獲取數據代碼如下：

lons?=?nj_data["經度"].values

lats?=?nj_data["緯度"].values

data?=?nj_data["PM2.5"].values

pykrige包計算插值結果

在經過上面的數據處理過程后，我們已經構建出符合pykrige包進行插值計算所需的全部參數數據，接下來，我們直接調用即可，具體操作代碼如下：

from?pykrige.ok?import?OrdinaryKriging

OK?=?OrdinaryKriging(lons,?lats,?data,?variogram_model='gaussian',nlags=6)

z1,?ss1?=?OK.execute('grid',?grid_lon,?grid_lat)

注意：

我們使用

OrdinaryKriging方法進行插值計算，此外，還有

UniversalKriging、RegressionKriging插值方法

variogram_model='gaussian'，我們設置高斯(gaussian)模型，其結果和一般的線性(linear)結果可能會有不同。

pykrige提供

linear, power, gaussian, spherical, exponential, hole-effect幾種variogram_model可供選擇，默認的為linear模型。

z1就是我們的插值結果，結果如下：

結果可以看出，其形狀為400*400(紅框中標出)，接下來我們進行插值結果的可視化繪制。

克里金(Kriging)插值結果可視化繪制

這里都是常用的方法了，我們直接給出代碼，大家不懂的可以查看之前的文章哈。

#轉換成網格

xgrid,?ygrid?=?np.meshgrid(grid_lon,?grid_lat)

#將插值網格數據整理

df_grid?=pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))

#添加插值結果

df_grid["Krig_gaussian"]?=?Krig_result

結果如下(部分)：

plotnine 可視化繪制

到這一步就很簡單了，我們直接給出繪圖代碼即可，這里我們先給出散點的可視化結果，方便對比插值結果。

「散點結果」

「克里金(Kriging)插值結果」繪圖代碼如下：

import?plotnine

from?plotnine?import?*

plotnine.options.figure_size?=?(5,?4.5)

Krig_inter_grid?=?(ggplot()?+

geom_tile(df_grid,aes(x='long',y='lat',fill='Krig_gaussian'),size=0.1)?+

geom_map(js,fill='none',color='gray',size=0.3)?+

scale_fill_cmap(cmap_name='Spectral_r',name='Krig_gaussian_result',

breaks=[30,40,60,80]

)+

scale_x_continuous(breaks=[117,118,119,120,121,122])+

labs(title="Map?Charts?in?Python?Exercise?02:?Map?point?interpolation",

)+

#添加文本信息

annotate('text',x=116.5,y=35.3,label="processed?map?charts?with?plotnine",ha="left",

size=10)+

annotate('text',x=120,y=30.6,label="Visualization?by?DataCharm",ha="left",size=9)+

theme(

text=element_text(family="Roboto?Condensed"),

#修改背景

panel_background=element_blank(),

axis_ticks_major_x=element_blank(),

axis_ticks_major_y=element_blank(),

axis_text=element_text(size=12),

axis_title?=?element_text(size=14),

panel_grid_major_x=element_line(color="gray",size=.5),

panel_grid_major_y=element_line(color="gray",size=.5),

))

可視化結果如下：

還是一樣，使用geopandas的clip()方法進行裁剪操作，唯一和上面繪圖不同的就是數據選取的不同，這里選擇裁剪后的數據。我們直接給出裁剪結果，如下：

Basemap的插值結果繪制

我們直接給出繪圖代碼及必要的可視化結果，具體不解的地方，可以查看之前的文章，代碼如下：

from?mpl_toolkits.basemap?import?Basemap

jiangsu_shp?=?r"F:\DataCharm\shpfile_data\JS\江蘇省_行政邊界"

fig,ax?=?plt.subplots(figsize=(6,4.5),dpi=130,facecolor="white")

map_base?=?Basemap(llcrnrlon=js_box[0],?urcrnrlon=js_box[2],?llcrnrlat=js_box[1],urcrnrlat=js_box[3],

projection="cyl",lon_0?=?119,lat_0?=?33,ax?=?ax)

#?lat?=?np.arange(30,36,1)

#?lon?=?np.arange(116,122,1)

map_base.drawparallels(np.arange(30,36,1),?labels=[1,0,0,0],fontsize=12,zorder=0)?#畫緯度線

map_base.drawmeridians(np.arange(116,122,1),?labels=[0,0,0,1],fontsize=12,zorder=0)?#畫經度線

map_base.readshapefile(shapefile?=?jiangsu_shp,?name?=?"Js",?default_encoding="ISO-8859-1",

drawbounds=True)

cp=map_base.pcolormesh(xgrid,?ygrid,?data=z1.data,cmap='Spectral_r')

colorbar?=?map_base.colorbar(cp,size='3%',pad="5%",label="Kriging_inter")

#設置colorbar

colorbar.outline.set_edgecolor('none')

for?spine?in?['top','left','right','bottom']:

ax.spines[spine].set_visible(None)?#隱去軸脊

ax.text(.5,1.1,"Map?Charts?in?Python?Exercise?02:Map?Kriging?Grid?line",transform?=?ax.transAxes,ha='center',

va='center',fontweight="bold",fontsize=14)

ax.text(.5,1.03,?"processed?map?charts?with?Basemap",

transform?=?ax.transAxes,ha='center',?va='center',fontsize?=?10,color='black')

ax.text(.83,-.06,'\nVisualization?by?DataCharm',transform?=?ax.transAxes,

ha='center',?va='center',fontsize?=?8,color='black')

可視化結果如下：

還可以通過：

ct=map_base.contour(xgrid,?ygrid,?data=z1.data,colors='w',linewidths=.7)

添加二維等值線，結果如下：

裁剪結果可視化如下：

添加等值線結果：

總結

到這里，Python的克里金(Kriging)插值計算方法及結果的可視化繪制就介紹完了，還是那句話，有現成的“輪子”可以用，大家盡量使用哈(當然，高度的定制化需求除外)，此外，懂得其計算原理也是很重要的哦。下一篇，我們將介紹使用R語言及其優秀的第三包進行克里金(Kriging)插值計算和插值結果可視化展示。

推薦閱讀

好文！必須在看

總結

以上是生活随笔為你收集整理的python 克里金空间插值_Python克里金(Kriging)插值计算及可视化绘制的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。