在 r 中添加多波段栅格

Question

我正在使用 Hantson & Chuvieco (2010) 提出的方法对 landsat 数据进行一些地形校正，该方法涉及一种基于特定 NDVI 值的覆盖分离方法，以区分植被和非植被区域。

所以，我从一个多波段栅格对象变为两个，分别处理，然后组合形成原始的地形校正场景。

问题是我在场景中有一些 NA 值，所以如果我不摆脱它们，我会得到一个空白图像。我使用sum 函数和定义为TRUE 的na.rm 选项来解决这个问题，但结果是单波段光栅对象而不是堆栈或光栅砖对象。

这是我的代码：

#topographic correction
VEGB1TC<- topocorr(VEGB1,slopeG,aspectG,sunelev,sunazimuth,method="ccorrection",na.value=NA)
VEGB2TC<- topocorr(VEGB2,slopeG,aspectG,sunelev,sunazimuth,method="ccorrection",na.value=NA)
....        
SDB1TC<- topocorr(SDB1,slopeG,aspectG,sunelev,sunazimuth,method="ccorrection",na.value=NA)
SDB2TC<- topocorr(SDB2,slopeG,aspectG,sunelev,sunazimuth,method="ccorrection",na.value=NA)
...

#rasterize corrected band and stack them 
VB1<-raster(VEGB1TC)
VB2<-raster(VEGB2TC)
....

ADEVEG<-stack(VB1, VB2,...)

SB1<-raster(SDB1TC)
SB2<-raster(SDB2TC)
....

ADESD<-stack(SB1, SB2,...)

#-----------combine images-------

ADE<-sum(stack(ADEVEG,ADESD),na.rm=T)

是否有可以添加带有na.rm 标签的多波段栅格的功能？我找到了do.call 函数，但在我看来，它将几个栅格相加成一个单波段对象。

感谢您的帮助。

ps：Hantson, S. Chuvieco, E. 2010。陆地卫星图像不同地形校正方法的评估。 International Journal of Applied Earth Observations and Geoinformation 13(2011):691-700p.

Answer 1

这是一个可重现的例子：

library(raster)
s1 <- stack(system.file("external/rlogo.grd", package="raster")) 
s2 <- flip(s1, 'x')
s2[1000:4000] <- NA  

现在你可以这样做了：

r1 <- s1 + s2

但这不考虑NA 值。那么您是否需要求助于：

r2 <- overlay(s1, s2, fun=function(x) sum(x, na.rm=TRUE))

或者这个：

s2r <- reclassify(s2, cbind(NA, 0))
r3 <- s1 + s2r

PS。我理解预期sum(s1, s2) 会返回一个多层对象，但我认为返回一层与base::sum 是一致的：

sum(1:4, 1:4)
# [1] 20