Indexing RasterStack in for loop

Question

I am trying to calculate canopy height. I have a RasterStack which has 8 rasters in it. The first 7 are DSM; the 8th is DTM. As you see in the function, I need to subtract every DSM with the same DTM. How can I specify x[[8]] this without saying 8?? or How can I use these rasters without stacking?? I want the code to be reproducible.

x in the function is rasterstack.

####Canopy Height Calculation####
CHM <- c() 
chmCalc <- function(x) { #class(x) must be a "RasterStack"
 for (i in (1:(nlayers(x)-1))){
 res_chm <- x[[i]]-x[[8]] #DTM data needed to be set
 CHM <- c(CHM,res_chm)
 }
 names(CHM) <- paste0("CHM_H", 1:length(CHM))
 return(CHM)
} ```

Answer 1

You may pass a second argument, a string to specify the name of the layer which is the DTM; in the following example I filter out the dtm layer from the layer names, because it may happen it's not the last layer

library(raster)
rs = stack() # build a stack to use for the example
for(i in 1:8) {
 m = matrix(runif(100), nrow = 10)
 r = raster(m)
 names(r) = paste0("r", i)
 rs = stack(rs, r)
}
names(rs)[8] = "dtm" # name this layer dtm
# include the stack into the function
chmCalc <- function(x, dtm) { # dtm argument included
 CHM <- stack() 
 for (i in names(x)[names(x) != "dtm"]){
 res_chm <- x[[i]] - x[[dtm]] # refer to dtm argument
 names(res_chm) = paste0("chm_", i) 
 CHM <- stack(CHM, res_chm)
 }
 return(CHM)
}
canopy = chmCalc(rs, "dtm")
plot(canopy)

EDIT: Notes on code performance

When using for loops it is advisable to preallocate your vectors; otherwise the whole object is copied twice on each iteration. You may as well use lapply:

# PREALLOCATED VECTOR
chmCalc_v <- function(x, dtm) { 
 CHM = vector(mode = "list", length = nlayers(x)-1) 
 names(CHM) = names(rs)[names(rs) != "dtm"]
 for (i in names(x)[names(x) != "dtm"]){
 CHM[[i]] = x[[i]] - x[[dtm]] # DTM data needed to be set
 }
 return(stack(CHM))
}
# LAPPLY
chmCalc_l <- function(x, dtm) { 
 CHM = lapply(names(x)[names(x) != "dtm"], 
 function(i) {x[[i]] - x[[dtm]] })
 names(CHM) = names(rs)[names(rs) != "dtm"]
 CHM = stack(CHM)
 }
# comparing non-preallocated for-loop; preallocated for-loop and lapply:
# with a 100*100*100 raster stack
microbenchmark::microbenchmark(times = 10, chmCalc(rs,"dtm"), chmCalc_v(rs, "dtm"), chmCalc_l(rs, "dtm"))
Unit: milliseconds
 expr min lq mean median uq max neval
 chmCalc(rs, "dtm") 4652.2381 4691.9770 4857.7803 4777.1149 4872.9359 5454.4707 10
 chmCalc_v(rs, "dtm") 541.3120 554.1429 582.8172 567.7251 604.2531 654.9167 10
 chmCalc_l(rs, "dtm") 545.6215 553.4041 630.5154 601.7907 669.8223 791.9517 10

• RS_girl: This is a super nice solution but I got DTM in result stack CHM which I don't want to keep it in calculation. As I see I have to specify DTM raster outside the stack beforehand running the function. @ElioDiaz