如何使用terra::tapp对连续两年的栅格数据批量应用自定义函数

我有一组逐日栅格数据，希望在连续两年的滑动时间窗口（如2001-01-01至2002-12-31、2002-01-01至2003-12-31）内应用自定义函数。目前用循环实现速度较慢，想改用terra包的tapp函数来优化效率。

核心思路 #

tapp函数的核心是通过分组索引对栅格栈中的图层进行批量聚合计算。针对滑动两年窗口，我们需要先为每个窗口生成对应的图层索引列表，再将自定义函数适配为可处理栅格子栈的形式，最后传入tapp执行计算。

步骤1：构建滑动窗口的分组索引 #

首先提取栅格栈的时间序列，然后定义所有需要的滑动两年窗口，再为每个窗口筛选出对应的图层位置：

library(terra)

# 生成示例数据（用户提供的代码）
rast_stack <- sapply(1:1826, function(i) rast(nrows=10, ncols=10, xmin=0, xmax=10, ymin=0, ymax=10))
rast_stack <- rast(rast_stack)
rast_stack[] <- rnorm(n=ncell(rast_stack))
terra::time(rast_stack) <- seq(as.Date("2001-01-01"), as.Date("2005-12-31"), by="day")

# 提取时间序列
dates <- time(rast_stack)

# 定义滑动两年窗口的起始年份（根据数据时间范围调整）
start_years <- 2001:2004

# 生成每个窗口的时间区间，并获取对应的图层索引
window_indices <- lapply(start_years, function(y) {
  start_date <- as.Date(paste0(y, "-01-01"))
  end_date <- as.Date(paste0(y+1, "-12-31"))
  which(dates >= start_date & dates <= end_date)
})

# 为每个窗口命名（方便后续识别结果）
names(window_indices) <- paste0(start_years, "-", start_years+1)

步骤2：适配自定义函数 #

tapp传入的函数，第一个参数是该组内所有图层组成的栅格子栈。根据你的Myfunction逻辑，它是对单个栅格图层做逐像元转换，因此我们可以先对每个图层应用Myfunction，再在窗口内做聚合（比如求和，你可以根据需求替换为mean、max等）：

# 原自定义函数（用户提供）
Myfunction <- function(TempRast, TDMIN, TDMAX, TCXSTOP){
  X <- TempRast
  X[is.na(TempRast)|is.nan(TempRast)] <- 0
  udevc <- TempRast
  udevc[] <- NA
  udevc[X <= TDMIN] <- 0
  udevc[TDMIN < X & X < TDMAX] <- X[TDMIN < X & X < TDMAX] - TDMIN
  udevc[TDMAX <= X & X < TCXSTOP] <- (TDMAX - TDMIN) / (TDMAX - TCXSTOP) * (X[TDMAX <= X & X < TCXSTOP] - TCXSTOP)
  udevc[X >= TCXSTOP] <- 0
  udevc[is.na(TempRast)|is.nan(TempRast)] <- NA
  return(udevc)
}

# 适配为适用于tapp的函数：先处理每个图层，再聚合
window_fun <- function(sub_stack, TDMIN, TDMAX, TCXSTOP) {
  # 对组内每个图层应用Myfunction
  processed <- app(sub_stack, function(x) Myfunction(x, TDMIN, TDMAX, TCXSTOP))
  # 对处理后的子栈做聚合（这里用求和，可替换为mean、max等）
  sum(processed, na.rm=TRUE)
}

如果你的需求是直接对窗口内的时间序列（每个像元的所有值）应用计算逻辑，可直接修改函数接受子栈，逐像元处理时间序列：

# 示例：直接对窗口内的时间序列做计算（根据你的需求调整）
window_fun_direct <- function(sub_stack, TDMIN, TDMAX, TCXSTOP) {
  app(sub_stack, function(x) {
    # x是单个像元在窗口内的所有值组成的向量
    processed <- sapply(x, function(val) {
      if (is.na(val)) return(NA)
      if (val <= TDMIN) return(0)
      if (TDMIN < val & val < TDMAX) return(val - TDMIN)
      if (TDMAX <= val & val < TCXSTOP) return((TDMAX - TDMIN)/(TDMAX - TCXSTOP)*(val - TCXSTOP))
      if (val >= TCXSTOP) return(0)
    })
    sum(processed, na.rm=TRUE)
  })
}

步骤3：用tapp执行计算 #

将分组索引和适配后的函数传入tapp，同时传入自定义函数的参数：

# 设置自定义函数的参数
TDMIN <- 0
TDMAX <- 20
TCXSTOP <- 30

# 执行计算
res <- tapp(rast_stack, index=window_indices, fun=window_fun, TDMIN=TDMIN, TDMAX=TDMAX, TCXSTOP=TCXSTOP)

# 查看结果
res

关键说明 #

• tapp的index参数支持列表格式，每个列表元素对应一个分组（滑动窗口）的图层索引，完美适配滑动窗口场景。
• 相比循环，tapp底层采用向量化和并行优化（可通过cores参数开启多线程），能显著提升计算速度。
• 如果你的自定义函数本身可以向量化处理（比如直接对向量操作而非单个值），可以进一步优化app内部的逻辑，减少循环开销。

内容的提问来源于stack exchange，提问作者Renan le roux