当前位置：首页 > news >正文

R绘图|3分钟复现瑞士“苏黎世大学”Nature全球地图——基于R包ggplot2+sf等

news 来源：原创 2025/5/5 6:56:39

一、引言

本文我们复现苏黎世大学团队François Keck等在Nature最新文章“The global human impact on biodiversity”中的全球地图。

之前的图纸是在平面坐标系里面进行绘制，本次我们在罗宾逊投影中进行绘制。整体代码逻辑非常简单，就是采样点坐标系的转换，地图坐标系转换图绘制(采用的sf包），ggplot绘制点图。

1.1 要复现的图纸

图源：Keck, F., Peller, T., Alther, R. et al. Nature (2025).

1.2 本代码绘图结果

可以说一模一样，这个图纸基本完美复现原文图片！！

二、示例数据和R代码

2.1 准备数据

world_shp 文件夹，全球地图数据信息集，用于绘制地图

data.json 数据信息总结，用于绘图的有经纬度列Longitude和Latitude，和alpha_diversity(地图上的点）

⚠️根据自己的数据进行整理后绘图，有问题可以留言获取帮助

2.2 R代码

⚠️该代码适用于上述数据，要完成自己的任务适当调整！ 创建一个地图，展示不同采样点的数据分布，特别是展示是否包含 alpha 多样性数据，并结合全球国界数据进行可视化。

2.2.1数据处理

读取世界地图数据:读取 Shapefile (TM_WORLD_BORDERS_SIMPL-0.3.shp) 以获取全球国界数据。

处理缺失值：删除包含缺失经纬度 (Longitude 和 Latitude) 的记录，保留有效数据。

空间数据转换：将数据从常规数据框转换为空间数据对象 (sf)，设置坐标系统为 WGS84 (EPSG:4326)。将坐标数据转换为罗宾逊投影坐标系统，以便于地图可视化。

数据随机化和标记：创建新列 has_alpha，根据是否包含 alpha 多样性数据标记为 "Community composition only" 或 "Local diversity and community composition"。随机打乱数据顺序，以便更好地展示。

关闭 S2 球面几何处理：禁用 S2 几何处理，以避免某些空间计算错误。

空间数据处理:

1.将多边形转换为单个多边形要素。

2.计算每个多边形的面积，并过滤掉面积过小的区域（如小岛，面积 < 10⁹ 平方米）。

3.通过 st_buffer(0) 修复可能的几何错误。

4.裁剪到标准经纬度范围（经度 -180° 至 180°，纬度 -90° 至 90°）。

5.转换为罗宾逊投影，优化全球地图展示效果。

2.2.2. 绘制地图

基础地图绘制 (map)

使用 ggplot2 绘制世界地图底图 (geom_sf)，填充颜色为灰色 (grey80)。

添加采样点 (map + geom_sf)

读取采样数据 (dat_rob)，并按 has_alpha 变量进行着色：

蓝色 (#00364eff)：仅包含群落组成数据。

青色 (#00b9ffff)：包含局部多样性数据。

随机分组采样点，避免绘图时数据点过度重叠。

设定点的形状（shape = 20，实心圆点）和大小（size = 0.5）。

优化地图样式 (theme_minimal())

移除标题、坐标轴标签和边距，保持地图简洁清晰。

调整图例位置至底部，并优化图例的文本大小和间距。

2.2.3 最终输出

保存高分辨率图片 (ggsave)

以 800 DPI 高分辨率将地图保存为 map.png，尺寸设定为 6×3 英寸，背景颜色为白色 (bg = "white")。

2.3 R包加载并加载数据

#清理所有对象
rm(list=ls())# 加载必要的库
library(dplyr)       # 数据操作和转换
library(ggplot2)     # 数据可视化
library(sf)          # 空间数据处理
library(forcats)     # 因子变量处理
library(rphylopic)   # 生物图标处理
library(scales)      # 图形比例尺调整
library(patchwork)   # 多图拼接# 读取JSON数据并转换为tibble格式
dat <- jsonlite::read_json("data.json", simplifyVector = TRUE) %>% as_tibble()

2.4 数据整理

## 处理数据
dat_rob <- dat %>% # 移除经纬度缺失的记录filter(!is.na(Longitude), !is.na(Latitude))dat_rob1 <- dat_rob %>% # 转换为sf对象(空间数据格式)，指定经纬度列和坐标系(WGS84)sf::st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>% # 转换到罗宾逊投影坐标系sf::st_transform(crs = sf::st_crs("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs")) dat_rob <- dat_rob1%>% # 创建新列标记是否包含alpha多样性数据mutate(has_alpha = ifelse(is.na(alpha_diversity), "Community composition only", "Local diversity and community composition")) %>% # 随机打乱数据顺序(prop=1表示全部数据)slice_sample(prop = 1)# 关闭S2球面几何处理(避免某些空间计算错误)
sf::sf_use_s2(FALSE)

2.5 绘图

## 创建地图
map <- sf::read_sf("World_shp/TM_WORLD_BORDERS_SIMPL-0.3.shp") %>% # 将多边形转换为单个多边形要素sf::st_cast("POLYGON") %>% # 计算每个多边形的面积mutate(area = sf::st_area(geometry)) %>% # 过滤掉面积过小的多边形(小于10^9平方米)filter(as.numeric(area) > 10^9) %>% # 修复可能的几何错误sf::st_buffer(0) %>% # 裁剪到标准经纬度范围sf::st_crop(sf::st_bbox(c(xmin = -180, xmax = 180, ymin = -90, ymax = 90))) %>% # 同样转换为罗宾逊投影sf::st_transform(crs = sf::st_crs("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs")) map <- map %>%# 开始ggplot绘图ggplot() +# 绘制世界地图底图geom_sf(fill = "grey80", colour = "grey80") +# 添加采样点数据，按has_alpha分组着色geom_sf(aes(color = has_alpha, group = cut(sample(seq(1, nrow(dat_rob))), 100)),  # 分组用于避免过度绘图shape = 20,       # 点形状(实心圆点)size = 0.5,       # 点大小data = dat_rob    # 使用的数据)    +# 设置颜色方案scale_color_manual(values = c("#00364eff", "#00b9ffff")) +# 设置x轴刻度(虽然地图投影后可能不直接使用)scale_x_continuous(breaks = seq(-180, 180, 45))+# 使用简约主题theme_minimal() +theme(plot.title = element_blank(),          # 无标题plot.margin = margin(t = 0, r = 0, b = 0, l = 0),  # 无边距legend.position = "bottom",            # 图例在底部legend.title = element_blank(),        # 无图例标题legend.text = element_text(size = 7),   # 图例文字大小legend.margin = margin(t = 0, r = 0, b = 0, l = 0), # 图例无边距legend.box.margin = margin(t = 0, r = 0, b = 0, l = 0), # 图例框无边距axis.text = element_blank(),           # 无坐标轴文字axis.title = element_blank()           # 无坐标轴标题)# 显示地图
mapggsave("map.png",plot=map,dpi=800,height =3,width=6,bg = "white")