I am seeking advice on how to build a workflow for individual tree segmentation (ITS) processing with the lidR package using catalog objects.

My current ITS workflow is parallel and leverages the lidR package, but it does not use LAScatalog objects. I am not quite sure how to make the LAScatalog object with a directory of DTMs and CSMs to use in a workflow like the one described here

All of the documentation for ITS in the wiki appears to assume that the DSMs and DTMs are prepared on the fly.

Perhaps we first create .vrt files for the DTMs and CSMs and pass to a custom ITS function?

My current workflow looks something like this

  • Iterate through tiles in parallel from 1 to N (parallel package)
    1. read in lidar with readLAS within a parallel node, also read in dtm and dsm using raster(...)
    2. use lastrees() with the CSM and watershed() to assign trees to points
    3. use lasfilter to remove points not associated with trees
    4. normalize tree heights and remove points below 2m with lasnormalize() and lasfilter()
    5. generate hulls using tree_hulls() and compute .stdmetrics
    6. replace values in PointSourceID with treeID*
    7. write lidar to new .laz file with only tree data
    8. write hulls to .shp file
    9. write hulls attributes to .csv file (slightly redundant)

Perhaps a simple improvement without using LAScatalogs would be to extract data from a catalog file instead of reading in las files directly. This would enable me to buffer tiles slightly and remove edge artifacts (additional workflow steps involved), but would likely add some additional processing time.

My current thought is to build a new function and supply it to catalog_apply(), something like what is below, although in version 3.0 of lidR apparently the segment_trees() function can also accept a catalog.

