How to fill pits for a DEM/CHM when using the LAScatalog engine?

Question

I'm trying to create a DEM/CHM using the LAScatalog engine. I've followed the steps in the lidRbook. In chapter 6 section 6.4 the post processing fills in areas no missing values for single las files as follows:

fill.na <- function(x, i=5) { if (is.na(x)[i]) { return(mean(x, na.rm = TRUE)) } else { return(x[i]) }}
w <- matrix(1, 3, 3)
chm <- rasterize_canopy(las, res = 0.5, algorithm = p2r(subcircle = 0.15), pkg = "terra")
filled <- terra::focal(chm, w, fun = fill.na)
smoothed <- terra::focal(chm, w, fun = mean, na.rm = TRUE)

I want produce the same filled in raster using the LAScatalog engine for multiple files. This is what I have done:

###create a DSM - supply a non-normalised height 
opt_output_files(ctg) <- paste0("~/dir/", "/{*}_dsm")
dsm <- rasterize_terrain(ctg, resolution, tin())

###classify ground returns
opt_output_files(ctg) <- paste0("~/dir/", "{*}_classified")
ws <- seq(3, 12, 3)
th <- seq(0.1, 1.5, length.out = length(ws))
classified_ctg <- classify_ground(ctg, algorithm = pmf(ws = ws, th = th))
gc()


###DTM
opt_output_files(classified_ctg) <- paste0("~/dir/", "/{*}_dtm")
dtm <- rasterize_terrain(classified_ctg, resolution, tin())
gc()

###normalise point cloud density
opt_output_files(classified_ctg) <-  paste0("~/dir/", "/{*}_norm_height")
norm_height_ctg <- normalize_height(classified_ctg, tin()) ###can also use tin or kriging

opt_output_files(norm_height_ctg) <-  paste0("~/dir/", "/{*}_dem")
dem <- rasterize_canopy(norm_height_ctg, resolution, algorithm = dsmtin(max_edge=12))
fill.na <- function(x, i=5) { if (is.na(x)[i]) { return(mean(x, na.rm = TRUE)) } else { return(x[i]) }}
w <- matrix(1, 3, 3)
filled <- NULL
opt_output_files(dem) <-  paste0("~/dir/", "/{*}_filled")
dem_filled <- terra::focal(dem, w, fun = fill.na)

But it outputs this error.

Error in opt_chunk_size(ctg) : no slot of name "chunk_options" for this object of class "SpatRaster" code here

Is this possible to do within the catalog engine? Or do I need to write a loop? I know this function is from the terra package but I was hoping there would be some functionality to allow for this from the LAScatalog engine.

Answer 1

What you are trying to do is to apply a custom function on the catalog. This should look like this (not tested)

chm_pit_filled = function(las, resolution)
{
  chm = rasterize_canopy(las, resolution, algorithm = dsmtin(max_edge=12))
  fill.na <- function(x, i=5) { if (is.na(x)[i]) { return(mean(x, na.rm = TRUE)) } else { return(x[i]) }}
  w <- matrix(1, 3, 3)
  dem_filled <- terra::focal(dem, w, fun = fill.na)
  smoothed <- terra::focal(chm, w, fun = mean, na.rm = TRUE)
  return(smooth)
}

opt_output_files(norm_height_ctg) <- paste0("~/dir/", "/{*}_dem")
opt_chunk_buffer(norm_height_ctg) <- 5 # need a buffer because of focal()
chm = catalog_map(norm_height_ctg, chm_pit_filled, resolution = 0.5)