如何基于四维栅格（多日期Sentinel-2数据）计算NDVI？

针对你提供的SpatRasterDataset对象（每个元素对应一个日期的Sentinel-2影像，包含B04红波段和B08近红外波段），可以通过以下方式批量计算每景影像的NDVI：

基础计算方案 #

直接遍历SpatRasterDataset中的每景影像，提取对应波段代入NDVI公式计算：

library(terra)

# 遍历每个日期的影像，计算NDVI
ndvi_sds <- sapply(s, function(img) {
  # 提取近红外(B08)和红波段(B04)
  nir <- img[["B08"]]
  red <- img[["B04"]]
  # 计算NDVI
  (nir - red) / (nir + red)
})

# 为NDVI数据集设置日期名称
names(ndvi_sds) <- names(s)

优化：处理除数为0的情况 #

如果影像中存在B08+B04=0的像素（比如完全黑的无效区域），直接计算会出现除以0的错误，可通过ifelse将这类像素的NDVI设为NA：

ndvi_sds <- sapply(s, function(img) {
  nir <- img[["B08"]]
  red <- img[["B04"]]
  denominator <- nir + red
  # 除数为0时返回NA，否则计算NDVI
  ifelse(denominator == 0, NA, (nir - red) / denominator)
})

names(ndvi_sds) <- names(s)

结果说明 #

• 最终的ndvi_sds是一个新的SpatRasterDataset对象，每个元素对应原日期的NDVI单波段影像
• 可通过plot(ndvi_sds)快速查看各日期的NDVI结果

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