tree_fn <- function(
 , vrt_dsm 
 , vrt_dtm 
 , th 
 , fn_metrics = .stdmetrics
 , dir_ply_out
  #something like this:
  dtm_all  <- raster::raster(vrt_dtm)
  csm_all  <- raster::raster(vrt_dsm)
  ht_chunk   <- lidR::normalize_height(las_chunk, dtm_all)
  algo_all <- lidR::watershed(csm_all, th = th)
  ht_ws_chunk  <- lidR::lastrees(ht_chunk, algo_all)
  trs_chunk <- lidR::lasfilter(ht_ws_chunk, !is.na(treeID))
  hulls_chunk  <- lidR::tree_hulls(trs_chunk, type = "concave", concavity = 2, func = fn_metrics)

  #these steps I don't know how to do, or perhaps return and let catalog_apply handle writing somehow?
  if(F) lidR::writeLAS(trs_chunk,file.path(dir_las_out,"??"))
  if(F) rgdall::writeOGR(hulls_chunk,file.path(dir_ply_out,"??"))
  if(F) write.csv(hulls_chunk@data,file.path(dir_csv_out,"??"))

  #or return both?

2 Answers 2


Your code is close to be functional. For clarity I won't answer to the question related to writing 3 outputs at a time. You can ask another focused question. So here we will assume you only want to write the shapefile of hulls. Also the recent release of v3.0.0 improved how individual tree segmentation is managed for a LAScatalog.

First it is a good start to make VRTs. This will enable to carry big raster in a lightweight manner. I'm assuming you have a CHM + DTM in a VRT format.

dtm_all  <- raster::raster(vrt_dtm_file)
csm_all  <- raster::raster(vrt_dsm_file)

tree_fn <- function(chunk, vrt_dsm, vrt_dtm, th)
  las_chunk = readLAS(chunk)
  if (is.empty(las_chunk)) return(NULL)

  ht_chunk <- lidR::normalize_height(las_chunk, vrt_dtm)
  algo_all <- lidR::watershed(vrt_dsm, th = th)
  ht_ws_chunk <- lidR::segment_trees(ht_chunk, algo_all, uniqueness = "xxx")
  trs_chunk <- lidR::filter_poi(ht_ws_chunk, !is.na(treeID))
  hulls_chunk <- lidR::delineate_crowns(trs_chunk, type = "concave", concavity = 2, func = .stdmetrics)

  # Removing the buffer is tricky on this one and
  # this is suboptimal. When used standalone with a
  # catalog delineate_crowns() does the job better than that
  hulls_chunk <- raster::crop(hulls_chunk, raster::extent(chunk))

opt_output_files(ctg) <- "templated/path/to/output/HULL_{XCENTER}_{YCENTER}"
opt_chunk_buffer(ctg) <- 40
out <- catalog_apply(ctg, tree_fn, vrt_dsm = csm_all, vrt_dtm = dtm_all, th = 2)

Replace xxx in uniqueness = "xxx" by what you prefer. See the documentation of this new parameter introduced in v3.0.0

From v3.0.0 segment_trees() will only run the segmentation within the bounding box of your LAS object ht_chunk. Not on the whole CHM provided. This is where the improvement of v3.0.0 change the game. Previously this would have crashed by performing the computation on the whole raster.

I didn't test it so you will probably have to fix few things but I'm sure this will give you a good starting points. Test it first on a small dataset.

  • Do you mean that perhaps a cleaner workflow would be to run delineate_crowns() directly on a lidar catalog composed of these same tiles? Would this require first making a duplicate copy of the lidar dataset with z values (elevations) replaced with height values? -thank you!
    – Jacob
    Commented Jun 20, 2020 at 6:33
  • It is the easiest option but the less efficient. You can also keep going with your code. It is up to you
    – JRR
    Commented Jun 20, 2020 at 10:56

Thanks for the help. With assistance, I came to a working solution. It is, however, very slow, even in parallel when reading off an SSD drive.

The previous approach to dealing with trees on interior edges (falling on the edge of two tiles) was to simply to cut them using

hulls_chunk <- raster::crop(hulls_chunk, raster::extent(chunk))

This works ok, but there are issues with topology for the tree objects, so it was necessary to first buffer the tree objects

rgeos::gBuffer(hulls_chunk, byid=TRUE, width=0)

The approach that I landed on, however, was to subset on centroid coordinates that fall within the tile extent. This eliminated most edge artifacts, especially trees split along edges, on interior tiles. The key is to buffer the tile by an amount that is larger than the great majority of tree crowns' radii.

Another important change was to enable reading of chunks to fail. I plan to process in parallel, and the chunk order is sequential (adjacent chunks). This means that for buffered tiles, two nodes may attempt to read the same file at the same time, causing processing failure.

A minor final tweak was to enable partial processing without a complete restart. I was initially motivated to enable this capability when I processed 10% of tiles, and then accidentally pressed the escape key while hovering on the RStudio session. I haven't checked for a "resume" feature of some kind in the lidR package.

dtm_all  <- raster::raster(vrt_dtm_file)
csm_all  <- raster::raster(vrt_dsm_file)

tree_fn <- function(chunk, vrt_dsm, vrt_dtm, th )

  #deal with file access clashes
  for(i in 1:attempt_max){
    las_chunk = try(readLAS(chunk))
    if(!class(las_chunk) == "try-error") i = attempt_max
    else Sys.sleep(5)
  if (is.empty(las_chunk)) return(NULL)

  ht_chunk <- lidR::normalize_height(las_chunk, vrt_dtm)
  algo_all <- lidR::watershed(vrt_dsm, th = th)
  ht_ws_chunk <- lidR::segment_trees(ht_chunk, algo_all, uniqueness = "bitmerge")
  trs_chunk <- lidR::filter_poi(ht_ws_chunk, !is.na(treeID))
  hulls_chunk <- lidR::delineate_crowns(trs_chunk, type = "concave", concavity = 2, func = lidR::.stdmetrics)
  hulls_chunk@data[,c("x","y")] = sp::coordinates(hulls_chunk)
  #remove trees with crowns outside extent
  dat_trs = hulls_chunk@data
  coordinates(dat_trs) = ~x+y
  tile0_ext = as(raster::extent(chunk),"SpatialPolygons")
  in_tile = rgeos::gIntersects(dat_trs, tile0_ext,byid=T)
  hulls_chunk1 = subset(hulls_chunk,subset=as.vector(in_tile))


#enable start / stop of processing
in_files = list.files(dir_in,pattern="[.]las",full.names=T)
out_files = list.files(dir_out,pattern="[.]gpkg",full.names=T)
out_exist = gsub("[.]las",".gpkg",basename(in_files)) %in% gsub("HULLS_","",basename(out_files))
ctg = lidR::readLAScatalog(in_files[!out_exist])

opt_output_files(ctg) <- "templated/path/to/output/HULL_{XCENTER}_{YCENTER}"
opt_chunk_buffer(ctg) <- 15
opt_chunk_size(ctg) = 0
ctg@output_options$drivers$Spatial$extension <- ".gpkg"

plan(multisession, workers = 8L)
out <- lidR::catalog_apply(ctg, tree_fn, vrt_dsm = csm_all, vrt_dtm = dtm_all, th = 2)

The result of this workflow looks very nice: no split trees, no duplicated trees, no overlap trees.

Canopy Hulls from adjacent lidar tiles

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.