Operationalizing the use of TLS in forest inventories: the R package FORTLS

Environmental Modelling and Software 150 (2022) 105337 Available online 28 January 2022 1364-8152/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Operationalizing the use of TLS in forest inventories: The R package FORTLS Juan Alberto Molina-Valero a , * , Adela Martínez-Calvo a , María Jos´ e Ginzo Villamayor b , Manuel Antonio Novo P´ erez c , Juan Gabriel ´ Alvarez-Gonz´ alez a , Fernando Montes d , C´ esar P´ erez-Cruzado e a Unidad de Gesti´ on Ambiental y Forestal Sostenible (UXAFORES), Departamento de Ingeniería Agroforestal, Escuela Polit´ ecnica Superior de Ingeniería, Universidade de Santiago de Compostela, Benigno Ledo s/n, Campus Terra, 27002 Lugo, Spain b Departamento de Estadística, An´ alisis Matem´ atico y Optimizaci´ on, Facultad de Matem´ aticas, Universidade de Santiago de Compostela, Rúa Lope G´ omez de Marzoa s/n, Campus Vida, 15782 Santiago de Compostela, Spain c Departamento de Matem´ aticas, Facultad de Inform´ atica, Universidade da Coru˜ na, Campus Elvi˜ na s/n, 15071 A Coru˜ na, Spain d Centro de Investigaci´ on Forestal (INIA, CSIC), Ctra. De la Coru˜ na km 7,5, 28040 Madrid, Spain e Proyectos y Planificaci´ on (PROEPLA), Departamento de Producci´ on Vegetal y Proyectos de Ingeniería, Escuela Polit´ ecnica Superior de Ingeniería, Universidade de Santiago de Compostela, Benigno Ledo s/n, Campus Terra, 27002 Lugo, Spain ARTICLE INFO Keywords: Forest monitoring Forest stands parameters LiDAR Precision forestry Remote sensing Terrestrial-based-technologies ABSTRACT Terrestrial Laser Scanning (TLS) devices show great potential for application in Forest Inventories (FIs) as they are capable of registering high resolution point clouds rapidly and automatically. Nevertheless, operational use of TLS for FI purposes has been hampered by the absence of algorithms for processing the acquired data, particularly in the single-scan mode, as occlusions result in loss of information. The R package FORTLS has been developed to overcome this obstacle, as it automates the processing of single-scan TLS point cloud data for forestry purposes and includes several features that deal with occlusions. FORTLS makes use of the main advantage of the single-scan scenario in FI, thus improving the efficiency of data acquisition and post-processing. All of these features of the FORTLS package are potentially valuable for the operational use of TLS in FIs, in combination with inference techniques derived from model-based and model-assisted approaches. 1. Introduction Information about forest resources is essential for sustainable forest management and development of forest policies, and forest inventories (FIs) are fundamental for estimating and monitoring the state and evolution of forest resources at both global and regional scales (Tomppo et al., 2010). FIs have improved since they were first introduced, owing to the continuous appearance of new technologies, especially in the last few decades since the emergence of remote and proximal sensing. In this technological context, light detection and ranging (LiDAR) systems provide 3-dimensional point clouds, which are suitable for estimating tree attributes and are very useful for many forestry applications (Dubayah and Drake, 2000). This technology has proved operationally viable for estimating essential FI variables at stand level, such as arithmetic mean height (h, m), basal area (G, m 2 ha −1 ) and volume (V, m 3 ha −1 ), with airborne laser scanning (ALS) devices (Wulder et al., 2012; White et al., 2016). Terrestrial laser scanning (TLS) devices, such as LiDAR devices with millimetric precision, are considered to show great potential for enhancing FIs (Dassot et al., 2011; White et al., 2016) and also forest ecology research (Calders et al., 2020; Danson et al., 2018). Apart from much higher spatial resolution under canopy, the main advantages of using TLS data rather than ALS data are better observation of near-ground vegetation (White et al., 2016) and thus better trunk coverage for estimating the woody component, which is one of the most important components in FIs. In fact, TLS-based approaches can provide very accurate estimates of the diameter at breast height (dbh, measured at 1.3 m from the ground) and the stem curve; with the single-scan approach yielding values of 1–4 and 1.3–6 cm respectively, * Corresponding author. E-mail addresses: [email protected] (J.A. Molina-Valero), [email protected] (A. Martínez-Calvo), [email protected] (M.J. Ginzo Villamayor), [email protected] (M.A. Novo P´ erez), [email protected] (J.G. ´ Alvarez-Gonz´ alez), [email protected] (F. Montes), cesar. [email protected] (C. P´ erez-Cruzado). Contents lists available at ScienceDirect Environmental Modelling and Software journal homepage: www.elsevier.com/locate/envsoft https://doi.org/10.1016/j.envsoft.2022.105337 Received 13 September 2021; Received in revised form 18 January 2022; Accepted 24 January 2022 Environmental Modelling and Software 150 (2022) 105337 2 depending on the stand conditions (Liang et al., 2018a), which are close to the values required in practical applications such as FIs. However, TLS devices have not been yet adopted in FIs, for several reasons: (i) difficulties in the automation of data processing to provide reliable measurements of important forest variables, (ii) high acquisition costs; (iii) limited software; and (iv) lack of trained personnel (Liang et al., 2016). Many researchers agree that affordability is the main key challenge to overcome, emphasizing that automation of point cloud processing with attainable and easy-to-use software able to extract information related to important forest attributes is essential (Dassot et al., 2011; Newnham et al., 2015; White et al., 2016; Liang et al., 2016, 2018a). As TLS data sets comprise millions of points, sophisticated methods for automatic processing are required. Many algorithms with a high level of automation and that are able to extract tree attributes, such as dbh, total height (h, m) and stem volume (v, m 3 ), have been developed in the last few decades (Cabo et al., 2018; Liang et al., 2012, 2018b; Olofsson et al., 2014; Olofsson and Holmgren, 2016; Zhang et al., 2019). Although most algorithms yield acceptable dbh and stem curve estimations according to FI requirements, stem detection and estimation of h cause bottlenecks in the process, especially with single scans (Krok et al., 2020; Liang et al., 2018a). Some of the algorithms developed have also been included in software applications, such as SimpleForest (Hackenberg et al., 2015), 3D Forest (Trochta et al., 2017) and AutoStem™ (Bienert et al., 2007), among others (Krok et al., 2020). However, these programs have some drawbacks for use in FIs: (i) they focus on single-tree rather than stand-level approaches (SimpleForest); (ii) they involve semi-automatic processing (3D Forest); and (iii) the software is only available commercially (AutoStem™) (i.e. it is not free or open source). Furthermore, the previous studies have mainly focused on replicating plot-based measurements, which does not extend conventional inventory approaches from a sampling perspective and thus limits the utility of TLS in FIs, the main purpose of which is to estimate important forest variables at larger scales (e.g. stand and regional) than tree and plot levels (Newnham et al., 2015; White et al., 2016). Methods that enable TLS to be used for FI purposes must therefore contemplate the use of different approaches (Newnham et al., 2015), and other procedures in which not all trees in the sample plots are measured may be feasible (Liang et al., 2018a). Thus, further research is required to address the challenges in the operational use of TLS (Liang et al., 2018a). Here we present FORTLS (Molina-Valero et al., 2021), an R package developed with the objective of automating TLS point cloud data processing and estimating variables for forestry purposes. To fulfil this objective, FORTLS enables (i) detection of trees and estimation of dbh and other tree attributes, (ii) estimation of some stand variables (e.g. h, G, V), (iii) computation of metrics related to important tree attributes estimated in FIs at stand-level, and (iv) optimization of plot design for combining TLS data and field measured data. The package also includes several features for correcting occlusion problems to improve the estimation of stand variables. The current version of FORTLS is based on single-scan TLS data, with the aim of facilitating operational use, and it has been designed as relatively easy-to-use open-source software aimed at use by both scientists and technical users. Relative to multi-scan and multi-single-scan approaches, the single-scan approach improves data acquisition, shortens the processing time and increases the sample size in a cost-efficient manner, mainly because it does not require pre-scanning tasks involving location of artificial reference objects (Holopainen et al., 2014) or automated post-processing matching methods (Liu et al., 2017). Finally, and as a case study, we have tested FORTLS for estimating common forestry variables in an experimental plot of 1 ha located in an even-aged Pinus sylvestris L. stand in northern Spain. These features of the FORTLS package may enable the operational use of TLS in FIs, in combination with model-based or model-assisted inference approaches. 2. Methods 2.1. Software design FORTLS (Molina-Valero et al., 2021) has been developed as an R package (R Core Team, 2021) because R is free statistical software which is accessible to any user interested in this tool. The initial stages of development of this package were outlined in Molina-Valero et al. (2020), although the first version of FORTLS was not available until March 2021. Currently, both the most recent stable version of the package and the most up-to-date version can be downloaded free of charge, from respectively the CRAN (https://CRAN.R-project.org/packa ge=FORTLS) and GitHub development (https://github.com/Mo lina-Valero/FORTLS/tree/devel) repositories. The R package FORTLS has been optimized by implementing C++ code in the most demanding computing processes by means of the Rcpp package (Eddelbuettel, 2013; Eddelbuettel and Balamuta, 2018; Eddelbuettel and François, 2011) and the RcppEigen package (Bates and Eddelbuettel, 2013), which enables integration of the Eigen C++library for specific matrix calculation. For operations with objects in spatial data classes, both the raster (Hijmans, 2020) and sp (Bivand et al., 2013; Pebesma and Bivand, 2005) packages have been used. For obtaining Voronoi polygons, the ggvoronoi package (Garrett et al., 2021) has been used in the simulations and metrics.variables functions. As TLS point clouds represent large data sets, FORTLS also imports the vroom package (Hester and Wickham, 2020) for accelerating loading and saving .txt files. We have used other important packages to generate and save interactive graphics, namely plotly (Sievert, 2020) and htmlwidgets (Vaidyanathan et al., 2020). Apart from the packages included in R base distribution and other accessory packages such as progress (Cs´ ardi and FitzJohn, 2019), scales (Wickham and Seidel, 2020) and tidyr (Wickham, 2021), the other external R packages used for more specific functions are mentioned below, with their respective functions. The functions and results compiled in this work are based on the stable version 1.0.6 of the FORTLS package available in CRAN. In the following sections, all steps involved in TLS point cloud data processing with FORTLS, as well as the most relevant algorithms, are described: (i) normalization; (ii) tree detection; and (iii) estimation of metrics and variables at stand-level. 2.1.1. Normalization The normalization process is a necessary first step in processing point cloud data, and it is implemented in the normalize function (Table 1), which for some processes uses the functions readLAS, clip_circle, classify_ground, grid_terrain and normalize_height included in the lidR package (Roussel et al., 2020; Roussel and Auty, 2020). Normalization involves obtaining the coordinates relative to plot centre for TLS point clouds supplied as .las or .laz files. The process includes the following steps: (i) classification of points as “ground”; (ii) generation of a digital terrain model (DTM); (iii) computation of coordinates relative to DTM (Cartesian, cylindrical and spherical); and (iv) reduction of point cloud density by the point cropping process (PCP). In the initial step, points are classified as “ground” or “not ground” with the Cloth Simulation Filter (CSF) algorithm (Zhang et al., 2016). The DTM is then generated by spatial interpolation of “ground” points. Two methods are available for executing this process: (i) spatial interpolation based on Delaunay triangulation (by default); and (ii) spatial interpolation using a k-nearest neighbour approach with inverse-distance weighting. The point cloud is then normalized by subtracting the DTM created. Once the point cloud has been normalized, Cartesian, cylindrical and spherical coordinates are calculated relative to the sampling point (TLS device establishment point). Finally, the normalize function applies the PCP algorithm developed by Molina-- Valero et al. (2019) to reduce the point density and thus produce a spatially homogeneous point cloud in which the distribution of points is proportional to the object size. During execution of the PCP, a selection J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 3 probability prob(p)is assigned to each cloud point p according to Eq. (1) (Fig. 1): prob(p) = rp rmax (1) where rp is the radial distance of the point p from the plot centre, and rmax is the radial distance of the farthest point from the plot centre. Finally, the point is selected if the selection probability is equal to or higher than a random value generated from the uniform distribution for the interval [0, 1]. 2.1.2. Tree detection Tree detection represents a very important and challenging step in estimating variables of interest in TLS-assisted FIs. This is partly due to occlusions, which are much more important when working with single scans. This process of tree detection is implemented in the tree.detection function (Table 1), which has been designed to detect as many trees as possible from point cloud data obtained as normalize function output. In the tree detection process, one or several horizontal slices from the original point cloud are extracted. Slices at heights of 1.0, 1.3 and 1.6 m ±5 cm over the terrain (all around 1.3 m as reference section to estimate dbh) are considered by default in tree.detection. The probability of detecting trees increases when more than one slice is considered. However, the tree.detection function can also extract (if specified in the arguments) slices at other heights when trees are not identifiable at these pre-established values. Each horizontal slice is processed by algorithms that are able to (i) remove branches and foliage points, (ii) detect point clusters corresponding to potential tree sections, and (iii) classify detected point clusters as tree sections, or not, according to several tests. The tree sections detected for all horizontal slices are then merged, and for each tree detected the tree.detection function estimates certain attributes by: (i) calculating the coordinates corresponding to tree normal section centre (1.3 m above ground level) and its horizontal distance from plot centre; (ii) estimating the dbh; (iii) classifying the tree as fully visible or partially occluded; and (iv) obtaining the number of points corresponding to normal section slice (1.3 m ±5 cm) for both original and reduced (by applying PCP) point clouds. The main steps involved in the previously mentioned algorithms for detecting tree sections for each horizontal slice and for estimating metrics and variables related to detected trees attributes are described in the following subsections and summarized in Fig. 2. 2.1.2.1. Removing branches and foliage points. For each horizontal slice, this first step aims to remove points corresponding to fine branches and foliage (e.g. leaves and shrubs) and mainly to retain stem points, for which we considered local surface variation, also known as the normal change rate (NCR). This is a quantitative measure of curvature feature useful for discerning some “noisy” points (Pauly et al., 2002) with higher values representing more curved surfaces (predictably fine branches and foliage). The NCR index is estimated at point level considering a local neighbourhood, as follows. Given a fixed radius r, the set of local Table 1 Main functions included in the FORTLS package. Function (arguments) Description correlations(simulations, variables =c ("N","G","V","d","dg","d.0","h","h.0"), method =c("pearson","spearman”), save.result =TRUE, dir.result =NULL) Correlation between field estimates and TLS metrics distance.sampling(tree.list.tls, id.plots =NULL, strata.attributes =NULL) Distance sampling methods for correcting occlusion effects estimation.plot.size(tree.list.tls, plot.parameters =list( radius.max =25, k.tree.max =50, BAF.max = 4), average =FALSE, all.plot.designs =FALSE) Assessment of consistency of metrics for simulated TLS plots metrics.variables(tree.list.tls, distance.sampling =NULL, plot.parameters, dir.data =NULL, save.result =TRUE, dir. result =NULL) Computation of metrics and variables for TLS plots normalize(las, x.center =NULL, y.center =NULL, max.dist =NULL, min.height =NULL, max. height =NULL, algorithm.dtm ="tin”, res.dtm =0.2, csf =list(cloth_resolution =0.5), id =NULL, file =NULL, dir.data =NULL, save.result =TRUE, dir. result =NULL) Production of relative coordinates and density reduction for TLS point clouds optimize.plot.design(correlations, variables =c ("N","G","V","d","dg","d.0","h","h.0"), dir.result =NULL) Optimization of plot design based on optimal correlations relative.bias(simulations, variables =c ("N","G","V","d","dg","d.0","h","h.0"), save.result =TRUE, dir.result =NULL) Relative bias between field estimates and TLS metrics simulations( tree.list.tls, distance.sampling =NULL, tree. list.field, plot.parameters =list( radius.max =25, k.tree.max =50, BAF.max = 4), dir.data =NULL, save.result =TRUE, dir. result =NULL) Computation of metrics and variables for simulated TLS and field plots tree.detection(data, dbh.min =7.5, dbh.max =200, ncr.threshold =0.1, tls.resolution =list(), breaks =c(1.0,1.3,1.6), plot.attributes =NULL, save.result =TRUE, dir.result =NULL) Tree detection and cross section estimation tree.detection.multiple(las.list, id =NULL, file =NULL, normalize.arguments =list( max.dist =NULL, min.height =NULL, max. height =NULL, algorithm.dtm ="tin”, res.dtm =0.2), tree.detection.arguments =list( dbh.min =7.5, dbh.max =200, ncr.threshold =0.1, tls.resolution =list(), breaks =c(1.0,1.3,1.6), plot.attributes =NULL), dir.data =NULL, save.result =TRUE, dir. result =NULL) Tree detection and cross section estimation for multiple plots Fig. 1. Selection probability of points corresponding to a single-scan according to the radial distance from the plot centre, implemented in the normalize function. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 4 neighbours of a point p is denoted by {pi}i∈Np,r, where Np,r is the index set of all cloud points satisfying the condition that d(p,pi)<r (where d is the Euclidean distance). The NCR for point p is then estimated by eigenanalysis of the 3 ×3 covariance matrix Cp of its local neighbourhood (Eq. (2)): Cp=1 Np,r∑ i∈Np,r(pi−p)(pi−p)T(2) where p is the centroid of the local neighbourhood of p. Obtaining the eigenvalues {λi}2 i=0 by singular value decomposition of Cp and assuming that λ0≤λ1≤λ2, λ0 describes the variation along the normal surface, the extent to which the points deviate from the tangent plane is estimated (Pauly et al., 2002). Hence, the NCR index for radius r at point p is defined as follows (Eq. (3)): NCRr(p)= λ0 λ0+λ1+λ2 (3) Once NCR is computed, p is retained as a stem point if its NCR value is lower than a pre-established threshold (Fig. 3a). In accordance with other studies, in the tree.detection function, the neighbourhood was established by a radius of 5 cm as suitable for calculating NCR for the stem separation in forests (Ma et al., 2016; Xia et al., 2015). The threshold value of NCR used to remove branches and foliage points was established as 0.1 by default, also according to other studies in which the index has already been used with the same objective and obtaining good results (e.g. Jin et al., 2016; Zhang et al., 2019). Nevertheless, other NCR thresholds can be specified by users in the corresponding argument of the tree.detection function. 2.1.2.2. Detection of point clusters. Once branches and foliage points have been removed from the horizontal slice, the next step is to detect point clusters corresponding to potential tree sections. This process involves several steps: (i) a clustering process is first applied to the horizontal projection of Cartesian coordinates of points; (ii) points corresponding to possible branches are then removed using surface density approaches; and (iii) clusters receiving fewer points than expected for a full visible stem are discarded. As mentioned above, potential tree sections are first detected through a clustering process applied to the horizontal projection of Cartesian coordinates. This clustering process is performed by the Density-Based Spatial Clustering of Applications based on the Noise (DBSCAN) method (Ester et al., 1996) and applied with the dbscan function in the R package dbscan (Hahsler et al., 2019). The size of the epsilon neighbourhood is established as the minimum distance between two consecutive points at the farthest distance from TLS in the respective horizontal slice, and a minimum of 5 points required in that epsilon neighbourhood. This algorithm detects as many as possible sections corresponding to trees (according to occlusion conditions), as well as other possible clusters generated by other items (branches, shrubs, etc.), and it has been used in previous studies of TLS with the same objective (Ferrara et al., 2018; Molina-Valero et al., 2019). For refining extracted stem points, we used a similar approach to that proposed by Zhang et al. (2019) in order to remove any branches remaining in the clusters, based on the principle that stems should generate more points than other parts of the tree (due to the fractal-size distribution according to West et al., 1999); in addition, these points should have a predominantly vertical distribution. Thus, if the point cloud is vertically projected and rasterized in a grid, cells over stems usually include more points than those located over branches and foliage. Each cluster is thus rasterized on the horizontal plane (Cartesian coordinates) with an adapted grid step size of twice the distance Fig. 2. FORTLS workflow for normalization (red) and tree detection (green) processes. All of the processes included here are described in sections 2.1.1 and 2.1.2. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 5 between two consecutive points at mean cluster distance from TLS, and those points included in cells with fewer points than the median value of the number of points per cell will be removed (Fig. 3b). At this point of the process, several tests were used to distinguish those clusters belonging to tree sections. Each cluster was first projected vertically with cylindrical coordinates (φ, z) and divided into regular strips bound by vertical scan resolution ( α v) (Fig. 3c). As stems completely visibly from TLS must generate the maximum possible number of points according to scan resolution and distance from the TLS instrument, we estimated the approximate number of points that each strip must include if the section corresponds to a fully visible tree. As normalized coordinates are projected on a horizontal plane, the slope effect was first corrected to calculate the vertical resolution (Δv) in coordinate z at cluster real distance from TLS (m), as follows (Eq. (4)): Δv =2×tan( α v/2) rcluster/cos(slope) (4) where α v is the vertical scan resolution (rad), rcluster is the mean radial distance in spherical coordinates from TLS to cluster (m), and slope is the mean slope of cluster according to DTM (rad). The number of points n that each strip must contain in fully visibility conditions is then computed as follows (Eq. (5)): n=Δz Δv (5) where Δz is the slice thickness (m), and Δv is the previously defined vertical resolution (Eq. (4)). Finally, we reduced the predicted number of points by 30% in very stepped terrains (>0.5 rad). Only those clusters with at least one strip containing the number of points n that every strip Fig. 3. Detection of potential tree sections implemented in the tree.detection function: a) points corresponding to branches and foliage (green) removed by applying the NCR index for a threshold of 0.025; b) refinement of extracted stem points where those points included in grid cells with fewer points than the median (red points) are removed; c) estimation of the approximate number of points (n) that each strip must contain if the section corresponds to a tree fully viewed by TLS, the green points represent those cells containing more or equal points than n; and d) estimated location of tree section centre (green point) and radius calculated (green wide arrow) as the average of all distances (red fine arrows) between cluster points and the grid intersection allocated as the tree centre. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 6 must contain in fully visibility conditions were selected. Once clusters have fulfilled all of the previous checks, the next step involves obtaining the centre of the potential tree section. With this aim, regular square grids of 1 cm were overlapped on each cluster selected. The tree section centre was then considered as the intersection grid point where the variance of the distances between this intersection and all the cluster points reaches the lowest value (Fig. 3d), considered to occur when the coefficient of variation of distances is smaller than 0.1, otherwise the cluster is dismissed. Finally, the radius of the tree section was computed as the average of all the distances between the estimated centre and remaining cluster points at the moment of algorithm processing. 2.1.2.3. Cluster classification. This step consists of checking multiple geometrics features and indices to verify which clusters finally correspond to tree sections. Most criteria used are based on those determined by Molina-Valero et al. (2019), and they were applied to each cluster. The first test involves checking whether the centre is located behind the cluster points relative to TLS. This is considered fulfilled when at least 95% of the cluster points have a lower cylindrical coordinate ρ (hence closer than plot centre) than the tree centre. The second test consists of checking for the absence of points behind the tree surface, which is considered true when at least 95% of the distances between cluster points and the tree centre are greater than half of the estimated radius. This may be visually checked by means of a distance histogram (Fig. 4a). The similarity between the cluster shape and the circumference arc is then assessed, checking that extreme points in cylindrical coordinate ρ are farther from TLS than central points (Fig. 4b). This is possible when trees are largely visible from TLS, but otherwise will not be possible due to partial occlusions. In such cases, the clusters were checked to determine whether they form a smaller arc of a circle (Fig. 4c). For this purpose, we calculated the Pearson coefficient correlations for φ values in increasing order and the correlative numbering (Fig. 4d). Those clusters with values below 0.995 were removed. 2.1.2.4. Estimating tree attributes. When several sections are identified at different heights, those corresponding to the same trees are joined using the DBSCAN algorithm on the horizontal projection, and some tree Fig. 4. Representation of some of the tests implemented in the tree.detection function to distinguish those clusters corresponding to tree sections: a) histogram of distances between cluster points and estimated tree centre to check for the absence of points behind tree surface; b) mean coordinates of points included in first and last φ coordinate percentiles (red points) and in the middle cluster φ ±TLS angle aperture (red triangle); c) tree partially occluded; and d) assessment of the Pearson’s correlation between φ values in increasing order and their corresponding correlative numbers. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 7 attributes are obtained: coordinates of normal section centre and horizontal distance from plot centre, estimated dbh, indicator of partial occlusion and number of points corresponding to normal section for original and reduced point clouds. When trees are exclusively detected at 1.3 m, the dbh is estimated directly as twice the radius estimated at this height. Conversely, when trees are detected from other section(s) (including 1.3 m or not), a linear taper equation is fitted with radius as the response variable radius and section height (hsech sec ) as the explanatory variable (Eq. (6)): radius =β0+β1⋅hsec (6) The radius at 1.3 m is then predicted, as follows (Eq. (7)):  radius1.3=radiusi+ β1⋅(1.3−hseci)(7) where radiusi and hseci are the estimated radius and the height corresponding to section i, and  β1 is the slope parameter fitted in the linear regression (Eq. (6)). Hence, dbh is computed as twice the averaged predicted radius. Finally, the number of points (num.points) corresponding to a normal section (+/- 5 cm) in the original point cloud is computed and also estimated for each detected tree j by using the dbh values previously computed as follows (Eq. (8)): num.points.estj=dbhj∑i∈I num.pointsi dbhi #I(8) where the index set I corresponds to all detected trees fully visible at 1.3 m, and #I is the number of trees fully visible at 1.3 m. The number of points and the estimated number of points for the point cloud reduced by PCP are obtained in a similar way. 2.1.3. Computing TLS metrics and variables at stand level Once normalization and tree detection processes for each point cloud are completed, TLS metrics and variables can be estimated at stand level for three different plot designs, all of which are included in the metrics. variables function (Table 1). These plot designs are circular fixed area, ktree and angle-count (Bitterlich, 1948) plots (Fig. 5), each of which is defined by a unique design parameter (radius, k and basal area factor (BAF), respectively) that must be specified in the function arguments. The metrics and variables computed for each plot design are summarized in Table 2, and more details are compiled below. The two approaches for optimizing plot design by means of FORTLS functions are also described. 2.1.3.1. Stand-level metrics. Metrics are computed using points directly from normalized point cloud in a relatively similar way as in the FUSION/LDV software for LiDAR data analysis and visualization (McGaughey, 2009). These are statistical descriptive measures such as percentiles or only number of points belonging to specific sections of point cloud. The following metrics are available: - Total number of points corresponding to the normal section (+/- 5 cm) of trees detected after removing some noisy points with NCR and refinement of extracted stem points processes (see section 2.1.2). These can be computed for raw point clouds (num.points) and reduced point clouds (num.points.hom). - Number of estimated points corresponding to the normal section (+/- 5 cm) of trees detected. The number of points corresponding to each tree can be estimated according to Eq. (8) for raw point clouds (num.points.est) and, analogously, for reduced point clouds (num. points.hom.est). - Percentiles of z coordinate (m). Computed percentiles are P01, P05, P10, P20, P25, P30, P40, P50, P60, P70, P75, P80, P90, P95 and P99. - Descriptive statistics of z coordinate distribution: mean, maximum (max), minimum (min), standard deviation (sd), variance (var), mode, kurtosis and skewness. - Percentage of points above mode (perc_on_mode) and mean (perc_on_mean) values of z coordinates. - Scale (weibull_b) and shape (weibull_c) parameters of a Weibull distribution fitted to z coordinates distribution. 2.1.3.2. Stand-level variables. Variables represent estimates based on the attributes of trees detected from TLS point cloud data, further aggregated at stand-level, and finally expanded to unit area (ha). The following variables are available: - Apparent stand density (N.tls, trees ha −1 ), which is estimated for trees detected in a similar procedure to that used in conventional inventories for circular fixed area and k-tree plots (Eq. (9)) and angle-count plots (Eq. (10)): N.tls =10000 π R2⋅n(9) N.tls =∑ n i=1 BAF gi (10) where R is the plot radius (m), n is the number of trees detected in the corresponding plot design, BAF is the basal area factor (m 2 ha −1 ), and gi is the basal area of the tree i (m 2 ). - Apparent stand basal area (G.tls, m 2 ha −1 ), which is estimated for trees detected in a similar procedure as that used in conventional inventories for circular fixed area and k-tree plots (Eq. (11)) and angle-count plots (Eq. (12)): G.tls =10000 π R2∑ n i=1 gi(11) Fig. 5. Example of the three plot designs considered in FORTLS for the same sampling point. Circular fixed area plot of 10 m radius (red) includes 7 trees; ktree plot for k =8 (green) yields a larger radius than circular fixed area; and angle-count plot (blue) with a particular BAF includes 7 trees. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 8 G.tls =BAF⋅n(12) - Apparent stand stem volume (V.tls, m 3 ha −1 ), which is estimated for trees detected by modelling stem profile as a paraboloid and calculating the volumes of revolution for fixed area and k-tree plots (Eq. (13)) and angle-count plots (Eq. (14)): V.tls =10000 π R2∑ n i=1 π ⋅h2 P99 i 2⋅(dbhi / 2 )2 (hP99i−1.3)2(13) V.tls =∑ n i=1 BAF gi ⋅ π ⋅h2 P99 i 2⋅(dbhi / 2 )2 (hP99i−1.3)2(14) where hP99 i and dbhi are the 99th percentile of points delimited by Voronoi polygons (m) (i.e. estimates of h) and dbh (m) for tree i, respectively. - Mean and dominant diameters (cm), which are estimated for arithmetic, quadratic, geometric and harmonic means. In the case of dominant diameters, only the n largest trees per ha (according to dbh) are considered. Although it can be specified in the arguments, the 100 largest trees per ha are considered by default. - Mean and dominant heights (m), which are estimated for arithmetic, quadratic, geometric and harmonic means. In the case of dominant heights, only the n largest trees per ha (according to dbh) are considered. Although the number of trees can be specified in the arguments, the 100 largest trees per ha are considered by default. In the previous calculations, for the k-tree design, the plot radius R is defined as the mean of horizontal distances of trees k and k+1 (Kleinn and Vilˇ cko, 2006). 2.1.3.3. Dealing with occlusions. All FORTLS functions used for estimating stand variables also include correction of occlusions approaches. In the case of angle-count plots, occlusion corrections are based on gap probability attenuation with distance from TLS depending on a Poisson distribution. In the case of circular fixed area and k-tree plots, distance sampling methods and shadowing effect correction are considered. In order to obtain occlusion corrections based on distance sampling methods, the distance.sampling function must be executed previously and the values obtained must be incorporated as an argument for functions which compute stand variables, whereas the other corrections are computed by default. A brief description of the implemented occlusion corrections is given below. 2.1.3.4. Poisson attenuation model. This method has been used in measurements with TLS (Strahler et al., 2008; Lovell et al., 2011) and optical (Montes et al., 2019) instruments to reduce the device-related bias in the relascope-based approach. It is based on geometric gap probability (Pgap), which decreases exponentially following a Poisson distribution (Eq. (15)): Pgap(λ,DE,R)=e(− λDER)(15) where λ is the number of trees per m 2 , DE is the effective dbh and R is the Table 2 FORTLS metrics and variables. Metrics and variables Description N a /N.tls Stand density (N, trees ha −1 ) N.hn b , N.hr b , N.hn.cov b , N.hr.cov b Stand density (N, trees ha −1 ) with occlusion corrections based on distance sampling methodologies N.sh b Stand density (N, trees ha −1 ) with correction of the shadowing effect N.pam c Stand density (N, trees ha −1 ) with occlusion correction based on a Poisson attenuation model G a /G.tls Stand basal area (G, m 2 ha −1 ) G.hn b , G.hr b , G.hn.cov b , G.hr.cov b Stand basal area (G, m 2 ha −1 ) with occlusion corrections based on distance sampling methodologies G.sh b Stand basal area (G, m 2 ha −1 ) with correction of the shadowing effect G.pam c Stand basal area (G, m 2 ha −1 ) with occlusion correction based on a Poisson attenuation model V a /V.tls Stand volume (V, m 3 ha −1 ) V.hn b , V.hr b , V.hn.cov b , V.hr.cov b Stand volume (V, m 3 ha −1 ) with occlusion corrections based on distance sampling methodologies V.sh b Stand volume (V, m 3 ha −1 ) with correction of the shadowing effect V.pam c Stand volume (V, m 3 ha −1 ) with occlusion correction based on a Poisson attenuation model num.points, num.points.est, num.points.hom, num. points.hom.est Number of points and estimated number of points corresponding to normal sections (+/- 5 cm) of trees in the original point cloud (num.points and num.points.est, respectively); and number of points and estimated number of points (points) corresponding to normal sections (+/- 5 cm) of trees in the reduced point cloud (num.points.hom and num.points.hom.est, respectively). d a , dg a , dgeom a , dharm a /d.tls, dg.tls, dgeom.tls, dharm.tls Stand mean dbh, using arithmetic mean (d, cm), quadratic mean (d g , cm), geometric mean (cm), and harmonic mean (cm), respectively. h a , hg a , hgeom a , hharm a /h.tls, hg.tls, hgeom.tls, hharm.tls Stand mean h, using arithmetic mean (h, m), quadratic mean (m), geometric mean (m), and harmonic mean (m), respectively. d.0 a , dg.0 a , dgeom.0 a , dharm.0 a /d.0.tls, dg.0.tls, dgeom.0.tls, dharm.0.tls Stand dominant mean dbh, using arithmetic mean (d 0 , cm), quadratic mean (cm), geometric mean (cm), and harmonic mean (cm), respectively. h.0 a , hg.0 a , hgeom.0 a , hharm.0 a /h.0.tls, hg.0.tls, hgeom.0.tls, hharm.0.tls Stand dominant mean h, using arithmetic mean (H 0 , m), quadratic mean (m), geometric mean (m), and harmonic mean (m), respectively. P01, P05, P10, P20, P25, P30, P40, P50, P60, P70, P75, P80, P90, P95, P99 Height percentiles (m) derived from z coordinates of TLS point clouds relative to ground level. mean d , max d , min d , sd d , var d , mode d , kurtosis d , skewness d Descriptive statistics of z coordinates distribution of TLS point clouds relative to ground level: mean, maximum, minimum, standard deviation, variance, mode, kurtosis and skewness, respectively. perc.on.mode d , perc.on.mean d Percentage of points above mode and mean values of z coordinates of TLS point clouds relative to ground level. weibull_b d , weibull_c d Scale and shape parameters, respectively, for Weibull distribution fitted for z coordinates of TLS point clouds relative to ground level. a Variables estimated from field measurement data. b Variables only estimated for circular fixed area and k-tree plot design. c Variables only estimated for angle-count plot design. d Variables not included in FORTLS example data because they were included in the previous version (1.0.3). J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 9 horizontal distance to the farthest tree. The corrected stand density (N. pam, trees ha −1 ) for angle-count plots is then obtained as follows (Eq. (16)): N.pam =N.tls /F(λDER)(16) where F is a function defined as F(t) = 2 t2(1−e−t(1+t)), with t=λDED 2 BAF √: see Strahler et al. (2008) and Lovell et al. (2011) for further details. Similarly, corrected stand basal area (G.pam, m 2 ha −1 ) and volume (V. pam, m 3 ha −1 ) are computed. 2.1.3.5. Point transect sampling. This approach is based on the point transects method from distance sampling methods (Buckland et al., 2001). These methods use detection functions g(r,θ)with variable r (distance from sampling point) and parameter θ, which describe how the probability of detection decreases as distance increases. We used Half-Normal (Eq. (17)) and Hazard-Rate (Eq. (18)) functions as these have been successfully used in measurements with TLS (Astrup et al., 2014) and optical (Montes et al., 2019) devices: g(r,θ)=e(−r2 2 σ 2)(17) g(r,θ)=1−e−(r σ )−b (18) Parameter θ includes the shape (b) (only in the Hazard-Rate function) and the scale ( σ ), which was also expanded with dbh as a covariate into an exponential function (Eq. (19)) according to Ducey and Astrup (2013) and Astrup et al. (2014): σ = α 0e( α 1dbh)(19) Parameter θ is estimated by maximum likelihood (Marques and Buckland, 2003; Miller and Thomas, 2015; Clark, 2016) with data left-truncated at 1 m according to Astrup et al. (2014). The fitting process is carried out by means of the ds function included in the R package Distance (Miller et al., 2019). Once the parameters of detection functions are estimated, the probability of tree detection (Pi) is estimated with Eq. (20), which is implemented in the expansion factor of circular fixed area and k-tree plots as in Eq. (21): Pi=2 R2∫R 0 rg(r, θ)(20) EF =∑ n i=1 10000 (Pi π R2)(21) where R is the plot radius (m), and EF is the expansion factor. Multiplying N.tls by these EFs yields corrected stand densities (N.hn, N.hr, N. hn.cov, N.hr.cov, trees ha −1 ). Similarly, corrected stand basal area (G. hn, G.hr, G.hn.cov, G.hr.cov, m 2 ha −1 ) and volume (V.hn, V.hr, V.hn. cov, V.hr.cov, m 3 ha −1 ) are computed. All of these estimates are obtained by previously executing the distance.sampling function (Table 1), which returns Pi and also values of parameter estimates of detection functions and the corresponding Akaike information criterion (AIC) of these fits. 2.1.3.6. Correcting the shadowing effect. This approach was developed by Seidel and Ammer (2014) for single-scan mode. These authors used the approach to correct the shadowing effect, which generates shaded unsampled areas (Eq. (22)), according to the shaded area percentage related to the total area sampled: Ashadow =[( π R2)−( π r2 tree) 360◦(dbh rtree)]−⎡ ⎢ ⎢ ⎢ ⎣ π (dbh / 2 )2 2 ⎤ ⎥ ⎥ ⎥ ⎦(22) where R is the radius of the plot (m), and rtree is the distance between the TLS instrument and the tree centre. This method is implemented for circular fixed area and k-tree plots and yields an expansion factor used to compute corrected estimates for stand density (N.sh, trees ha −1 ), basal area (G.sh, m 2 ha −1 ) and volume (V.sh, m 3 ha −1 ). 2.1.3.7. Optimizing the plot design. Two different approaches can be used to find the best possible plot design depending on whether validation field data are available or not. The approaches are represented in the FORTLS workflow (Fig. 6) and detailed below. 2.1.3.8. Analysis of estimation stability. The function estimation.plot. size estimates both apparent tree density (N.tls, trees ha −1 ) and apparent basal area (G.tls, m 2 ha −1 ) for all of the aforementioned plot designs. In the case of circular fixed area plot design, concentric plots, in regular increments of 0.1 m radius (by default) to the maximum radius specified in the arguments, are simulated for computing N.tls and G.tls. As a result, line charts with estimates through plot size are obtained. For ktree design, all possible plots are defined by k ={1, 2, …, n}, where 1 is the nearest tree and n the farthest tree considered in the argument k.tree. max (or the farthest detected/existing tree if the argument is not specified). Finally, for the angle-count design, variables will be estimated for regular BAF increments comprised from 0.1 to the BAF.max specified in the arguments. All of these line charts were inspired by Fig. 3 in Brunner and Gizachew (2014). 2.1.3.9. Validation with field measurements. For cases when field data are available, we designed a set of interconnected functions able to assess the performance of processed TLS data relative to the corresponding field data: simulations, relative.bias and correlations. The field data necessary to conduct the analysis described hereinafter are tree dimensions (dbh and h) and positions relative to the TLS scanner. Analysis of the performance is based on comparisons between these two data sources for the different plot designs and sizes. The first function is simulations, which computes (in a similar way as estimation. plot.size) all of the metrics and variables aforementioned for TLS data 2.1.3and the corresponding variables based on field data (Table 2). The relative.bias function was designed for direct comparison of TLSbased estimates and field based measurements, by means of relative bias (Eq. (23)). Relative bias = 1 n∑n i=1yi−1 n∑n i=1xi 1 n∑n i=1xi (23) where xi and yi are the values of the field estimate and its TLS counterpart, respectively, corresponding to plot i for i =1, …, n. Relative bias is assessed for all the simulations to find the best possible plot design for each variable of interest. For other possible approaches apart from direct variables estimations, the package has other functions that assess the best possible plot designs according to the correlations between variable estimates from field data and metrics/variables derived from TLS. The correlations function computes both Pearson and Spearman correlation coefficients for common set of plots and all simulations and plot designs considered. For each variable of interest, this function produces the optimum correlations for all simulations. The optimize.plot.design function then produces a graphical representation of the strongest correlations for all variables of interest. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 16 best performance in all cases, h was always underestimated, reaching similar relative bias only for circular fixed area plots of around 18 m radius (S2). This again shows that height variables are systematically underestimated (Liang et al., 2016, 2018a; Krok et al., 2020). On the other hand, AutoStem™ was assessed for a stand planted with Picea sitchensis (Bong.) Carr. with a current density of 600 trees per ha (Mengesha et al., 2015). This was measured in nine randomly located plots of 15 m radius, in which all trees were measured by conventional methods and single TLS scans were made from the plot centre with a resolution of 6.28 mm at 10 m. The overall difference between TLS-derived and conventional volume estimates was 5.6% when occluded trees (not detected by TLS) were excluded from the analysis of both sources, i.e. TLS and field data, and 10.2% when estimates based on TLS were corrected by simple correction factors. FORTLS yielded much lower relative bias in volume estimates, especially for the k-tree plot design and estimates corrected with distance sampling methods (G.hr. cov), which yielded values of around 0% for 12–14 trees (S2). In any case, comparison with these findings should be done with caution, as we used paraboloid function for estimating volume, in contrast to Mengesha et al. (2015), who applied methods based on retrieved stem profiles. Nevertheless, this is the most interesting comparison because the aforementioned authors used a very similar sampling methodology, with single randomly located scans (systematic in our case) and the corresponding field plot measurements. All plot designs considered (circular fixed area, k-tree and anglecount plots) yielded stable estimates in N.tls and G.tls for certain plot size ranges (Fig. 9), as occurred for G estimates using a 2D TLS device (Fig. 3, Brunner and Gizachew, 2014). However, the variables were slightly underestimated, unlike in Brunner and Gizachew (2014), in which unbiased G estimates between real and scan basal area were observed for 5–10 m radius plot size (only for one of the studied stands). However, our findings are generally consistent with the most recent findings, which indicate underestimation of N and G due to deficit in tree detection caused by occlusions, especially for single scans (Liang et al., 2016, 2018a; Krok et al., 2020). In our study, uncorrected estimates (G. tls) yielded lower bias for smaller plot sizes, due to lower occlusion rates, which is consistent with the findings of Corona et al. (2019), who concluded that under easy-to-measure stand conditions, plots of 10 m radius in which occlusion corrections are not considered may be good enough for estimating G. To overcome underestimates derived from occlusion effects, we incorporated several methods applied in TLS single scans in other studies, which improved estimations in larger plot sizes, where G.hn yielded the best estimates for 17.5–19 m radius (Fig. 10). This finding is consistent with those of Astrup et al. (2014) who reported that a larger detection radius seems to improve estimates when distance sampling methods are applied. The results presented here indicate that correction occlusion methods can improve estimates as observed for G in angle-count (Strahler et al., 2008; Lovell et al., 2011) and circular fixed area plots (Seidel and Ammer, 2014; Astrup et al., 2014). Because these methods can be assessed in FORTLS in terms of relative bias through continuous plot size increment and different plot designs, this represents an advantage for determining the best possible plot design in execution of a single workflow. FORTLS also evaluates correlations between variables of interest and TLS-derived metrics and variables. This contribution implies a new perspective, enabling selection of the best possible plot design according to statistical correlation measures instead of measurement accuracy. This approach may be considered for estimating forestry variables assisted by or based on models, in a similar way as ABA inference developed for ALS devices (Næsset, 2002). In some cases, simple linear regressions may be fitted when strong relationships are observed for only one TLS metrics and/or variables, as between h.0 and P95 for the study case, with correlations above 0.94 (Fig. 11). In this case the plot design considered was an angle-count plot for a BAF value of 1.3. In addition, this concept provides an opportunity to solve the systematic problem regarding underestimation of height variables derived from TLS measurements (Liang et al., 2016, 2018a; Krok et al., 2020). Different groups of variables yielded the highest correlations at different plot sizes (Fig. 12); when the correlations for height variables were strongest in larger plots, diameter variables retained approximate stable correlations, and the highest correlations for other variables (N, G and V) were reached in smaller plots. Thus, the plot design could be adapted to our stand conditions and target variables for more efficient sampling. Here we have demonstrated the utility of the R package FORTLS in FIs in a case study. As FORTLS works with single scan data, coregistration of point clouds in specific software and placement of targets for field measurements are not required. This improves data acquisition and shortens the processing time, as well as enabling the sample size to be increased in a cost-efficient manner, which is one of the most desirable features of TLS in FIs (Liang et al., 2016). Further research to consolidate FORTLS for the approaches mentioned here should encompass the following: (i) larger and more complex study cases; (ii) consideration of more metrics and variables with high potential for correlation with other forest attributes (leaf area index, species, etc.); (iii) exploration of the possibility of making inferences assisted by models, by developing an adequate sampling methodology; and (iv) improvement of the computation process as much as possible, in relation to both algorithms and computing time. 5. Conclusions The R package FORTLS is useful software for processing TLS data for forestry purposes. It has the advantage of working with single scans and conducting automatic data processing, which may overcome the major challenge of affordability in data acquisition and data processing. It has yielded good results for conventional variables, based on a preliminary case study with direct estimates, as well as good correlations between field-derived variables and TLS-derived metrics and variables. However, its potential for producing model-assisted inferences from metrics and/ or variables has not yet been demonstrated. In addition, one of the most valuable features of the software is its flexibility to adapt to the best possible plot design for each variable, enabling multiple plot designs to be used in a single sampling design. Further research considering larger and more complex case studies is necessary to consolidate FORTLS as an operational tool in FIs, as well as to develop new metrics and variables. 6. Software availability Name of software: FORTLS 1.0.6 Developers: Juan Alberto Molina-Valero, María Jos´ e Ginzo Villamayor, Manuel Antonio Novo P´ erez, Adela Martínez-Calvo, Juan Gabriel ´ Alvarez-Gonz´ alez, Fernando Montes, C´ esar P´ erez-Cruzado Contact Address: Unit for Sustainable Environmental and Forest Management (UXAFORES), Department of Agroforestry Engineering, Higher Polytechnic Engineering School, Universidade de Santiago de Compostela, Benigno Ledo s/n, Campus Terra, 27002 Lugo, Spain. Email: [email protected] Software required: R ≥3.5.0 First available: March 2, 2021 Availability: https://CRAN.R-project.org/package=FORTLS Installation in R: install.packages(‘FORTLS’) Program languages: R and C++ License: GPL 3 Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 17 Acknowledgements This work was supported by the Spanish Ministry of Science and Innovation [AGL2016-76769-C2-2-R; PID2020-119204RB-C22] and Galician Regional Government [2020-CP031; ED431F 2020/02]; JAMV was supported by the Spanish Ministry of Science, Innovation and Universities through the FPU program [FPU16/03057]; AMC was supported by Galician Regional Government within the framework of the agreement “Development of the Galician continuous forest inventory” [2020-CP031]; CPC was supported by the Spanish Ministry of Science and Innovation [RYC2018-024939-I]. The authors thank Diego Lombardero Barrera, Joel Rodríguez Ruiz, Mario L´ opez Fern´ andez and ´ Oscar L´ opez ´ Alvarez for help with fieldwork. We are also grateful for the comments made by anonymous reviewers, which helped us to improve the quality of the paper. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.envsoft.2022.105337. References Astrup, R., Ducey, M.J., Granhus, A., Ritter, T., von Lüpke, N., 2014. Approaches for estimating stand-level volume using terrestrial laser scanning in a single-scan mode. Canadian J. For. Res. 44 (6), 666–676. https://doiorg/10.1139/cjfr-2013-0535. Bates, D., Eddelbuettel, D., 2013. Fast and elegant numerical linear algebra using the RcppEigen package. J. Stat. Software 52 (5), 1–24. https://doi.org/10.18637/jss. v052.i05. Bienert, A., Scheller, S., Keane, E., Mohan, F., Nugent, C., 2007. Tree detection and diameter estimations by analysis of forest terrestrial laserscanner point clouds. In: ISPRS Workshop on Laser Scanning, vol. 36, pp. 50–55 (Espoo, Finland). Bitterlich, W., 1948. Die winkelz¨ ahlprobe. Allg. Forst-und holzwirt. Zeitung Wien., 59, pp. 4–5 (Winkelz¨ ahlprobe =angle count sampling). Bivand, R.S., Pebesma, E., Gomez-Rubio, V., 2013. Applied Spatial Data Analysis with R, vol. 2. Springer, New York. Brunner, A., Gizachew, B., 2014. Rapid detection of stand density, tree positions, and tree diameter with a 2D terrestrial laser scanner. Eur. J. For. Res. 133, 819–831. https://doi.org/10.1007/s10342-014-0799-1. Buckland, S.T., Anderson, D.R., Burnham, K.P., Laake, J.L., Borchers, D.L., Thomas, L., 2001. Introduction to Distance Sampling: Estimating Abundance of Biological Populations. Oxford University Press. Cabo, C., Ord´ o˜ nez, C., L´ opez-S´ anchez, C.A., Armesto, J., 2018. Automatic dendrometry: tree detection, tree height and diameter estimation using terrestrial laser scanning. Int. J. Appl. Earth Obs. Geoinf. 69, 164–174. https://doi.org/10.1016/j. jag.2018.01.011. Calders, K., Adams, J., Armston, J., Bartholomeus, H., Bauwens, S., Bentley, L.P., Chave, J., Danson, F.M., Demol, M., Disney, M., Gaulton, R., Krishna Moorthy, S.M., Levick, S.R., Saarinen, N., Schaaf, C., Stovall, A., Terryn, L., Wilkes, P., Verbeeck, H., 2020. Terrestrial laser scanning in forest ecology: expanding the horizon. Remote Sens. Environ. 251, 112102 https://doi.org/10.1016/j.rse.2020.112102. Clark, R.G., 2016. Statistical efficiency in distance sampling. PLoS One 11 (3), e0149298. https://doi.org/10.1371/journal.pone.0149298. Corona, P., Di Biase, R.M., Fattorini, L., D’Amati, M., 2019. A Monte Carlo appraisal of tree abundance and stand basal area estimation in forest inventories based on terrestrial laser scanning. Canadian J. For. Res. 49 (1), 41–52. https://doi.org/ 10.1139/cjfr-2017-0462. Cs´ ardi, G., FitzJohn, R., 2019. Progress: Terminal Progress Bars. R package version 1.2.2. https://CRAN.R-project.org/package=progress. Danson, F.M., Disney, M.I., Gaulton, R., Schaaf, C., Strahler, A., 2018. The terrestrial laser scanning revolution in forest ecology. Interface Focus 8 (2), 20180001. https:// doi.org/10.1098/rsfs.2018.0001. Dassot, M., Constant, T., Fournier, M., 2011. The use of terrestrial LiDAR technology in forest science: application fields, benefits and challenges. Annals For. Sci. 68, 959–974. https://doi.org/10.1007/s13595-011-0102-2. Du, S., Lindenbergh, R., Ledoux, H., Stoter, J., Nan, L., 2019. AdTree: accurate, detailed, and automatic modelling of laser-scanned trees. Rem. Sens. 11 (18), 2074. https:// doi.org/10.3390/rs11182074. Dubayah, R.O., Drake, J.B., 2000. Lidar remote sensing for forestry. J. For. 98 (6), 44–46. https://doi.org/10.1093/jof/98.6.44. Ducey, M.J., Astrup, R., 2013. Adjusting for nondetection in forest inventories derived from terrestrial laser scanning. Can. J. Rem. Sens. 39 (5), 410–425. Eddelbuettel, D., 2013. Seamless R and C++ Integration with Rcpp. Springer, New York, ISBN 978-1-4614-6867-7. Eddelbuettel, D., Balamuta, J.J., 2018. Extending R with C++: a brief introduction to Rcpp. Am. Statistician 72 (1), 28–36. https://doi.org/10.1080/ 00031305.2017.1375990. Eddelbuettel, D., François, R., 2011. Rcpp: seamless R and C++ integration. J. Stat. Software 40 (8), 1–18. https://doi.org/10.18637/jss.v040.i08. Ester, M., Kriegel, H.P., Sander, J., Xu, X., 1996. A density-based algorithm for discovering clusters in large spatial databases with noise. In Kdd 96 (34), 226–231. Ferrara, R., Virdis, S.G.P., Ventura, A., Ghisu, T., Duce, P., Pellizzaro, G., 2018. An automated approach for wood-leaf separation from terrestrial LIDAR point clouds using the density based clustering algorithm DBSCAN. Agric. For. Meteorology 262, 434–444. https://doi.org/10.1016/j.agrformet.2018.04.008. Garrett, R.C., Nar, A., Fisher, T.J., 2021. Ggvoronoi: Voronoi Diagrams and Heatmaps with ’ggplot2. R package version 0.8.4. https://CRAN.R-project.org/packa ge=ggvoronoi. Hackenberg, J., Spiecker, H., Calders, K., Disney, M., Raumonen, P., 2015. SimpleTree—an efficient open source tool to build tree models from TLS clouds. Forests 6 (11), 4245–4294. https://doi.org/10.3390/f6114245. Hahsler, M., Piekenbrock, M., Doran, D., 2019. Dbscan: fast density-based clustering with R. J. Stat. Software 91 (1), 1–30. https://doi.org/10.18637/jss.v091.i01. Hester, J., Wickham, H., 2020. Vroom: Read and Write Rectangular Text Data Quickly. R package version 1.3.2. https://CRAN.R-project.org/package=vroom. Hijmans, R.J., 2020. Raster: Geographic Data Analysis and Modeling. R package version 3.4-5. https://CRAN.R-project.org/package=raster. Holopainen, M., Vastaranta, M., Hyypp¨ a, J., 2014. Outlook for the next generation’s precision forestry in Finland. Forests 5 (7), 1682–1694. https://doi.org/10.3390/ f5071682. Jin, S., Tamura, M., Susaki, J., 2016. A new approach to retrieve leaf normal distribution using terrestrial laser scanners. J. For. Res. 27, 631–638. https://doi.org/10.1007/ s11676-015-0204-z. Kleinn, C., Vilˇ cko, F., 2006. A new empirical approach for estimation in k-tree sampling. For. Ecol. Manag. 237 (1), 522–533. https://doi.org/10.1016/j.foreco.2006.09.072. Krok, G., Kraszewski, B., Stere´ nczak, K., 2020. Application of terrestrial laser scanning in forest inventory–an overview of selected issues. For. Res. Pap. 81 (4), 175–194. https://doi.org/10.2478/frp-2020-0021. Liang, X., Litkey, P., Hyyppa, J., Kaartinen, H., Vastaranta, M., Holopainen, M., 2012. Automatic stem mapping using single-scan terrestrial laser scanning. IEEE Trans. Geosci. Rem. Sens. 50 (2), 661–670. https://doi.org/10.1109/TGRS.2011.2161613. Liang, X., Kankare, V., Hyypp¨ a, J., Wang, Y., Kukko, A., Haggr´ en, H., Yu, X., Kaartinen, H., Jaakkola, A., Guan, F., Holopainen, M., Vastaranta, M., 2016. Terrestrial laser scanning in forest inventories. ISPRS J. Photogrammetry Remote Sens. 115, 63–77. https://doi.org/10.1016/j.isprsjprs.2016.01.006. Liang, X., Hyypp¨ a, J., Kaartinen, H., Lehtom¨ aki, M., Py¨ or¨ al¨ a, J., Pfeifer, N., Holopainen, M., Brolly, G., Francesco, P., Hackenberg, J., Huang, H., Jo, H., Katoh, M., Liu, L., Makros, M., Morel, J., Olofsson, K., Poveda-Lopez, J., Trochta, J., Wang, D., Wang, J., Xi, Z., Yang, B., Zheng, G., Kankare, V., Louma, V., Yu, X., Chen, L., Vastaranta, M., Saarinen, N., Wang, Y., 2018a. International benchmarking of terrestrial laser scanning approaches for forest inventories. ISPRS J. Photogrammetry Remote Sens. 144, 137–179. https://doi.org/10.1016/j. isprsjprs.2018.06.021. Liang, X., Kukko, A., Hyypp¨ a, J., Lehtom¨ aki, M., Py¨ or¨ al¨ a, J., Yu, X., Kaartinen, H., Jaakkola, A., Wang, Y., 2018b. In-situ measurements from mobile platforms: an emerging approach to address the old challenges associated with forest inventories. ISPRS J. Photogrammetry Remote Sens. 143, 97–107. https://doi.org/10.1016/j. isprsjprs.2018.04.019. Liu, J., Liang, X., Hyypp¨ a, J., Yu, X., Lehtom¨ aki, M., Py¨ or¨ al¨ a, J., Zhu, L., Wang, Y., Chen, R., 2017. Automated matching of multiple terrestrial laser scans for stem mapping without the use of artificial references. Int. J. Appl. Earth Obs. Geoinf. 56, 13–23. https://doi.org/10.1016/j.jag.2016.11.003. Lovell, J.L., Jupp, D.L.B., Newnham, G.J., Culvenor, D.S., 2011. Measuring tree stem diameters using intensity profiles from ground-based scanning lidar from a fixed viewpoint. ISPRS J. Photogrammetry Remote Sens. 66 (1), 46–55. https://doi.org/ 10.1016/j.isprsjprs.2010.08.006. Ma, L., Zheng, G., Eitel, J.U., Moskal, L.M., He, W., Huang, H., 2016. Improved salient feature-based approach for automatically separating photosynthetic and non photosynthetic components within terrestrial lidar point cloud data of forest canopies. IEEE Trans. Geosci. Rem. Sens. 54 (2), 679–696. https://doi.org/10.1109/ TGRS.2015.2459716. Marques, F.F., Buckland, S.T., 2003. Incorporating covariates into standard line transect analyses. Biometrics 59 (4), 924–935. https://doi.org/10.1111/j.0006341X.2003.00107.x. McGaughey, R.J., 2009. FUSION/LDV: software for LIDAR data analysis and visualization. In: US Department of Agriculture, Forest Service, Pacific Northwest Research Station: Seattle, WA, USA, vol. 123, 2. Mengesha, T., Hawkins, M., Nieuwenhuis, M., 2015. Validation of terrestrial laser scanning data using conventional forest inventory methods. Eur. J. For. Res. 134 (2), 211–222. https://doi.org/10.1007/s10342-014-0844-0. Miller, D.L., Thomas, L., 2015. Mixture models for distance sampling detection functions. PLoS One 10 (3), e0118726. https://doi.org/10.1371/journal.pone.0118726. Miller, D.L., Rexstad, E., Thomas, L., Marshall, L., Laake, J.L., 2019. Distance sampling in R. J. Stat. Software 89 (1), 1–28. https://doi.org/10.18637/jss.v089.i01. Molina Valero, J.A., Ginzo Villamayor, M.J., Novo P´ erez, M.A., ´ Alvarez-Gonz´ alez, J.G., P´ erez-Cruzado, C., 2019. Estimaci´ on del ´ area basim´ etrica en masas maduras de Pinus sylvestris en base a una única medici´ on del esc´ aner l´ aser terrestre (TLS). Cuad. Soc. Esp. Cienc. For. 45 (3), 97–116. https://doi.org/10.31167/csecfv0i45.19887. Molina-Valero, J.A., Ginzo Villamayor, M.J., Novo P´ erez, M.A., ´ Alvarez-Gonz´ alez, J.G., Montes, F., Martínez-Calvo, A., P´ erez-Cruzado, C., 2020. FORTLS: an R package for processing TLS data and estimating stand variables in forest inventories. Environ. Sci. Proc. 3, 38. https://doi.org/10.3390/IECF2020-08066. Molina-Valero, J.A., Ginzo Villamayor, M.J., Novo P´ erez, A.M., Martínez-Calvo, A., ´ Alvarez-Gonz´ alez, J.G., Montes, F., P´ erez-Cruzado, C., 2021. FORTLS: Automatic J.A. Molina-Valero et al. Environmental Modelling and Software 150 (2022) 105337 18 Processing of TLS Point Cloud Data for Forestry Purposes. R package version 1.0.6. https://CRAN.R-project.org/package=FORTLS. Montes, F., Rubio-Cuadrado, A., S´ anchez-Gonz´ alez, M., Aull´ o-Maestro, I., Cabrera, M., G´ omez, C., 2019. Occlusion probability in operational forest inventory field sampling with ForeStereo. Photogramm. Eng. Rem. Sens. 85 (7), 493–508. https:// doi.org/10.14358/PERS.85.7.493. Næsset, E., 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sens. Environ. 80 (1), 88–99. https://doi.org/10.1016/S0034-4257(01)00290-5. Newnham, G.J., Armston, J.D., Calders, K., Disney, M.I., Lovell, J.L., Schaaf, C.B., Sthrahler, A.H., Danson, F.M., 2015. Terrestrial laser scanning for plot-scale forest measurement. Curr. Forestry Rep. 1, 239–251. https://doi.org/10.1007/s40725015-0025-5. Olofsson, K., Holmgren, J., 2016. Single tree stem profile detection using terrestrial laser scanner data, flatness saliency features and curvature properties. Forests 7 (9), 207. https://doi.org/10.3390/f7090207. Olofsson, K., Holmgren, J., Olsson, H., 2014. Tree stem and height measurements using terrestrial laser scanning and the RANSAC algorithm. Rem. Sens. 6 (5), 4323–4344. https://doi.org/10.3390/rs6054323. Othmani, A., Voon, L.F.C.L.Y., Stolz, C., Piboule, A., 2013. Single tree species classification from terrestrial laser scanning data for forest inventory. Pattern Recogn. Lett. 34 (16), 2144–2150. https://doi.org/10.1016/j.patrec.2013.08.004. Pauly, M., Gross, M., Kobbelt, L.P., 2002. Efficient simplification of point-sampled surfaces. In: IEEE Conference on Visualization, pp. 163–170. https://doi.org/ 10.1109/VISUAL.2002.1183771. Boston, USA. Pebesma, E.J., Bivand, R.S., 2005. Classes and methods for spatial data in R. R. News 5 (2), 9–13. https://CRAN.R-project.org/doc/Rnews/. R Core Team, 2021. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/. Raumonen, P., Kaasalainen, M., Åkerblom, M., Kaasalainen, S., Kaartinen, H., Vastaranta, M., Holopainen, M., Disney, M., Lewis, P., 2013. Fast automatic precision tree models from terrestrial laser scanner data. Rem. Sens. 5 (2), 491–520. https://doi.org/10.3390/rs5020491. Roussel, J.R., Auty, D., 2020. Airborne LiDAR Data Manipulation and Visualization for Forestry Applications. R package version 3.0.4. https://cran.r-project.org/ package=lidR. Roussel, J.R., Auty, D., Coops, N.C., Tompalski, P., Goodbody, T.R.H., Sanchez Meador, A., Bourdon, J.F., de Boissieu, F., Achim, A., 2020. lidR: an R package for analysis of Airborne Laser Scanning (ALS) data. Remote Sens. Environ. 251, 112061 https://doi.org/10.1016/j.rse.2020.112061. Seidel, D., Ammer, C., 2014. Efficient measurements of basal area in short rotation forests based on terrestrial laser scanning under special consideration of shadowing. iFor. Biogeosci. For. 7 (4), 227–232. https://doi.org/10.3832/ifor1084-007. Sievert, C., 2020. Interactive Web-Based Data Visualization with R, Plotly, and Shiny. CRC Press, ISBN 9781138331457. https://plotly-r.com. Strahler, A.H., Jupp, D.L.B., Woodcock, C.E., Schaaf, C.B., Yao, T., Zhao, F., Yang, X., Lovell, J., Culvenor, D., Newnham, G., Ni-Miester, W., Boykin-Morris, W., 2008. Retrieval of forest structural parameters using a ground-based lidar instrument (Echidna®). Can. J. Rem. Sens. 34 (Suppl. 2), S426–S440. https://doi.org/10.5589/ m08-046. Tomppo, E., Gschwantner, T., Lawrence, M., McRoberts, R.E., Gabler, K., Schadauer, K., Cienciala, E., 2010. National forest inventories. Pathways for common reporting. European Science Foundation 1, 541–553. https://doi.org/10.1007/978-90-4813233-1. Trochta, J., Krůˇ cek, M., Vrˇ ska, T., Kr´ al, K., 2017. 3D Forest: an application for descriptions of three-dimensional forest structures using terrestrial LiDAR. PLoS One 12 (5), e0176871. https://doi.org/10.1371/journal.pone.0176871. Vaidyanathan, R., Xie, Y., Allaire, J.J., Cheng, J., Sievert, C., Russell, K., 2020. Htmlwidgets: HTML Widgets for R. R package version 1.5.3. https://CRAN.R-project .org/package=htmlwidgets. West, G.B., Brown, J.H., Enquist, B.J., 1999. The fourth dimension of life: fractal geometry and allometric scaling of organisms. Sci 284 (5420), 1677–1679. https:// doi.org/10.1126/science.284.5420.1677. White, J.C., Coops, N.C., Wulder, M.A., Vastaranta, M., Hilker, T., Tompalski, P., 2016. Remote sensing technologies for enhancing forest inventories: a review. Can. J. Rem. Sens. 42 (5), 619–641. https://doi.org/10.1080/07038992.2016.1207484. Wickham, H., 2021. Tidyr: Tidy Messy Data. R package version 1.1.3. https://CRAN.R-p roject.org/package=tidyr. Wickham, H., Seidel, D., 2020. Scales: Scale Functions for Visualization. R package version 1.1.1. https://CRAN.R-project.org/package=scales. Wickham, H., Hester, J., Chang, W., 2021. Devtools: Tools to Make Developing R Packages Easier. R package version 2.4.2. https://CRAN.R-project.org/pack age=devtools. Wulder, M.A., White, J.C., Nelson, R.F., Næsset, E., Ørka, H.O., Coops, N.C., Hilker, T., Bater, C.W., Gobakken, T., 2012. Lidar sampling for large-area forest characterization: a review. Remote Sens. Environ. 121, 196–209. https://doi.org/ 10.1016/j.rse.2012.02.001. Xia, S., Wang, C., Pan, F., Xi, X., Zeng, H., Liu, H., 2015. Detecting stems in dense and homogeneous forest using single-scan TLS. Forests 6 (11), 3923–3945. https://doi. org/10.3390/f6113923. Zhang, W., Qi, J., Wan, P., Wang, H., Xie, D., Wang, X., Yan, G., 2016. An easy-to-use airborne LiDAR data filtering method based on Cloth simulation. Rem. Sens. 8 (6), 501. https://doi.org/10.3390/rs8060501. Zhang, W., Wan, P., Wang, T., Cai, S., Chen, Y., Jin, X., Yan, G., 2019. A novel approach for the detection of standing tree stems from plot-level terrestrial laser scanning data. Rem. Sens. 11 (2), 211. https://doi.org/10.3390/rs11020211. J.A. Molina-Valero et al.