如何先栅格化多面要素计算1km网格物种丰富度（R语言实现）

依赖包安装 #

使用sf处理矢量数据、terra处理栅格数据，运行效率优于传统的raster包：

# 安装包（仅首次运行需要）
install.packages(c("sf", "terra"), dependencies = TRUE)
# 加载包
library(sf)
library(terra)

步骤1：读取物种分布矢量数据 #

# 读取裁剪后的IUCN物种分布shp文件，替换为你的本地文件路径
sp_dist <- st_read("your_file_path/species_dist.shp", stringsAsFactors = FALSE)
# 可选：查看属性表结构，确认物种名字段（通常为scientific_name）、类群字段（通常为class）
head(sp_dist)

步骤2：坐标系统一与目标栅格创建 #

注意：必须使用米为单位的投影坐标系（如研究区对应的UTM投影），不能使用WGS84等地理坐标系（单位为度），否则1km分辨率设置无效

# 替换为你研究区对应的UTM投影EPSG编码，如东亚北纬区域可使用EPSG:32650
target_crs <- "EPSG:32650"
# 矢量数据转投影
sp_dist_proj <- st_transform(sp_dist, crs = target_crs)
# 获取研究区范围
study_ext <- ext(sp_dist_proj)
# 创建1km分辨率的空白栅格
r_target <- rast(study_ext, resolution = 1000, crs = crs(sp_dist_proj))

步骤3：计算物种丰富度 #

场景1：计算所有类群总丰富度 #

# 获取所有物种的唯一学名列表
sp_unique <- unique(sp_dist_proj$scientific_name)
# 初始化栅格栈存储单个物种的分布栅格
sp_rast_stack <- rast()
# 循环处理每个物种
for (sp_name in sp_unique) {
  # 提取单个物种的分布范围
  single_sp <- sp_dist_proj[sp_dist_proj$scientific_name == sp_name, ]
  # 栅格化：物种存在的像元赋值为1，不存在为0
  single_sp_rast <- rasterize(single_sp, r_target, field = 1, background = 0)
  # 加入栅格栈
  add(sp_rast_stack) <- single_sp_rast
}
# 栅格叠加求和，得到每个像元的物种总数
total_richness <- sum(sp_rast_stack)

场景2：计算单一类群（如爬行动物）丰富度 #

仅需在上述流程前增加类群筛选步骤即可：

# 筛选爬行动物数据，class为属性表中存储类群信息的字段名，根据实际情况调整
reptile_dist <- sp_dist_proj[sp_dist_proj$class == "reptiles", ]
# 后续步骤和全类群计算一致，将sp_dist_proj替换为reptile_dist即可

步骤4：结果导出与可视化 #

# 导出丰富度栅格为tif格式
writeRaster(total_richness, "1km_species_richness.tif", overwrite = TRUE)
# 可视化结果
plot(total_richness, col = hcl.colors(100, "viridis"), main = "1km网格物种丰富度")

优化提示 #

• 若物种数量较多运行慢，可提前合并同一物种的多个分布斑块：

  sp_dist_proj <- aggregate(sp_dist_proj, by = "scientific_name", FUN = head)
  

• 若内存不足，可调整terra包的内存限制参数，或分块处理栅格。

内容的提问来源于stack exchange，提问作者user_3264