R-ggplot2 标准中国地图制作

ggplot2 + sf 中国标准地图绘制

cowplot 包实现地图子图插入

ggplot2 + sf 中国标准地图绘制

本期的绘图我们是按照如下图表进行绘制的(该图来自朋友公号)：

其实我很早就想试着用R进行此类图表的绘制了，但由于强大的Arcgis和绘图进度等原因，导致迟迟没有进行，这期就进行此类图表的绘制。注：图中散点的经纬度信息已从朋友那获知，但涉及具体数据则是虚构的。

地图数据、代码下载，后台回复“sf”获得链接

- *.geojson，地图数据文件

- *.R/Rmd R脚本文件

地图读取和散点标注

这里的中国地图格式为 geojson 格式，还是使用sf包进行灵活读取，具体代码如下：

china_shp <- "中国省级地图GS（2019）1719号.geojson"
nine <- "九段线GS（2019）1719号.geojson"
china <- sf::read_sf(china_shp)
nine_line <- sf::read_sf(nine)

数据预览如下(部分)：

接下来我们进行探索性绘图，主要使用geom_sf() 进行绘制，代码如下：

library(ggspatial)
ggplot() + 
  geom_sf(data = china,fill="NA",size=1,color="black") + 
  geom_sf(data = nine_line) + 
  coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+
  annotation_scale(location = "bl") +
      # spatial-aware automagic north arrow
       annotation_north_arrow(location = "tl", which_north = "false",
                             style = north_arrow_fancy_orienteering)

这里可以看到，我们没有使用常规的WGS84坐标，而是采用 +proj=laea +lat_0=40 +lon_0=104 进行绘制，可视化结果如下：

比例尺和指北针的添加，则是依靠ggspatial包。

散点数据的读取

由于只获得数据经纬度信息，其他数据均为虚构，数据预览如下：

这里注意点，要想使数据绘制在地图上，我们需使用st_as_sf()对其进行转换,处理：

scatter_df_tro <- st_as_sf(scatter_df,coords = c("lon", "lat"),crs = 4326)

结果如下：注意红色框中的内容

接下来将数据在地图中显示，具体代码如下：

ggplot() + 
  geom_sf(data = china,fill="NA",size=.5,color="black") + 
  geom_sf(data = nine_line) + 
  geom_sf(data = scatter_df_tro,aes(fill=class,size=data),shape=21,colour='black',stroke=.25)+
  coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+
  scale_fill_manual(values = c("#E21C21","#3A7CB5","#51AE4F"))+
  scale_size(range = c(1,5))+
  annotation_scale(location = "bl") +
      # spatial-aware automagic north arrow
       annotation_north_arrow(location = "tl", which_north = "false",
                             style = north_arrow_fancy_orienteering)+
  theme_linedraw()+
  theme(text = element_text(family = 'Times_New_Roman',size = 12,face = "bold"),
        panel.background = element_rect(fill = NA),
        panel.grid.major = element_line(colour = "grey50"),
        #axis.line = element_line(size = 1, colour = "black"),
        panel.ontop = FALSE
        )

可视化结果如下：

这里还存在如下问题：

(1) 尽量将南海部分放置在右下角(虽然这样也可以)

(2) 图例只是默认的，没有进行定制化操作，影响美观

cowplot包子图添加

cowplot包可是实现对绘图结果的任意组合，很适合绘制地图子图的绘制，这里直接给出绘图代码，如下：

library(cowplot)

map <- ggplot() + 
  geom_sf(data = china,fill=NA) + 
  geom_sf(data = nine_line,color='gray50',size=.8)+
  geom_sf(data = scatter_df_tro,aes(fill=class,size=data),shape=21)+
  coord_sf(ylim = c(-2387082,1654989),crs="+proj=laea +lat_0=40 +lon_0=104")+
  scale_fill_manual(values = c("#E21C21","#3A7CB5","#51AE4F"))+
  scale_size(range = c(1,5))+
  annotation_scale(location = "bl",text_face = "bold",
                   text_family = "Times_New_Roman") +
      # spatial-aware automagic north arrow
       annotation_north_arrow(location = "tl", which_north = "false",
                             style = north_arrow_fancy_orienteering,
                            )+
  #定制化图例：这一步可以设计出亮眼的图例哦
  guides(fill = guide_legend(override.aes = list(size = 3),
                             title = "",
                             label.position = "right",
                             ncol=3,
                             ),
         size = guide_legend(
                             title = "",
                             label.position = "right",
                             ncol=5)) +
 labs(
          caption = 'Visualization by DataCharm')+
 #theme_bw()+
 theme(
     text = element_text(family = 'Times_New_Roman',size = 18,face = "bold"),
     
     panel.background = element_rect(fill = NA),
     panel.grid.major = element_line(colour = "grey80",size=.2),
     legend.key = element_rect(fill = "white"),
     legend.position = "bottom",
 )

nine_map <- ggplot() +
  geom_sf(data = china,fill='NA') + 
  geom_sf(data = nine_line,color='gray70',size=1.)+
  #geom_sf(data = scatter_df_tro,aes(fill=class,size=data),shape=21,colour='black',stroke=.25)+
  coord_sf(ylim = c(-4028017,-1877844),xlim = c(117131.4,2115095),crs="+proj=laea +lat_0=40 +lon_0=104")+
  theme(
    #aspect.ratio = 1.25, #调节长宽比
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_blank(),
    panel.border = element_rect(fill=NA,color="grey10",linetype=1,size=1.),
    plot.margin=unit(c(0,0,0,0),"mm"))


gg_inset_map = ggdraw() +
  draw_plot(map) +
  draw_plot(nine_map, x = 0.8, y = 0.15, width = 0.1, height = 0.3)

可视化效果如下：

这里需注意的是xlim和ylim的设置，具体值都是在 +proj=laea +lat_0=40 +lon_0=104 投影坐标系下的值，而常规坐标转不同坐标下的操作也很简单，可使用st_transform进行操作。如下：

df2 <- tibble(lon = c(105,125),
       lat = c(3,25),
       ) 

df2 %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) -> df2_sf

df2_sf_pre <- sf::st_transform(df2_sf,crs="+proj=laea +lat_0=40 +lon_0=104")
df2_sf_pre

即可获取对应值，进而进行范围设置。

猜你喜欢

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

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

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

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

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

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

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

编程模板: Shell  R Perl

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

写在后面

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

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

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