Python pykrige package Kriging interpolation calculation and visual rendering

DataCharm 2021-02-22 23:08:09
python pykrige package kriging interpolation

In the first two tweets, we introduced the use of Python and R Conduct IDW( Inverse distance weighting method ) The calculation of interpolation and the visualization of the results , The details are as follows :

This issue tweets , We will show you how to use Python Make Kriging (Kriging) Interpolation calculation and visualization of interpolation results . The main knowledge points involved are as follows :

  • Kriging (Kriging) Introduction to interpolation
  • Python-pykrige Application of kukriging interpolation
  • Kriging (Kriging) The interpolation results are visualized

Kriging (Kriging) Introduction to interpolation

Kriging method (Kriging) Is based on covariance function for random processes / Spatial modeling and prediction of random fields ( interpolation ) The regression algorithm of . In a particular random process , For example, in an inherently stationary process , Kriging method can give the optimal linear unbiased estimation (Best Linear Unbiased Prediction, BLUP), Therefore, it is also called spatial optimal unbiased estimator in geostatistics (spatial BLUP)( The above definition comes from the network ). still IDW Interpolation is the same as , Let's get rid of the tedious process of formula derivation , The schematic diagram is as follows :

(Kriging Interpolation diagram )

While using Python Conduct Kriging Interpolation calculation does not need to customize complex functions , Here we call directly pykrige Package progress Kriging Interpolation calculation , And all we have to do is figure out pykrige The package plug-in can calculate the required parameter data .

Interpolation grid making

Whether it's custom or call package , We all need to create a mesh of our interpolated regions (grid), The method is also very simple , First of all, according to the map file (js) Get its latitude and longitude range , Here we use geopandas Read geojson Map file , And get the total_bounds attribute , The specific code is as follows :

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)

Here we still set up 400*400 The grid of , Be careful np.linspace() Methods and the last issue R Of seq() The use of is different . besides , We also need to get the latitude and longitude information of known sites (lons、lats) And the corresponding value (data), Here's a preview of the point data , as follows :

The data code is as follows :

lons = nj_data[" longitude "].values
lats = nj_data[" latitude "].values
data = nj_data["PM2.5"].values

pykrige Packet computing interpolation results

After the above data processing process , We've built a system that fits pykrige Package all parameter data needed for interpolation calculation , Next , We can call it directly , The specific operation code is as follows :

from pykrige.ok import OrdinaryKriging
OK = OrdinaryKriging(lons, lats, data, variogram_model='gaussian',nlags=6)
z1, ss1 = OK.execute('grid', grid_lon, grid_lat)

Be careful :

  • We use OrdinaryKriging Method for interpolation calculation , Besides , also UniversalKriging、RegressionKriging Interpolation method
  • variogram_model='gaussian', We set up Gauss (gaussian) Model , The result is linear (linear) The result may be different .
  • pykrige Provide linear, power, gaussian, spherical, exponential, hole-effect several variogram_model To choose from , The default is linear Model .

z1 It's our interpolation result , give the result as follows :

It turns out that , Its shape is 400*400( It's marked in the red box ), Next, we can visualize the interpolation results .

Kriging (Kriging) The interpolation results are visualized

Here are the common methods , We give the code directly , You don't understand can check the previous article ha .

# Convert to grid
xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
# Arrange the interpolation grid data
df_grid =pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))
# Add interpolation results
df_grid["Krig_gaussian"] = Krig_result

give the result as follows ( part ):

plotnine Visual rendering

It's easy to get to this point , We can give the drawing code directly , Here we first give the visualization results of scatter points , It is convenient to compare interpolation results .

「 Scatter results 」

