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 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 .
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,js_box,400) grid_lat = np.linspace(js_box,js_box,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
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 ：
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 .
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 )：
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 ：
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, urcrnrlon=js_box, llcrnrlat=js_box,urcrnrlat=js_box, 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') colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="Kriging_inter") # 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 ：
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 email@example.com Delete .
Original publication time ： 2020-12-21
Participation of this paper Tencent cloud media sharing plan , You are welcome to join us , share .