I am attempting to create a raster of pixel metrics for a large collection of lidar data using the lidR:: package. I want to first remove any outlier points in the point cloud, normalize the point cloud to a digital terrain model, and finally, calculate the standard z pixel metrics on a 20 m X 20 m grid. I followed the guidance on the lidR:: package's book and its vignettes for using the catalog_apply() engine. I have created a "low level API" function that first has a conditional to check if the input is a LAScatalog, and then runs the function through catalog_apply , then checks if the input is a LAScluster, and then runs the function directly and clips the chunk buffers from the output, and then finally checks if the input is a LAS, and then explicitly runs the function. I am struggling with getting the function to run properly on a LAScatalog. When I run the function on a LAS file, it works with out error, however, when I run it on a LAScatalog, all chunks show an error on the plot, and when the routine finishes, it throws this error:
Error in any_list[[1]] : subscript out of bounds
In addition: There were 15 warnings (use warnings() to see them)
This error makes me think that I am missing some sort of catalog_apply engine or SpatRaster driver option that tells the function how to merge the output chunks back together to form the final output, but I am not sure which option that would be, and I haven't been able to find any answers on the lidR:: wiki page, vignettes, or book, nor can I find a similar issue here on Stackoverflow. Any advice would be much appreciated. Below is my reproducible example:
##Loading Necessary Packages##
library(lidR)
library(future)
#Reading in the data##
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
ctg <- readLAScatalog(LASfile) # As LAScatalog
las<-readLAS(LASfile) # As LAS
####Custom function####
raster_metrics<-function(las, dtm_ras, grid_size, sensitivity){#start function
if(is(las, "LAScatalog")){#Start first conditional for LASCatalog
options <-list(automerge=TRUE, need_buffer=TRUE)
output<-catalog_apply(las, raster_metrics, grid_size=grid_size, sensitivity=sensitivity, .options=options)
return(output)
} else { #end first condition start first else
if (is(las, "LAScluster")){ #start second conditional for LAScluster
las<-readLAS(las)
if (is.empty(las)){return(NULL)}#Conditional for empty chunk (self contained)
output_tmp<-raster_metrics(las, dtm_ras, grid_size, sensitivity)
bbox<-sf::st_bbox(las)
output<-st_crop(output_tmp, bbox)
return(output)
} else {# End second conditional begin second else
if (is(las, "LAS")){
p95 <- pixel_metrics(las, ~quantile(Z, probs = 0.95), grid_size)
las <- merge_spatial(las, p95, "p95")
las <- filter_poi(las, Z < p95*sensitivity)
las$p95 <- NULL
norm<-las - dtm_ras
output<-pixel_metrics(norm, .stdmetrics_z, grid_size)
return(output)
}else { #end final conditional begin final else
stop("This type is not supported.")
}#end final else
} #end second else
} #end first else
} #end function
##Creating a rasterized dtm to feed to the function##
dtm_ras<-rasterize_terrain(ctg, algorithm = knnidw())
##Defining some function and engine option setttings##
grid_size<-20.0
sensitivity<-1.2
chunk_size<-grid_size*50
chunk_buffer<-grid_size*2
##Setting driver and engine option parameters##
opt_output_files(ctg)<-paste0(tempdir(), "/{XCENTER}_{YCENTER}_{ID}_Norm_Height")
opt_chunk_size(ctg)<-chunk_size
opt_chunk_buffer(ctg)<-chunk_buffer
opt_wall_to_wall(ctg)<-TRUE
opt_stop_early(ctg)<-FALSE
opt_filter(ctg)<-"-drop_withheld"
opt_select(ctg)<-"xyz"
ctg@output_options$drivers$SpatRaster$param$overwrite<-TRUE
##Setting up parallel processing##
plan(multisession, workers = nbrOfWorkers()-1)
set_lidr_threads(nbrOfWorkers()-1)
##Running the function##
example1<-raster_metrics(las=ctg, dtm_ras = dtm_ras, grid_size = grid_size, sensitivity = sensitivity)#Throws error
example2<-raster_metrics(las=las, dtm_ras = dtm_ras, grid_size = grid_size, sensitivity = sensitivity)#Works without error
UPDATE 2/3/2023
Doing a little digging on my own, it appears that this error gets thrown by the internal lidR::: function engine_merge(), which has an argument any_list=. This makes me think that somehow my function violates one of the template rules of catalog_apply(), but I copied the template verbatim from the vignette. Hoping this elucidates the source of my error.
dtm_rasoutput<-catalog_apply(las, raster_metrics, dtm_ras = dtm_ras, grid_size=grid_size, sensitivity=sensitivity, .options=options)
bbox <-terra::ext(las)
output<-terra::crop(output_tmp, bbox)
With the following function it works in sequential mode
raster_metrics<-function(las, dtm_ras, grid_size, sensitivity)
{
if(is(las, "LAScatalog"))
{
options <-list(automerge=FALSE, need_buffer=TRUE)
output<-catalog_apply(las, raster_metrics, dtm_ras = dtm_ras, grid_size=grid_size, sensitivity=sensitivity, .options=options)
return(output)
}
else if (is(las, "LAScluster"))
{
las<-readLAS(las)
if (is.empty(las)){return(NULL)}
output_tmp <- raster_metrics(las, dtm_ras, grid_size, sensitivity)
bbox <-terra::ext(las)
output<-terra::crop(output_tmp, bbox)
return(output)
}
else if (is(las, "LAS"))
{
p95 <- pixel_metrics(las, ~quantile(Z, probs = 0.95), grid_size)
las <- merge_spatial(las, p95, "p95")
las <- filter_poi(las, Z < p95*sensitivity)
las$p95 <- NULL
norm <- las - dtm_ras
output<-pixel_metrics(norm, .stdmetrics_z, grid_size)
return(output)
}
else
{
stop("This type is not supported.")
}
}
However it does not work in parallel because terra's SpatRaster are not serializable. To say it simple, when the dtm_ras is sent to each worker, it no longer exists. This is not an issue with lidR it is an issue with terra. In lidR functions, I use an internal workarounds to deal with SpatRaster by converting them to raster.
On your side the simplest option is to use a RasterLayer from raster.