--- title: "R workflow for forest structure metrics computation from ALS data" author: "Jean-Matthieu Monnet, with contributions from B. Reineking and A. Glad" date: "`r Sys.Date()`" output: pdf_document: default html_document: default papersize: a4 bibliography: "../bib/bibliography.bib" --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) # Set so that long lines in R will be wrapped: knitr::opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE) knitr::opts_chunk$set(fig.align = "center") ``` --- Licence: GNU GPLv3 / [Source page](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/forest.structure.metrics.Rmd) 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: * 1D height metrics + 2D metrics of the canopy height model (CHM) + tree segmentation metrics (see also [tree segmentation tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/tree.detection.Rmd)) + forest gaps and edges metrics (see also [gaps and edges detection tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/gaps.edges.detection.Rmd)) 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*)(@Glad20). 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. ## Parameters Load `foreach` package and set the number of cores to use for parallel computing ```{r include = TRUE} 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. ```{r parameters, include = TRUE} # output metrics map resolution (m) resolution <- 10 # canopy height model resolution (m) resCHM <- 0.5 # buffer size (m) for tile processing (20 m is better for gaps metrics, 10 m is enough for # tree metrics) b.size <- 20 # height threshold (m) for high points removal (points above this threshold are considered # as outliers) h.points <- 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 l.sigma <- as.list(sigma) # height breaks for penetration ratio and density breaksH <- c(-Inf, 0, 0.5, 1, 2, 5, 10, 20, 30, 60, Inf) # percentiles of height distribution percent <- c(0.10, 0.25, 0.5, 0.75, 0.9) # surface breaks for gap size (m2) breaksGapSurface <- c(4, 16, 64, 256, 1024, 4096, Inf) # # gap surface names n.breaksG <- gsub("-","",paste("G.s", breaksGapSurface[c(-length(breaksGapSurface))],"_", breaksGapSurface[c(-1)],sep="")) # height bin names n.breaksH <- gsub("-","",paste("H.nb", breaksH[c(-length(breaksH))],"_", breaksH[c(-1)], sep="")) ``` 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. ```{r lascatalog, include = TRUE, fig.dim = c(3.5, 2.5), out.width='40%', warning=FALSE} # create catalog of LAS files cata <- lidR::catalog("../data/forest.structure.metrics") # set coordinate system sp::proj4string(cata) <- sp::CRS(SRS_string = "EPSG:2154") # set sensor type lidR::sensor(cata) <- "ALS" # disable display of catalog processing lidR::opt_progress(cata) <- FALSE # display sp::plot(cata, main="Bounding box of las files") ``` ## Workflow exemplified on one tile ### Load point cloud The tile is loaded, plus a buffer on adjacent tiles. Buffer points have a column "buffer" filled with 1. ```{r loadPointCloud, include = TRUE, message=FALSE, warnings=FALSE, fig.dim = c(4.5, 3.5), out.width='60%', } # 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]-b.size, b.box[2]-b.size, b.box[3]+b.size, b.box[4]+b.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 ``` 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. ```{r filterPointCloud, include = TRUE, message=FALSE, warnings=FALSE} # remove unwanted point classes, and points higher than height threshold a <- lidR::filter_poi(a, is.element(Classification, class.points) & Z<=h.points) # set negative heights to 0 a@data$Z[a@data$Z<0] <- 0 # summary(a@data) ``` ### Canopy height model The next step is to compute the canopy height model (CHM). It will be used to derive: * 2D canopy height metrics related to multi-scale vertical heterogeneity (mean and standard deviation of CHM, smoothed at different scales) + tree metrics from tree top segmentation + gaps and edges metrics from gap segmentation 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. ```{r computeCHM, include = TRUE, out.width = '50%', fig.dim=c(5, 4.5), message=FALSE} # compute chm chm <- lidaRtRee::points2DSM(a,res=resCHM) # replace NA, low and high values chm[is.na(chm) | chm<0 | chm>h.points] <- 0 # display CHM raster::plot(chm, asp=1, main="Canopy height model") ``` ## Metrics computation ### 2D CHM metrics 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. ```{r 2dchmMetrics, include = TRUE, fig.dim = c(9, 4.5), out.width='100%', message=FALSE} # for each value in list of sigma, apply filtering to chm and store result in list st <- lapply(l.sigma, FUN = function(x) { lidaRtRee::demFiltering(chm, nlFilter = "Closing", nlSize = 5, sigmap = x)[[2]] }) # 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. ```{r 2dmetrics2points, include = TRUE, message=FALSE, warning=FALSE} # 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/resCHM)^2 # modify layer names names(st) <- paste0("layer.", 1:6) # compute raster metrics metrics.2dchm <- lidaRtRee::rasterMetrics(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 ``` 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. ```{r 2dmetricsOutputs, include = TRUE, fig.dim = c(9, 5), out.width='100%', warning=FALSE} raster::plot(metrics.2dchm[[c("CHM0.sd", "CHM2.sd", "CHM8.sd", "CHM.mean", "CHM.PercInf1", "CHM.PercSup10")]]) ``` ### Gaps and edges metrics Gaps are computed with the function `gapDetection`. ```{r gapMetrics, include = TRUE, fig.dim = c(9, 3), out.width='100%', warning=FALSE, message=FALSE} # compute gaps gaps <- lidaRtRee::gapDetection(chm, ratio=2, gap.max.height=1, min.gap.surface=min(breaksGapSurface), closing.height.bin=1, nlFilter="Median", nlsize=3, gapReconstruct=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. ```{r gapStats, include = TRUE, warning=FALSE, message=FALSE} # 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::rasterMetrics(gaps.surface, res=resolution, fun=function(x) { hist(x$layer,breaks=breaksGapSurface, plot=F)$counts*(resCHM/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.breaksG, paste("G.s", min(breaksGapSurface), "_Inf",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 ``` Edges are detected with the function `edgeDetection` as the outer envelope of the previously delineated gaps. ```{r edgeDetection, include = TRUE, fig.dim = c(9, 3), out.width='100%', warning=FALSE, message=FALSE} # Perform edge detection by erosion edges <- lidaRtRee::edgeDetection(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. ```{r edgeStats, include = TRUE, warning=FALSE, message=FALSE} # 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::rasterMetrics(edges, res=resolution, fun=function(x) { sum(x$layer)*(resCHM/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 ``` 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. ```{r gapDisplay, include = TRUE, fig.dim = c(9, 3.2), out.width='100%', warning=FALSE, message=FALSE} par(mfrow=c(1,3)) # raster::plot(metrics.gaps[["G.s4_16"]], main="Prop of gaps 4-16m2") raster::plot(metrics.gaps[["G.s64_256"]], main="Prop of gaps 64-256m2") raster::plot(metrics.gaps[["G.s4_Inf"]], main="Prop of gaps >4m2") raster::plot(metrics.edges, main="Proportion of edges") ``` ### Tree metrics Tree tops are detected with the function `treeSegmentation` and then extracted with `treeExtraction`. ```{r treeMetrics, include = TRUE, fig.dim = c(9, 4.5), out.width='100%', , warning=FALSE, message=FALSE} # tree top detection (default parameters) segms <- lidaRtRee::treeSegmentation(chm, hmin=5) # extraction to data.frame trees <- lidaRtRee::treeExtraction(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 `stdTreeMetrics` 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. ```{r cropTreeMetrics, include = TRUE, warning=FALSE, message=FALSE} # 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::rasterMetrics(trees[,-1], res=resolution, fun=function(x) { lidaRtRee::stdTreeMetrics(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.treechm <- segms$filled.dem # set chm to NA in non segment area r.treechm[segms$segments.id==0] <- NA # compute raster metrics dummy <- lidaRtRee::rasterMetrics(raster::crop(r.treechm, 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/resCHM)^2, mean(x$filled.dem,na.rm=T))}, output="raster") names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot") # dummy <- raster::extend(dummy, metrics.trees) metrics.trees <- raster::addLayer(metrics.trees, dummy) # metrics.trees ``` Some outputs are displayed hereafter. ```{r displayTreeMetrics, include = TRUE, fig.width = 9, fig.height = 3.2, warning=FALSE, message=FALSE} 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)") ``` ### Point cloud (1D) metrics Buffer points are first removed from the point cloud, then metrics are computed from all points and from vegetation points only. ```{r 1dmetrics, include = TRUE} # 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("H.p",percent*100,sep=""),"nbSol","I.mean.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=breaksH, right=F, plot=F)$counts) ),res = resolution) names(dummy) <- c("H.mean","H.max","H.sd","nbVeg",n.breaksH) # 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) ```{r 1dmetricsAdditional, include = TRUE, warnings=FALSE} # Simpson index metrics.1d$H.simpson <- raster::stackApply(metrics.1d[[n.breaksH[c(-1,-length(n.breaksH))]]], 1, function(x, ...){vegan::diversity(x, index="simpson")}) # Relative density for (i in n.breaksH[c(-1,-length(n.breaksH))]) { metrics.1d[[paste0(i,"relative_density")]] <- metrics.1d[[i]] / (metrics.1d$nbVeg + metrics.1d$nbSol) } # Penetration ratio # cumulative sum cumu.sum <- metrics.1d[["nbSol"]] for (i in n.breaksH) { cumu.sum <- raster::addLayer(cumu.sum, cumu.sum[[dim(cumu.sum)[3]]] + metrics.1d[[i]]) } names(cumu.sum) <- c("nbSol", n.breaksH) # 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) <- paste0(names(cumu.sum)[-1], "ratio") # merge rasterstacks metrics.1d <- raster::addLayer(metrics.1d, intercep.ratio) # rm(cumu.sum, intercep.ratio) # metrics.1d ``` Some outputs are displayed hereafter. ```{r display1dmetrics, include = TRUE, fig.width = 9, fig.height = 3.2, warning=FALSE, message=FALSE} # display results par(mfrow=c(1,3)) raster::plot(metrics.1d$I.mean.1stpulse, main="Mean intensity (1st return)") raster::plot(metrics.1d$H.p50, main="Percentile 50 of heights") raster::plot(metrics.1d$H.nb2_5ratio, main="Interception ratio 2-5 m") ``` ## Batch processing 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. ```{r batchProcessing, include = TRUE, eval=TRUE, warning=FALSE, message=FALSE, fig.width = 9, fig.height = 2.6} # 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]-b.size, b.box[2]-b.size, b.box[3]+b.size, b.box[4]+b.size)) # if (class(a)=="try-error") # check if points are successfully loaded { # else NULL output temp <- NULL } else { # 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<=h.points) # check number of remaining points if (length(a)==0) { # output NULL temp <- NULL } else { # set negative heights to 0 a@data$Z[a@data$Z<0] <- 0 # # compute chm chm <- lidaRtRee::points2DSM(a, res=resCHM) # replace NA, low and high values chm[is.na(chm) | chm<0 | chm>h.points] <- 0 # ########################## # compute 2D CHM metrics # for each value in list of sigma, apply filtering to chm and store result in list st <- lapply(l.sigma, FUN = function(x) { lidaRtRee::demFiltering(chm, nlFilter = "Closing", nlSize = 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/resCHM)^2 # compute raster metrics metrics.2dchm <- lidaRtRee::rasterMetrics( 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::gapDetection(chm, ratio=2, gap.max.height=1, min.gap.surface=min(breaksGapSurface), closing.height.bin=1, nlFilter="Median", nlsize=3, gapReconstruct=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::rasterMetrics(gaps.surface, res=resolution, fun=function(x){ hist(x$layer, breaks=breaksGapSurface, plot=F)$counts*(resCHM/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.breaksG, paste("G.s", min(breaksGapSurface), "_Inf",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::edgeDetection(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::rasterMetrics(edges, res=resolution, fun=function(x) { sum(x$layer)*(resCHM/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::treeSegmentation(chm, hmin=5) # extraction to data.frame trees <- lidaRtRee::treeExtraction(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::rasterMetrics(trees[,-1], res=resolution, fun=function(x){ lidaRtRee::stdTreeMetrics(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.treechm <- segms$filled.dem # set chm to NA in non segment area r.treechm[segms$segments.id==0] <- NA # compute raster metrics dummy <- lidaRtRee::rasterMetrics( raster::crop(r.treechm, 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/resCHM)^2, mean(x$filled.dem,na.rm=T)) } , output="raster") names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot") 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("H.p",percent*100,sep=""),"nbSol","I.mean.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=breaksH, right=F, plot=F)$counts) ),res = resolution) names(dummy) <- c("H.mean","H.max","H.sd","nbVeg",n.breaksH) # 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 } # end of buffered area points presence check } # success of file reading 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$H.simpson <- raster::stackApply(metrics[[n.breaksH[c(-1,-length(n.breaksH))]]], 1, function(x, ...){vegan::diversity(x, index="simpson")}) # # Relative density for (i in n.breaksH[c(-1,-length(n.breaksH))]) { metrics[[paste0(i,"relative_density")]] <- metrics[[i]] / (metrics$nbVeg + metrics$nbSol) } # Penetration ratio # compute cumulative sum cum.sum <- metrics[["nbSol"]] for (i in n.breaksH) { cum.sum <- raster::addLayer(cum.sum, cum.sum[[dim(cum.sum)[3]]] + metrics[[i]]) } names(cum.sum) <- c("nbSol", n.breaksH) # interception ratio intercep.ratio <- cum.sum[[-1]] for (i in 2:dim(cum.sum)[3]) { intercep.ratio[[i-1]] <- 1 - cum.sum[[i-1]]/cum.sum[[i]] } names(intercep.ratio) <- paste0(names(cum.sum)[-1], "ratio") # 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("I.mean.1stpulse","H.simpson")]], main=c("Mean intensity of 1st pulse","Simpson indice of point heights")) ``` ## References