# Indexing RasterStack in for loop

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)
} ```
``````

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

``````
• 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 Mar 5 at 23:22