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

最新推荐文章于 2025-02-26 14:49:03 发布

数据魅力(DataCharm)

最新推荐文章于 2025-02-26 14:49:03 发布

阅读量1.8k

点赞数 27

分类专栏：数据可视化文章标签： r语言开发语言数据可视化 python

本文链接：https://blog.csdn.net/qq_40483688/article/details/135776458

版权

数据可视化专栏收录该内容

39 篇文章

订阅专栏

本文介绍了如何在R中利用R-gstat包进行IDW插值计算，配合ggplot2进行可视化，包括数据预处理、网格构建、IDW计算以及结果的可视化展示。作者还比较了Python和R在空间插值方面的应用差异。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

上一篇文章，我们使用了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的两个版本哦，大家敬请期待！