Licence: GNU GPLv3 / Source page
Many thanks to Pascal Obstétar for checking code and improvement suggestions.
The R code below presents a forest structure metrics computation workflow from Airborne Laser Scanning (ALS) data. Workflow is based on functions from packages lidaRtRee and lidR, packages vegan and foreach are also required. Metrics are computed for each cell of a grid defined by a resolution. Those metrics are designed to describe the 3D structure of forest. Different types of metrics are computed:
The forest structure metrics derived from airborne laser scanning can be used for habitat suitability modelling and mapping. This workflow has been applied to compute the metrics used in the modeling and mapping of the habitat of Capercaillie (Tetrao urogallus)(Glad et al. (2020)). For more information about tree segmentation and gaps detection, please refer to the corresponding tutorials.
The workflow processes normalized point clouds provided as las/laz tiles of rectangular extent. Parallelization is used for faster processing, packages foreach, future and doFuture are used. A buffer is loaded around each tile to prevent border effects in tree segmentation and CHM processing.
Load foreach package and set the number of cores to use for parallel computing
library(foreach)
# create parallel frontend, specify to use two parallel sessions
doFuture::registerDoFuture()
future::plan("multisession", workers = 2L)
Numerous parameters have to be set for processing.
# output metrics map resolution (m)
resolution <- 10
# canopy height model resolution (m)
res_chm <- 0.5
# buffer size (m) for tile processing (20 m is better for gaps metrics, 10 m is
# enough for tree metrics)
buffer_size <- 20
# height threshold (m) for high points removal (points above this threshold are
# considered as outliers)
points_max_h <- 60
# points classes to retain for analysis (vegetation + ground)
class_points <- c(0, 1, 2, 3, 4, 5)
# ground class
class_ground <- 2
# Gaussian filter sigma values for multi-scale smoothing
sigma <- c(0, 1, 2, 4, 8, 16)
# convert vector of sigma values in list
sigma_l <- as.list(sigma)
# height breaks for penetration ratio and density
breaks_h <- c(-Inf, 0, 0.5, 1, 2, 5, 10, 20, 30, 60, Inf)
# percentiles of height distribution
percent <- c(0.1, 0.25, 0.5, 0.75, 0.9)
# surface breaks for gap size (m2)
breaks_gap_surface <- c(4, 16, 64, 256, 1024, 4096, Inf)
# gap surface names
n_breaks_gap <- gsub("-", "", paste0("G_s", breaks_gap_surface[c(-length(breaks_gap_surface))],
"to", breaks_gap_surface[c(-1)]))
# height bin names
n_breaks_h <- gsub("-", "", paste0("nb_H", breaks_h[c(-length(breaks_h))], "to",
breaks_h[c(-1)]))
#
The first step is to create a catalog of LAS files (should be normalized, non-overlapping rectangular tiles). Preferably, tiles should be aligned on a multiple of resolution, and points should not lie on the northern or eastern border when such borders are common with adjacent tiles.
# create catalog of LAS files
cata <- lidR::readALSLAScatalog("../data/forest.structure.metrics")
# set coordinate system
lidR::projection(cata) <- 2154
# disable display of catalog processing
lidR::opt_progress(cata) <- FALSE
# display
sp::plot(cata, main = "Bounding box of las files")
The tile is loaded, plus a buffer on adjacent tiles. Buffer points have a column “buffer” filled with 1.
# tile to process
i <- 1
# tile extent
b_box <- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
a <- lidR::clip_rectangle(cata, b_box[1] - buffer_size, b_box[2] - buffer_size, b_box[3] +
buffer_size, b_box[4] + buffer_size)
# add 'buffer' flag to points in buffer with TRUE value in this new attribute
a <- lidR::add_attribute(a, a$X < b_box[1] | a$Y < b_box[2] | a$X >= b_box[3] | a$Y >=
b_box[4], "buffer")
a
## class : LAS (v1.1 format 1)
## memory : 26.7 Mb
## extent : 944800, 944920, 6638800, 6638900 (xmin, xmax, ymin, ymax)
## coord. ref. : RGF93 / Lambert-93
## area : 11996.4 m²
## points : 333.3 thousand points
## density : 27.78 points/m²
Only points from desired classes are retained. In case some mis-classified, high or low points remain in the data set, the point cloud is filtered and negative heights are replaced by 0.
# remove unwanted point classes, and points higher than height threshold
a <- lidR::filter_poi(a, is.element(Classification, class_points) & Z <= points_max_h)
# set negative heights to 0
a@data$Z[a@data$Z < 0] <- 0
#
summary(a@data)
## X Y Z gpstime
## Min. :944800 Min. :6638800 Min. : 0.00 Min. :1.52e+08
## 1st Qu.:944828 1st Qu.:6638825 1st Qu.: 2.33 1st Qu.:1.52e+08
## Median :944860 Median :6638850 Median : 9.94 Median :1.52e+08
## Mean :944860 Mean :6638849 Mean :10.90 Mean :1.52e+08
## 3rd Qu.:944891 3rd Qu.:6638872 3rd Qu.:18.14 3rd Qu.:1.52e+08
## Max. :944920 Max. :6638900 Max. :37.45 Max. :1.52e+08
## Intensity ReturnNumber NumberOfReturns ScanDirectionFlag
## Min. : 2.0 Min. :1.000 Min. :1.000 Min. :0
## 1st Qu.: 27.0 1st Qu.:1.000 1st Qu.:2.000 1st Qu.:0
## Median : 67.0 Median :1.000 Median :2.000 Median :0
## Mean :101.7 Mean :1.713 Mean :2.435 Mean :0
## 3rd Qu.:135.0 3rd Qu.:2.000 3rd Qu.:3.000 3rd Qu.:0
## Max. :821.0 Max. :7.000 Max. :7.000 Max. :0
## EdgeOfFlightline Classification Synthetic_flag Keypoint_flag
## Min. :0 Min. :0.000 Mode :logical Mode :logical
## 1st Qu.:0 1st Qu.:5.000 FALSE:333229 FALSE:333229
## Median :0 Median :5.000
## Mean :0 Mean :4.122
## 3rd Qu.:0 3rd Qu.:5.000
## Max. :0 Max. :5.000
## Withheld_flag ScanAngleRank UserData PointSourceID buffer
## Mode :logical Min. :-2.00 Min. :0 Min. :0 Mode :logical
## FALSE:333229 1st Qu.: 1.00 1st Qu.:0 1st Qu.:0 FALSE:278011
## Median :86.00 Median :0 Median :0 TRUE :55218
## Mean :44.68 Mean :0 Mean :0
## 3rd Qu.:87.00 3rd Qu.:0 3rd Qu.:0
## Max. :87.00 Max. :0 Max. :0
The next step is to compute the canopy height model (CHM). It will be used to derive:
The CHM is computed and NA values are replaced by 0. A check is performed to make sure low or high points are not present.
# compute chm
chm <- lidaRtRee::points2DSM(a, res = res_chm)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
# replace NA, low and high values
chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
# display CHM
raster::plot(chm, asp = 1, main = "Canopy height model")
The CHM is smoothed with a Gaussian filter with different sigma values. Smoothed results are stored in a list and then converted to a raster stack.
# for each value in list of sigma, apply filtering to chm and store result in
# list
st <- lapply(sigma_l, FUN = function(x) {
lidaRtRee::dem_filtering(chm, nl_filter = "Closing", nl_size = 5, sigmap = x)$smoothed_image
})
# convert to raster
st <- raster::stack(st)
# display
raster::plot(st)
The raster stack is then converted to points. Points outside the tile extent are removed and summary statistics (mean and standard deviation) of smoothed CHM heights are computed for each pixel at the output metrics resolution. Canopy cover in different height layers are also computed for the initial CHM after applying a non-linear filter designed to fill holes.
# crop to bbox
st <- raster::crop(st, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
# multiplying factor to compute proportion
mf <- 100/(resolution/res_chm)^2
# modify layer names
names(st) <- paste0("layer.", 1:6)
# compute raster metrics
metrics_2dchm <- lidaRtRee::raster_metrics(st, res = resolution, fun = function(x) {
c(sd(x$layer.1), sd(x$layer.2), sd(x$layer.3), sd(x$layer.4), sd(x$layer.5),
sd(x$layer.6), mean(x$layer.1), sum(x$layer.1 < 0.5) * mf, sum(x$layer.1 <
1) * mf, sum(x$layer.1 > 5) * mf, sum(x$layer.1 > 10) * mf, sum(x$layer.1 >
20) * mf, (sum(x$layer.1 < 5) - sum(x$layer.1 < 1)) * mf, (sum(x$layer.1 <
5) - sum(x$layer.1 < 2)) * mf)
}, output = "raster")
names(metrics_2dchm) <- c("CHM0_sd", "CHM1_sd", "CHM2_sd", "CHM4_sd", "CHM8_sd",
"CHM16_sd", "CHM_mean", "CHM_PercInf0_5", "CHM_PercInf1", "CHM_PercSup5", "CHM_PercSup10",
"CHM_PercSup20", "CHM_Perc1_5", "CHM_Perc2_5")
#
metrics_2dchm
## class : RasterBrick
## dimensions : 10, 10, 100, 14 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 944800, 944900, 6638800, 6638900 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## source : memory
## names : CHM0_sd, CHM1_sd, CHM2_sd, CHM4_sd, CHM8_sd, CHM16_sd, CHM_mean, CHM_PercInf0_5, CHM_PercInf1, CHM_PercSup5, CHM_PercSup10, CHM_PercSup20, CHM_Perc1_5, CHM_Perc2_5
## min values : 0.09293137, 0.08238094, 0.12682748, 0.21107244, 0.42359988, 0.17968139, 0.08870000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000
## max values : 11.990108, 10.747519, 9.463671, 7.492695, 4.271958, 2.322395, 28.713500, 99.250000, 100.000000, 100.000000, 100.000000, 100.000000, 55.250000, 41.250000
Some outputs are displayed hereafter. On the first line are the standard deviation of CHM smoothed with different sigma values. On the second line are the CHM mean and percentages of CHM values below 1 m and above 10 m.
raster::plot(metrics_2dchm[[c("CHM0_sd", "CHM2_sd", "CHM8_sd", "CHM_mean", "CHM_PercInf1",
"CHM_PercSup10")]])
Gaps are computed with the function gap_detection.
# compute gaps
gaps <- lidaRtRee::gap_detection(chm, ratio = 2, gap_max_height = 1, min_gap_surface = min(breaks_gap_surface),
closing_height_bin = 1, nl_filter = "Median", nl_size = 3, gap_reconstruct = TRUE)
# display maps
par(mfrow = c(1, 3))
raster::plot(chm, main = "CHM")
raster::plot(log(gaps$gap_surface), main = "Log of gap surface")
# display gaps
dummy <- gaps$gap_id
dummy[dummy == 0] <- NA
raster::plot(dummy%%8, asp = 1, main = "Gap id (random colors)", col = rainbow(8),
legend = FALSE)
Summary statistics are then computed: pixel surface occupied by gaps in different size classes. Results are first cropped to the tile extent.
# crop results to tile size
gaps_surface <- raster::crop(gaps$gap_surface, raster::extent(b_box[1], b_box[3],
b_box[2], b_box[4]))
# compute gap statistics at final metrics resolution, in case gaps exist
if (!all(is.na(raster::values(gaps_surface)))) {
metrics_gaps <- lidaRtRee::raster_metrics(gaps_surface, res = resolution, fun = function(x) {
hist(x$layer, breaks = breaks_gap_surface, plot = F)$counts * (res_chm/resolution)^2
}, output = "raster")
# compute total gap proportion
metrics_gaps$sum <- raster::stackApply(metrics_gaps, rep(1, length(names(metrics_gaps))),
fun = sum)
# if gaps are present
names(metrics_gaps) <- c(n_breaks_gap, paste("G_s", min(breaks_gap_surface),
"toInf", sep = ""))
# replace possible NA values by 0 (NA values are present in lines of the
# target raster where no values >0 are present in the input raster)
metrics_gaps[is.na(metrics_gaps)] <- 0
} else {
metrics_gaps <- NULL
}
#
metrics_gaps
## class : RasterBrick
## dimensions : 10, 10, 100, 7 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 944800, 944900, 6638800, 6638900 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## source : memory
## names : G_s4to16, G_s16to64, G_s64to256, G_s256to1024, G_s1024to4096, G_s4096toInf, G_s4toInf
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 0.0750, 0.3200, 0.5175, 0.0000, 1.0000, 0.0000, 1.0000
Edges are detected with the function edge_detection as the outer envelope of the previously delineated gaps.
# Perform edge detection by erosion
edges <- lidaRtRee::edge_detection(gaps$gap_id > 0)
# display maps
par(mfrow = c(1, 3))
raster::plot(chm, main = "CHM")
raster::plot(gaps$gap_id > 0, main = "Gaps", legend = FALSE)
raster::plot(edges, main = "Edges", legend = FALSE)
The percentage of surface occupied by edges is then computed as summary statistic. Results are first cropped to the tile extent.
# crop results to tile size
edges <- raster::crop(edges, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
# compute gap statistics at final metrics resolution, in case gaps exist
if (!all(is.na(raster::values(edges)))) {
metrics_edges <- lidaRtRee::raster_metrics(edges, res = resolution, fun = function(x) {
sum(x$layer) * (res_chm/resolution)^2
}, output = "raster")
names(metrics_edges) <- "edges.proportion"
# replace possible NA values by 0 (NA values are present in lines of the
# target raster where no values >0 are present in the input raster)
metrics_edges[is.na(metrics_edges)] <- 0
} else {
metrics_edges <- NULL
}
#
metrics_edges
## class : RasterLayer
## dimensions : 10, 10, 100 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 944800, 944900, 6638800, 6638900 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## source : memory
## names : edges.proportion
## values : 0, 0.1225 (min, max)
Some outputs are displayed hereafter. The proportion of surface occupied by gaps of various sizes is depicted as well as the proportion of surface occupied by edges.
par(mfrow = c(1, 3))
# raster::plot(metrics_gaps[['G_s4to16']], main='Prop of gaps 4-16m2')
raster::plot(metrics_gaps[["G_s64to256"]], main = "Prop of gaps 64-256m2")
raster::plot(metrics_gaps[["G_s4toInf"]], main = "Prop of gaps >4m2")
raster::plot(metrics_edges, main = "Proportion of edges")
Tree tops are detected with the function treeSegmentation and then extracted with treeExtraction.
# tree top detection (default parameters)
segms <- lidaRtRee::tree_segmentation(chm, hmin = 5)
# extraction to data.frame
trees <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
# display maps
par(mfrow = c(1, 2))
raster::plot(chm, main = "CHM and tree tops")
points(trees$x, trees$y, pch = 4, cex = 0.3)
# display segments, except ground segment
dummy <- segms$segments_id
dummy[dummy == 0] <- NA
raster::plot(dummy%%8, asp = 1, main = "Segments (random colors)", col = rainbow(8),
legend = FALSE)
points(trees$x, trees$y, pch = 4, cex = 0.3)
Results are cropped to the tile extent. Summary statistics from std_tree_metrics are then computed for each pixel. Canopy cover in trees and mean canopy height in trees are computed in a second step because they can not be computed from the data of the extracted trees which crowns may span several output pixels. They are calculated from the images obtained during the previous tree segmentation.
# remove trees outside of tile
trees <- trees[trees$x >= b_box[1] & trees$x < b_box[3] & trees$y >= b_box[2] & trees$y <
b_box[4], ]
# compute raster metrics
metrics_trees <- lidaRtRee::raster_metrics(trees[, -1], res = resolution, fun = function(x) {
lidaRtRee::std_tree_metrics(x, resolution^2/10000)
}, output = "raster")
# remove NAs and NaNs
metrics_trees[!is.finite(metrics_trees)] <- 0
# compute canopy cover in trees and canopy mean height in trees in region of
# interest, because it is not in previous step.
r_tree_chm <- segms$filled_dem
# set chm to NA in non segment area
r_tree_chm[segms$segments_id == 0] <- NA
# compute raster metrics
dummy <- lidaRtRee::raster_metrics(raster::crop(r_tree_chm, raster::extent(b_box[1],
b_box[3], b_box[2], b_box[4])), res = resolution, fun = function(x) {
c(sum(!is.na(x$filled_dem))/(resolution/res_chm)^2, mean(x$filled_dem, na.rm = T))
}, output = "raster")
names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
#
dummy <- raster::extend(dummy, metrics_trees)
metrics_trees <- raster::addLayer(metrics_trees, dummy)
#
metrics_trees
## class : RasterStack
## dimensions : 10, 10, 100, 13 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 944800, 944900, 6638800, 6638900 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## names : Tree_meanH, Tree_sdH, Tree_giniH, Tree_density, TreeInf10_density, TreeSup10_density, TreeSup20_density, TreeSup30_density, Tree_meanCrownSurface, Tree_meanCrownVolume, TreeCanopy_meanH, TreeCanopy_cover_in_plot, TreeCanopy_meanH_in_plot
## min values : 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 6.113102
## max values : 37.1800000, 15.1957247, 0.2050181, 500.0000000, 200.0000000, 500.0000000, 400.0000000, 200.0000000, 159.0000000, 3136.9250000, 28.4389711, 1.0000000, 28.7135000
Some outputs are displayed hereafter.
par(mfrow = c(1, 3))
raster::plot(metrics_trees$Tree_meanH, main = "Mean (detected) tree height (m)")
points(trees$x, trees$y, cex = trees$h/40)
raster::plot(metrics_trees$TreeSup20_density, main = "Density of (detected) trees > 20m (/ha)")
raster::plot(metrics_trees$Tree_meanCrownVolume, main = "Mean crown volume of detected trees (m3)")
Buffer points are first removed from the point cloud, then metrics are computed from all points and from vegetation points only.
# remove buffer points
a <- lidR::filter_poi(a, buffer == 0)
# all points metrics
metrics_1d <- lidR::grid_metrics(a, as.list(c(as.vector(quantile(Z, probs = percent)),
sum(Classification == 2), mean(Intensity[ReturnNumber == 1]))), res = resolution)
names(metrics_1d) <- c(paste("Hp", percent * 100, sep = ""), "nb_ground", "Imean_1stpulse")
# vegetation-only metrics
a <- lidR::filter_poi(a, Classification != class_ground)
dummy <- lidR::grid_metrics(a, as.list(c(mean(Z), max(Z), sd(Z), length(Z), hist(Z,
breaks = breaks_h, right = F, plot = F)$counts)), res = resolution)
names(dummy) <- c("Hmean", "Hmax", "Hsd", "nb_veg", n_breaks_h)
# merge rasterstacks
dummy <- raster::extend(dummy, metrics_1d)
metrics_1d <- raster::addLayer(metrics_1d, dummy)
rm(dummy)
# crop to tile extent
metrics_1d <- raster::crop(metrics_1d, raster::extent(b_box[1], b_box[3], b_box[2],
b_box[4]))
Based on those metrics, additional metrics are computed (Simpson index, relative density in height bins and penetration ratio in height bins)
# Simpson index
metrics_1d$Hsimpson <- raster::stackApply(metrics_1d[[n_breaks_h[c(-1, -length(n_breaks_h))]]],
1, function(x, ...) {
vegan::diversity(x, index = "simpson")
})
# Relative density
for (i in n_breaks_h[c(-1, -length(n_breaks_h))]) {
metrics_1d[[gsub("nb", "relativeDensity", i)]] <- metrics_1d[[i]]/(metrics_1d$nb_veg +
metrics_1d$nb_ground)
}
# Penetration ratio cumulative sum
cumu_sum <- metrics_1d[["nb_ground"]]
for (i in n_breaks_h) {
cumu_sum <- raster::addLayer(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics_1d[[i]])
}
names(cumu_sum) <- c("nb_ground", n_breaks_h)
# compute interception ratio of each layer
intercep_ratio <- cumu_sum[[-1]]
for (i in 2:dim(cumu_sum)[3]) {
intercep_ratio[[i - 1]] <- 1 - cumu_sum[[i - 1]]/cumu_sum[[i]]
}
names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
# merge rasterstacks
metrics_1d <- raster::addLayer(metrics_1d, intercep_ratio)
# rm(cumu_sum, intercep_ratio)
metrics_1d
## class : RasterStack
## dimensions : 10, 10, 100, 40 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 944800, 944900, 6638800, 6638900 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## names : Hp10, Hp25, Hp50, Hp75, Hp90, nb_ground, Imean_1stpulse, Hmean, Hmax, Hsd, nb_veg, nb_HInfto0, nb_H0to0.5, nb_H0.5to1, nb_H1to2, ...
## min values : 0.000000e+00, 0.000000e+00, 1.000000e-02, 3.000000e-02, 6.000000e-02, 5.200000e+01, 7.795880e+01, 5.387324e-02, 7.300000e-01, 8.174446e-02, 6.760000e+02, 0.000000e+00, 1.400000e+01, 0.000000e+00, 0.000000e+00, ...
## max values : 1.762000e+01, 2.269000e+01, 2.572500e+01, 2.791000e+01, 3.167000e+01, 9.820000e+02, 3.193222e+02, 2.488985e+01, 3.718000e+01, 1.161968e+01, 3.850000e+03, 0.000000e+00, 7.060000e+02, 1.510000e+02, 2.240000e+02, ...
Some outputs are displayed hereafter.
# display results
par(mfrow = c(1, 3))
raster::plot(metrics_1d$Imean_1stpulse, main = "Mean intensity (1st return)")
raster::plot(metrics_1d$Hp50, main = "Percentile 50 of heights")
raster::plot(metrics_1d$intercepRatio_H2to5, main = "Interception ratio 2-5 m")
The following code uses parallel processing to handle multiple files of a catalog, and arranges all metrics in a raster stack. Code in the “parameters” section has to be run beforehand.
# processing by tile
metrics <- foreach::foreach(i = 1:nrow(cata), .errorhandling = "remove") %dopar%
{
# tile extent
b_box <- as.numeric(cata@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")])
# load tile extent plus buffer
a <- try(lidR::clip_rectangle(cata, b_box[1] - buffer_size, b_box[2] - buffer_size,
b_box[3] + buffer_size, b_box[4] + buffer_size))
# check if points are successfully loaded
if (class(a) == "try-error") {
return(NULL)
}
# add 'buffer' flag to points in buffer with TRUE value in this new
# attribute
a <- lidR::add_attribute(a, a$X < b_box[1] | a$Y < b_box[2] | a$X >= b_box[3] |
a$Y >= b_box[4], "buffer")
# remove unwanted point classes, and points higher than height
# threshold
a <- lidR::filter_poi(a, is.element(Classification, class_points) & Z <=
points_max_h)
# check number of remaining points
if (length(a) == 0) {
return(NULL)
}
# set negative heights to 0
a@data$Z[a@data$Z < 0] <- 0
# compute chm
chm <- lidaRtRee::points2DSM(a, res = res_chm)
# replace NA, low and high values
chm[is.na(chm) | chm < 0 | chm > points_max_h] <- 0
#-----------------------
# compute 2D CHM metrics for each value in list of sigma, apply
# filtering to chm and store result in list
st <- lapply(sigma_l, FUN = function(x) {
lidaRtRee::dem_filtering(chm, nl_filter = "Closing", nl_size = 5, sigmap = x)[[2]]
})
# convert to raster
st <- raster::stack(st)
# crop to bbox
st <- raster::crop(st, raster::extent(b_box[1], b_box[3], b_box[2], b_box[4]))
# multiplying factor to compute proportion
mf <- 100/(resolution/res_chm)^2
# compute raster metrics
metrics_2dchm <- lidaRtRee::raster_metrics(st, res = resolution, fun = function(x) {
data.frame(CHM0.sd = sd(x$non_linear_image), CHM1.sd = sd(x$smoothed_image.1),
CHM2.sd = sd(x$smoothed_image.2), CHM4.sd = sd(x$smoothed_image.4),
CHM8.sd = sd(x$smoothed_image.8), CHM16.sd = sd(x$smoothed_image.16),
CHM.mean = mean(x$non_linear_image), CHM.PercInf0.5 = sum(x$non_linear_image <
0.5) * mf, CHM.PercInf1 = sum(x$non_linear_image < 1) * mf, CHM.PercSup5 = sum(x$non_linear_image >
5) * mf, CHM.PercSup10 = sum(x$non_linear_image > 10) * mf, CHM.PercSup20 = sum(x$non_linear_image >
20) * mf, CHM.Perc1_5 = (sum(x$non_linear_image < 5) - sum(x$non_linear_image <
1)) * mf, CHM.Perc2_5 = (sum(x$non_linear_image < 5) - sum(x$non_linear_image <
2)) * mf)
}, output = "raster")
# ------------------- compute gap metrics
gaps <- lidaRtRee::gap_detection(chm, ratio = 2, gap_max_height = 1, min_gap_surface = min(breaks_gap_surface),
closing_height_bin = 1, nl_filter = "Median", nl_size = 3, gap_reconstruct = TRUE)
# crop results to tile size
gaps_surface <- raster::crop(gaps$gap_surface, raster::extent(b_box[1], b_box[3],
b_box[2], b_box[4]))
# compute gap statistics at final metrics resolution, in case gaps
# exist
if (!all(is.na(raster::values(gaps_surface)))) {
metrics_gaps <- lidaRtRee::raster_metrics(gaps_surface, res = resolution,
fun = function(x) {
hist(x$layer, breaks = breaks_gap_surface, plot = F)$counts * (res_chm/resolution)^2
}, output = "raster")
# compute total gap proportion
metrics_gaps$sum <- raster::stackApply(metrics_gaps, rep(1, length(names(metrics_gaps))),
fun = sum)
# if gaps are present
names(metrics_gaps) <- c(n_breaks_gap, paste("G_s", min(breaks_gap_surface),
"toInf", sep = ""))
# replace possible NA values by 0 (NA values are present in lines
# of the target raster where no values >0 are present in the input
# raster)
metrics_gaps[is.na(metrics_gaps)] <- 0
} else {
metrics_gaps <- NULL
}
# -------------------- compute edge metrics Perform edge detection by
# erosion
edges <- lidaRtRee::edge_detection(gaps$gap_id > 0)
# crop results to tile size
edges <- raster::crop(edges, raster::extent(b_box[1], b_box[3], b_box[2],
b_box[4]))
# compute edges statistics at final metrics resolution, in case edges
# exist
if (!all(is.na(raster::values(edges)))) {
metrics_edges <- lidaRtRee::raster_metrics(edges, res = resolution, fun = function(x) {
sum(x$layer) * (res_chm/resolution)^2
}, output = "raster")
names(metrics_edges) <- "edges.proportion"
# replace possible NA values by 0 (NA values are present in lines
# of the target raster where no values >0 are present in the input
# raster)
metrics_edges[is.na(metrics_edges)] <- 0
} else {
metrics_edges <- NULL
}
# -------------------- compute tree metrics tree top detection (default
# parameters)
segms <- lidaRtRee::tree_segmentation(chm, hmin = 5)
# extraction to data.frame
trees <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima,
segms$segments_id)
# remove trees outside of tile
trees <- trees[trees$x >= b_box[1] & trees$x < b_box[3] & trees$y >= b_box[2] &
trees$y < b_box[4], ]
# compute raster metrics
metrics_trees <- lidaRtRee::raster_metrics(trees[, -1], res = resolution,
fun = function(x) {
lidaRtRee::std_tree_metrics(x, resolution^2/10000)
}, output = "raster")
# remove NAs and NaNs
metrics_trees[!is.finite(metrics_trees)] <- 0
# extend layer in case there are areas without trees
metrics_trees <- raster::extend(metrics_trees, metrics_2dchm, values = 0)
# compute canopy cover in trees and canopy mean height in trees in
# region of interest, because it is not in previous step.
r_tree_chm <- segms$filled_dem
# set chm to NA in non segment area
r_tree_chm[segms$segments_id == 0] <- NA
# compute raster metrics
dummy <- lidaRtRee::raster_metrics(raster::crop(r_tree_chm, raster::extent(b_box[1],
b_box[3], b_box[2], b_box[4])), res = resolution, fun = function(x) {
c(sum(!is.na(x$filled_dem))/(resolution/res_chm)^2, mean(x$filled_dem,
na.rm = T))
}, output = "raster")
names(dummy) <- c("TreeCanopy_cover_in_plot", "TreeCanopy_meanH_in_plot")
dummy <- raster::extend(dummy, metrics_2dchm, values = 0)
#
metrics_trees <- raster::addLayer(metrics_trees, dummy)
# ------------------------- compute 1D height metrics remove buffer
# points
a <- lidR::filter_poi(a, buffer == 0)
#
if (length(a) == 0) {
metrics_1d <- NULL
} else {
# all points metrics
metrics_1d <- lidR::grid_metrics(a, as.list(c(as.vector(quantile(Z, probs = percent)),
sum(Classification == 2), mean(Intensity[ReturnNumber == 1]))), res = resolution)
names(metrics_1d) <- c(paste("Hp", percent * 100, sep = ""), "nb_ground",
"Imean_1stpulse")
# vegetation-only metrics
a <- lidR::filter_poi(a, Classification != class_ground)
if (length(a) != 0)
{
dummy <- lidR::grid_metrics(a, as.list(c(mean(Z), max(Z), sd(Z),
length(Z), hist(Z, breaks = breaks_h, right = F, plot = F)$counts)),
res = resolution)
names(dummy) <- c("Hmean", "Hmax", "Hsd", "nb_veg", n_breaks_h)
# merge rasterstacks
dummy <- raster::extend(dummy, metrics_1d)
metrics_1d <- raster::addLayer(metrics_1d, dummy)
rm(dummy)
# ------------- merge metrics
metrics_1d <- raster::extend(metrics_1d, metrics_2dchm, values = 0)
metrics_1d <- raster::crop(metrics_1d, metrics_2dchm)
metrics_gaps <- raster::extend(metrics_gaps, metrics_2dchm, values = 0)
metrics_gaps <- raster::crop(metrics_gaps, metrics_2dchm)
metrics_edges <- raster::extend(metrics_edges, metrics_2dchm, values = 0)
metrics_edges <- raster::crop(metrics_edges, metrics_2dchm)
metrics_trees <- raster::crop(metrics_trees, metrics_2dchm)
#
temp <- raster::addLayer(metrics_1d, metrics_2dchm, metrics_gaps,
metrics_edges, metrics_trees)
temp[is.na(temp)] <- 0
} # end of vegetation-only points presence check
} # end of buffer-less points presence check
return(temp)
}
# -------------- mosaic rasters
names_metrics <- names(metrics[[1]])
metrics <- do.call(raster::merge, metrics)
names(metrics) <- names_metrics
# -------------------------- compute additional metrics Simpson index
metrics$Hsimpson <- raster::stackApply(metrics[[n_breaks_h[c(-1, -length(n_breaks_h))]]],
1, function(x, ...) {
vegan::diversity(x, index = "simpson")
})
# Relative density
for (i in n_breaks_h[c(-1, -length(n_breaks_h))]) {
metrics[[gsub("nb", "relativeDensity", i)]] <- metrics[[i]]/(metrics$nb_veg +
metrics$nb_ground)
}
# Penetration ratio compute cumulative sum
cumu_sum <- metrics[["nb_ground"]]
for (i in n_breaks_h) {
cumu_sum <- raster::addLayer(cumu_sum, cumu_sum[[dim(cumu_sum)[3]]] + metrics[[i]])
}
names(cumu_sum) <- c("nb_ground", n_breaks_h)
# interception ratio
intercep_ratio <- cumu_sum[[-1]]
for (i in 2:dim(cumu_sum)[3]) {
intercep_ratio[[i - 1]] <- 1 - cumu_sum[[i - 1]]/cumu_sum[[i]]
}
names(intercep_ratio) <- gsub("nb", "intercepRatio", names(cumu_sum)[-1])
# merge rasterstacks
metrics <- raster::addLayer(metrics, intercep_ratio)
# ---------------------- export as raster files for (i in 1:names(metrics)) {
# print(i)
# writeRaster(metrics[[i]],file=paste('output/raster_',i,'_',resolution,'m.tif',sep='')
# } ------- display
raster::plot(metrics[[c("Imean_1stpulse", "Hsimpson")]], main = c("Mean intensity of 1st pulse",
"Simpson indice of point heights"))