上篇推文我们介绍了使用Python-pykrige包实现了克里金(Kriging)插值计算及对应的可视化结果绘制，详细内容点击下方链接：Python-pykrige包-克里金(Kriging)插值计算及可视化绘制，相信你也感受到了Python的简单方便性。本期推文，我们就推出使用R-gstat包实现克里金(Kriging)插值的计算及对应结果的可视化绘制，主要知识点如下：
- gstat.krige()实现克里金插值计算
- 插值结果的可视化绘制
gstat.krige()实现克里金插值计算
model选择
我们之前使用gstat包进行IDW插值计算，本期的推文全部在上次操作的基础之上(可能有些地方大家不是很明白，这个系列结束，我会分享完整的源码、文档和数据的)，大家不明白的地方，也可以参考之前的文档(R-gstat-ggplot2 IDW计算及空间插值可视化绘制)。在进行krige计算之前，需要对数据进行“查看”，根据数据的分布情况选择合适的核函数进行拟合计算。代码如下：
semivariog<-variogram(PM2.5~1, locations=scatter_df, data=scatter_df)
plot(semivariog)
这里的scatter_df 都是经sp包转换后的地理类型数据，可得数据距离分布情况如下(根据距离分布选择合适的model)：
我们根据数据分布选择 model="Exp",使用如下代码进行拟合线的绘制：
model.variog<-vgm(psill=125, model="Exp", nugget=45, range=.6)
fit.variog<-fit.variogram(semivariog, model.variog)
plot(semivariog, fit.variog)
结果如下：
由于上篇Python的结果我们使用高斯函数，这里我们也选择相同的model,代码如下：
model.variog<-vgm(psill=125, model="Gau", nugget=45, range=.6)
fit.variog<-fit.variogram(semivariog, model.variog)
plot(semivariog, fit.variog)
结果如下：
注意：
- 以上psill=125、nugget=45, range=.6等参数则是根据数据分布情况进行合理设置。接下来我们就选择对应model进行克里金插值计算。
- 使用vgm()函数即可查看gstat包支持的model种类。
krige 计算
「model="Exp"」
首先，我们对model="Exp" 进行krige 计算，主要代码如下：
model.variog<-vgm(psill=125, model="Exp", nugget=45, range=.6)
#插值计算
krig<-krige(formula=PM2.5 ~ 1, locations=scatter_df, newdata=grid, model=model.variog)
#将结果转成df 数据类型
krig_output=as.data.frame(krig)
names(krig_output)[1:3]<-c("long","lat","OK_pred")
head(krig_output)
得到的插值结果如下(部分)：
「model="Gau"」
model.variog<-vgm(psill=125, model="Gau", nugget=45, range=.6)
#插值计算
krig_Gau<-krige(formula=PM2.5 ~ 1, locations=scatter_df, newdata=grid, model=model.variog)
krig_output_Gau=as.data.frame(krig_Gau)
names(krig_output_Gau)[1:3]<-c("long","lat","OK_Gau_pred")
head(krig_output_Gau)
结果如下(部分)：
接下来我们就这两种情况进行可视化绘制。
插值结果的可视化绘制
我们有了规整好的df类型数据，这就可以方便的使用ggplot2进行可视化绘制。
针对model="Exp"的结果
首先我们绘制网格数据可视化结果，代码如下：
library(sf) library(tidyverse) library(ggspatial) library(RColorBrewer) library(ggtext) library(hrbrthemes)
#自定义颜色
my_colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
OK_Map_title <- ggplot() +
geom_tile(data = krig_output,aes(x=long,y=lat,fill=OK_pred)) +
geom_sf(data = jiangsu,fill="NA",size=.5,color="gray40") +
#geom_sf(data = scatter_df_tro,aes(fill=PM2.5),shape=21,size=4,show.legend = FALSE)+
scale_fill_gradientn(colours = my_colormap)+
annotation_scale(location = "bl") +
# spatial-aware automagic north arrow
annotation_north_arrow(location = "tr", which_north = "false",
style = north_arrow_fancy_orienteering) +
labs(x="",y="",
title = "Map Charts in R Exercise 02: <span style='color:#D20F26'>Map OK Interpolation</span>",
subtitle = "processed map charts with <span style='color:#1A73E8'>geom_tile()</span>",
caption = "Visualization by <span style='color:#DD6449'>DataCharm</span>") +
theme_ipsum(base_family = "Roboto Condensed") +
theme(
plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
size = 20, margin = margin(t = 1, b = 12)),
plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
plot.caption = element_markdown(face = 'bold',size = 12),
)
可视化结果如下：
当然，还可以通过*geom_contour()*绘制二维等值线：
geom_contour(data = krig_output,aes(x=long,y=lat,z=OK_pred),colour="white")
可视化结果如下：
「裁剪操作」
这步骤说了很多遍了，这里直接给出的代码：
#需要对数据进行投影转换
OK_output_sf <- st_as_sf(krig_output,coords = c("long", "lat"),crs = 4326)
OK_mask_result <- sf::st_intersection(OK_output_sf,jiangsu,tolerance=.01)
#可视化绘制
OK_mask_map <- ggplot() +
geom_sf(data = OK_mask_result,aes(color=OK_pred)) +
geom_sf(data = jiangsu,fill="NA",size=.5,color="gray40",alpha=.4) +
scale_color_gradientn(colours = my_colormap)+
annotation_scale(location = "bl") +
# spatial-aware automagic north arrow
annotation_north_arrow(location = "tr", which_north = "false",
style = north_arrow_fancy_orienteering) +
labs(
title = "Map Charts in R Exercise 02: <span style='color:#D20F26'>Map OK Interpolation Mask</span>",
subtitle = "processed map charts with <span style='color:#1A73E8'>geom_sf()</span>",
caption = "Visualization by <span style='color:#DD6449'>DataCharm</span>") +
theme_ipsum(base_family = "Roboto Condensed") +
theme(
plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
size = 20, margin = margin(t = 1, b = 12)),
plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
plot.caption = element_markdown(face = 'bold',size = 12),
)
可视化结果如下：
注意：这里的裁剪方法存在一些问题：在面对较大面积时，裁剪处不够明确(出现过多的范围，这个问题在本系列结束时给出的完整文档中会给出解决方法)。
针对model="Gau"的结果
由于步骤类似，我们这里直接给出代码和绘图结果即可：
library(sf)
library(tidyverse)
library(ggspatial)
library(RColorBrewer)
library(ggtext)
library(hrbrthemes)
#自定义颜色
my_colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
OK_Map_title_Gau <- ggplot() +
geom_tile(data = krig_output_Gau,aes(x=long,y=lat,fill=OK_Gau_pred)) +
geom_sf(data = jiangsu,fill="NA",size=.5,color="gray40") +
geom_contour(data = krig_output_Gau,aes(x=long,y=lat,z=OK_Gau_pred),colour="white")+
#geom_sf(data = scatter_df_tro,aes(fill=PM2.5),shape=21,size=4,show.legend = FALSE)+
scale_fill_gradientn(colours = my_colormap)+
annotation_scale(location = "bl") +
# spatial-aware automagic north arrow
annotation_north_arrow(location = "tr", which_north = "false",
style = north_arrow_fancy_orienteering) +
labs(x="",y="",
title = "Map Charts in R Exercise 02: <span style='color:#D20F26'>Map OK Interpolation</span>",
subtitle = "processed map charts with <span style='color:#1A73E8'>geom_tile()</span>",
caption = "Visualization by <span style='color:#DD6449'>DataCharm</span>") +
theme_ipsum(base_family = "Roboto Condensed") +
theme(
plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
size = 20, margin = margin(t = 1, b = 12)),
plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
plot.caption = element_markdown(face = 'bold',size = 12),
)
可视化结果如下：
总结
到这里，R版本的克里金(Kriging)插值计算结果及可视化绘制就完成了，相比于Python-pykrige包计算的结果，由于计算及部分参数设置的不同，导致结果有所偏差，大家可以根据自己的实际情况进行选择(个人建议选择R版本的)。大家不理解的地方可以查看gstat官网。目前小编在制作类别空间插值可视化绘制(Categorical spatial interpolation),希望可以同时制作Python和R两个版本的，可能还会涉及到机器学习的内容，大家尽请期待哦！(有小伙伴问我是GIS方面的学生吗？为啥做空间可视化多一点？我的回答是”我只是数据可视化设计爱好者，当然，研究生期间也涉及空间GIS方法，最重要的是空间可视化作品比较炫！！“)
