lidR: Writing output from pixel_metrics() from "multiple-output-function" to separate raster files

Question

First some background from the original question:

I'm using the lidR package to make rasters of metrics calculated on point clouds. I have a custom function that calculates 9 metrics outputed as a named list and want to produce one raster per metric.

Outputed rasters should keep the filename of the LAS file but get the name of the metric attached (pasted) to the name (e.g. LASfilename_metric.tif) by using the bylayer and suffix arguments of the writeRaster function in the raster package.

I also want to use a slightly customized driver function when writing the raster-files to make sure output rasters had the right CRS. These "adjustments" are introduced by modifying the catalog output_options (write and param) as per inspiration from https://github.com/r-lidar/lidR/wiki/Modify-the-LAScatalog-drivers.

This have all worked fine using the "old" grid_metrics() function in the lidR package but has stopped working when updating the lidR package (where pixel_metrics have replaced grid_metrics).

An example that reproduce the problem:

library("lidR")
library("raster")
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")

myMetrics <- function(z){
  list(max=max(z), mean=mean(z))
}

Raster_MS <- function(x, filename, crs, ...){
  suppressWarnings(crs(x) <- paste0('EPSG:', crs))
  suppressWarnings(raster::writeRaster(x, filename, ...))
}

ctg<-readLAScatalog(LASfile)
opt_chunk_buffer(ctg) <- 0
ctg@output_options$drivers$Raster$write <- Raster_MS
ctg@output_options$drivers$Raster$param$bylayer <- T
ctg@output_options$drivers$Raster$param$suffix <- quote(names(x))
ctg@output_options$drivers$Raster$param$crs <- 26917

opt_output_files(ctg) <- paste0("path/to/output/folder/", '{ORIGINALFILENAME}')

metricsrast1 <- grid_metrics(ctg, ~myMetrics(Z), res = 10)

What happens is that I get the following error:

Processing [===========================================================================================================] 100% (1/1) eta:  0s
Error: formal argument "sources" matched by multiple actual arguments

I have seen this related question where the solution was to write multilayer rasters but I would really prefer to have several (2 in my reproducible example) single layer rasters produced (and named accordingly, i.e. Megaplot_max.tif and Megaplot_mean.tif).

Answer 1

I'm suprised your code used to work in first place. It should have not. Let me explain: in lidR, for each chunk a file is written. The path is retained and, for rasters, all the path to the files are used to build a virtual raster mosaic. The problem with your code is that the file that is supposed to be written (e.g. /tmp/Megaplot.tif) is NOT actually written and instead what is written is /tmp/Megaplot_mean.tif and /tmp/Megaplot_max.tif. But internally the engine does not know it and was not designed like that. Internally it records that the written file is /tmp/Megaplot.tif. Thus at the end it should fail when building the virtual raster because the files path recorded by the engine is not the one written on disk. There is no version of the package were your use case is working as far as I can tell. If it used to work it is an happy accident and I'm really wondering how is it possible?!?

Anyway what you want to do is doable. The rasters will actually be written on disk. The failure happens when merging the output of each chunk. grid_metric() will send you a message of merging failure and return a list of strings (path to raster written on disk) instead of a properly merged raster. The list of string itself will be inccorect since it contains stuff like /tmp/Megaplot.tif instead of /tmp/Megaplot_mean.tif and /tmp/Megaplot_max.tif but at the end, on the disk your rasters are ok. However you must combine them yourself as lidR was not design to work this way.

So here my solution. I'm using raster::writeRaster because terra does not have the by layer option. I removed your custom write function that is not required here. As expected the function fails at the end but the rasters are properly written on disk.

library("lidR")
library("raster")
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")

myMetrics <- function(z){
  list(max=max(z), mean=mean(z))
}

ctg <-readLAScatalog(LASfile)
st_crs(ctg) <- 26917
opt_chunk_size(ctg) <- 150 # to test with 4 chunks (remove)
opt_chunk_alignment(ctg) <- c(0, 70) # to test with 4 chunks (remove)
opt_chunk_buffer(ctg) <- 0
ctg@output_options$drivers$Raster$write <- raster::writeRaster
ctg@output_options$drivers$Raster$param$bylayer <- TRUE
ctg@output_options$drivers$Raster$param$suffix <- quote(names(x))

opt_output_files(ctg) <- paste0(tempdir(), '/{ID}_test')

metricsrast1 <- grid_metrics(ctg, ~myMetrics(Z), res = 10)

max = terra::rast(paste0(tempdir(), '/1_test_max.tif'))
max
terra::plot(max)

To avoid error message you can disable final merging

opt_merge(ctg) = FALSE

But anyway the output of grid_metrics() is incorrect in this specific case since the paths returned do not exists