## 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) +
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",
)+
# Add text message
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"),
# Change the background
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),
))
```

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",
drawbounds=True)
cp=map_base.pcolormesh(xgrid, ygrid, data=z1.data,cmap='Spectral_r')
# Set up colorbar
colorbar.outline.set_edgecolor('none')
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',
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')
```

The visualization results are as follows ：

You can also use ：

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

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

The tailoring results are visualized as follows ：

Add isoline results ：

## summary

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 .