基于草地占比≥60%的栅格条件提取多边形平均SOC值

背景数据 #

• 30×30m分辨率的土地覆盖栅格文件
• 250×250m分辨率的土壤有机碳（SOC）栅格文件
• 包含不同大小连续多边形的矢量文件

核心需求 #

针对每个多边形，计算其中草地占比≥60%的像素的平均SOC值

现有方案的问题 #

之前的代码仅通过ifel(lulc == 3, 1, NA)或subst(landuse, 3, 1, others=NA)判断单元格是否为草地，但这种逻辑仅识别优势土地利用类型——即使单元格内草地占比极低（如仅10%），只要优势类型是草地就会被纳入计算，导致SOC值与实际高占比草地区域的数值偏差较大。

已完成栅格分辨率统一操作：

g <- terra::resample(carbon, landuse, method = 'near')

修正后的解决方案代码 #

library(terra)

# 替换为你的实际文件路径
landcover <- rast("土地覆盖栅格路径.tif")
soc <- rast("SOC栅格路径.tif")
pols <- vect("矢量多边形路径.shp")

# 统一栅格分辨率（将SOC重采样到土地覆盖的30m分辨率，后续再聚合）
soc_resampled <- resample(soc, landcover, method = 'near')

# 标记草地像素（假设草地编码为3，替换为你的实际编码）
grass_pixels <- ifel(landcover == 3, 1, 0)

# 聚合到250m分辨率，计算每个大像素内草地的占比
# fact为聚合因子，250/30≈8.33，可根据栅格对齐情况调整
grass_ratio <- aggregate(grass_pixels, fact = round(250/30), fun = mean, na.rm = TRUE)

# 筛选出草地占比≥60%的像素
grass_qualified <- ifel(grass_ratio >= 0.6, 1, NA)

# 提取符合条件的SOC值
soc_qualified <- mask(soc, grass_qualified)

# 计算每个多边形内的平均SOC值
avg_soc <- extract(soc_qualified, pols, mean, na.rm = TRUE, ID = FALSE)

# 将结果关联到矢量属性表
pols$avg_soc_grass_60 <- avg_soc

# 可选：保存结果矢量文件
writeVector(pols, "结果多边形文件.shp")

代码说明 #

• 像素标记：将草地像素标记为1，非草地为0，便于后续计算占比
• 聚合计算占比：把30m的土地覆盖数据聚合到250m分辨率，得到每个大像素内草地的比例
• 筛选达标像素：仅保留草地占比≥60%的像素，其余设为NA
• 掩码提取SOC：通过掩码过滤掉不符合条件的SOC值
• 多边形统计：对每个多边形内的达标SOC像素计算平均值

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