绘图提高篇 | R-gstat-ggplot2 IDW计算及空间插值可视化

上一篇文章，我们使用了Python 自定义IDW插值函数进行了IDW空间插值及可视化的plotnine、Basemap的绘制方法，本期推文我们将使用R-gstat进行IDW插值计算和使用ggplot2进行可视化绘制，主要涉及的知识点如下：

• R-gstat包IDW插值计算

• R-ggplot2 IDW插值结果可视化绘制

R-gstat包IDW插值计算

得益于优秀且丰富的R语言第三方包，我们可以直接使用空间统计计算的R-gstat包实现包括IDW在内的多种插值方法，使用R-sf包完美绘制空间可视化绘制。还是老样子，我们对所需数据(散点值+地图数据）的基本情况进行预览，结果如下：

绘图数据预览

• 散点情况(scatter_df)

  

• 地图文件(jiangsu)

  

  散点在地图上的可视化结果如下(之前推文已有过操作方法，不明白的小伙伴可以参考下)：

  

gstat-IDW计算

接下来，我们使用gstat包进行IDW计算，在计算之前，我们需使用sp包对数据进行相关处理，具体操作如下：

将数据准换成空间数据格式：

sp::coordinates(scatter_df)?=?~经度+纬度

数据格式准换如下：

「制作网格数据」： 更具地图文件的经纬度范围信息，我们进行网格(grid)的构建(我们需要插值的区域)，代码如下：

sf::st_bbox(jiangsu)
#生成400*400的网格
grid?<-?expand.grid(x=seq(from?=?st_bbox(jiangsu)[1],to?=?st_bbox(jiangsu)[3],length.out?=?400),
????????????????????y=seq(from?=?st_bbox(jiangsu)[2],to?=?st_bbox(jiangsu)[4],length.out?=?400))

这样就可以完成了400x400的网格点构建，接下来要将构建的网格点转换成空间数据格式，还是使用sp包操作，代码如下：

sp::coordinates(grid)?<-?~x+y
sp::gridded(grid)?<-?TRUE

以上操作我们就完成了gstat包所需的所有数据处理操作。

「IDW 插值操作」

由于有现成的函数可以调用，这里我们直接使用，代码如下：

idw?<-?idw(formula?=?PM2.5?~?1,locations?=?scatter_df,?newdata=grid)

解释如下：

• formula = PM2.5~ 1 为固定样式，用于需要进行插值的纬度数据。

• locations = scatter_df 为sp包处理过的样例点位置信息。

• newdata=grid 为需要插值的网格大小，sp对象。 由于计算的idw结果较多(400*400),这里将其转换成df数据类型，同时对列进行重命名，也为以后的绘图提供方便，转换代码如下：

idw_output?<-?as.data.frame(idw)
names(idw_output)[1:3]<-c("long","lat","IDW_Result")

结果如下(部分)：

ggplot2 可视化IDW插值结果

经过上面的数据规整，我们直接可以进行可视化操作，代码如下：

library(sf)
library(tidyverse)
library(ggspatial)
library(RColorBrewer)
library(ggtext)
library(hrbrthemes)

#自定义颜色
my_colormap?<-?colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
IDW_Map_title?<-?ggplot()?+?
??geom_tile(data?=?idw_output,aes(x=long,y=lat,fill=IDW_Result))?+
??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(
???????title?=?"Map?Charts?in?R?Exercise?02:?<span style='color:#D20F26'>Map?IDW?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),
??)

注意：这里我们将散点绘制在插值结果之上，也为了查看插值的效果，最终的可视化结果如下：

sf包裁剪操作

上面的可视化结果只是将网格插值结果全部绘制出来，没有将目标区域进行单独绘制(地图文件)，这里使用sf::st_intersection() 函数进行实现“裁剪”操作，这里不再赘述，不明白的可以查看我之前的推文。最终的可视化结果如下：

注意： 小伙伴们可能也发现了，这样裁剪的结果不是完全的按照地图文件进行裁剪的，会有部分“溢出”，特别是在绘制较大范围的空间图表的时候，这里可以转换成栅格数据，然后再使用*mask()*方法也是可以操作的，具体其他的，后面我会专门出一期推文介绍这两种"裁剪"操作的不同，这里先不做介绍。

总结

继上期我们推出Python 版本的IDW 空间插值之后，本期我们又推出了R版本的，大家可以对比下两种插值的结果(可能会存在些许的不同)。还是那句话，在绘制空间图表时，R因其完整的绘图体系及优秀的第三方包，可以较好的完成绘图需求(各种空间绘图元素的添加)，但Python因其简单好学，也具有一定优势，大家可以选择适合自己的方法进行学习，至于对比两种语言绘图不同，就交给小编来做吧。下期，我们继续空间插值(克里金：Kriging)的计算及可视化绘制，还是Python和R的两个版本哦，大家敬请期待！