python 克里金空间插值_Python克里金(Kriging)插值计算及可视化绘制
前面两篇推文我们分别介绍了使用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)插值计算及可视化绘制相关推荐
- js插值计算_Python IDW插值计算及可视化绘制
前面几篇推文我们分辨介绍了使用Python和R绘制了二维核密度空间插值方法,并使用了Python可视化库plotnine.Basemap以及R的ggplot2完成了相关可视化教程的绘制推文,详细内容如 ...
- 【数据可视化应用】Python反距离权重(IDW)插值计算及可视化绘制
本文我们将介绍IDW(反距离加权法(Inverse Distance Weighted)) 插值的Python计算方法及插值结果的可视化绘制过程.主要涉及的知识点如下: IDW简介 自定义Python ...
- python多行字符串变单行_Python 正则表达式里的单行s和多行m模式
Python 的 re 模块内置函数几乎都有一个 flags 参数,规定了正则匹配时的各种策略模式,其中有两个模式:单行(re.DOTALL, 或者re.S)和多行(re.MULTILINE, 或者r ...
- python中gt是什么意思_python代码里出现gt;gt;gt;是啥意思
一种开发高质量软件的方法是为每一个函数开发测试代码,并且在开发过程中经常运行这些测试代码. doctest模块提供一个工具,这个工具可以扫描一个模块并验证确认内嵌到程序中的文档字符串测试代码.测试构造 ...
- python 网格数据插值_python – 网格数据的快速插值
当然!有两个选项可以做不同的事情,但是既能利用原始数据的定期网格性质. 第一个是scipy.ndimage.zoom.如果你只想通过内插原始数据生成一个更加密集的规则网格,那就是要走的路. 第二个是s ...
- 用python做一个数据查询软件_Python实现功能简单的数据查询及可视化系统
欢迎点击右上角关注小编,除了分享技术文章之外还有很多福利,私信学习资料可以领取包括不限于Python实战演练.PDF电子文档.面试集锦.学习资料等. image.png 前言 数据时代,数据的多源集成 ...
- python 怎么得到图像深度图 软件_Python/OpenCV:从立体图像计算深度图
我有两个立体图像要用来计算深度图.虽然我不幸不知道C/C++,但我知道Python--所以当我发现this tutorial时,我是乐观的. 不幸的是,教程似乎有些过时了.它不仅需要调整以运行(将&q ...
- python中delta是什么意思_python – 根据dataframe中的值计算delta
我有这个DataFrame(这只是一个例子,而不是真实的数据): In [1]: import pandas as pd my_data = [{'client_id' : '001', 'items ...
- Kriging插值类毕业论文文献有哪些?
本文是为大家整理的Kriging插值主题相关的10篇毕业论文文献,包括5篇期刊论文和5篇学位论文,为Kriging插值选题相关人员撰写毕业论文提供参考. 1.[期刊论文]基于Kriging插值的校园内 ...
最新文章
- linux 查看强制位,linux强制位与冒险位
- pyqt5必须和python对应_python 使用PyQt5
- ssh(Spring+Spring mvc+hibernate)——EmpController
- 前端性能优化:Add Expires headers
- mysql 拼接符是什么_mysql 字符串拼接
- ubuntu下c 调用java_ubuntu下使用JNI Java调用C++的例子
- .Net报文请求转义
- MYSQL 引擎的情况
- php判断ie的内核,js判断浏览器版本以及浏览器内核的方法_javascript技巧
- Windows设置程序开机自启动的几种方法(整理发布)
- Javaeve博客教你怎么发带图片的博客,非其他网络连接图片
- docker 安装RabbitMQ(镜像安装)
- 关于Linux系统重启过慢问题解决方案
- java 进阶笔记线程与并发之ForkJoinPool简析
- 夜里走了很多路,醒来还是在床上
- Python3 调用谷歌翻译
- 怎么办理质量管理体系认证证书ISO9001?
- PCB板框的绘制——AD19
- Mysql数据库练习题之商品库
- 【Unity3D】3D 视图操作 ( 视图基本元素 | 导航器 | 栅格 | 天空盒 | 3D 视图操作 | 视图旋转 | 视图缩放 | 视图平移 | 导航器操作 | 恢复方向 | 顶、右、前视图 )