使用terra包ifel函数处理NASA黑大理石NTL影像坏像素问题

我正在使用NASA的黑大理石月度夜间灯光（NTL）产品VNP46A3，数据为.h5格式，除辐射产品外还附带质量保证（QA）数据。选用All_Angle_Composite_Snow_Free辐射NTL影像进行分析，已通过以下代码提取辐射影像及QA栅格：

library(terra)

wd <- "path/"

s <- sds(paste0(wd, "VNP46A3.A2018182.h06v05.001.2021125183820.h5"))

# extract single image as spatraster from h5
r <- s[5] # s[] are the radiance and QA rasters

writeRaster(r, paste0(wd, "NTL.tif"))

接下来需要对NTL影像做预处理剔除坏质量像素，现有两个QA栅格：

• Quality为二进制栅格，优质像素值为0
• Num为生成月度栅格所用的观测次数，需大于0

期望实现：当Quality=0且Num>0时保留NTL栅格值，否则设为缺失值。但尝试编写ifel语句时出现报错，错误代码如下：

NTL <- rast(paste0(wd, "NTL.tif"))
Num <- rast(paste0(wd, "Num.tif"))
Quality <- rast(paste0(wd, "Quality.tif"))

r <- ifel(Quality = 0, NTL, NULL)
# Error in .local(test, ...) : argument "no" is missing, with no default

r1 <- ifel(Num > 0, NTL, NULL)
# Error in .local(test, ...) : argument "no" is missing, with no default

报错原因有两点：一是ifel函数不支持NULL作为替代值，栅格数据的缺失值需用NA表示；二是条件判断的语法错误（Quality = 0是赋值而非比较，应该用Quality == 0）。可以直接将两个条件合并为一个逻辑表达式，一次完成坏像素剔除：

library(terra)

wd <- "path/"

# 加载所有栅格数据
NTL <- rast(paste0(wd, "NTL.tif"))
Num <- rast(paste0(wd, "Num.tif"))
Quality <- rast(paste0(wd, "Quality.tif"))

# 组合条件：仅保留Quality为0且Num大于0的像素值，其余设为NA
cleaned_NTL <- ifel(Quality == 0 & Num > 0, NTL, NA)

# 可选：将处理后的结果写入文件
writeRaster(cleaned_NTL, paste0(wd, "Cleaned_NTL.tif"), overwrite = TRUE)

关键说明 #

• 栅格逐像素逻辑运算需用向量化运算符&，而非用于单值判断的&&
• 栅格缺失值必须用NA，NULL无法被Terra包的栅格函数识别
• 合并多条件可以避免多次运算，提升处理效率

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