rlidr

Subscript out of bounds error with custom function in lidR::catalog_apply()


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.


Solution

    1. You missed to propagate dtm_ras
    output<-catalog_apply(las, raster_metrics,  dtm_ras = dtm_ras, grid_size=grid_size, sensitivity=sensitivity, .options=options)
    
    1. You used incorrect package to crop
    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.