3

I am porting my existing code to lidR version 2.0 and would like someone to let me know if I am coding efficiently. My objective is:

  • to process a directory of LAS files
  • normalize las files
  • calculate stdmetrics using grid_metrics
  • write out results in a data frame that I can use for subsequent analysis

To perform all these step at once I'm trying to use the function catalog_apply to apply user-defined function. What I did so far is:

Load a LAScatalog and set some processing options

las.catalogue <- catalog("F://Work//LASTest")
opt_cores(las.catalogue) <- 4L
opt_chunk_size(las.catalogue) <- 260
opt_chunk_buffer(las.catalogue) <- 10

Define my function that is expected to read sequentially some chunks of the point cloud, normalize this point cloud chunk and computed metrics in a area based approach.

stat.function <- function(las)
{
  las = lasnormalize(las, tin())
  metrics <- grid_metrics(las, .stdmetrics_z, res = 20)
  # I think I am supposed to remove buffers but need advice on how to do
  return(metrics)
}

Then with catalog_apply I sequentially run my function on my catalog

metrics <- catalog_apply(las.catalogue, stat.function, .options = opt)

To finish I merge the list of outputs for each chunk into a single RasterBrick. Because I want a data.frame I convert the RasterBrickinto a data.frame

metric.merge  <- do.call(raster::merge, metrics)
metrics <- raster::as.data.frame(metric.merge,xy=TRUE, na.rm = TRUE)

I'm not yet comfortable with the new design of the lidRpackage. Am I on the right track to perform my computation?

  • 2
    Have you written all this code in hope, or have you tested it? Does it work? Does it not work? – Spacedman Jan 16 at 15:25
  • 1
    Welcome to GIS SE. Please edit your post to include a single, focused question. Currently, this is bordering on “opinion based”. – Aaron Jan 16 at 17:49
2

Below I give you a more comprehensive answer. But here I suggest a better option in my opinion. In your code, you try to normalize the dataset each time you want to perform the computation. This have a very high cost in term of computation. I recommend to have an already normalized dataset. Normalize your point cloud with LAStools (recommended if you have a license) or with lidR (less efficient). Then, you can simply write:

las.catalogue <- catalog("F://Work//LASTest")
metrics <- grid_metrics(las.catalogue, .stdmetrics_z, res = 20)
metrics <- raster::as.data.frame(metric, xy=TRUE, na.rm = TRUE)

That being said if you want to use catalog_apply, you can! But this function is dedicated to advanced users because it is nothing else than the core engine of the package.

First if you want to use catalog_applysafely you must respect the template given in the documentation of catalog_apply

stat.function <- function(cluster)
{
   las = readLAS(cluster)
   if (is.empty(las)) return(NULL)
   # Do something
   # return(something)
}

The statement if (is.empty(las)) return(NULL) is important because with a chunk_size of 260 some chunks may be empty. This may happen in a middle of a lake removed by the data provider or in non rectangular dataset for examples.

To remove extra pixels of your output raster that belong in the buffer you can use the bounding box of the processed chunk:

metrics <- raster::crop(metrics, raster::extent(cluster))

Also you probably misunderstood the role of the .options argument in catalog_apply. There is a dedicated section in the documentation of catalog apply. So if I understand your code what you tried looks like that:

las.catalogue <- catalog("F://Work//LASTest")
opt_cores(las.catalogue) <- 4L
opt_chunk_size(las.catalogue) <- 260
opt_progress(las.catalogue) <- TRUE
opt_chunk_buffer(las.catalogue) <- 10

stat.function <- function(chunk)
{
  las = readLAS(chunk)
  if (is.empty(las)) return(NULL
  las = lasnormalize(las, tin())
  metrics <- grid_metrics(las, .stdmetrics_z, res = 20)
  metrics <- raster::crop(metrics, raster::extent(cluster)
  return(metrics)
}

opt     <- list(raster_alignment = 20)
metrics <- catalog_apply(las.catalogue, stat.function, .options = opt)
metrics <- do.call(raster::merge, metrics)
metrics <- raster::as.data.frame(metric, xy=TRUE, na.rm = TRUE)

I also recommend you to read the vignettes dedicated to LAScatalog

  • Thank you. I was headed in this direction and had already realized it was best to separate the task of normalizing the LAS files. I have a license to LASTools and so I have done so using it. Thank you for your assistance. – adam.dick Jan 16 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.