R-ggplot2+sf 核密度空间插值可视化绘制

上篇推文我们介绍了使用Python的plotnine、Basemap包对空间kde插值结果进行了可视化绘制，当然也包括了具体的插值过程，详细内容大家可以点击下方链接查看：Python-plotnine 核密度空间插值可视化绘制 、Python-Basemap核密度空间插值可视化绘制。

本期推文我们就介绍下使用R进行核密度估计、空间插值计算以及ggplot2+sf的可视化绘制操作。涉及的主要知识点如下：

• R-sm包计算核密度估计结果
• R-SP包转换网格插值结果
• R-ggplot2+sf包绘制网格插值结果
• R-sf包实现完美“裁剪”

R-sm包计算核密度估计结果

sf包散点位置可视化

首先，我们需要读取的数据还是地图文件(jiangsu)和点数据(scatter_df),散点数据预览如下：在计算核密度估计之前，我们先使用sf包进行散点的可视化绘制。首先我们使用sf包将其转换成sf地图对象，操作代码如下：

scatter_df_tro <- st_as_sf(scatter_df,coords = c("经度", "纬度"),crs = 4326)

结果如下(部分)：接下来，我们将使用sf包进行点位置的可视化绘制，代码如下：

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

#自定义颜色
my_colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
Map_point <- ggplot() + 
  geom_sf(data = jiangsu,fill="NA",size=.6,color="black") + 
  #coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+ #中国地图常用投影
  geom_sf(data = scatter_df_tro,aes(fill=PM2.5),shape=21,size=5) +
  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 01: <span style='color:#D20F26'>Map point</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),
  )

可视化结果如下：

sm包计算核密度估计结果

在上述可视化结果之后，我们需要根据已有的点进行核密度估计，在R中，ks、gss、KernSmooth以及sm包都可以实现核密度估计操作，在考虑定制化设置上，我们最终选择sm包进行空间核密度计算，具体代码如下：

point_dens<- sm.density(data.frame(scatter_df$`经度`, scatter_df$`纬度`),
                        display= "image", ngrid=400,
                        ylim=st_bbox(jiangsu)[c(2,4)],
                        xlim=st_bbox(jiangsu)[c(1,3)])

library(sp)

Density_result<-SpatialPoints(expand.grid(x=point_dens$eval.points[,1], 
                                          y=point_dens$eval.points[,2]))

Density_result<-SpatialPixelsDataFrame(Density_result, data.frame(kde_value = array(point_dens$estimate, 
                                                                  length(point_dens$estimate))))
#准换成df格式
df_density<-data.frame(Density_result)
head(df_density)

结果如下(部分)：知识点：

1. sm.density()根据散点经纬度计算核密度估计结果

point_dens<- sm.density(data.frame(scatter_df$`经度`, scatter_df$`纬度`),
                        display= "image", ngrid=400,
                        ylim=st_bbox(jiangsu)[c(2,4)],
                        xlim=st_bbox(jiangsu)[c(1,3)])

其中：2. sp包SpatialPoints、SpatialPixelsDataFrame方法将插值结果转网格再转空间df数据。

1. 使用data.frame()将结果转成data.frame()类型便于ggplot2和sf包绘制。

R-ggplot2+sf包绘制网格插值结果

接下来，我们将上方的核密度估计结果进行可视化绘制，首先，我们绘制插值的网格结果：

#自定义颜色
my_colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
Map_point_kde_nomask <- ggplot() + 
  geom_sf(data = jiangsu,fill="NA",size=.8,color="black") + 
  geom_tile(data = df_density,aes(x=x,y=y,fill=kde_value))+
  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 01: <span style='color:#D20F26'>Map kde nomask</span>",
       subtitle = "processed map charts with <span style='color:#1A73E8'>geom_title()</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 = 5)),
        plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
        plot.caption = element_markdown(face = 'bold',size = 12),
  )

可视化结果如下：

可以看到，还是出现了和我们使用Python绘制的结果一样，都是没有对感兴趣区域(地图文件)进行裁剪出来，接下里我们使用sf包进行 “裁剪” 操作。

R-sf包实现完美“裁剪”

st_intersection() 裁剪操作

sf包提供了多种方便、灵活、简单的空间数据操作函数(抽时间可以系统进行主要函数的详细讲解)，在这里，我们就使用了sf::st_intersection() 将我们感兴趣的区域进行获取，首先，我们将kde插值结果转换成sf对象类型，代码如下：

df_density_df <-  st_as_sf(df_density,coords = c("x", "y"),crs = 4326)
head(df_density_df)

结果如下(部分)：再结合之前的地图文件(jiangsu,sf包读取的geojson文件格式)，使用*st_intersection()*方法即可，如下：

mask_result <- sf::st_intersection(df_density_df,jiangsu)

geom_sf()裁剪结果可视化绘制

在获取裁剪的结果之后，我们就可以使用geom_sf()方法进行绘制了，这里注意:aes(color=kde_value)操作，具体代码如下：

Map_point_kde <- ggplot() + 
  geom_sf(data = mask_result,aes(color=kde_value)) +
  geom_sf(data = jiangsu,fill="NA",size=.5,color="gray40") +
  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 01: <span style='color:#D20F26'>Map point kde</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),
  )

最终的可视化结果如下：

总结

这一篇推文我们详细介绍了R核密度估计、空间网格数据以及裁剪之后的可视化绘制结果，我们可以看出，R在操作空间数据上较Python 还是灵活下，特别是功能较为强大的sf包，此外，R在绘制地图可视化作品时，其元素添加也较为丰富(比例尺、指北针等)。接下里，我将继续使用R和Python(两个版本)， 探索空间插值应用较为广泛的方法及对应的可视化结果，感受空间可视化带给我们的视觉盛宴！希望小伙伴们能够喜欢。

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读