0

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).

4
  • Please make a minimal, reproducible example and explain your problem straight to the point. Sounds like a very interesting question, but sadly there are too much information in your message, most being irrelevant to the problem and I did not read everything.
    – JRR
    Jan 18 at 14:21
  • I hope I've succeeded in making an example and the question more succinct Jan 18 at 16:53
  • So you confirm that the code you provided was working? I'm just asking because what you did is clever but also surprising. I just want to be sure that your example is supposed to work.
    – JRR
    Jan 18 at 18:15
  • Yes, the code in my original question has been used to produce many tif-files as desired. Unfortunately, I don't remember what version of the lidR package I used nor the version of R but I produced such files in december 2022 and I tend to work with updated versions of both packages and R. Jan 18 at 20:25

1 Answer 1

0

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

7
  • I don't really understand why the chunks need to be merged (which was the cause of failure in grid_metrics). I was under the impression that when having chunk size = 0 (as I have in my catalog) each whole LAS-file is a chunk. But I will read up on the lidR documentation again to find out how my code could work in the first place. And of course try out your solution. Jan 19 at 0:19
  • By the way, I did try the solution of writing multilayer rasters according to the related question referensed above and it worked but I noticed that a vrt-file was written as well. However, I have no recollection that such a vrt-file was written previously when my code worked. Is this a new feature? Could this be part of the problem/reason why it worked before? Is there any way to turn off the writing of such a vrt-file or is it a necessary feature? Jan 19 at 0:33
  • Chunk is a generic term to talk about a processable sub-area of the dataset. Your dataset is already physically chunked in files, so usually a chunk is a file. The output is merged at the end to provide a useful and manageable output to the users. For rasters the best options is to create a vrt and this feature is in lidR since version 2.0.0 at least. You can disable final merging. See my edit.
    – JRR
    Jan 19 at 2:18
  • I have now tried the suggested solution (with or without removing the alterations to chunk_size and chunk_alignment and with or without opt_merge(ctg) <- FALSE), The problem persist. The error is still "formal argument "sources" matched by multiple actual arguments" and it brakes before anything gets written to disk (the tempdir() never even gets created) and there is no output from grid_metrics. Jan 19 at 13:54
  • 1
    I'm now really embarrassed. I did as you suggested and updated raster and terra. My original code works! I really thought I had updated versions of all packages (I recently reinstalled everything...). While I'm sorry to have wasted your time, I can now (happily) say that my (original) code works, even if it shouldn't. I can produce separate tif-files for each metric for each LAS-file where the filenames follow the pattern "LASfilename_metric.tif". I know I shouldn't, but thanks! Jan 19 at 20:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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