2

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

1 Answer 1

2

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

1
  • 1
    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
    – RS_girl
    Mar 5, 2021 at 23:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.