1

I have to create several canopy height model (CHM) rasters from DTM & DSM rasters that are stored in different folders. The problem is that I have to repeat the process for more than 300 rasters.

Is there a way to loop the code for applying to one entire folder??

Here is my code:

#Packages needed
 library("raster")
 library("rgdal")

# load the DTM & DSM rasters
DTM <- raster("H:/ALS/DTM/33-122-117.tif")
DSM <- raster("H:/ALS/DOM/33-122-117.tif")

# CHM-overlay
CHM_ov<- overlay(DSM,
            DTM,
              fun=function(r1, r2){return(r1-r2)})

# export CHM object to new GeotIFF
writeRaster(CHM_ov, "H:/ALS/nDSMs/33-122-117.tiff",
        format="GTiff",  
        overwrite=TRUE, 
        NAflag=-9999)
1
  • R has for for loops, and list.files for getting the contents of a folder, and paste for constructing strings like the filenames for input and output. Where are you stuck?
    – Spacedman
    Oct 19 '17 at 13:15
5

300 rasters, huh? Sounds like a problem for parallel! In R, it is more common to use apply functions than your typical for loop. Using apply functions you can run processes in parallel pretty easily as well. If you are running into performance issues then this might be the way to go. Here is my solution using data.table as a raster storage container and clusterMap from the parallel library. It helps to have a lot of RAM here.

## Change to the parent directory that stores the "DTM" and "DSM" folders.
setwd("myParentDirectory")

## Packages needed (You can use non-standard evaluation here, i.e. leave out the quotes)
library(raster)
library(rgdal)
library(data.table)
library(parallel)

## get the the names of the files as separate vectors
dtmFile <- list.files("./DTM")
dsmFile <- list.files("./DSM")

## Use data.table as a storage container for your rasters.
## This is a super convenient way to run raster processes in parallel.
## Essentially, we are storing your rasters in columns of a table.
## This code block reads your rasters into columns.
rastTab <- data.table(dtmFile, dsmFile)
rastTab[,dtmRast:=lapply(file.path("./DTM", dtmFile), raster),]
rastTab[,dsmRast:=lapply(file.path("./DSM", dsmFile), raster),]

## CHM-overlay
## Here we define your function and what it will return to the table.
overFun <- function(DSM, DTM){
    require(raster)
    return(overlay(DSM, DTM, fun=function(r1, r2){return(r1-r2)}))
}

## Now we get to run your analysis in parallel using as many cores as your
## machine has minus 1 using clusterMap. You can modify this as you see fit. 
## This will take each row element for our two raster columns and apply the
## overFun function to those elements. Hopefully this makes sense.
cl <- makeCluster(detectCores() - 1)
rastTab[,CHM_ov:=clusterMap(cl, overFun, dsmRast, dtmRast),]
stopCluster(cl)

## export CHM raster object to a new GeoTIFF. This will export each new raster
## to the parent directory where you are storing your input rasters.
## You could potentially run this in parallel as well swapping mapply for clusterMap.
## see `?clusterMap` for more details.
exportFun <- function(outRast, name){
    writeRaster(outRast,
                name,
                format = "GTiff",
                overwrite = T,
                NAflag = -9999)
}

mapply(exportFun,
       rastTab$CHM_ov,
       paste0(gsub("\\..*", "", rastTab$dsmFile), gsub("\\..*", "", rastTab$dtmFile)))
2

I've create a working directory with a DOM and a DTM folder with ten test rasters in, and an empty nDSMs folder for the results.

This gets me a list of file names in DTM - note it only gives the file names and not the path, which is useful...

files = list.files("./DTM/")
files
 [1] "1-2.tif"   "10-11.tif" "2-3.tif"   "3-4.tif"   "4-5.tif"   "5-6.tif"  
 [7] "6-7.tif"   "7-8.tif"   "8-9.tif"   "9-10.tif" 

Now we want to process each pair of tifs. So loop over the file name, construct the path to each raster in the DTM and DOM folder, operate on them to get the output raster, and save it, constructing a path in the nDSMs folder:

 > for(f in files){
  DTM = raster(file.path("DTM",f))
  DSM = raster(file.path("DOM",f))
  out = foo(DTM,DSM) # whatever
  writeRaster(out, file.path("nDSMs",f))
 }

Produces ten rasters in the nDSMs folder.

Note this will fail if there's not one tif in DOM for every tif in DTM.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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