如何在R中基于多要素土地覆盖shp按TYPE批量统计NDVI影像均值

问题根源 #

你用terra::extract得到的是每个矢量要素（对应CODE1字段）的NDVI均值，而不是按TYPE土地覆盖类型分组的结果——因为extract默认对矢量的每一行（单个要素）计算统计值，不会自动按指定字段聚合。

单张影像的正确处理方法 #

有两种思路可以实现按TYPE统计：

方法1：先聚合矢量再提取（推荐，效率更高）

先把相同TYPE的矢量要素合并成单个面，再对合并后的面提取统计值，直接得到分组结果：

library(terra)

# 读取土地覆盖矢量
u <- vect("D:/muusVutm/type_49N.shp")
# 按TYPE字段聚合矢量，合并同类型要素
u_type_agg <- aggregate(u, by = "TYPE")

# 读取单张NDVI影像
ndvi_rast <- rast("D:/ndvi2001.tif")
# 提取聚合后每个TYPE的平均NDVI
type_stats <- extract(ndvi_rast, u_type_agg, fun = mean, na.rm = TRUE)
# 合并TYPE字段与统计结果
final_stats <- cbind(TYPE = u_type_agg$TYPE, mean_NDVI = type_stats[[1]])

方法2：先提取要素均值再分组汇总

如果不想修改矢量，可以先提取所有要素的NDVI均值，再用dplyr按TYPE分组计算：

library(terra)
library(dplyr)

u <- vect("D:/muusVutm/type_49N.shp")
ndvi_rast <- rast("D:/ndvi2001.tif")

# 提取每个矢量要素的NDVI均值
element_stats <- extract(ndvi_rast, u, fun = mean, na.rm = TRUE)
# 合并TYPE字段，按TYPE分组计算均值
type_stats <- cbind(TYPE = u$TYPE, element_stats) %>%
  as.data.frame() %>%
  group_by(TYPE) %>%
  summarise(mean_NDVI = mean(ndvi2001, na.rm = TRUE))

数千张影像的批量处理方案 #

通过循环或lapply遍历所有NDVI文件，批量生成统计结果并合并：

library(terra)
library(dplyr)
library(tools)

# 1. 准备工作：读取并聚合矢量
u <- vect("D:/muusVutm/type_49N.shp")
u_type_agg <- aggregate(u, by = "TYPE")

# 2. 获取所有NDVI影像路径（替换为你的NDVI文件夹路径）
ndvi_file_paths <- list.files(path = "D:/ndvi_all", 
                              pattern = "\\.tif$", 
                              full.names = TRUE)

# 3. 定义单张影像处理函数
process_single_ndvi <- function(file_path) {
  # 读取影像
  ndvi_rast <- rast(file_path)
  # 获取影像文件名（比如年份，作为结果列名）
  img_label <- file_path_sans_ext(basename(file_path))
  # 提取TYPE对应的均值
  stats <- extract(ndvi_rast, u_type_agg, fun = mean, na.rm = TRUE)
  # 返回整理后的结果：TYPE + 当前影像的均值列
  return(data.frame(TYPE = u_type_agg$TYPE, 
                    !!img_label := stats[[1]]))
}

# 4. 批量处理所有影像
all_stats_list <- lapply(ndvi_file_paths, process_single_ndvi)

# 5. 合并所有结果（按TYPE列匹配）
final_all_stats <- Reduce(function(x, y) merge(x, y, by = "TYPE"), all_stats_list)

# 6. 保存结果到CSV文件
write.csv(final_all_stats, "D:/ndvi_zonal_stats_all.csv", row.names = FALSE)

注意事项 #

• 如果影像文件过大，可分批次处理（比如每100张一组），避免内存溢出；
• 确保TYPE字段在矢量中无缺失值，否则聚合或分组会出错；
• 如果需要统计其他指标（如中位数、最大值），只需修改extract中的fun参数即可。

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