您现在的位置是:首页 > 技术教程 正文

空间绘图 | Python-pykrige包-克里金(Kriging)插值计算及可视化绘制

admin 阅读: 2024-03-24
后台-插件-广告管理-内容页头部广告(手机)

前面两篇推文我们分别介绍了使用Python和R进行IDW(反距离加权法) 插值的计算及结果的可视化过程,详细内容可见如下:

  • Python - IDW插值计算及可视化绘制

  • R-gstat-ggplot2 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属性,具体代码如下:

  1. js_box = js.geometry.total_bounds
  2. grid_lon = np.linspace(js_box[0],js_box[2],400)
  3. grid_lat = np.linspace(js_box[1],js_box[3],400)

这里我们还是设置400*400的网格,注意np.linspace()方法和上期中 R的seq() 的使用不同。除此之外,我们还需要获取已知站点的经纬度信息(lons、lats)和对应值(data),这里给出点数据预览,如下:

获取数据代码如下:

  1. lons = nj_data["经度"].values
  2. lats = nj_data["纬度"].values
  3. data = nj_data["PM2.5"].values

pykrige包计算插值结果

在进过上面的数据处理过程后,我们已经构建出符合pykrige包进行插值计算所需的全部参数数据,接下来,我们直接调用即可,具体操作代码如下:

  1. from pykrige.ok import OrdinaryKriging
  2. OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian',nlags=6)
  3. 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)插值结果可视化绘制

这里都是常用的方法了,我们直接给出代码,大家不懂的可以查看之前的文章哈。

  1. #转换成网格
  2. xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
  3. #将插值网格数据整理
  4. df_grid =pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))
  5. #添加插值结果
  6. df_grid["Krig_gaussian"= Krig_result

结果如下(部分):

plotnine 可视化绘制

到这一步就很简单了,我们直接给出绘图代码即可,这里我们先给出散点的可视化结果,方便对比插值结果。

「散点结果」

「克里金(Kriging)插值结果」 绘图代码如下:

  1. import plotnine
  2. from plotnine import *
  3. plotnine.options.figure_size = (54.5)
  4. Krig_inter_grid = (ggplot() + 
  5.            geom_tile(df_grid,aes(x='long',y='lat',fill='Krig_gaussian'),size=0.1+
  6.            geom_map(js,fill='none',color='gray',size=0.3+ 
  7.            scale_fill_cmap(cmap_name='Spectral_r',name='Krig_gaussian_result',
  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 point interpolation",
  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),
  25.                panel_grid_major_x=element_line(color="gray",size=.5),
  26.                panel_grid_major_y=element_line(color="gray",size=.5),
  27.             ))

可视化结果如下:

还是一样,使用geopandas的clip()方法进行裁剪操作,唯一和上面绘图不同的就是数据选取的不同,这里选择裁剪后的数据。我们直接给出裁剪结果,如下:

Basemap的插值结果绘制

我们直接给出绘图代码及必要的可视化结果,具体不解的地方,可以查看之前的文章,代码如下:

  1. from mpl_toolkits.basemap import Basemap
  2. jiangsu_shp = r"F:\DataCharm\shpfile_data\JS\江苏省_行政边界"
  3. fig,ax = plt.subplots(figsize=(6,4.5),dpi=130,facecolor="white")
  4. map_base = Basemap(llcrnrlon=js_box[0], urcrnrlon=js_box[2], llcrnrlat=js_box[1],urcrnrlat=js_box[3],
  5.                   projection="cyl",lon_0 = 119,lat_0 = 33,ax = ax)
  6. # lat = np.arange(30,36,1)
  7. # lon = np.arange(116,122,1)
  8. map_base.drawparallels(np.arange(30,36,1), labels=[1,0,0,0],fontsize=12,zorder=0) #画纬度线
  9. map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=0) #画经度线
  10. map_base.readshapefile(shapefile = jiangsu_shp, name = "Js"default_encoding="ISO-8859-1",
  11.                        drawbounds=True)
  12. cp=map_base.pcolormesh(xgrid, ygrid, data=z1.data,cmap='Spectral_r')  
  13. colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="Kriging_inter")
  14. #设置colorbar
  15. colorbar.outline.set_edgecolor('none')
  16. for spine in ['top','left','right','bottom']:
  17.     ax.spines[spine].set_visible(None) #隐去轴脊
  18. ax.text(.5,1.1,"Map Charts in Python Exercise 02:Map Kriging Grid line",transform = ax.transAxes,ha='center'
  19.         va='center',fontweight="bold",fontsize=14)
  20. ax.text(.5,1.03"processed map charts with Basemap",
  21.         transform = ax.transAxes,ha='center', va='center',fontsize = 10,color='black')
  22. ax.text(.83,-.06,'\nVisualization by DataCharm',transform = ax.transAxes,
  23.         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)插值计算和插值结果可视化展示。

标签:
声明

1.本站遵循行业规范,任何转载的稿件都会明确标注作者和来源;2.本站的原创文章,请转载时务必注明文章作者和来源,不尊重原创的行为我们将追究责任;3.作者投稿可能会经我们编辑修改或补充。

在线投稿:投稿 站长QQ:1888636

后台-插件-广告管理-内容页尾部广告(手机)
关注我们

扫一扫关注我们,了解最新精彩内容

搜索
排行榜