Export more than one object with catalog_apply

Question

I am looking for advice on how to export more than one type of object when using catalog_apply() from the lidR package.

For example, in the code below the custom function fn_bufferTile() computes a lidar tile extent for a buffered lidar tile. I would like to return or write both the buffered original las tile and the tile extent SpatialPolygons object. I am not sure what mechanism to use to pass back multiple objects, or how to simply write outputs directly from within a chunk?

fn_bufferTile = function(chunk)
{
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)

  ext_plyi = as(extent(chunk), "SpatialPolygons")

  # perhaps there is a way to generate custom paths here?
  #writeOGR(ext_plyi, "somePath1.shp")
  #writeLAS(las, "somePath2.las")

  # this works but write only the shapefiles
  return(ext_plyi)

  # doesn't work and fail because writing list is not supported
  #return(list(las, ext_plyi))
}

# set things up
ctg = readLAScatalog("folder/")
opt_output_files(ctg) <-  "folder/extent_{ORIGINALFILENAME}"
opt_chunk_buffer(ctg) <- 100
opt_chunk_size(ctg) <- 0
ctg@output_options$drivers$Spatial$extension <- ".gpkg"

#run function
out <- catalog_apply(ctg,  fn_bufferTile)

Not terribly important, but related, I am also wondering if the exported extent polygons can be written as separate layers to the same geopackage (e.g. the sf GeoPackage writing approach with append = T) instead of having to write to separate .gpkg files for every exported spatial object?

Answer 1

What you want to do is doable but definitively not standard. I'll try to give you a comprehensive answer so it can be used as a kind of documentation.

First you can't write from within the function. Well, technically you can. But this would invalidate the options that gives the choice to write in file or not with opt_output_file(). So don't do it.

We have a functions that returns a list

fn_bufferTile = function(chunk)
{
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  ext_plyi <- as(extent(chunk), "SpatialPolygons")
  return(list(ext_plyi, las))
}

We then define a special write function capable of writing two outputs. This function will be fed with two objects. The list you want to write and the path where to write which is the parsed template string meaning that if the template was {XCENTER} this template has been replaced by the adequate value. But here you need another supplementary parsing step to generate two filenames. We will add @@@ in the path to detect the new replacement.

mySpecialWrite = function(output_list, file)
{
  extent = output_list[[1]]
  las = output_list[[2]]
  
  path1 = gsub("@@@","extent", file)
  path2 = gsub("@@@","points", file)
  path1 = paste0(path1, ".shp")
  path2 = paste0(path2, ".las")
  
  shapefile(extent, path1, overwrite = TRUE)
  writeLAS(las, path2)
}

We then create a driver like this

mydriver = list(
  write = mySpecialWrite,
  extension = "",
  object = "output_list",
  path = "file",
  param = list())

We can now create our catalog

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile)
opt_output_files(ctg) <-  paste0(tempdir(), "/@@@_{ID}")
opt_chunk_buffer(ctg) <- 20
opt_chunk_size(ctg) <- 150

And register the driver. The class of the output is list so it is the list driver

ctg@output_options$drivers$list <- mydriver

It may be interesting to look at the summary

summary(ctg)
#> class       : LAScatalog (v1.2 format 1)
#> extent      : 684766.4, 684993.3, 5017773, 5018007 (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=utm +zone=17 +datum=NAD83 +units=m +no_defs 
#> area        : 53133.17 m²
#> points      : 81.6 thousand points
#> density     : 1.5 points/m²
#> num. files  : 1 
#> proc. opt.  : buffer: 20 | chunk: 150
#> input opt.  : select: * | filter: 
#> output opt. : on disk | w2w guaranteed | merging enabled
#> drivers     :
#>  - Raster : format = GTiff  NAflag = -999999  
#>  - LAS : no parameter
#>  - Spatial : overwrite = FALSE  
#>  - SimpleFeature : quiet = TRUE  
#>  - DataFrame : no parameter
#>  - list : no parameter

We can now use catalog_apply()

out <- catalog_sapply(ctg,  fn_bufferTile)

And you can see that you have now 2 files per chunk. One is extent_i.shp the other one is point_i.las. The only problem is that out does not contains the correct path because the engine was designed in such a way that your case is not well handled. This might be improved in the future if requested I guess.

out
#> "/tmp/RtmpVBgvAw/@@@_1" "/tmp/RtmpVBgvAw/@@@_2" "/tmp/RtmpVBgvAw/@@@_3" "/tmp/RtmpVBgvAw/@@@_4" "/tmp/RtmpVBgvAw/@@@_5" "/tmp/RtmpVBgvAw/@@@_6"

---

That being said, if what you want is really the buffered file + shapefile you already have it. LAScatalog is a shapefile (use as.spatial) and you can use catalog_retile() to buffer your las file which is much more efficient.