catalog engine used in batch processing of tree detection

R/gaps.edges.detection.Rmd

@@ -170,7 +170,7 @@ gaps1r <- lidaRtRee::gap_detection(chm,
 # display detected gaps
 par(mfrow = c(1, 3))
 # plot gaps
-terra::plot(chm, main = "Canopy heightmodel")
+terra::plot(chm, main = "Canopy height model")
 terra::plot(gaps1$gap_id > 0, main = "Gaps", legend = FALSE)
 terra::plot(gaps1r$gap_id > 0, main = "Gaps extended by reconstruction", legend = FALSE)
 ```

R/tree.detection.Rmd

@@ -570,13 +570,6 @@ lidaRtRee::plot_tree_inventory(tree_metrics_h
 # write.csv(tree_metrics_h, "./data/output/correct_detections.csv", row.names = FALSE)
 # write.csv(tree_inventory_chablais3, "./data/output/tree_inventory_chablais.csv", row.names = FALSE)
 ```
-```{r test, include=FALSE, fig.dim=c(9, 9)}
-# add attributes of apices to polygons by merging
-# !!!! for some reason the segment_id is messed up during the process
-# segments_v@data <- merge(segments_v@data, metrics, by.x = "segments_id", by.y = "seg_id")
-# remove ground segment
-# segments_v <- segments_v[segments_v$segments_id != 0,]
-```
 
 ## Display point cloud
 
@@ -699,56 +692,55 @@ for (i in 1:nrow(apices))
 
 ## Batch processing
 
-The following code exemplifies how to process numerous LAS files and extract apices for the whole area with parallel processing. The processing runs faster if data is provided as chunks to the segmentation algorithm, and results are then aggregated, rather than running on the full coverage of the data. 
-In order to avoid border effects, chunks are provided to the algorithm as overlapping tiles. Results are cropped to prevent the same tree from appearing twice in the final results. Tile size and buffer size are important parameters :
+The following code exemplifies how use the catalog engine of package `lidR` to process numerous LAS files and extract apices for the whole area with parallel processing. The catalog engine provides the data as chunks to the segmentation algorithm, and results are aggregated afterwards. 
+In order to avoid border effects, a buffer zone is added to chunks. Chunk results are cropped to prevent the same tree from appearing twice in the final results. Chunk size and buffer size are important parameters :
 
-- tile size is a trade-off between the number of chunks to process and the amount of RAM required to process a single tile ;
-- buffer size is a trade-off between redundant processing of the overlap area, and assuring that a tree is which treetop is located at the border of a tile has its full crown within the buffer size.
+- `chunk_size` is a trade-off between the number of chunks to process and the amount of RAM required to process a single chunk ;
+- `chunk_buffer` size is a trade-off between redundant processing of the overlap area, and assuring that a tree which treetop is located at the border of a chunk has its full crown within the buffer size.
 
-Tiles can be processed with parallel computing, within limits of the cluster's RAM and number of cores.
+Chunks can be processed with parallel computing, within limits of the cluster's RAM and number of cores.
 
 The steps for processing a batch of las/laz files are :
 
-- build catalog of files
-- provide tiling parameters
+- build catalog of files and set chunk processing options
 - provide segmentation parameters
 - provide output parameters
 - set cluster options for parallel processing
-- compute the X and Y coordinates of tiles
-- parallelize the processing with the package `future`, by sending to the parallel sessions the coordinates of tiles to process, and a function with the instructions to proceed :
-    - load tile corresponding to the coordinates
+- use the catalog engine to proceed with the following steps for each chunk:
     - compute CHM
     - segment and extract apices
+    - remove apices inside buffer zone
 - aggregate list of results
 
 Be aware that in case tree segments are vectorized, some obtained polygons might overlap. The segmentation algorithm might not be deterministic and borders are sometimes not consistent when adjacent polygons are pasted from different tiles.
 
 ```{r batch, include=TRUE, warning = FALSE}
 rm(list = ls())
-# BUILD CATALOG OF FILES
+# BUILD CATALOG OF FILES AND SET CATALOG PROCESSING OPTIONS
 # folder containing the files
 lazdir <- "../data/forest.structure.metrics"
-# build catalog
-cata <- lidR::readALSLAScatalog(lazdir)
-# disable display of catalog processing
-lidR::opt_progress(cata) <- FALSE
-# set coordinate system
-lidR::projection(cata) <- 2154
-# option to read only xyzc attributes (coordinates, intensity, echo order and classification) from height files
-lidR::opt_select(cata) <- "xyz"
-#
-# BATCH PROCESSING PARAMETERS
-# tile size to split area into chunks to process
+# build catalog and set options
+# - progress: disable display of catalog processing
+# - select: read only xyzc attributes (coordinates, intensity,
+# echo order and classification) from height files
+# - chunk_size: tile size to split area into chunks to process,
 # trade-off between RAM capacity VS total number of tiles to process
-tile_size <- 70 # here 70 for example purpose with small area
-# buffer size: around each tile to avoid border effects on segmentation results
+# here 70 for example purpose with small area
+# - buffer_size: around each tile to avoid border effects on segmentation results
 # trade-off between making sure a whole tree crown is processed in case its top
 # is on the border VS duplicate processing
-buffer_size <- 10 # 5 m is minimum, 10 is probably better depending on tree size
+# 5 m is minimum, 10 is probably better depending on tree size
+cata <- lidR::readALSLAScatalog(lazdir,
+                                progress = FALSE, 
+                                select = "xyz",
+                                chunk_size = 70,
+                                chunk_buffer = 10)
+# set coordinate system
+lidR::projection(cata) <- 2154
 #
 # TREE SEGMENTATION PARAMETERS
 # set raster resolution
-resolution <- 1
+res <- 1
 #
 # OUTPUT PARAMETERS
 # option to vectorize crowns (set to FALSE if too slow)
@@ -764,112 +756,79 @@ future::plan("multisession", workers = 2L)
 # remove warning when using random numbers in parallel sessions
 options(future.rng.onMisuse = "ignore")
 #
-# COORDINATES OF TILES
-bbox <- sf::st_bbox(cata)
-# low left corner
-x <- seq(
-  from = floor(bbox$xmin / tile_size) * tile_size,
-  by = tile_size,
-  to = ceiling(bbox$xmax / tile_size - 1) * tile_size
-) # [1:2]
-length_x <- length(x)
-y <- seq(
-  from = floor(bbox$ymin / tile_size) * tile_size,
-  by = tile_size,
-  to = ceiling(bbox$ymax / tile_size - 1) * tile_size
-) # [1:2]
-length_y <- length(y)
-# repeat coordinates for tiling while area
-x <- as.list(rep(x, length_y))
-y <- as.list(rep(y, each = length_x))
-#
-# PARALLEL PROCESSING
-# function takes coordinates of tile as arguments
-resultats <- future.apply::future_mapply(
-  # i and j are the coordinated of the low left corner of the tile to process
-  FUN = function(i, j) {
-    # initialize output
-    output <- list()
-    output$name <- paste0(i, "_", j)
-    # extract tile plus buffer
-    las_tile <- lidR::clip_rectangle(
-      cata,
-      i - buffer_size,
-      j - buffer_size,
-      i + tile_size + buffer_size,
-      j + tile_size + buffer_size
+# FUNCTION TO APPLY TO CATALOG
+routine <- function(chunk, resolution = res)
+{
+  # empty list for output
+  output <- list()
+  # bounding box of region, without buffer
+  bbox <- sf::st_bbox(chunk)
+  # read chunk of data
+  las <- lidR::readLAS(chunk)
+  # return NULL if empty
+  if (lidR::is.empty(las)) return(NULL)
+  # normalization if required
+  # las <- lidR::normalize_height(las, lidR::tin())
+  # in this example LAS tiles are already normalized
+  # compute canopy height model
+  chm <- lidR::rasterize_canopy(las, resolution, algorithm = lidR::p2r(), pkg = "terra")
+  #
+  # check all chm is not NA
+  if (all(is.na(terra::values(chm)))) return(NULL)
+  #
+  # tree detection (default parameters)
+  segms <- lidaRtRee::tree_segmentation(chm)
+  # tree extraction
+  apices <- lidaRtRee::tree_extraction(
+    segms$filled_dem,
+    segms$local_maxima,
+    segms$segments_id
+  )
+  # remove apices in buffer area
+  apices <- sf::st_crop(apices, sf::st_bbox(chunk))
+  # add tile id
+  apices$tile <- paste0(bbox$xmin, "_", bbox$ymin)
+  # put apices into output slot
+  output$apices <- apices
+  #
+  # crop chm if requested output (global variable)
+  if (out_chm | out_save_chm)
+  {
+    chm <- terra::crop(chm, terra::ext(bbox$xmin, bbox$xmax, bbox$ymin, bbox$ymax))
+  }
+  # output chm if asked for 
+  if (out_chm)
+  {
+    output$chm <- terra::wrap(chm)
+  }
+  # if save on disk
+  if (out_save_chm) {
+    terra::writeRaster(chm,
+    file = paste0("chm_", bbox$xmin, "_", bbox$ymin, ".tif"),
+    overwrite = TRUE
     )
-    # check if points are present
-    if (nrow(las_tile) > 0) {
-      # normalization if required
-      # las_tile <- lidR::normalize_height(las_tile, lidR::tin())
-      # in this example LAS tiles are already normalized
-      # compute canopy height model
-      chm <- lidR::rasterize_canopy(las_tile, resolution, algorithm = lidR::p2r(), pkg = "terra") 
-      # check all chm is not NA
-      if (!all(is.na(terra::values(chm)))) {
-        # output chm if asked for
-        if (out_chm) {
-          output$chm <- terra::wrap(terra::crop(chm, terra::ext(
-            i, i + tile_size,
-            j, j + tile_size
-          )))
-        }
-        # save on disk
-        if (out_save_chm) {
-          terra::writeRaster(terra::crop(
-            chm,
-            terra::ext(
-              i, i + tile_size,
-              j, j + tile_size
-            )
-          ),
-          file = paste0("chm_", i, "_", j, ".tif"),
-          overwrite = TRUE
-          )
-        }
-        #
-        # tree detection
-        segms <- lidaRtRee::tree_segmentation(chm)
-        # tree extraction
-        apices <- lidaRtRee::tree_extraction(
-          segms$filled_dem,
-          segms$local_maxima,
-          segms$segments_id
-        )
-        # remove apices in buffer area
-        # copy coordinates
-        apices <- cbind(apices, sf::st_coordinates(apices))
-        apices <- apices[apices$X >= i & apices$X < i + tile_size &
-          apices$Y >= j & apices$Y < j + tile_size, ]
-        # add tile id to apices to avoid duplicates in final file
-        apices$tile <- paste0(i, "_", j)
-        # convert to vectors if option is TRUE
-        if (out_vectorize_apices) {
-          # vectorize
-          apices_v <- terra::as.polygons(segms$segments_id)
-          # remove polygons which treetop is in buffer
-          apices_v <- apices_v[is.element(apices_v$segments_id, apices$id), ]
-          # convert to sf
-          apices_v <- sf::st_as_sf(apices_v)
-          names(apices_v)[1] <- "id"
-          # add attributes
-          apices_v <- merge(apices_v, sf::st_drop_geometry(apices), all.x = TRUE)
-          # save in list
-          output$apices_v <- apices_v
-        }
-        # save apices in list
-        output$apices <- apices
-      } # end of raster is not all NAs check
-    } # end of nrow LAS check
-    output
-  }, x, y, SIMPLIFY = FALSE
-) # function applied to the lists of coordinates (x and y)
+  }
+  # convert to vectors if option is TRUE
+  if (out_vectorize_apices) {
+    # vectorize
+    apices_v <- terra::as.polygons(segms$segments_id)
+    # remove polygons which treetop is in buffer
+    apices_v <- apices_v[is.element(apices_v$segments_id, apices$id), ]
+    # convert to sf
+    apices_v <- sf::st_as_sf(apices_v)
+    names(apices_v)[1] <- "id"
+    # add attributes
+    apices_v <- merge(apices_v, sf::st_drop_geometry(apices), all.x = TRUE)
+    # save in list
+    output$apices_v <- apices_v
+  }
+  output
+} # end of routine function
 #
-# RESULTS AGGREGATION
-# extract results from nested list into separate lists and then bind data
-id <- unlist(lapply(resultats, function(x) x[["name"]]))
+# APPLY FUNCTION TO CATALOG
+resultats <- lidR::catalog_apply(cata, routine)
 #
+# RESULTS AGGREGATION
 # apices
 apices <- lapply(resultats, function(x) x[["apices"]])
 # remove NULL elements
@@ -884,7 +843,6 @@ if (out_chm) {
   # merge chm
   # no names in list otherwise do.call returns an error
   chm_all <- do.call(terra::merge, chm)
-  names(chm) <- id
 }
 # apices_v
 if (out_vectorize_apices) {