Tree segmentation
Jean-Matthieu Monnet
tree.detection.Rmd
Source: vignettes/tree.detection.Rmd
The code below presents a tree segmentation workflow from Airborne
Laser Scanning (lidar remote sensing) 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), and it includes the following steps:
- treetop detection,
- crown segmentation,
- accuracy assessment with field inventory,
- species classification
Steps 1 and 3 are documented in (Monnet et al. 2010; Monnet 2011). The detection performance of this algorithm was evaluated in a benchmark (Eysn et al. 2015), which part of the dataset is published in open acces (Monnet and Eysn 2023).
Material
Field inventory
The field inventory corresponds to a 50 m x 50 m plot located in the
Chablais mountain (France) (Monnet 2011,
34). All trees with a diameter at breast height above 7.5 cm are
inventoried. The data is available in package
lidaRtRee
.
# load dataset from package (default)
data(tree_inventory_chablais3, package = "lidaRtRee")
Otherwise you can load your own data provided positions and heights are measured.
# import field inventory
fichier <- "chablais3_listeR.csv"
tree_inventory_chablais3 <- read.csv(file = fichier, sep = ";", header = F,
stringsAsFactors = TRUE)
names(tree_inventory_chablais3) <- c("x", "y", "d", "h", "n", "s", "e", "t")
# save as rda for later access
# save(tree_inventory_chablais3,file="tree_inventory_chablais3.rda")
Attributes are:
-
x
easting coordinate in Lambert 93 (epsg:2154) -
y
northing coordinate in Lambert 93 -
d
diameter at breast height (cm) -
h
tree height (m) -
n
tree number -
s
species abreviated as GESP (GEnus SPecies) -
e
appearance (0: missing or lying, 1: normal, 2: broken treetop, 3: dead with branches, 4: snag) -
t
tilted (0: no, 1: yes)
## x y d h n s e t
## 1 974353.3 6581643 37.6 23.6 1 PIAB 1 0
## 2 974351.0 6581648 15.7 13.9 2 PIAB 1 0
## 3 974348.5 6581650 34.2 23.0 3 PIAB 1 0
Function plot_tree_inventory
is designed to plot forest
inventory data. Main species are European beech (Fagus
sylvatica), Norway spruce (Picea abies) and silver fir
(Abies alba).
# display inventoried trees
lidaRtRee::plot_tree_inventory(tree_inventory_chablais3
[, c("x", "y")],
tree_inventory_chablais3$h,
species = as.character(tree_inventory_chablais3$s)
)
The ggplot2
package also provides nice outputs.
# use table of species of package lidaRtRee to always use the same color for a
# given species
plot.species <- lidaRtRee::species_color()[levels(tree_inventory_chablais3$s),
"col"]
library(ggplot2)
ggplot(tree_inventory_chablais3, aes(x = x, y = y, group = s)) +
geom_point(aes(color = s, size = d)) +
coord_sf(datum = 2154) +
scale_color_manual(values = plot.species) +
scale_radius(name = "Diameter") +
geom_text(aes(label = n, size = 20), hjust = 0, vjust = 1) +
labs(color = "Species") # titre de la légende
We define the region of interest (ROI) to crop ALS data to
corresponding extent before further processing. ROI is set on the extent
of tree inventory, plus a 10 meter buffer. The tree table is converted
to a sf
object to make spatial processing easier in the
following steps.
# duplicate coordinates to ensure they remain in the data.frame
tree_inventory_chablais3[, c("X", "Y")] <-
tree_inventory_chablais3[, c("x", "y")]
# convert to spatial sf object
tree_inventory_chablais3 <- sf::st_as_sf(tree_inventory_chablais3,
coords = c("X", "Y"), crs = 2154)
# buffer to apply around ROI (meters)
ROI_buffer <- 10
# ROI limits: bounding box of trees
ROI_range <- round(sf::st_bbox(tree_inventory_chablais3))
Airborne Laser Scanning data
In this tutorial, ALS data available in the lidaRtRee
package is used.
# load data in package lidaRtRee (default)
LASfile <- system.file("extdata", "las_chablais3.laz", package="lidaRtRee")
las_chablais3 <- lidR::readLAS(LASfile)
# set projection
lidR::projection(las_chablais3) <- 2154
Otherwise, ALS data is loaded with functions of package
lidR
. First a catalog of files containing ALS data is
built. Then points located inside our ROI are loaded.
# directory for laz files
lazdir <- "./data/tree.detection"
# build catalog of files
# specifying ALS data
cata <- lidR::readALSLAScatalog(lazdir)
# set coordinate system
lidR::projection(cata) <- 2154
# extract points in ROI plus additional 5m buffer
las_chablais3 <- lidR::clip_roi(
cata,
ROI_range + (ROI_buffer + 5) * c(-1, -1, 1, 1)
)
# save as rda for easier access:
# save(las_chablais3, file="las_chablais3.rda", compress = "bzip2")
Data preparation
Digital Elevation Models
From the ALS point cloud, digital elevation models are computed (Monnet 2011, 43–46). The Digital Terrain Model (DTM) represents the ground surface, it is computed by bilinear interpolation of points classified as ground. The Digital Surface Model (DSM) represents the canopy surface, it is computed by retaining in each raster’s pixel the value of the highest point contained in that pixel. The Canopy Height Model (CHM) is the normalized height of the canopy. It is computed by subtracting the DTM to the DSM.
# define extent and resolution of raster
output_raster <- terra::rast(resolution = 0.5,
xmin = ROI_range$xmin - ROI_buffer,
xmax = ROI_range$xmax + ROI_buffer,
ymin = ROI_range$ymin - ROI_buffer,
ymax = ROI_range$ymax + ROI_buffer,
crs = sf::st_crs(las_chablais3)$wkt
)
# terrain model computed from points classified as ground
dtm <- lidR::rasterize_terrain(las_chablais3, output_raster, lidR::tin())
# surface model
dsm <- lidR::rasterize_canopy(las_chablais3, output_raster, lidR::p2r())
# canopy height model
chm <- dsm - dtm
# save for later use
# chm_chablais3 <- terra::wrap(chm)
# save(chm_chablais3, file = "~/R/lidaRtRee/data/chm_chablais3.rda")
Visual comparison of field inventory and ALS data
A plot mask is computed from the inventoried positions, using a height-dependent buffer. Indeed, tree tops are not necessarily located verticaly above the trunk positions. In order to compare detected tree tops with inventoried trunks, buffers have to be applied to trunk positions to account for the non-verticality of trees.
# plot mask computation based on inventoried positions
# convex hull of union of points geometry
mask_chull <-
sf::st_convex_hull(sf::st_union(sf::st_geometry(tree_inventory_chablais3)))
# union of buffers around points geometry
mask_tree <-
sf::st_union(sf::st_buffer(sf::st_geometry(tree_inventory_chablais3),
2.1 + 0.14 * tree_inventory_chablais3$h))
# union of convex hull and tree buffers
mask_plot_v <- sf::st_union(mask_chull, mask_tree)
# rasterize mask
mask_plot <- terra::rasterize(terra::vect(mask_plot_v), dsm)
Displaying inventoried trees on the CHM shows a pretty good agreement of crowns visible in the CHM with trunk locations and sizes.
# display CHM
terra::plot(chm,
col = gray(seq(0, 1, 1 / 255)),
main = "Canopy Height Model and tree positions"
)
# add inventoried trees
lidaRtRee::plot_tree_inventory(tree_inventory_chablais3[, c("x", "y")],
tree_inventory_chablais3$h,
species = as.character(tree_inventory_chablais3$s), add = TRUE
)
# display plot mask
terra::plot(mask_plot_v, border = "red", add = TRUE)
Tree delineation
Tree delineation consists in two steps, which are exemplified in the following paragraph:
- image processing for local maxima detection, selection and image
segmentation (function
lidaRtRee::tree_segmentation()
) - extraction of tree information (function
lidaRtRee::tree_extraction()
).
The function lidaRtRee::tree_detection()
is a shortcut
that combines those two steps, and allows to process a catalog of files.
It is presented in the batch processing paragraph.
Segmentation
Tree segmentation is performed on the Canopy Height Model by using
the function lidaRtRee::tree_segmentation()
which consists
in the following steps:
- image pre-processing (non-linear filtering and smoothing for noise removal),
- local maxima filtering and selection for apex (local maximum) detection,
- image segmentation with a watershed algorithm for crown delineation.
The first two steps are documented in Monnet (2011), pp. 47-52.
# tree detection (default settings), applied on canopy height model
segms <- lidaRtRee::tree_segmentation(chm)
#
par(mfrow = c(1, 3))
# display pre-processed chm
terra::plot(segms$smoothed_dem, main = "Pre-processed CHM")
# display selected local maxima
terra::plot(segms$local_maxima, main = "Selected local maxima")
# display segments, except ground segment
dummy <- segms$segments_id
dummy[dummy == 0] <- NA
terra::plot(dummy, main = "Segments (random colors)",
col = rainbow(8), type = "classes", legend = FALSE)
Extraction of apices positions and attributes
A data.frame
of detected apices located in the plot mask
is then extracted with the function
lidaRtRee::tree_extraction()
. Crown polygons can be
vectorized from the segments. Attributes are :
-
id
: apex id -
x
: easting coordinate of apex -
y
: northing coordinate of apex -
h
: height of apex -
dom_radius
: distance of apex to nearest higher pixel of CHM -
s
: crown surface -
v
: crown volume -
sp
(if plot mask is provided): crown surface inside plot -
vp
(if plot mask is provided): crown volume inside plot -
crown
(optional): 2D crown polygon in WKT format
# tree extraction only inside plot mask for subsequent comparison
apices <- lidaRtRee::tree_extraction(segms, r_mask = mask_plot, crown = TRUE)
# convert WKT field to polygons
crowns <- sf::st_as_sf(sf::st_drop_geometry(apices), wkt = "crown")
# remove WKT field from apices
apices <- apices[, -which(names(apices)=="crown")]
head(apices, n = 3L)
## Simple feature collection with 3 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 974347.2 ymin: 6581649 xmax: 974380.8 ymax: 6581672
## Projected CRS: RGF93 v1 / Lambert-93
## id x y h dom_radius s v sp vp
## 1 3 974349.2 6581672 24.99 6.0 75.50 1191.5775 71.75 NaN
## 2 6 974347.2 6581649 21.03 4.5 50.25 724.0775 50.25 724.0775
## 3 7 974380.8 6581662 15.91 6.0 32.25 385.1500 32.25 385.1500
## geometry
## 1 POINT (974349.2 6581672)
## 2 POINT (974347.2 6581649)
## 3 POINT (974380.8 6581662)
#
# display initial image
terra::plot(chm, col = gray(seq(0, 1, 1 / 255)),
main = "CHM and detected positions")
# display segments border
terra::plot(sf::st_geometry(crowns), border = "white", add = T, col = NA)
# display plot mask
terra::plot(mask_plot_v, border = "red", add = T)
# display detected apices
plot(apices["h"], col = "blue", cex = apices$h / 20, pch = 2, add = TRUE)
Detection evaluation
Tree matching
To assess detection accuracy, reference (field) trees should be
linked to detected apices. Despite the possibility of error, automated
matching has the advantage of making the comparison of results from
different algorithms and settings reproducible and objective. The
algorithm presented below is based on the 3D distance between detected
treetops and inventory positions and heights (Monnet et al. 2010). It returns an object with
the pairs of matched reference trees and detected apices. The function
lidaRtRee::plot_matched()
then plots the results.
# match detected apices with field trees based on relative distance of apices
matched <- lidaRtRee::tree_matching(
tree_inventory_chablais3[, c("x", "y", "h")],
apices[, c("x", "y", "h")]
)
# display matching results
lidaRtRee::plot_matched(
tree_inventory_chablais3[, c("x", "y", "h")],
apices[, c("x", "y", "h")], matched, chm, mask_plot_v
)
Detection accuracy
Detection accuracy is evaluated thanks to the number of correct
detections (53), false detections (4) and omissions (57). In
heterogeneous stands, detection accuracy is height-dependent, it is
informative to display those categories on a height histogram with
lidaRtRee::hist_detection()
.
# height histogram of detections
detection_stats <- lidaRtRee::hist_detection(
tree_inventory_chablais3[, c("x", "y", "h")],
apices[, c("x", "y", "h")], matched
)
Height estimation accuracy
For true detections, estimated heights can be compared to
field-measured heights. Here, the root mean square error is 0.8m, while
the bias is -0.2m. The linear regression is displayed hereafter
(lidaRtRee::height_regression()
).
# linear regression between reference height and estimated height
height_reg <- lidaRtRee::height_regression(
tree_inventory_chablais3[, c("x", "y", "h")],
apices[, c("x", "y", "h")],
matched,
species = tree_inventory_chablais3$s
)
Species classification
The previous steps also output a segmentation of the CHM, i.e. each detected apex is associated to a crown segment in 2D. Next we aim at checking for relationships between the point cloud contained in the segment, and the species of the largest field tree contained in the segment. The processing steps are:
- calculate statistical indices describing the point cloud for each segment,
- extract the highest reference tree from each segment,
- exploratory analysis of the relationship between the indices and the species.
Points in segments
Before computation of point cloud metrics in each segment, the whole point cloud is normalized to convert point altitude to height above ground. Points are then labeled with the id of the segment they belong to. A list of LAS objects corresponding to the segments is then created.
# normalize point cloud
lasn <- lidR::normalize_height(las_chablais3, lidR::tin())
# add segment id in LAS object
lasn <- lidR::merge_spatial(lasn, segms$segments_id, "seg_id")
# put all seg_id values in ordered list
list_seg_id <- sort(unique(lasn$seg_id))
# set names of list equal to values
names(list_seg_id) <- list_seg_id
# extract point cloud for each segment id in a list
las_l <- lapply(list_seg_id, function(x) {lidR::filter_poi(lasn, seg_id == x)})
Metrics computation
Basic point cloud metrics are computed in each segment (maximum,
mean, standard deviation of height, mean and standard deviation of
intensity), with the function lidaRtRee::clouds_metrics()
which applies a function to all point clouds in a list, by passing it to
lidR::cloud_metrics()
. Please refer to the help of this
last function for the expected syntax of the user-defined function to
compute metrics.
# compute basic las metrics in each segment
metrics <- lidaRtRee::clouds_metrics(las_l, func = ~ list(
maxZ = max(Z), meanZ = mean(Z),
sdZ = sd(Z), meanI = mean(Intensity),
sdI = sd(Intensity)
))
# add segment id attribute
metrics$seg_id <- row.names(metrics)
head(metrics, n = 3L)
## maxZ meanZ sdZ meanI sdI seg_id
## 0 8.68 0.6028999 1.299651 84.81373 61.59432 0
## 1 26.43 14.4630361 6.079366 46.28273 38.56279 1
## 3 24.97 11.9435605 6.580607 43.45585 34.96201 3
Merge with reference trees and detected apices
Computed metrics are merged with reference trees and detected apices, thanks to the segment id.
# associate each reference tree with the segment that contains its trunk.
dummy <- terra::extract(
segms$segments_id,
sf::st_coordinates(tree_inventory_chablais3)
)
tree_inventory_chablais3$seg_id <- dummy$segments_id
# create new data.frame by merging metrics and inventoried trees
# (without geometry)
# based on segment id
tree_metrics <- base::merge(sf::st_drop_geometry(tree_inventory_chablais3),
metrics)
# remove non-tree segment
tree_metrics <- tree_metrics[tree_metrics$seg_id != 0, ]
# add metrics to extracted apices data.frame
apices <- base::merge(apices, metrics, by.x = "id", by.y = "seg_id", all.x = T)
Plotting the reference trees with the mean intensity of lidar points in the segments shows that when they are dominant, broadleaved trees tend to have higher mean intensity than coniferous trees.
# create raster of segment' mean intensity
r_mean_intensity_segm <- segms$segments_id
# match segment id with id in metrics data.frame
dummy <- match(terra::values(r_mean_intensity_segm), apices$id)
# replace segment id with corresponding mean intensity
terra::values(r_mean_intensity_segm) <- apices$meanI[dummy]
# display tree inventory with mean intensity in segment background
terra::plot(r_mean_intensity_segm, main = "Mean intensity in segment")
# display tree inventory
lidaRtRee::plot_tree_inventory(tree_inventory_chablais3[, c("x", "y")],
tree_inventory_chablais3$h,
species = as.character(tree_inventory_chablais3$s), add = T
)
Exploratory analysis
A boxplot of mean intensity in segments per species shows that mean intensity distribution is different between species. The analysis can be run by considering only the highest inventoried trees in each segment, as smaller trees are likely to be suppressed and have smaller contribution to the point cloud.
# create new variable orderer by segment id then decreasing height
tree_metrics_h <- tree_metrics[order(tree_metrics$seg_id, -tree_metrics$h), ]
i <- 2
# leave only first (highest) tree per segment id
while (i < nrow(tree_metrics_h)) {
if (tree_metrics_h$seg_id[i] == tree_metrics_h$seg_id[i + 1]) {
tree_metrics_h <- tree_metrics_h[-(i + 1), ]
} else {
i <- i + 1
}
}
tree_metrics_h$s <- factor(tree_metrics_h$s)
par(mfrow = c(1, 2))
boxplot(meanI ~ s,
data = tree_metrics[, c("s", "maxZ", "meanZ", "sdZ", "meanI", "sdI")],
ylab = "Mean intensity in segment", xlab = "Specie",
main = "All inventoried trees", las = 2, varwidth = TRUE
)
boxplot(meanI ~ s,
data = tree_metrics_h, ylab = "Mean intensity in segment",
xlab = "Specie", main = "Highest inventoried tree in segment", las = 2,
varwidth = TRUE
)
A linear discriminant analysis shows that it might be possible to discriminate between spruce, fir and the deciduous species, thanks to a combination of height and intensity variables.
# principal component analysis
pca <-
ade4::dudi.pca(tree_metrics_h[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")],
scannf = F)
# linear discriminant analysis
lda <- ade4::discrimin(pca, tree_metrics_h$s, scannf = F, nf = 2)
plot(lda)
The following line creates a new factor variable Groups
with three groups: the species PIAB, ABAL and all others together.
tree_metrics_h$Groups <- tree_metrics_h$s
levels(tree_metrics_h$Groups)[!is.element(levels(tree_metrics_h$Groups),
c("ABAL", "PIAB"))] <- "Other"
A linear discriminant analysis is performed on this new variable.
Classification
The function MASS::lda
is used for prediction purposes.
First, the cross-validation option (CV = TRUE
) is chosen in
order to produce a confusion matrix in a prediction case.
lda_MASS <-
MASS::lda(tree_metrics_h[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")],
tree_metrics_h$Groups, CV = TRUE)
# confusion matrix
matrix_confusion <- table(tree_metrics_h$Groups, lda_MASS$class)
matrix_confusion
##
## ABAL Other PIAB
## ABAL 8 4 4
## Other 2 19 0
## PIAB 5 1 9
# percentage of good classification
round(sum(diag(matrix_confusion)) / sum(matrix_confusion) * 100, 1)
## [1] 69.2
# confidence interval of the percentage
binom.test(sum(diag(matrix_confusion)), sum(matrix_confusion))$conf.int
## [1] 0.5489760 0.8128266
## attr(,"conf.level")
## [1] 0.95
Then the model is fitted and applied to all segments of the map.
# build model
lda_MASS <-
MASS::lda(tree_metrics_h[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")],
tree_metrics_h$Groups)
# apply model to metrics
metrics$predicted_s <-
predict(lda_MASS, metrics[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")])$class
# apply model to trees
tree_metrics_h$predicted_s <- predict(lda_MASS,
tree_metrics_h[, c("maxZ", "meanZ", "sdZ",
"meanI", "sdI")])$class
To display the classification results, an image of the classified segments is created and the highest field trees in each segment are displayed on top of it. Results are correct when:
- purple dots are on purple segments (correct ABAL classification),
- blue dots are on blue segments (correct PIAB classification),
- other colors are on green segments (correct “others” classification).
# create image of predicted species
species <- terra::deepcopy(segms$segments_id)
# replace segment id by id of predicted species in image
terra::values(species) <-
as.numeric(metrics$predicted_s)[match(terra::values(segms$segments_id),
metrics$seg_id)]
# remove ground segment
species[segms$segments_id == 0] <- NA
# build raster attribute table rat
rat <- data.frame(id = 1:length(levels(metrics$predicted_s)),
Species = levels(metrics$predicted_s))
# retrieve reference colors
rat$col <- lidaRtRee::species_color()[rat$Species, "col"]
# set NA color to green
rat$col[is.na(rat$col)] <- "green"
# convert to factor (add RAT in SpatRaster)
levels(species) <- rat
# display results
terra::plot(species, col = rat$col)
terra::plot(sf::st_geometry(crowns), add = TRUE, border = "black")
lidaRtRee::plot_tree_inventory(tree_metrics_h
[, c("x", "y")],
tree_metrics_h$h,
bg = lidaRtRee::species_color()[as.character(tree_metrics_h$s), "col"],
col = "black",
pch = 21, add = TRUE
)
Display point cloud
The point cloud can be displayed colored by segment, with poles at the location of inventoried trees.
rgl::par3d(mouseMode = "trackball") # parameters for interaction with mouse
rgl::bg3d(color = "white")
# select segment points and offset them to avoid truncated coordinates in 3d plot
points.seg <- lasn[which(lasn$seg_id != 0), c("X", "Y", "Z", "seg_id")]
points.seg$X <- points.seg$X - 974300
points.seg$Y <- points.seg$Y - 6581600
# plot point cloud
rgl::plot3d(points.seg@data[, c("X", "Y", "Z")],
col = points.seg$seg_id %% 10 + 1, aspect = FALSE)
#
# add inventoried trees
tree_inventory_chablais3$z <- 0
for (i in 1:nrow(tree_inventory_chablais3))
{
rgl::lines3d(
rbind(
tree_inventory_chablais3$x[i] - 974300,
tree_inventory_chablais3$x[i] - 974300
),
rbind(
tree_inventory_chablais3$y[i] - 6581600,
tree_inventory_chablais3$y[i] - 6581600
),
rbind(
tree_inventory_chablais3$z[i],
tree_inventory_chablais3$z[i] + tree_inventory_chablais3$h[i]
)
)
}
Batch processing
The function lidaRtRee::tree_detection()
encompasses
both the segmentation and extraction steps, and can take as input:
- a Canopy Height Model provided as
SpatRaster
object, - a LAS object,
- a LAS-catalog object.
If the LAS object or files are not normalized, the function has an
argument (normalization = TRUE
) to do it on-the-fly.
A region of interest (ROI) can be specified in order to output only trees contained in the ROI, as in the example below which reproduces the results of previous paragraphs.
# perform tree detection
apices <- lidaRtRee::tree_detection(chm, ROI = mask_plot_v)
# display initial image
terra::plot(chm, col = gray(seq(0, 1, 1 / 255)),
main = "CHM and detected positions")
# display plot mask
terra::plot(mask_plot_v, border = "red", add = T)
# display detected apices
plot(apices["h"], col = "blue", cex = apices$h / 20, pch = 2, add = TRUE)
In case of a LAS-catalog object, the lidR
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.
-
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 tree which treetop is located at the border of a chunk has its full crown within the buffer area.
Chunks can be processed with parallel computing, within limits of the cluster’s RAM and number of cores. If a ROI is provided, only files intersecting it are processed.
The steps for processing a batch of las/laz files are :
- build catalog of files and set catalog processing options
- provide segmentation parameters
- provide output parameters
- set cluster options for parallel processing
- apply
lidaRtRee::tree_detection()
function to catalog - the Canopy Height Model is not an output of the function, it can be
computed separately by applying
lidR::rasterize_canopy()
to the catalog (if LAS files are normalized; more information on this step in the preprocessing tutorial).
Be aware that in case tree segments are vectorized, some obtained crowns might overlap. The segmentation algorithm is not be deterministic and borders might not be consistent when adjacent polygons are pasted from different chunks.
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 = 70,
chunk_buffer = 10)
## Be careful, a chunk size smaller than 250 is likely to be irrelevant.
# set coordinate system
lidR::projection(cata) <- 2154
# proceed with processing in case one chunk returns errors
cata@processing_options$stop_early <- FALSE
#
# 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
apices <- lidaRtRee::tree_detection(cata, res = resolution, crown = TRUE)
# create crown polygons object
crowns <- sf::st_as_sf(sf::st_drop_geometry(apices), wkt = "crown",
crs = sf::st_crs(apices))
# compute canopy height model (apply lidR::rasterize_canopy to catalog,
# with same resolution)
chm <- lidR::rasterize_canopy(cata, res = resolution)
The following image displays the results for the whole area.
# threshold outsiders in chm
chm[chm > 40] <- 40
chm[chm < 0] <- 0
# display chm
terra::plot(chm,
main = "Canopy Height Model and segments"
)
# display segments border
plot(sf::st_geometry(crowns), border = "black", add = T)
# add apices
plot(sf::st_geometry(apices), cex = apices$h / 40, add = TRUE, pch = 2)
The following lines save outputs to files (and overwrite existing files).
# merged chm
terra::writeRaster(chm, file = "chm.tif", overwrite = TRUE)
# apices
sf::st_write(apices, "apices.gpkg", delete_dsn = TRUE)
# crowns
sf::st_write(crowns, "crowns.gpkg", delete_dsn = TRUE)