在R中实现类GIS分区统计：多栅格多边形均值提取及函数疑问

嘿，我来帮你理清这两个问题：

1. zonal与extract的区别仅在于支持的几何类型吗？ #

否。这俩的核心差异不止几何类型支持：

• zonal() 是专门为面状区域设计的批量统计工具，直接输出你要的均值、总和这类预定义统计结果，不会返回单个像元的数据，处理大区域+多栅格的场景效率很高，毕竟它底层是针对区域统计优化的。
• extract() 更偏向“原始值提取”——它会把每个多边形覆盖的所有栅格像元值都给你，你得自己再用mean()这类函数算统计量。好处是灵活（比如你想同时算均值和标准差），但处理2500+多边形的大场景时，速度可能不如zonal()，而且会生成大量中间数据。

2. 提取多栅格均值并按文件名设列名的实用方案 #

我给你一套用terra包（现在替代旧raster包的主流工具，性能拉满）实现的代码，完美适配你的需求：

代码步骤拆解 #

1. 先加载工具包，读取你的多边形矢量文件
2. 批量抓取文件夹里的所有栅格
3. 循环计算每个栅格的多边形均值，用栅格文件名当列名
4. 把所有结果合并到一个数据框里，方便后续分析

完整代码示例 #

# 加载必备工具包
library(terra)
library(sf)

# 1. 读取你的多边形shapefile，替换成自己的文件路径
poly_data <- st_read("你的多边形文件路径.shp")
# 转成terra兼容的矢量格式，方便栅格计算
poly_vect <- vect(poly_data)

# 2. 目标栅格文件夹路径，替换成你的路径
raster_folder <- "你的栅格文件夹路径"
# 抓取文件夹里所有tif格式的栅格（后缀可以改成你需要的，比如.img）
raster_list <- list.files(raster_folder, pattern = "\\.tif$", full.names = TRUE)

# 3. 初始化结果表，先保留多边形原有的属性列
final_result <- as.data.frame(poly_vect)

# 4. 循环处理每个栅格
for (raster_path in raster_list) {
  # 读取当前栅格
  current_raster <- rast(raster_path)
  # 提取文件名（去掉路径和后缀）作为结果列名
  column_name <- tools::file_path_sans_ext(basename(raster_path))
  # 计算每个多边形的均值，na.rm=TRUE是忽略栅格里的NA值
  poly_mean <- zonal(current_raster, poly_vect, fun = "mean", na.rm = TRUE)[, 2]
  # 把当前栅格的均值结果添加到结果表里
  final_result[[column_name]] <- poly_mean
}

# 看看处理后的结果
head(final_result)

踩坑提示 #

• 要是运行报错，先检查栅格和多边形的坐标系是否一致，不一致的话用project(poly_vect, crs(current_raster))转投影。
• 如果用旧的raster包，逻辑差不多，但terra处理2500+多边形的速度会快很多，更推荐用这个。
• 要是你的栅格文件名有特殊字符（比如空格、中文），记得提前处理下，避免列名出问题。

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