Constructing a function to pass to lidR::grid_metrics

Question

I am attempting to write an R function, as part of a package, that takes a lidar tile (LAS object) and a vector of cut-points defining vertical layers, and returns a RasterStack where each layer has a layer point count as the cell value. The function makes use of the grid_metrics function in the lidR package. Here is a simplified version of the code to illustrate:

# Main package function
get_layer_counts <- function(las, breaks, res = 10) {
  f <- function(Z) layer_count_metrics(Z, breaks)
  stats <- lidR::grid_metrics(las, func = f(Z), res = res)
  lidR:::as.raster.lasmetrics(stats)
}

# Metrics function to count points within height intervals
layer_count_metrics <- function(Z, breaks) {
  zcat <- findInterval(Z, breaks, left.open = TRUE)
  as.list(table(zcat))
}

Unfortunately this gives an error: "Error in eval(call) : object 'f' not found" which comes from the lidR:::lasaggregate function.

It seems (I think) that unless the function passed to grid_metrics is defined in the global environment it will not be found. Is that correct and, if so, can anyone suggest a work-around?

Answer 1

grid_metrics is designed in a way that your expression is evaluated within the frame of the data.table that contains the point cloud. Your function is found within the loaded R packages. Thus your code cannot work for two reasons:

1. f is defined within the frame of get_layer_counts
2. breaks does not exists within the frame of the point cloud

This version works:

# Main package function
get_layer_counts <- function(las, res = 10) 
{
  stats <- lidR::grid_metrics(las, func = layer_count_metrics(Z, breaks = seq(0,20,2)), res = res)
  lidR:::as.raster.lasmetrics(stats)
}

# Metrics function to count points within height intervals
layer_count_metrics <- function(Z, breaks) 
{
  zcat <- findInterval(Z, breaks, left.open = TRUE)
  tcount = as.list(table(factor(zcat, levels = 1:length(breaks))))
  names(tcount) <- paste0("layer", 1:length(breaks))
  tcount
}

library(lidR)

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

stats = get_layer_counts(las)

You can notice one important drawback. You can't propagate the argument breaks through the get_layer_count function because it looks for breaks within the frame of the point cloud. When explicitly passed there is no name to look for and thus it just works. I hope to improve that in future releases.

Also notice that I added factor before table. Indeed the function must returns always the same number of metrics for each pixel. Otherwise the pixel aggregation fails.