「 Kriging (Kriging) Interpolation results 」 The drawing code is as follows :

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) +
labs(title="Map Charts in Python Exercise 02: Map point interpolation",
# Add text message
annotate('text',x=116.5,y=35.3,label="processed map charts with plotnine",ha="left",
annotate('text',x=120,y=30.6,label="Visualization by DataCharm",ha="left",size=9)+
text=element_text(family="Roboto Condensed"),
# Change the background
axis_title = element_text(size=14),

The visualization results are as follows :

Still the same , Use geopandas Of clip() Method for clipping operation , The only difference from the above plot is the data selection , Here we choose the cropped data . We're going to give you the tailoring directly , as follows :

Basemap The result of interpolation is drawn

We give the drawing code and necessary visualization results directly , What's really puzzling is , You can see the previous article , The code is as follows :

from mpl_toolkits.basemap import Basemap
jiangsu_shp = r"F:\DataCharm\shpfile_data\JS\ Jiangsu Province _ Administrative boundaries "
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) # Draw the latitude line
map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=0) # Draw a longitude line
map_base.readshapefile(shapefile = jiangsu_shp, name = "Js", default_encoding="ISO-8859-1",
cp=map_base.pcolormesh(xgrid, ygrid,,cmap='Spectral_r')
colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="Kriging_inter")
# Set up colorbar
for spine in ['top','left','right','bottom']:
ax.spines[spine].set_visible(None) # Remove the axial ridge
ax.text(.5,1.1,"Map Charts in Python Exercise 02:Map Kriging Grid line",transform = ax.transAxes,ha='center',
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')

The visualization results are as follows :

You can also use :

ct=map_base.contour(xgrid, ygrid,,colors='w',linewidths=.7)

Add two-dimensional contours , give the result as follows :

The tailoring results are visualized as follows :

Add isoline results :


Come here ,Python Kriging (Kriging) The interpolation calculation method and the visualization of the results are introduced , Or that sentence , Have a ready-made “ wheel ” It can be used , Use it as much as possible ( Of course , Except for highly customized requirements ), Besides , It is also very important to understand the calculation principle . Next , We will introduce the use of R Language and its excellent third package carry Kriging (Kriging) Interpolation calculation and visualization of interpolation results .

This article is from WeChat official account. - DataCharm(shujumeili) , author : Ning Haitao

The source and reprint of the original text are detailed in the text , If there is any infringement , Please contact the Delete .

Original publication time : 2020-12-21

Participation of this paper Tencent cloud media sharing plan , You are welcome to join us , share .


  1. 使用Python开发DeFi项目
  2. python 函数详解
  3. Python工程师是做什么的?前景如何?
  4. Python - zip() 函数
  5. 30 周年生日,Python 先驱是怎么评价这门语言的?
  6. python将excel自适应导入数据库
  7. 从小白到大师,这里有一份Pandas入门指南
  8. [Python] 茎叶图和复合饼图的画法
  9. [Python interface automation] - regular use case parameterization
  10. Translation: practical Python Programming 02_ 02_ Containers
  11. Two years of Java, to write Python and go
  12. Translation: practical Python Programming 02_ 02_ Containers
  13. Two years of Java, to write Python and go
  14. Python-geoplot 空间核密度估计图绘制
  15. Python-seaborn 经济学人经典图表仿制
  16. python空间绘图- regionmask掩膜操作示例
  17. Python 空间绘图 - Cartopy 经纬度添加
  18. Python-pykrige包-克里金(Kriging)插值计算及可视化绘制
  19. Python 批量重采样、掩膜、坡度提取
  20. python - 多种交通方式可达圈分析
  21. Python 空间绘图 - 房价气泡图绘制
  22. Translation: practical Python Programming 02_ 02_ Containers
  23. Research on Portfolio Optimization Based on particle swarm optimization
  24. Ubuntu deploying Django project
  25. Two years of Java, write Python and go without byte beating
  26. Translation: practical Python Programming 02_ 02_ Containers
  27. So learn python, grandfather learned! Introduction to super simple Python
  28. python3 多线程 与 mongo亿级消费日志数据 新鲜demo 【优化第一版】
  29. Summary of Chinese word segmentation based on Jieba
  30. I've heard it n times, but I'm not impressed. After reading this, you'll understand
  31. Summary of Chinese word segmentation based on Jieba
  32. From movie art to Python code to realize God's reverse thinking mode
  33. Summary of Chinese word segmentation based on Jieba
  34. ARIMA模型预测CO2浓度时间序列-python实现
  35. Python belongs to back-end development or front-end development? Introduction to Python!
  36. python isinstance()
  37. I've heard it n times, but I'm not impressed. After reading this, you'll understand
  38. This article will familiarize you with the transformation process of Python - & gt; cafe - & gt; om model
  39. 如何用Python一键修改上万个文件名
  40. One day quick start to Python
  41. Python 学习笔记: List
  42. 翻译:《实用的Python编程》02_03_Formatting
  43. Is there any age requirement for learning Python? Is 30 OK?
  44. Professor Tsinghua! The most complete Python tutorial in 12 hours (free sharing at the end of the article)
  45. Using Python to develop defi project
  46. Detailed explanation of Python function
  47. Python 可变类型作为函数默认参数时的副作用
  48. What do Python engineers do? What's their future?
  49. 这是我见过最好的Python教程:十分钟带你认识Python
  50. Python欢喜冤家:爬虫与反爬虫带着处理方案来给大家拜年了
  51. Python - zip() function
  52. 写Python会遇到如下的错误:ModuleNotFoundError: No module named 'email.mime'; 'email' is not a package
  53. Python类的调用以及私有和公有属性方法的调用
  54. Python类的专有方法
  55. Python基础之:数字字符串和列表
  56. How did Python pioneers evaluate this language on their 30th birthday?
  57. Python基础之:数字字符串和列表
  58. Python基础之:数字字符串和列表
  59. 窥探未来不是梦,python数据分析轻松实现
  60. This article will familiarize you with the transformation process of Python - & gt; cafe - & gt; om model