基于起止日栅格创建仿stackApply的栅格栈像元区间求和函数

原来的逐像元循环虽然能实现需求，但面对大栅格时效率实在拉胯。我们可以写一个专门的函数，利用Raster包的特性和向量化操作来提速，同时保持用法简洁。

函数实现（基础版） #

这个版本兼容性好，适合大多数场景：

stackIntervalSum <- function(raster_stack, start_raster, end_raster) {
  # 先做输入合法性检查
  if (!inherits(raster_stack, c("RasterStack", "RasterBrick"))) {
    stop("raster_stack必须是RasterStack或RasterBrick类型哦")
  }
  if (!inherits(start_raster, "RasterLayer") || !inherits(end_raster, "RasterLayer")) {
    stop("start_raster和end_raster都得是RasterLayer对象哈")
  }
  
  # 提取关键信息
  n_cells <- ncell(start_raster)
  n_layers <- nlayers(raster_stack)
  start_days <- getValues(start_raster)
  end_days <- getValues(end_raster)
  
  # 筛选出起始/结束日期都非NA的有效像元
  valid_cells <- which(!is.na(start_days) & !is.na(end_days))
  
  # 初始化结果向量
  result_vals <- rep(NA, n_cells)
  
  # 对每个有效像元计算区间和
  result_vals[valid_cells] <- sapply(valid_cells, function(cell) {
    start <- max(start_days[cell], 1)  # 防止起始日小于1
    end <- min(end_days[cell], n_layers)  # 防止结束日超过总层数
    if (start > end) return(NA)  # 区间无效时返回NA
    sum(getValues(raster_stack[[start:end]], cell))
  })
  
  # 把结果转成栅格返回
  result_raster <- start_raster
  values(result_raster) <- result_vals
  return(result_raster)
}

快速版（适合大栅格） #

如果你的栅格数据量很大，上面的sapply还是有点慢，可以试试这个基于矩阵的版本，内存访问效率更高：

stackIntervalSum_fast <- function(raster_stack, start_raster, end_raster) {
  # 输入检查
  if (!inherits(raster_stack, c("RasterStack", "RasterBrick"))) {
    stop("raster_stack必须是RasterStack或RasterBrick类型哦")
  }
  if (!inherits(start_raster, "RasterLayer") || !inherits(end_raster, "RasterLayer")) {
    stop("start_raster和end_raster都得是RasterLayer对象哈")
  }
  
  # 把栅格栈转成矩阵（行=像元，列=日期层）
  stack_matrix <- as.matrix(raster_stack)
  start_days <- getValues(start_raster)
  end_days <- getValues(end_raster)
  
  valid_cells <- which(!is.na(start_days) & !is.na(end_days))
  result_vals <- rep(NA, nrow(stack_matrix))
  
  # 循环计算（这里的循环是操作内存矩阵，比操作栅格对象快很多）
  for (cell in valid_cells) {
    start <- max(start_days[cell], 1)
    end <- min(end_days[cell], ncol(stack_matrix))
    if (start <= end) {
      result_vals[cell] <- sum(stack_matrix[cell, start:end])
    }
  }
  
  result_raster <- start_raster
  values(result_raster) <- result_vals
  return(result_raster)
}

用法示例 #

调用起来超级简单，直接传入你的三个栅格对象就行：

# 读取输入数据（和你原来的代码一致）
s <- stack(list.files("~/dailyraster", full.names=TRUE))
r_start <- raster("~/stackSumSTART.asc")
r_end <- raster("~/stackSumEND.asc")

# 使用自定义函数计算区间和
x <- stackIntervalSum(s, r_start, r_end)
# 或者用快速版
x <- stackIntervalSum_fast(s, r_start, r_end)

小提示 #

• 两个函数都做了边界检查，比如起始日小于1会自动设为1，结束日超过总层数会设为最后一层，避免索引报错
• 如果你的栅格有投影或者其他属性，结果栅格会自动继承start_raster的所有属性，不用额外处理

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