multi outputs with LAScatalog processing engine

[wiesehahn]

How would I write multiple outputs when processing a LAScatalog with the processing engine?

using opt_output_files(ctg) <- "path/name" I can define one name per chunk. As it is given without extension it might be possible to output e.g. one raster and one vector data set per chunk which are written when they are returned from the function (I have not tried that).
Is there a way to write multiple files of the same driver per chunk? How would I write a DTM ({ORIGINALFILENAME}_dtm.tif) and corresponding hillshade ({ORIGINALFILENAME}_hillshade.tif) given the example from #569?

[Jean-Romain]

https://gis.stackexchange.com/questions/387853/is-it-possible-to-return-2-objects-from-a-function-when-running-catalog-apply/387863#387863

[wiesehahn]

I did not thought it would be that complicated. Tried it anyway. For LAScatalog it seems to work but not for LAS. It kind of feels like a mess.

library(lidR)

hillshade <- function(las, resolution = 1, storage_path = here::here())
{
  write_two_rast = function(two_rast, file)
  {
    file1 = file2 = tools::file_path_sans_ext(file)
    file1 = paste0(file1, "_dtm.tif")
    file2 = paste0(file2, "_hillshade.tif")
    terra::writeRaster(two_rast[[1]], file1)
    terra::writeRaster(two_rast[[2]], file2)
  }
  
  create_hillshade <- function(chunk, resolution)
  {
    las <- readLAS(chunk)
    if (is.empty(las)) return(NULL)
    
    dtm <- rasterize_terrain(las, resolution, tin(), pkg = "terra")
    
    dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
    dtm_hillshade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
    
    dtm <- terra::crop(dtm, ext(chunk))
    dtm_hillshade <- terra::crop(dtm_hillshade, ext(chunk))
    
    output <- list(dtm, dtm_hillshade)
    return(output)
  }
  
  
  if (is(las, "LAScatalog"))  {
    opt_output_files(las) <- paste0(storage_path, "/{ORIGINALFILENAME}")
    opt_chunk_buffer(las) <- 30
    opt_chunk_size(las) <- 0
    opt_select(las) <- "xyzc"
    opt_filter(las) <- "-keep_class 2"
    opt_merge(las) <- FALSE
    
    las@output_options$drivers$list = list(
      write = write_two_rast,
      object = "two_rast",
      path = "file",
      extension = "",
      param = list())
    
    
    options <- list(automerge = FALSE, need_buffer = TRUE)
    catalog_apply(las, create_hillshade, resolution=resolution, .options = options)
  }
  else if (is(las, "LAScluster")) {
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    create_hillshade(las, resolution)
  }
  else if (is(las, "LAS")) {
    dtm <- rasterize_terrain(las, resolution, tin(), pkg = "terra")
    dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
    dtm_hillshade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
    
    output <- list(dtm, dtm_hillshade)
    return(output)
  }
  else {
    stop("This type is not supported.")
  }
}

[Jean-Romain]

Rule of thumb: even if it is technically valid, if you define a function inside a function you have likely a design flaw.

Yet your code looks valid at first glance but the following is better (didn't try)

library(lidR)

write_two_rast = function(two_rast, file)
{
  file1 = file2 = tools::file_path_sans_ext(file)
  file1 = paste0(file1, "_dtm.tif")
  file2 = paste0(file2, "_hillshade.tif")
  terra::writeRaster(two_rast[[1]], file1)
  terra::writeRaster(two_rast[[2]], file2)
}

hillshade <- function(las, resolution = 1)
{
  if (is(las, "LAScatalog")) {
    las@output_options$drivers$list = list(
      write = write_two_rast,
      object = "two_rast",
      path = "file",
      extension = "",
      param = list())
    options <- list(automerge = FALSE, need_buffer = TRUE)
    opt_select(las) <- "xyzc"
    opt_filter(las) <- "-keep_class 2"
    return(catalog_apply(las, hillshade, resolution=resolution, .options = options))
  }
  else if (is(las, "LAScluster")) {
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    res = hillshade(las, resolution)
    dtm <- terra::crop(res$dtm, ext(chunk))
    dtm_hillshade <- terra::crop(res$hillshade, ext(chunk))
    output <- list(dtm = dtm, hillshade = dtm_hillshade)
    return(output)
  }
  else if (is(las, "LAS")) {
    dtm <- rasterize_terrain(las, resolution, tin(), pkg = "terra")
    dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
    dtm_hillshade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
    output <- list(dtm = dtm, hillshade = dtm_hillshade)
    return(output)
  }
  else {
    stop("This type is not supported.")
  }
}

las = readLAScatalog()
opt_output_files(las) <- paste0(storage_path, "/{ORIGINALFILENAME}")
opt_chunk_buffer(las) <- 30
opt_chunk_size(las) <- 0
opt_merge(las) <- FALSE
hillshade(las)

Is it complex? Yes certainly! But this is the price to pay to be able to do anything even stuff not originally intended by the developer (me). When I designed the LAScatalog engine I tried to design it in such a way that users can do stuff I did not even thought about. So far I never had to change the core of the engine even for edge cases like yours. But yes it is a bit complicated. But once it is done you have a function that is easy to use and infinitly reusable by others 😉