如何基于已知年份计算SpatRaster的dNBR等火烧后光谱恢复指标？

数据背景 #

我有研究区（AOI）内38年的拟合NBR值SpatRaster（由LandTrendR导出），以及单场火灾的火烧年份：

# class       : SpatRaster 
# dimensions  : 967, 1193, 38  (nrow, ncol, nlyr)
# resolution  : 8.983153e-05, 8.983153e-05  (x, y)
# extent      : -115.5009, -115.3937, 51.15043, 51.2373  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# source      : NBR_FTV_fit_neg_stack.tif 
# names       : 1984, 1985, 1986, 1987, 1988, 1989, ...

year_of_fire <- 2014

需求说明 #

需要编写通用函数，无需手动传入火烧前后年份，仅基于year_of_fire对SpatRaster图层索引，计算以下火烧后光谱恢复指标：

• dNBR：公式为 dNBR = NBR_pre - NBR_post（默认取火烧前1年NBR作为pre，火烧后1年作为post，可按需调整）
• R80P（80%恢复率）：公式为 R80P = max(火烧后第4/5年NBR) / (火烧前3年NBR平均值 * 0.8)
• 函数需支持批量处理数十场火灾

解决方案 #

核心思路 #

利用terra包的图层名称（年份）匹配year_of_fire，提取对应时序窗口的图层，再逐像元计算指标。无需依赖tapp，直接通过图层索引和原生函数实现更直观。

1. 辅助工具函数：提取指定年份窗口的图层 #

先编写辅助函数，快速提取火烧前后目标时序的NBR图层：

library(terra)

extract_nbr_window <- function(nbr_stack, fire_year, pre_years = -3:-1, post_years = 1:5) {
  # 获取所有图层的年份
  layer_years <- as.integer(names(nbr_stack))
  # 提取火烧前/后窗口的图层索引
  pre_idx <- which(layer_years %in% (fire_year + pre_years))
  post_idx <- which(layer_years %in% (fire_year + post_years))
  
  # 返回对应的SpatRaster子集
  list(
    pre_nbr = nbr_stack[[pre_idx]],
    post_nbr = nbr_stack[[post_idx]]
  )
}

2. 计算dNBR #

calc_dnbr <- function(nbr_stack, fire_year) {
  layer_years <- as.integer(names(nbr_stack))
  # 匹配火烧前1年和后1年的图层索引
  pre_idx <- which(layer_years == fire_year - 1)
  post_idx <- which(layer_years == fire_year + 1)
  
  # 校验年份是否存在
  if(length(pre_idx) == 0 || length(post_idx) == 0) {
    stop("火烧前后目标年份不在NBR堆栈中")
  }
  
  # 计算dNBR并命名
  dnbr <- nbr_stack[[pre_idx]] - nbr_stack[[post_idx]]
  names(dnbr) <- paste0("dNBR_", fire_year)
  return(dnbr)
}

# 调用示例
dnbr_result <- calc_dnbr(nbr_stack = your_nbr_raster, fire_year = 2014)

3. 计算R80P #

calc_r80p <- function(nbr_stack, fire_year) {
  layer_years <- as.integer(names(nbr_stack))
  # 匹配火烧前3年和后4-5年的图层索引
  pre_idx <- which(layer_years %in% (fire_year - 3):(fire_year - 1))
  post_idx <- which(layer_years %in% (fire_year + 4):(fire_year + 5))
  
  # 校验年份窗口完整性
  if(length(pre_idx) < 3 || length(post_idx) < 2) {
    stop("火烧前后目标年份窗口不完整")
  }
  
  # 分步计算指标
  pre_mean <- mean(nbr_stack[[pre_idx]])
  post_max <- max(nbr_stack[[post_idx]])
  r80p <- post_max / (pre_mean * 0.8)
  
  names(r80p) <- paste0("R80P_", fire_year)
  return(r80p)
}

# 调用示例
r80p_result <- calc_r80p(nbr_stack = your_nbr_raster, fire_year = 2014)

4. 批量处理多场火灾 #

遍历火灾年份列表，批量计算并合并结果：

# 假设有多个火灾年份
fire_years <- c(2010, 2014, 2018)

# 批量计算dNBR
batch_dnbr <- lapply(fire_years, function(y) calc_dnbr(your_nbr_raster, y))
batch_dnbr <- rast(batch_dnbr) # 合并为多图层SpatRaster

# 批量计算R80P同理
batch_r80p <- lapply(fire_years, function(y) calc_r80p(your_nbr_raster, y))
batch_r80p <- rast(batch_r80p)

替代方案（针对像元级火烧年份栅格） #

如果year_of_fire是每个像元对应不同火烧年份的SpatRaster，可使用ifel逐像元计算：

calc_dnbr_per_pixel <- function(nbr_stack, fire_year_raster) {
  layer_years <- as.integer(names(nbr_stack))
  dnbr <- rast(fire_year_raster) # 创建空结果栅格
  
  # 遍历所有唯一火烧年份
  unique_fire_years <- unique(values(fire_year_raster))
  unique_fire_years <- unique_fire_years[!is.na(unique_fire_years)]
  
  for(y in unique_fire_years) {
    pre_idx <- which(layer_years == y - 1)
    post_idx <- which(layer_years == y + 1)
    if(length(pre_idx) == 0 || length(post_idx) == 0) next
    
    # 对该火烧年份的像元赋值dNBR
    dnbr <- ifel(fire_year_raster == y, nbr_stack[[pre_idx]] - nbr_stack[[post_idx]], dnbr)
  }
  
  names(dnbr) <- "dNBR_per_pixel"
  return(dnbr)
}

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