Is it possible to return 2 objects from a function when running catalog_apply?

Question

I am normalizing a set of LAS files. In the function after normalize_height, I filter out vegetation with heights below 0 using filter_poi, and returning this as input to catalog_apply. However, I would also like to save vegetation with heights below 0 as a check.

It doesn't appear possible to return 2 objects (veg_ge_0 and veg_lt_0). It is also not possible to use opt_filter outside of the function since this seems to work on the input to the function and not on the output.

Is there a more efficient way (like returning 2 objects from the function to use in a second catalog_apply) other than creating a new function combined with catalog_apply for filtering the negatives?

library(lidR)


setwd("D:/Projects2021/p2021_0101_lidR_vs_LasTools/")
dirLASfiles_Raw <- "LasFiles__original/"
dirLASfiles_Norm <- "LasFiles_normalizedWith_lidR_wTIN/"

interpolationMethod <- tin('', extrapolate = knnidw(k = 10, p = 2, rmax = 50))


ctg <- readLAScatalog(dirLASfiles_Raw,
                      progress = TRUE,
                      filter = "-drop_withheld -drop_overlap -drop_z_below -1")


groundWaterVeg_filterAndNormalization_Function = function(chunk) {
  lasChunk = readLAS(chunk)
  if (is.empty(lasChunk)) return(NULL)
  
  ground <- filter_poi(lasChunk, Classification == 2L | Classification == 9L)
  if (is.empty(ground))
    return(NULL)
  
  vegetation <- filter_poi(lasChunk, Classification == 3L | Classification == 4L | Classification == 5L)
  if (is.empty(vegetation))
    return(NULL)

  lasNorm <- normalize_height(lasChunk, algorithm = interpolationMethod, na.rm = TRUE, use_class = c(2L, 9L))

  #lasNormFltrNegVeg <- filter_poi(lasNorm, ((Classification == 3L | Classification == 4L | Classification ==5L) & (Z < 0)), buffer == 0)
  lasNormFltr <- filter_poi(lasNorm, ((Classification == 3L | Classification == 4L | Classification ==5L) & (Z >= 0 & Z <= 30)) | (Classification == 2L | Classification == 9L), buffer == 0)  #! LASFilter is depracated. use Filter_POi
  return(lasNormFltr) #,lasNormFltrNegVeg)
}

opt_chunk_buffer(ctg) <- 25
opt_output_files(ctg) <- paste(dirLASfiles_Norm,"{ORIGINALFILENAME}_normTIN_vegge0PlusGround", sep = "")
mylasCatalogueNorm <- catalog_apply(ctg, groundWaterVeg_filterAndNormalization_Function)

Answer 1

What you want to do is definitively complex. I had to read my own documentation to achieve it. The LAScatalog processing engine allows to do almost everything but is consequently hard to use. Everything is documented but it is definitively hard to find and put all pieces together. So we will go step by step and this answer will serve as reference.

There are several way to code it. I'm using my own style following this vignette but you can do it differently.

Reproducible example

First I'm generating data to make a reproducible example.

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile)
las$Classification[las$Classification == LASUNCLASSIFIED] <- LASLOWVEGETATION
f = tempfile(fileext = ".las")
writeLAS(las, f)
ctg = readLAScatalog(f)

groundWaterVeg_filterAndNormalization_Function = function(las)
{
  UseMethod("groundWaterVeg_filterAndNormalization_Function", las)
}

Modification of your function

I modified your function so it is more readable (imo) and more importantly it works for a LAS object and returns 2 outputs in a list

groundWaterVeg_filterAndNormalization_Function.LAS = function(las) 
{
  ground <- las$Classification %in% c(LASGROUND, LASWATER)
  if (!any(ground)) return(NULL)
  
  vegetation <- las$Classification %in% c(LASLOWVEGETATION, LASMEDIUMVEGETATION, LASHIGHVEGETATION)
  if (!any(vegetation)) return(NULL)
  
  interpolationMethod <- tin('', extrapolate = knnidw(k = 10, p = 2, rmax = 50))
  lasNorm <- normalize_height(las, algorithm = interpolationMethod, na.rm = TRUE, use_class = c(2L, 9L))
  
  lasNormFltrNegVeg <- filter_poi(lasNorm, (Classification %in% c(LASLOWVEGETATION, LASMEDIUMVEGETATION, LASHIGHVEGETATION) & (Z < 0)))
  lasNormFltr <- filter_poi(lasNorm, (Classification %in% c(LASLOWVEGETATION, LASMEDIUMVEGETATION, LASHIGHVEGETATION) & (Z >= 0 & Z <= 30)) |  Classification %in% c(LASGROUND, LASWATER))
  output <- list(lasNormFltr, lasNormFltrNegVeg)
  return(output)
}

We can try it to check if it works. It returns a list with two LAS object

out = groundWaterVeg_filterAndNormalization_Function(las)

Your function for LAScluster

We have a function that works for a LAS now we need a version compatible with catalog_apply(). In your example you only created the latest but I always split in two parts to make versatile functions that can be easily tested. Here the key is that you must clear the buffer of the two objects.

groundWaterVeg_filterAndNormalization_Function.LAScluster = function(las)
{
  las = readLAS(las)
  if (is.empty(las)) return(NULL)
  
  out = groundWaterVeg_filterAndNormalization_Function(las)
  if (is.null(out)) return(NULL)
  
  out[[1]] = filter_poi(out[[1]], buffer == 0L)
  out[[2]] = filter_poi(out[[2]], buffer == 0L)
  
  return(out)
}

The example file being tiny we can test it as is without redirecting the output into disk files. The output is a list with 4 elements. Each element is a list of two LAS objects

opt_chunk_size(ctg) <- 250
out = catalog_apply(ctg, groundWaterVeg_filterAndNormalization_Function)
out

Redirection in files

This is the trickiest part. The output of each file is a list of LAS and the catalog engine does not know how to write that in files. First we create a function capable of writing two LAS objects using a single path as input

write_two_las = function(two_las, file)
{
  file1 = file2 = tools::file_path_sans_ext(file)
  file1 = paste0(file1, "_NormFltr.las")
  file2 = paste0(file2, "_NormFltrNegVeg.las")
  writeLAS(two_las[[1]], file1)
  writeLAS(two_las[[2]], file2)
}

Now we must tell the catalog how to write a list because it does not know how to do it. We will create a new driver for this LAScatalog. This is documented in ?lidR-LAScatalog-drivers and in this vignette

ctg@output_options$drivers$list = list(
  write = write_two_las,
  object = "two_las",
  path = "file",
  extension = "",
  param = list())

We have now a driver that knows how to write list of 2 LAS so, we can use opt_output_file()

opt_filter(ctg) <- "-drop_withheld -drop_overlap -drop_z_below -1"
opt_output_files(ctg) <- paste0(tempdir(),"/file_{ID}")
out = catalog_sapply(ctg, groundWaterVeg_filterAndNormalization_Function)
out
#> [1] "/tmp/RtmpgBhEaH/file_1" "/tmp/RtmpgBhEaH/file_2" "/tmp/RtmpgBhEaH/file_3" "/tmp/RtmpgBhEaH/file_4"

Notice that it returns the name of the files as defined by the engine not the one you actually used by overwriting the path generated by the engine.

Final output

ctgNormFltr = readLAScatalog(paste0(out, "_NormFltr.las"))
ctgNormFltrNegVeg = readLAScatalog(paste0(out, "_NormFltrNegVeg.las"))