0

In the case of a complex function, I would like to try the autocrop option, instead of crop inside the function, to resize output from catalog_apply() to the chunk size. In the code below, the commented-out code works, but the autocrop option doesn't.

# ---------------------------------------------------------------------------------
# How to crop raster chunks?
# ---------------------------------------------------------------------------------
library(terra)
library(lidR)

setwd("D:/Temp")
interpolationMethod <- tin('', extrapolate = knnidw(k = 10, p = 2, rmax = 50))

sampleLAS <- readLAS(system.file("extdata", "MixedConifer.laz", package="lidR"))
sampleLAS$Classification[sampleLAS$Classification == 1L] <- 4L # veg was misclassified as 1's. change to 4
# plot(sampleLAS, color = "Classification")

sampleLASNorm <- normalize_height(sampleLAS, algorithm = interpolationMethod, na.rm = TRUE, use_class = c(2L, 9L))
writeLAS(sampleLASNorm, "MixedConiferNorm.las", index = TRUE)
ctgNorm <- readLAScatalog("MixedConiferNorm.las")

ctgNorm@output_options$drivers$Raster$param$overwrite <- TRUE

makeCHM_function = function(chunk) {
  las = readLAS(chunk)
  if (is.empty(las))
    return(NULL)
  vegetation <- filter_poi(las, (Classification == 3L | Classification == 4L | Classification == 5L) & Z > 0)
  if (is.empty(vegetation))
    return(NULL)
  if (npoints(vegetation) < 3)
    return(NULL)
  chmPitFree <- rasterize_canopy(las, 1, pitfree(thresholds = c(0,2,5,10,15), max_edge = c(3, 1.5), subcircle = 0.15))
  chmP2R <- rasterize_canopy(las, res = 1, algorithm = p2r()) # fill voids with p2r interpolations
  chmP2R[is.na(chmP2R[])] <- 0 # still some NAs so replacing NA's by zero
  chmPitFree_voidsP2R <- terra::cover(chmPitFree, chmP2R, values = NA) # put them together

  # Can I use AUTOCROP = TRUE in options below to replace the following 3 lines of code. I.e., crop to chunk extent.
  # rast(chmPitFree_voidsP2R)
  # bbox <- lidR::extent(chunk)
  # chmPitFree_voidsP2R  <- terra::crop(chmPitFree_voidsP2R, bbox)
  
  return(chmPitFree_voidsP2R)
}

opt_chunk_size(ctgNorm) <- 20
opt_chunk_buffer(ctgNorm) <- 2
opt_progress(ctgNorm) <- TRUE
opt_output_files(ctgNorm) <- "D:/Temp/chmNorm_tile_{XLEFT}_{YBOTTOM}"

# Here the autocrop isn't doing what I'm expecting (reducing extent to the chunk extent)
opt <- list(autocrop = TRUE)
myCHM <- catalog_apply(ctgNorm, makeCHM_function, .options = opt)
# ---------------------------------------------------------------------------------

1 Answer 1

0

Your code looks correct and works fine. I can't see any problem. You could add opt <- list(autocrop = TRUE, automerge = TRUE)

That being said your function is still very simple (it returns a raster) and you can use catalog_map which simplifies the process by automatically reading and cropping.

makeCHM_function = function(las) 
{
  vegetation <- filter_poi(las, (Classification == 3L | Classification == 4L | Classification == 5L) & Z > 0)
  if (is.empty(vegetation))
    return(NULL)
  if (npoints(vegetation) < 3)
    return(NULL)
  chmPitFree <- rasterize_canopy(las, 1, pitfree(thresholds = c(0,2,5,10,15), max_edge = c(3, 1.5), subcircle = 0.15))
  chmP2R <- rasterize_canopy(las, res = 1, algorithm = p2r()) # fill voids with p2r interpolations
  chmP2R[is.na(chmP2R[])] <- 0 # still some NAs so replacing NA's by zero
  chmPitFree_voidsP2R <- terra::cover(chmPitFree, chmP2R, values = NA) # put them together
  
  return(chmPitFree_voidsP2R)
}

opt_chunk_size(ctgNorm) <- 20
opt_chunk_buffer(ctgNorm) <- 2
opt_progress(ctgNorm) <- TRUE
opt_output_files(ctgNorm) <- "D:/Temp/chmNorm_tile_{XLEFT}_{YBOTTOM}"

myCHM <- catalog_map(ctgNorm, makeCHM_function)
4
  • This code in the answer with catalog_map() works well so I'll go with that. But, for me, my code with autocrop = TRUE still leaves the overlap ... maybe it's my system?
    – Ray J
    Sep 27, 2023 at 11:20
  • I don't think so, Why do you think it leaves overlap. In this stage with filename as output it is not easy to see it. So I added automerge and the output looks good.
    – JRR
    Sep 27, 2023 at 11:28
  • automerge does work, but my output needs to be tiles. Ideally those tiles are 1 km by 1 km without the buffer. The overlap can be seen in a GIS like ArcPro
    – Ray J
    Sep 27, 2023 at 11:39
  • I'll have a look
    – JRR
    Sep 27, 2023 at 12:10

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.