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

你家的邻居

前面两篇推文我们分别介绍了使用Python和R进行IDW(反距离加权法)插值的计算及结果的可视化过程，详细内容可见如下：Python-IDW插值计算及可视化绘制R-gstat-ggplot2IDW计算及空间插值可视化绘制（需修改链接）本期推文，我们将介绍如何使用Python进行克里金(Kriging)插值计算及插值结果的可视化绘制。主要涉及的知识点如下：克里金(Kriging)插值简介Python-pykrige库克里金插值应用克里金(Kriging)插值结果可视化绘制克里金(Kriging)插值简介克里金法（Kriging）是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计（BestLinearUnbiasedPrediction,BLUP），因此在地统计学中也被称为空间最优无偏估计器(spatialBLUP)(以上定义来自于网络)。还是IDW插值介绍一样，我们省去繁琐的公式推导过程，示意图如下：(Kriging插值示意图)而使用Python进行Kriging插值计算无需自定义复杂的函数，这里我们直接调用pykrige包进行Kriging插值计算，二我们所要作的就是将计算出pykrige包插件计算所需要的参数数据即可。插值网格制作无论十自定义还是调用包，我们都需要制作出我们插值区域的网格(grid)，方法也十分简单，首先根据地图文件(js)获取其经纬度范围，这里我们使用geopandas读取geojson地图文件，并获取total_bounds属性，具体代码如下：js_box = js.geometry.total_boundsgrid_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["经度"].valueslats = nj_data["纬度"].valuesdata = nj_data["PM2.5"].valuespykrige包计算插值结果在进过上面的数据处理过程后，我们已经构建出符合pykrige包进行插值计算所需的全部参数数据，接下来，我们直接调用即可，具体操作代码如下：from pykrige.ok import OrdinaryKrigingOK = 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 plotninefrom 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 ython 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 Basemapjiangsu_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")#设置colorbarcolorbar.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 ython 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)插值计算和插值结果可视化展示。