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