Gaps and edges detection • lidaRtRee

Source: vignettes/gaps.edges.detection.Rmd

This tutorial presents R code for forest gaps and edges detection from Airborne Laser Scanning (ALS) data. The workflow is based on functions from R packages lidaRtRee (tested with version 4.0.7) and lidR (tested with version 4.0.4).

Study area and ALS data

The study area is a 500 m x 500 m forest located in the Jura mountain (France). The following code shows how to extract the ALS data from a folder containing normalized LAS files, and compute the Canopy Height Model (CHM) at 1 m resolution. For processing a set of several LAS/LAZ files, please see the Batch processing section.

# build catalog from folder with normalized LAS/LAZ files
# option to read only xyzc attributes (coordinates, intensity, echo order and 
# classification) from height files
# option to drop points with height above 50
cata <- lidR::readALSLAScatalog("/las.folder/",
                                progress = FALSE,
                                select = "xyzirnc",
                                filter = "-drop_z_above 50"
                                )
# "/media/reseau/lessem/ProjetsCommuns/Lidar/data/Franche-Comte/norm.laz/"
# set coordinate system
lidR::projection(cata) <- 2154
#
# define region of interest where points should be extracted
centre <- c(944750, 6638750)
size <- 250
# extract data on square area centered on centre
las <- lidR::clip_rectangle(
  cata, centre[1] - size, centre[2] - size, centre[1] + size,
  centre[2] + size
)
# compute canopy height model from normalized point cloud
chm <- lidR::rasterize_canopy(
  las, 1, terra::ext(centre[1] - size, centre[1] + size, centre[2] - size,
  centre[2] + size
)
# save chm in Rdata file
# save(chm, file="chm.rda")

If the CHM has been previously calculated.

# load dataset
load(file = "./data/gap.detection/chm.rda")
chm <- terra::rast(chm)

Missing, low and high values are replaced for better visualization and processing.

# replace NA by 0
chm[is.na(chm)] <- 0
# replace negative values by 0
chm[chm < 0] <- 0
# display CHM and areas where vegetation height is lower than 1m
par(mfrow = c(1, 2))
terra::plot(chm, main = "CHM (m)")
terra::plot(chm < 1, main = "CHM < 1m")

Gap detection

Gaps will be detected according to two different sets of criteria:

Vegetation height lower than 1 m, distance to surrounding vegetation larger than 0.5 times vegetation height and minimum gap surface of 100 m².
Vegetation height lower than 1 m, and no restrictions on distance to surrounding vegetation and surface.

The procedure of gap detection is exemplified in Glad et al. (2020).

# Perform gap detection with gap width criterion
gaps1 <- lidaRtRee::gap_detection(chm, ratio = 2, gap_max_height = 1,
                                  min_gap_surface = 100)
# Perform gap detection without gap width nor min gap surface criteria
gaps2 <- lidaRtRee::gap_detection(chm, ratio = NULL, gap_max_height = 1,
                                  min_gap_surface = 0)

The following figure displays the obtained gaps, colored by size.

# display CHM and areas where vegetation height is lower than 1m
par(mfrow = c(1, 2))
# max value of surface (log)
z_lim <- log(max(terra::values(gaps1$gap_surface), na.rm = TRUE))
# plot gap surface
terra::plot(log(gaps1$gap_surface),
  main = "log(Gap surface) (1)", range = c(0, ceiling(z_lim)),
  col = rainbow(255, end = 5 / 6)
)
terra::plot(log(gaps2$gap_surface),
  main = "log(Gap surface) (2)", range = c(0, ceiling(z_lim)),
  col = rainbow(255, end = 5 / 6)
)

Gaps can be vectorized for display over the original CHM.

# convert to vector
gaps1_v <- terra::as.polygons(gaps1$gap_id)
# display gaps over original CHM
terra::plot(chm, main = "Gaps over CHM (1)")
terra::plot(gaps1_v, add = T)

Statistics on gaps (number, size) can be derived from the raster objects. Example for gaps (2)

# compute gaps surface
gaps_df <- data.frame(t(table(terra::values(gaps2$gap_id))))[, -1]
# change column names
names(gaps_df) <- c("id", "surface")
# convert surface in pixels to square meters
gaps_df$surface <- gaps_df$surface * terra::xres(chm)^2
# data.frame
head(gaps_df, n = 3)

##   id surface
## 1  1   18226
## 2  6      26
## 3  7     113

# cumulative distribution of gap surface by surface
Hmisc::Ecdf(log(gaps_df$surface), weights = gaps_df$surface,
            xlab = "log(surface)")

Histogram can be used to compare the distribution of surface in different classes of gap surface. For better visualization of the differences between the set of criteria, gaps bigger than 500 m² are not taken into account.

# extract vector of total gap surface of pixel
surface1 <- terra::values(gaps1$gap_surface)
surface2 <- terra::values(gaps2$gap_surface)
# remove large gaps
surface1[surface1 > 500] <- NA
surface2[surface2 > 500] <- NA
# display histogram of surface
hist(surface2, col = "blue", xlab = "Gap surface classes",
     main = "Histogram of gap surface")
hist(surface1, border = "red", col = "red", density = 30, add = TRUE)
legend("topright", c("1", "2"), fill = c("red", "blue"))

Edges detection

In the previous part, a pixel that fulfills the maximum height criterion but does not fulfill the distance ratio criterion (i.e. it is too close to surrounding vegetation based on the distance / vegetation height ratio) is removed from gap surface. For edge detection, it seems more appropriate to integrate those pixels in the gap surface, otherwise some edges would be detected inside flat areas. The difference is exemplified in the following plots. The second image exhibits more gaps, because some gaps which are not reconstructed do not reach the minimum surface. The gaps in the second image are also larger, because gaps extend to all neighboring pixels complying with the height threshold, even if they do not comply with the distance criterion.

# load canopy height model
chm <- terra::rast(lidaRtRee::chm_chablais3)
chm[is.na(chm)] <- 0
# Perform gap detection with distance ratio criterion
gaps1 <- lidaRtRee::gap_detection(chm,
  ratio = 2, gap_max_height = 1, min_gap_surface = 50,
  nl_filter = "None"
)
# Perform gap detection with distance ratio criterion and gap reconstruction
gaps1r <- lidaRtRee::gap_detection(chm,
  ratio = 2, gap_max_height = 1, min_gap_surface = 50,
  gap_reconstruct = TRUE, nl_filter = "None"
)
# display detected gaps
par(mfrow = c(1, 3))
# plot gaps
terra::plot(chm, main = "Canopy height model")
terra::plot(gaps1$gap_id > 0, main = "Gaps", col = "green", legend = FALSE)
terra::plot(gaps1r$gap_id > 0, main = "Gaps extended by reconstruction",
            col = "green", legend = FALSE)

Edge detection is performed by extraction of the difference between a binary image of gaps and the result of a morphological erosion or dilation applied to the same image.

# Perform edge detection by erosion
edge_inside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0)
# Perform edge detection by dilation
edge_outside <- lidaRtRee::edge_detection(gaps1r$gap_id > 0, inside = FALSE)
# display zoom on edges
par(mfrow = c(1, 2))
terra::plot(edge_inside, main = "Edges (detection by erosion)", legend = FALSE)
terra::plot(edge_outside, main = "Edges (detection by dilation)",
            legend = FALSE)

