R语言栅格与物种缓冲点重叠计算问题求助

先把数据预处理捋顺 #

• 确保物种点数据带物种ID字段（比如命名为species），转成sf点对象后，先检查和可可种植区栅格的投影是否一致——投影不匹配是空间分析的常见坑，用st_crs(buffer_sf) == st_crs(cocoa_rast)验证，不一致就用st_transform()转成栅格的投影。
• 生成缓冲时别用地理坐标系（WGS84，单位是度），转成UTM这类投影坐标系再操作，不然缓冲距离没有实际意义。
• 先确认可可栅格的有效值：用unique(values(cocoa_rast))查看，确认是否是1代表种植区、0/NA代表非种植区，后续计算都要基于这个规则。

处理exact_extract返回的大列表 #

exact_extract的列表每个元素对应一个缓冲多边形的提取结果，直接转成数据框就能关联原数据：

library(exactextractr)
library(dplyr)

# 假设缓冲后的sf对象是buffer_sf，可可种植区栅格是cocoa_rast
extract_list <- exact_extract(cocoa_rast, buffer_sf, include_cols = "species")

# 转成数据框，计算单个缓冲的重叠面积
overlap_df <- bind_rows(extract_list) %>%
  filter(value == 1) %>% # 只筛选种植区的栅格单元
  group_by(species, polygon_id) %>% # polygon_id是exact_extract自动生成的，对应每个缓冲
  summarise(overlap_area = sum(coverage_fraction * cell_area)) %>% # 覆盖比例×栅格面积求和
  ungroup()

coverage_fraction是每个栅格被缓冲覆盖的比例，cell_area是栅格单元的实际面积，相乘再求和就是该缓冲与可可种植区的重叠面积。

计算每个物种的总重叠面积 #

基于上面的overlap_df，直接按物种分组求和即可：

species_total_overlap <- overlap_df %>%
  group_by(species) %>%
  summarise(total_overlap = sum(overlap_area)) %>%
  ungroup()

这就是每个物种所有缓冲点与可可种植区的总重叠面积。

关于terra::extract的警告（可选替代方案） #

之前的警告大概率是投影不匹配或缓冲与栅格范围不兼容导致的，试试转成terra矢量格式再提取：

library(terra)

# sf对象转terra矢量
buffer_vect <- vect(buffer_sf)
# 提取每个缓冲内种植区栅格的数量，乘以单个栅格面积得到重叠面积
extract_result <- extract(cocoa_rast, buffer_vect, fun = sum, na.rm = TRUE)
buffer_sf$overlap_area <- extract_result[,2] * cellSize(cocoa_rast)[1]

如果栅格是地理坐标系，cellSize返回的是平方度，要转成平方公里的话，用area(buffer_vect)的单位换算，或者直接把栅格转成投影坐标系后再处理。

验证结果正确性 #

挑几个缓冲样本手动计算交集面积，和工具提取的结果对比：

# 取第一个缓冲作为样本
sample_buf <- buffer_sf[1,]
# 把可可栅格转成种植区矢量
cocoa_sf <- as.polygons(cocoa_rast, values = TRUE) %>%
  st_as_sf() %>%
  filter(value == 1)
# 手动计算交集面积
manual_overlap <- st_intersection(sample_buf, cocoa_sf) %>%
  st_area() %>%
  sum()
# 和exact_extract的结果对比
cat("手动计算面积：", manual_overlap, "\n")
cat("工具提取面积：", overlap_df$overlap_area[1], "\n")

两者误差在栅格精度范围内就说明提取结果正确。

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