前面两篇推文我们分别介绍了使用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_dat