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

DataCharm 2021-02-22 13:02:30
Python 克里 python-pykrige pykrige 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)插值计算和插值结果可视化展示。

本文分享自微信公众号 - DataCharm(shujumeili) ,作者:宁海涛

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间: 2020-12-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

版权声明
本文为[DataCharm]所创,转载请带上原文链接,感谢
https://cloud.tencent.com/developer/article/1790265

  1. Python 3 entry, see this is enough
  2. 华为大佬打造的400集Python视频学起来,学完万物皆可爬
  3. 400 episodes of Python video created by Huawei boss
  4. django之csrf_exempt解决跨域请求的问题
  5. CSRF of Django_ Exempt solves the problem of cross domain requests
  6. 1.7 万 Star!一个简单实用的 Python 进度条库
  7. 17000 stars! A simple and practical Python progress bar library
  8. Python爬虫:设置Cookie解决网站拦截并爬取蚂蚁短租
  9. Python crawler: setting cookie to solve website interception and crawling ant short rent
  10. Python-Net编程
  11. Python net programming
  12. 学习Python数学英语基础重要吗?Python教程!
  13. Is it important to learn the basics of math and English in Python!
  14. Python数据分析常用库有哪些?Python学习!
  15. What are the common libraries for Python data analysis? Learn Python!
  16. win 创建python虚拟环境
  17. Creating Python virtual environment with win
  18. In order to automatically collect B station barrage, I developed a tool in Python
  19. 用Python编程语言来实现阿姆斯特朗数的检查
  20. Using python programming language to check Armstrong number
  21. Python中的解决中文字符编码的问题
  22. Solving the problem of Chinese character coding in Python
  23. Translation: practical Python Programming 02_ 01_ Datatypes
  24. Installation and use of Python and tensorflow in win10 environment (Python version 3.6, tensorflow version 1.6)
  25. Python series 46
  26. Linux安装Python3
  27. 【python接口自动化】- 正则用例参数化
  28. Python RestFul Api 设计
  29. filecmp --- 文件及目录的比较│Python标准库
  30. Installing python3 on Linux
  31. [Python] Matplotlib 圖表的繪製和美化技巧
  32. (資料科學學習手札108)Python+Dash快速web應用開發——靜態部件篇(上)
  33. 翻譯:《實用的Python程式設計》02_01_Datatypes
  34. 【python接口自动化】- 正则用例参数化
  35. 翻译:《实用的Python编程》02_02_Containers
  36. 两年Java,去字节跳动写Python和Go
  37. [Python interface automation] - regular use case parameterization
  38. Python restful API design
  39. 翻译:《实用的Python编程》02_02_Containers
  40. 两年Java,去字节跳动写Python和Go
  41. 翻译:《实用的Python编程》02_02_Containers
  42. Python基于粒子群优化的投资组合优化研究
  43. ubuntu部署django项目
  44. 兩年Java,去位元組跳動寫Python和Go
  45. 翻譯:《實用的Python程式設計》02_02_Containers
  46. 这样学习Python,爷爷都学会了!超简单Python入门
  47. [Python] 基于 jieba 的中文分词总结
  48. 【python】递归听了N次也没印象,读完这篇你就懂了
  49. [Python] 基于 jieba 的中文分词总结
  50. 人理解迭代,神则体会递归,从电影艺术到Python代码实现神的逆向思维模式
  51. [Python] 基於 jieba 的中文分詞總結
  52. Python属于后端开发还是前端开发?Python入门!
  53. 【python】递归听了N次也没印象,读完这篇你就懂了
  54. 一天快速入门python
  55. 学习Python对年龄有没有要求?30岁可以吗?
  56. 清华教授!12小时整理的最全Python教程(文末无偿分享)
  57. Filecmp -- comparison of files and directories
  58. Drawing and beautifying skills of [Python] Matplotlib chart
  59. Python + dash rapid web application development static components
  60. Translation: practical Python Programming 02_ 01_ Datatypes