Percentage of edges can be calculated as the ratio between edge pixels surface and total surface, multiplied by 100. The result is dependent on the raster resolution. Obtained percentage for the “erosion” method is 4.1, whereas it is 3.9 for method “dilation”.

Batch processing

In case of a LAS-catalog object referencing LAS/LAZ files (which should be normalized; if not, please refer to the preprocessing tutorial, the lidR catalog engine provides the data as chunks to the detection algorithm, and results are aggregated afterwards. In order to avoid border effects, a buffer zone is added to chunks.

chunk_size is a trade-off between the number of chunks to process and the amount of RAM required to process a single chunk. It is advisable to set its value to 0 (processing by file) in case files are non-overlapping tiles and they can be loaded in RAM.
chunk_buffer size is a trade-off between redundant processing of the overlap area, and making sure that a gap at the border of a chunk is fully contained in the buffer and hence total surface is correctly estimated.

In case a gap spans over two files:

if it is included in the buffer, then the gap portion in each tile has a different id, but the total surface of the gap is correct.
if it is partially included in the buffer, then the total surface is under-estimated.

Chunks can be processed with parallel computing, within limits of the cluster’s RAM and number of cores. I

The steps for processing a batch of las/laz files are :

build catalog of files and set catalog processing options
provide gap detection parameters
set cluster options for parallel processing
apply lidaRtRee::gap_detection() function to catalog.

rm(list = ls())
# BUILD CATALOG OF FILES AND SET CATALOG PROCESSING OPTIONS
# folder containing the files
lazdir <- "./data/forest.structure.metrics"
# build catalog and set options
# - progress: disable or not 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
# here 70 for example purpose with small area
# - chunk_buffer: 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 5 m is minimum, 
# 10 is probably better depending on tree size
cata <- lidR::readALSLAScatalog(lazdir,
                                progress = FALSE, 
                                select = "xyzc",
                                chunk_size = 0,
                                chunk_buffer = 20)
# set coordinate system
lidR::projection(cata) <- 2154
#
# CLUSTER PARAMETERS
# specify to use two parallel sessions
future::plan("multisession", workers = 2L)
# remove warning when using random numbers in parallel sessions
options(future.rng.onMisuse = "ignore")
#
# DETECTION PARAMETERS
resolution <- 1
#
# PROCESSING
# perform tree detection with crown polygon output
gaps <- lidaRtRee::gap_detection(cata, res = resolution)

The figures below show that the large gap at the top-left had by chance the same id in both chunks, but its surface is underestimated when computed from the right tile.

# display CHM and areas where vegetation height is lower than 1m
par(mfrow = c(1, 2))
# plot gap surface
terra::plot(gaps$gap_id,
  main = "Gap id", type = "classes",
  col = rainbow(8)
)
terra::plot(log(gaps$gap_surface),
  main = "log(Gap surface)", type = "continuous",
  col = rainbow(255, end = 5 / 6)
)

In case the spatial extent is large, it might not be possible to hold the result for all files in memory, so that results for each chunk should be save on disk (see corresponding documentation of lidR package).

References

Glad, Anouk, Björn Reineking, Marc Montadert, Alexandra Depraz, and Jean-Matthieu Monnet. 2020. “Assessing the Performance of Object-Oriented LiDAR Predictors for Forest Bird Habitat Suitability Modeling.” Remote Sensing in Ecology and Conservation 6 (1): 5–19. https://doi.org/10.1002/rse2.117.