python 克里金空间插值_Python-pykrige包-克里金(Kriging)插值计算及可视化绘制
前面两篇推文我们分别介绍了使用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属性,具体代码如下: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包进行插值计算所需的全部参数数据,接下来,我们直接调用即可,具体操作代码如下:frompykrige.okimportOrdinaryKriging
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)插值结果」绘图代码如下:importplotnine
fromplotnineimport*
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=
"MapChartsinPythonExercise02:Mappointinterpolation",
)+
#添加文本信息
annotate(
"text",x=116.5,y=35.3,label=
"processedmapchartswithplotnine",ha=
"left",
size=10)+
annotate(
"text",x=120,y=30.6,label=
"VisualizationbyDataCharm",ha=
"left",size=9)+
theme(
text=element_text(family=
"RobotoCondensed"),
#修改背景
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的插值结果绘制
我们直接给出绘图代码及必要的可视化结果,具体不解的地方,可以查看之前的文章,代码如下:frommpl_toolkits.basemapimportBasemap
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")
forspine
in[
"top",
"left",
"right",
"bottom"]:
ax.spines[spine].set_visible(None)
#隐去轴脊
ax.text(.5,1.1,
"MapChartsinPythonExercise02:MapKrigingGridline",transform=ax.transAxes,ha=
"center",
va=
"center",fontweight=
"bold",fontsize=14)
ax.text(.5,1.03,
"processedmapchartswithBasemap",
transform=ax.transAxes,ha=
"center",va=
"center",fontsize=10,color=
"black")
ax.text(.83,-.06,
"\nVisualizationbyDataCharm",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-pykrige包-克里金(Kriging)插值计算及可视化绘制相关推荐
- python 克里金空间插值_Python克里金(Kriging)插值计算及可视化绘制
前面两篇推文我们分别介绍了使用Python和R进行IDW(反距离加权法) 插值的计算及结果的可视化过程,详细内容可见如下: 本期推文,我们将介绍如何使用Python进行克里金(Kriging)插值计算 ...
- js插值计算_Python IDW插值计算及可视化绘制
前面几篇推文我们分辨介绍了使用Python和R绘制了二维核密度空间插值方法,并使用了Python可视化库plotnine.Basemap以及R的ggplot2完成了相关可视化教程的绘制推文,详细内容如 ...
- 【数据可视化应用】Python反距离权重(IDW)插值计算及可视化绘制
本文我们将介绍IDW(反距离加权法(Inverse Distance Weighted)) 插值的Python计算方法及插值结果的可视化绘制过程.主要涉及的知识点如下: IDW简介 自定义Python ...
- python判断是否包含某数字_python如何判断数组里是否有某个数字
Python 3语言开发教程.Python 3语言in操作符使用教程.Python 3语言如何判断数组内是否存在某一个元素? 工具/原料 电脑 Editplus 方法/步骤 1 在Python语言中, ...
- python爬虫beautifulsoup爬当当网_Python爬虫包 BeautifulSoup 递归抓取实例详解_python_脚本之家...
Python爬虫包 BeautifulSoup 递归抓取实例详解 概要: 爬虫的主要目的就是为了沿着网络抓取需要的内容.它们的本质是一种递归的过程.它们首先需要获得网页的内容,然后分析页面内容并找到 ...
- python导入自己写的模块_Python:包、模块和导入
南京著名风景区--牛首山 预计阅读时间--5分钟 Python的好处在于,你不需要懂很多概念,你就有机会投入工作,同样,问题也有机会随时发生. 包.模块 foo --bar.py --sim.py - ...
- python用循环打出阶梯图形_Python制图你真的会吗?一文学会如何绘制漂亮的阶梯图...
说到Python制图就不得不提matplotlib这个最为常用的库,matplotlib库作为Python经典的二维绘图库,在Python的数据可视化方面是最为常用的,今天呢,咱们接着上次和大家所探讨 ...
- python 编译成exe vmp加密_Python vmp包_程序模块 - PyPI - Python中文网
vmpy是评估typical的工具箱. 骑行性能指标来自骑行数据,如功率.心率.速度, 梯度,节奏流. 包中的所有函数都遵循惯例,其中输入/输出 格式要么是传统的python内置数据结构 或者是nd数 ...
- python爬取历史天气查询_python爬虫爬取各个城市历史天气及数据可视化
import asyncio import aiohttp from lxml import etree import re from collections import namedtuple Ar ...
- ArcGIS中使用协同克里金插值(co-kriging interplotation )对气象数据插值
ArcGIS中如何使用协同克里金插值(co-kriging interplotation )对气象数据插值 ANUSPLIN气象站点数据插值局限性 百度搜索ArcGIS 克里金插值 搭建梯子搜索Arc ...
最新文章
- r语言读取csv文件赋值gamma_tidyfst vs pandas(1):csv文件读写
- Linux环境下进入MySQL环境报权限问题:Access denied for user 'root@localhost' (using password:YSE)...
- nofollow标签_如何Nofollow外链
- 华为nove计算机在哪里,华为Nova手机备忘录怎么同步到电脑
- 游戏开发者怎么做出以假乱真的画面效果?大气散射渲染了解一下
- python执行gradle脚本
- 收件服务器主机名未响应,邮箱收件服务器主机名是什么
- Spring Boot Logback 配置详解
- matplotlib--python的数据可视化入门
- 【技术认证题库】SCCA理论aDesk-2考试【初级】
- 轻量级的java HTTP Server——NanoHttpd
- 计算机控制手机源码,Total Control电脑控制手机助手
- 风吹雨PHP多应用授权系统【开源】
- HP惠普服务器做RAID
- Chloe Orm的使用(一)
- Spring Cloud Stream报错:Invalid bean definition with name:bean definition with this name already exist
- 第五次作业:《国际贸易学》—WTO及区域经济一体化
- Cognos BI交流群共享资料整理
- php下对中国内地身份证进行验证
- 学生上课睡觉班主任怎么处理_学生上课睡觉老师该怎么办:高级教师教你一招...