在R的Terra中计算1km网格内草本覆盖超10%的30m栅格细胞占比

问题描述 #

处理30m分辨率的草本覆盖百分比栅格数据，需计算每个1km×1km网格（存储为含cellid和geometry列的sf数据框）内，栅格值大于10（草本覆盖超10%）的细胞占比。现有一组与sf网格同投影的年度SpatRasters栈，尝试通过循环实现逐年计算，预期输出各年度数据框（含grid cellid及对应占比列），但循环运行极慢，且遍历完整tortgrid_1km_strata对象时触发报错。

报错信息 #

Error in data.frame(layer = i, aggregate(w, list(x), sum, na.rm = FALSE)) : arguments imply differing number of rows: 1, 0

原始代码 #

herb_10_LIST <- vector("list",length(stack))

for (i in 1:nlyr(stack)){ ## Just running for first year raster
  herb_10 <- c() ## Making an empty vector to store the proportion values with >10% herb cover for each grid cell
  for (j in 1:nrow(tortgrid_1km_strata)){ 
    cell <- terra::vect(tortgrid_1km_strata[j,]) # 1 km grid cell
    ext <- terra::extract(x = stack[[i]], y = cell, fun=table, weights=TRUE, exact=FALSE) # Extract 30m cells under 1km grid and find vlaues and weights
    vals <- colnames(ext) # For some reason the way the table comes out, the raster values  are the column names...
    ## If there are no raster cells under the 1km grid, paste NA
      if (length(vals) <=2){ ## The first two columns are not actual raster values!
        herb_10[j] <- NA } 
    ## Otherwise,
    else {
      vals <- as.numeric(vals[3:length(vals)])
      perc <- ext[1,] ## The first row of the output table is the cell weights
      perc <- perc[,c(3:ncol(perc))] %>% as.numeric()
      tab <- as.data.frame(cbind(vals,perc)) # Making into a df of values and weights for each 30 cell
      length <- as.numeric(nrow(tab)) # Number of cells under the grid
      perc_10 <- filter(tab, vals >= 10) # Filtering for cells > 10% grass cover
      ## If there is at least one cell with > 10% grass cover
        if (nrow(perc_10) >= 1) { 
          prop_10 <- as.numeric(unname(colSums(perc_10)))
          prop_10 <- prop_10[2]/ length } # Finding prop of all cells that had > 10% cover
      ## Otherwise, assign prop as 0
        else { prop_10 <- 0 }
      ## Put the proportion into the vector
      herb_10[j] <- prop_10 
        }
  } # End of inner loop
  herb_10_LIST[[i]] <- as.data.frame(cbind(cellids,herb_10)) ## Binding to cellids (object I  created that represents just the gridcell ids

} # End of outer loop

优化方案与报错修复 #

核心问题 #

嵌套循环逐个处理网格单元格的方式效率极低，且extract+table的组合容易因空网格（无栅格覆盖）触发行数不匹配的报错。改用terra::zonal函数进行批量分区统计，可大幅提升效率并自动处理空值场景。

优化代码 #

library(terra)
library(dplyr)

# 将sf网格转换为SpatVector（与栅格投影一致）
grid_vector <- vect(tortgrid_1km_strata)

# 初始化结果列表
herb_10_LIST <- vector("list", nlyr(stack))

# 逐年处理栅格
for (i in 1:nlyr(stack)) {
  # 获取当前年度栅格
  annual_raster <- stack[[i]]
  
  # 创建二值栅格：草本覆盖>10%标记为1，其余为0（NA保留）
  binary_raster <- ifel(annual_raster > 10, 1, 0)
  
  # 分区计算每个网格内的加权均值（即目标占比），weights=TRUE考虑部分覆盖的栅格单元
  zonal_stats <- zonal(binary_raster, grid_vector, fun = "mean", weights = TRUE, na.rm = TRUE)
  
  # 合并cellid与计算结果，生成年度数据框
  annual_result <- tortgrid_1km_strata %>%
    select(cellid) %>%
    mutate(herb_cover_prop = zonal_stats$mean)
  
  # 存入结果列表
  herb_10_LIST[[i]] <- annual_result
}

关键说明 #

1. 效率提升：zonal函数批量处理所有网格，避免了逐个单元格的循环，计算速度可提升数倍至数十倍。
2. 自动处理空网格：无栅格覆盖的网格会返回NA，无需手动编写空值判断逻辑，彻底解决原代码的报错问题。
3. 结果准确性：加权均值直接对应栅格值>10的面积占比，自动处理栅格与网格的部分重叠场景，无需手动计算权重求和。

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