Fine-scale three-dimensional modeling of boreal forest plots to improve forest characterization with remote sensing

Remote Sensing of Environment 219 (2018) 99–114

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Fine-scale three-dimensional modeling of boreal forest plots to improve forest characterization with remote sensing

T

⁎

Jean-François Côtéa, , Richard A. Fournierb, Joan E. Lutherc, Olivier R. van Lierd a

Natural Resources Canada, Canadian Forest Service, Canadian Wood Fibre Centre, Quebec, QC G1V 4C7, Canada Department of Applied Geomatics, Centre d'applications et de recherche en télédétection, Université de Sherbrooke, Sherbrooke, QC J1K 2R1, Canada c Natural Resources Canada, Canadian Forest Service, Atlantic Forestry Centre, Corner Brook, NL A2H 6J3, Canada d Natural Resources Canada, Canadian Forest Service, Canadian Wood Fibre Centre, Corner Brook, NL A2H 6J3, Canada b

A R T I C LE I N FO

A B S T R A C T

Keywords: 3D architectural model Terrestrial LiDAR Tree structure Surrogate forest plots Airborne LiDAR

Improving the quality of information that can be obtained from forest inventories can enhance planning for the best use of forest resources. In this study, we demonstrate the capability to improve the characterization of forest inventory attributes using terrestrial laser scanner (TLS) data, a fine-scale architectural model (L-Architect), and airborne laser scanner (ALS) data. Terrestrial laser scanning provides detailed and accurate three-dimensional data and has the potential to characterize forest plots with comprehensive structural information. We use TLS data and in situ measurements as input to L-Architect to create reference plots. The use of L-Architect for modeling was validated by comparing selected attributes of the reference plots with validation plots produced using simulated TLS data, with normalized root-mean square error (NMRSE) values below 17%. Surrogate plots were then created using a library of tree models where individual trees were selected according to three attributes—tree height, diameter at breast height, and crown projected area—either measured from in situ plots or derived from ALS data. The accuracy of the surrogate plots was assessed by comparing several key forest attributes from the reference plots, including branching structure (e.g., number of whorls, knot surface), crown shape and size (e.g., base height, asymmetry), heterogeneity (e.g., lacunarity, fractal dimension), tree volume, and the spatial distribution of material (e.g., Weibull fit, leaf area index). Overall, the surrogate plots reproduced the attributes of the reference plots with NRMSE mean value of 17% (R2 = 0.68) using in situ ground measurements and 24% (R2 = 0.51) using inputs estimated with ALS. Some attributes, such as leaf area index, knot surface, and fractal dimension, were well predicted (R2 > 0.80), whereas others, like crown asymmetry and lacunarity, had weak correspondence (R2 < 0.16). The ability to create surrogate forest plots with L-Architect makes it possible to estimate detailed structural attributes that are difficult to measure with conventional forest mensuration techniques and that can be used for model calibration with above-canopy remote-sensing data sets.

1. Introduction Accurate and up-to-date information regarding the state, health, and development of forests is needed to quantify and predict, among other things: (i) the canopy–atmosphere exchange of material and

carbon (C) fixation, (ii) the spatial dynamics of cover from forested to nonforested resulting from disturbance or fire, (iii) the radiation regime within the canopy, (iv) the biophysical and biochemical properties of forests, (v) the dynamics of biodiversity and habitats, and (vi) the available resources for the timber industry to support sustainable

Abbreviations: ASYM, crown asymmetry; BR, branchiness ratio of the largest branch diameter to the diameter at breast height; Btot, total number of branches; CBH, crown base height; CHM, canopy height model; CPA, maximum crown projected area; CRMSE, centered root mean square error; CSA, canopy surface area; CSV, canopy surface volume; CW, maximum crown width; DBH, diameter at breast height; Dδ, fractal dimension; FC, fraction cover; HCPA, height at maximum crown projected area; HGT, total height; HLC, height to live crown; Ksurf, total knot surface on the main stem; LA, total leaf area; LAI, leaf area index; NRMSE, normalized root mean square error; PAI, plant area index; R, correlation coefficient; R2, coefficient of determination; RMSE, root mean square error; Tcoef, stem taper coefficient; VOL, crown volume; Wtot, total number of whorls; Wαleaf, scale parameter of the leaf area Weibull distribution; Wαwood, scale parameter of the wood volume Weibull distribution; Wβleaf, shape parameter of the leaf area Weibull distribution; Wβwood, shape parameter of the wood volume Weibull distribution; Λ, lacunarity index of tree and plot canopy structure; ΛCHM, lacunarity index of the canopy height model; σCHM, standard deviation of the canopy height model ⁎ Corresponding author. E-mail addresses: [email protected] (J.-F. Côté), [email protected] (R.A. Fournier), [email protected] (J.E. Luther), [email protected] (O.R. van Lier). https://doi.org/10.1016/j.rse.2018.09.026 Received 1 May 2017; Received in revised form 20 September 2018; Accepted 30 September 2018 0034-4257/ Crown Copyright © 2018 Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/).

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Moreover, it is not easy to distinguish between wood and foliage material from the return signal with current instruments (Béland et al., 2014). Thus, even when scans from different points of view are available, identification of structural elements and construction of a topologically and geometrically correct 3D structure relies on architectural quantitative modeling (Côté et al., 2011; Hackenberg et al., 2014; Raumonen et al., 2013). As an exact representation of the forest structure is not possible in practice, 3D architectural models provide convenient simplifications of forest structure at tree, branch (Cescatti, 1997; Livny et al., 2011), and down to the individual conifer shoot or leaf level (Pirk et al., 2012; Potapov et al., 2016; Widlowski et al., 2014; Xu et al., 2007). Representing tree structure is particularly challenging when dealing with mature trees in complex environments (Palubicki et al., 2009; Runions et al., 2007). In such cases, 3D architectural models have to deal with a high level of competition (between and/or intra species) for space, light and nutrient resources, which increases the number of unknown variables required to reproduce realistic forest scenarios. Several mathematical formalisms are available to express tree architecture (Godin, 2000; Godin and Caraglio, 1998; Prusinkiewicz and Lindenmayer, 1990). Among them, is the Open L-Systems formalism proposed by Měch and Prusinkiewicz (1996) to control branch growth and foliage addition within the tree crown. Based on Open L-Systems, the L-Architect (LiDAR to tree Architecture) model was developed to address limitations inherent in the TLS point cloud to extract the stem and main branches and then to construct the fine branching structure with the addition of foliage (Côté et al., 2009, 2011, 2012). L-Architect was further improved and validated using structural measurements for two coniferous species found in Newfoundland (Canada), namely, Abies balsamea (L.) Mill. (balsam fir or hereafter referred to as fir) and Picea mariana (Mill.) B.S.P. (black spruce or hereafter referred to as spruce) (Côté et al., 2013). L-Architect is thus capable of simulating fir and spruce trees, but procedures to simulate “surrogate” plots remain to be developed. Surrogate plot simulation with L-Architect provides fine-scale structural attributes of all the trees of the plot. This method has the potential to improve on the development of empirical relationships (Mahoney et al., 2018), the calibration of growth models (Falkowski et al., 2010) and the retrieval of biophysical properties (GastelluEtchegorry et al., 2015) with remote sensing observations in support of large-area mapping of forest attributes. Consequently, the general objective of this study was to develop and validate procedures for creating surrogate plots with L-Architect. Specific objectives were to (i) create a series of “reference plots” using trees scanned with TLS, (ii) validate LArchitect's ability to simulate plot-level attributes using an independent set of “validation plots”, (iii) create “surrogate ground plots” with a library of trees that are not necessarily collected at the plot and (iv) expand the application of L-Architect to generate “surrogate ALS plots” from above-canopy remote-sensing data. For clarity, the reference plots – created using trees scanned with TLS – provide the base for evaluating the capability to generate detailed tree- and plot-level attributes with LArchitect. The validation plots – generated with simulated TLS data – serve to validate that L-Architect produces results close to reality. The surrogate plots – generated with just a few key attributes measured in situ or simulated with ALS data – generalize spatially the simulations from either the ground plots or the ALS data and provide the fine-scale calibration data.

management. Remote sensing has proven to be a useful tool to monitor forest environments by providing repetitive measurements within a wide range of spatial, spectral, radiometric, and temporal resolutions. Interpretation of remote-sensing data requires the use of numerical methods and quantitative models. These models explain the nature of the measured physical signal (Chen and Leblanc, 1997; Disney et al., 2000; Gastellu-Etchegorry et al., 1996; Govaerts and Verstraete, 1998; Widlowski et al., 2006), characterize the state of the system under observation (Ceccato et al., 2002; Myneni et al., 2002; Pinty et al., 2006; Verstraete et al., 1996) or quantify empirical relationships between the variables of interest and the remote measurements (Koch, 2010; Luther et al., 2014; van Leeuwen and Nieuwenhuis, 2010; Wulder et al., 2004, 2012). The complexity of the three-dimensional (3D) forest structure has largely contributed to the use of empirical models because of the difficulty in defining an exact estimate of the state variables—such as structural attributes—of the observed system. Progress towards the development of models linking remote-sensing data with structural attributes or other derived variables largely depends on the ability to characterize the structure of forest stands with proper ground data. Although some forest structural attributes can be measured directly in the field (e.g., stem diameters and tree heights), practical constraints prohibit the measurement of fine-scale structure. The number of trees or plots that can be measured is often limited by transportation issues for remote locations, practical aspects of working in harsh environments, and required financial/logistical resources. Conventional forest inventories focus on measuring several key structural attributes at the tree and plot levels, mostly related to wood volume (Avery and Burkhart, 2002; Kangas, 2010). Measured plot-level attributes usually consist of species composition, an approximation of height and forest cover, characterization of the soil, drainage, site potential, and in some cases, a description of the understory. Measured tree-level attributes usually consist of the species, stem diameter at breast height (DBH), height, and qualitative estimates of crown quality, vitality, and health. These attributes are useful for volume estimation but they are limited in their ability to characterize fine-level structures. Generally, the measurement of fine-scale structural elements of forest canopies requires destructive sampling with intensive field work (Champion et al., 2001; Landry et al., 1997). Measuring fine-scale structural attributes requires such intensive measurement that it is rarely undertaken, and if so, it only includes a few plots because of resource limitations (Calders et al., 2015; Gonzalez de Tanago et al., 2017; Stovall et al., 2017). Therefore, alternative methods are required to estimate fine-scale tree- and plotlevel structural attributes. Light detection and ranging (LiDAR) systems provide 3D returns from canopy elements. In particular, airborne laser scanners (ALS) are well suited to characterize the shape of the canopy surface. Stand biomass, height, and the vertical distribution of forest structure are now routinely estimated using ALS data (Bouvier et al., 2015; Kankare et al., 2013; van Leeuwen and Nieuwenhuis, 2010; Wulder et al., 2012). These estimates of structural attributes characterize forest ecosystem structure, diversity, and function at finer spatial and temporal scales than previously possible. Terrestrial laser scanners (TLS) collect large amounts of 3D data on the fine-scale structure of trees and forest stands (Dassot et al., 2011). The detailed 3D tree and canopy structure information obtained from TLS is invaluable for the validation of abovecanopy remote sensing measurement and derived products (Calders et al., 2018; Widlowski et al., 2015). For each TLS scan, millions of returned signals are registered, providing positional information on all elements surrounding the instrument. Despite the potential for TLS to generate high quality data, these data can be unreliable when acquired in natural forest environments. The main problem for a complete measurement of a specific scene is the occlusion (often referred to as shadowing) of the laser pulse. Groups of clustered or opaque objects prevent detection of elements in occluded areas (Hopkinson et al., 2004; van der Zande et al., 2006, 2008; Watt and Donoghue, 2005).

2. Material 2.1. In situ measurements The study made use of a network of permanent sample plots (PSPs) located in Newfoundland, Canada (Fig. 1) and representing the eastern extent of the boreal forest region of North America (Rowe, 1972). Permanent sample plots are routinely measured by the provincial forest management service of Newfoundland and Labrador (Newfoundland 100

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Fig. 1. Location of permanent sample plots created with the architectural model.

Forest Service, 2011). A sample of 42 fir- and 35 spruce-dominated plots was selected to represent the range of species, height, crown density, and site quality classes present for the merchantable forest area of the island, whereby the dominant species made up > 25% of the basal area of the plot. By design, the operational plots were of variable sizes according to stem density for immature and semimature stand types and fixed at 0.04 ha for mature and over-mature types. All plots were rectangular in shape, with a width of 14 m. The center and four corner locations of each plot were measured using a Trimble® GeoXH™ global positioning system (GPS). The GPS locations were postprocessed to provide plot locations with submeter accuracy. At each plot, the species and DBH (measured at 1.3 m) were recorded for each tree above a specified size according to plot maturity. Height (HGT) was measured for a sample of trees, and species-specific relationships between DBH and HGT were developed according to Damman type (Meades and Moores, 1994), and used to predict heights for all trees. Maps of tree locations were created using, first, the precise GPS location of the PSP

central point. Then, the distance and azimuth of all trees in the plot were measured from the PSP center using a Trimble® LaserAce™ 1000 rangefinder. The (x,y,z) position of each measured tree was calculated using the GPS location, distance, and azimuth information. 2.2. Terrestrial laser scanner (TLS) data TLS data were collected as part of a coincident study on wood quality within the Newfoundland Fibre Project (Blanchette et al., 2015; Lessard et al., 2014; Luther et al., 2014). A total of 229 individual trees from the 77 plots were scanned with the Ilris-3D (Optech Inc., Vaughan, Ontario) in 2009 and the Imager 5006i (Zoller & Fröhlich, Wangen, hereafter called Z+F) in 2010. The Ilris-3D emits light at 1500 nm and scans with a maximum viewing window of 40° × 40°, which can be adjusted to fit target objects. This system uses the time-of-flight to measure distance by calculating time delay between the emission and the reception of signal backscattering intensity in either first or last 101

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

based on shadow propagation (Palubicki et al., 2009) computes a coarse estimate of light exposure for each bud. Light availability serves to determine which buds could produce new shoots using a priority resource allocation model that operates at the level of the entire tree (Palubicki et al., 2009). The total number of shoots is made to fit the total leaf area used as input. Branch shedding is simulated by calculating the total amount of resources gathered by a branch and the number of internodes it supports to decide if the branch is considered a liability for the tree. Finally, branch diameters are estimated following a pipe modeling approach depending on the number of children a parent branch supported (Xu et al., 2007). For this study, conifer shoot models were created using a generic description of the shape and distribution of the needles as well as field photographs for the selected species in a similar environment (for fir: Morris (1951); Piene (1983), and for spruce: Fraser (1966); Colombo (1986); Weng and Jackson (2000)). For fir, each shoot contained 150 needles of 2 cm length and 1 mm of diameter distributed over a 10-cm twig. The needle angle for insertion in the twig was set to 60°. The fascicle angle between pairs of needles was set to 60°. The needle angle distribution around the axis followed a Fibonacci sequence with a divergence angle of 180°. For spruce, each shoot contained 200 needles of 1.5 cm length and 1 mm of diameter distributed over a twig of 10 cm. The needle angle for insertion in the twig was set to 45°. The fascicle angle between pairs of needles was set to 45°. The needle angle distribution around the axis followed a Fibonacci sequence with a divergence angle of 8/13*2π. For both species, the twig diameter varied with the age of the shoot up to a maximum of 5 years of needle retention. Twig diameters were assigned at 3, 4.6, 7.2, 11.2 and 17.6 mm for age 1, 2, 3, 4, and 5 years respectively. A pipe model assuming three children per shoot with coefficient of 2.49 was used to calculate the twig width. Plot creation was accomplished using the L-Architect model with in situ measurements of DBH and HGT (Fig. 2). For each tree in the reference plot, the scanned tree that best fit the measurements was selected by: (1) matching recorded characteristics, namely the same

return mode. To reduce the negative impact of signal occlusion, scans were taken with the Ilris-3D from two to three viewpoint locations around a target tree using the first-return mode. The trees were scanned at a distance ranging from 10 to 15 m, giving a point spacing (footprint) of 5.0 to 7.5 mm (6.85 to 7.27 mm) respectively. The Z+F emits light at 680 nm and scans with a near complete spherical viewing window of 310° (vertical) × 360° (horizontal). The system uses the phase-shift technology to measure distance for a maximum range of 80 m. The very wide field of view captured a minimum of 10 live merchantable trees and their surrounding environment. The filtering procedure available in the LaserControl software from Z+F was applied to all point cloud data acquired by the Z+F TLS to remove points resulting from noise or erroneous returns. Typically, a total of five scans were performed within and outside the tree group area to minimize occlusion: four scans located outside (~ North, South, East and West cardinal points) and one at the center. The spatial resolution of the scans corresponded approximately to one point every 6.3 mm and a footprint of 4.1 mm on a sphere at 10 m. The 3D point clouds taken either by the Ilris-3D or the Z +F from different viewpoints were aligned into one geometric coordinate system with the software Pointstream (3DImageSuite®) and Z +F LaserControl respectively. Pointstream was then used to process all coregistered point clouds to select and extract the individual trees used as input for L-Architect. 3. Methods This study uses L-Architect to model forest plots. Detailed descriptions of L-Architect are available in Côté et al. (2011, 2013). First, reference plots were created with L-Architect using in situ measurements including TLS scans of individual trees. Then, a variety of structural attributes were calculated from the reference plots to be used as a fully controlled data set for validating the modeling approach. The modeling approach was validated on these attributes using independent (validation) plots generated with L-Architect using simulated TLS data. Finally, L-Architect was used to produce surrogate plots with just a few key attributes estimated in situ or with simulated ALS data. L-Architect was chosen for its capacity to reproduce the two coniferous species, fir and spruce, in the mature boreal forests of Newfoundland. However, the generic framework of the proposed approach could be easily adapted for any quantitative model of trees validated for a specific forest ecosystem. 3.1. Creation of reference plots with L-Architect L-Architect reproduces realistic tree structures based on TLS data and allometric relationships that define the total amount of foliage following two main steps: branch growth and addition of foliage within the crown. The Open L-Systems formalism (Měch and Prusinkiewicz, 1996) controls the two construction steps. Specific rules for simulating the tree development are implemented at the bud level. Two processes are implemented within an Open L-System framework to construct models of realistic trees in natural forests: (i) determining space availability for branch growth, and (ii) estimating light exposure for supporting foliage addition with a resource allocation model. The branching structure uses specific rules for tree development in conjunction with a space colonization algorithm (Runions et al., 2007) to take into account the close environment of the branch. The algorithm assumes that each bud is surrounded by a spherical occupancy zone and a conical perception volume. The set of points from the TLS coregistered scans, called attractors, represents the space available for tree growth. At every step, the algorithm calculates the optimal growth direction of a bud according to the position of the attractors in its perception volume, and deletes attractors in its occupancy zone. This step is repeated to generate the branching structure of the tree. Foliage addition starts at the end of the branching structure development if the tree is considered alive (no foliage otherwise). A method

Fig. 2. Method to produce reference plots using L-Architect to simulate detailed tree- and plot-level attributes. 102

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

species, site, and crown density class and (2) minimizing the following constraint between the scanned tree (subscript ‘s’) and the measured tree (subscript ‘m’):

plant and dividing by the plot area. The vertical distributions of leaf area and wood volume were characterized by fitting a two-parameter Weibull function:

|DBHm − DBHs | |HGTm − HGTs | + . DBHm HGTm

W (x ; α, β ) =

β

(1)

Dδ = lim

ln N (δ )

δ→0

3.2. Calculation of structural attributes

Tcoef

∑ i=1

di + 1 − di , hi + 1 − hi

1

ln δ

, (4)

where δ is the size of the voxels and N is the number of voxels (of dimension δ) intercepted. Dδ was estimated from the slope of the re1 gression between lnN(δ) and ln δ . The lacunarity index of the volume distribution was calculated with a modified gliding box algorithm for 3D gray images (Frazer et al., 2005; Kirkpatrick and Weishampel, 2005; Plotnick et al., 1996), i.e., the 3D arrays of voxels containing values of the total volume of material. Lacunarity defines a scale-dependent deviation of a geometrical object from translational invariance or homogeneity. Geometrical patterns with low lacunarity are finely textured and homogeneous with respect to the size distribution and spatial arrangement of gaps. In contrast, patterns with high lacunarity show substantial heterogeneity in the size distribution and spatial arrangement of gaps. The lacunarity Λ(δ) was calculated for a range of scales, or sizes, δ identical to the ones used for the fractal dimension. The box of size δ slides iteratively one (10 cm) voxel unit at a time through the 3D array and each box value was set to the total volume of material V within the volume δ3. With the first and second statistical moments, Z(1)(δ) and Z(2)(δ), and the probability distribution Q(V, δ), the lacunarity at box of size δ was defined as:

Three types of structural attributes were calculated from the reference plots: tree-level attributes, plot-level attributes, and indicators of heterogeneity. Calculating such a wide variety of attributes was made possible by the three formats adopted for spatial representation of the tree elements, namely (1) the GEOM modeling language (GML; Boudon et al., 2003), (2) the multiscale tree graph (MTG; Godin and Caraglio, 1998), and (3) the volume elements (voxels) representation. The first two formats are supported by the PlantGL open-source graphic toolkit for the creation, simulation, and analysis of 3D virtual plants (Pradal et al., 2009). The third format implied the spatialization of the geometric representation of the tree elements within a 3D array of 10 cm × 10 cm × 10 cm voxels using a slicing based method for computation of volume fractions in multimaterial solids (Sen and Srikanth, 2008). The tree-level attributes were calculated by traversing the MTG structure. It included the total number of branches (Btot) and whorls (Wtot), the maximum crown width (CW), the DBH, the total height (HGT), the height to live crown (HLC), the total knot surface (Ksurf) on the main stem, the total leaf area (LA), the branchiness ratio of the largest branch diameter to the DBH (BR) and the stem taper coefficient (Tcoef) calculated as (Lenz et al., 2012): K

(3)

where α and β are the scale and shape parameter. W(x; α, β) is zero for x < 0. Wαleaf, Wβleaf, Wαwood, Wβwood therefore characterize the vertical distribution of leaf area and wood volume, respectively. Fractal dimension and lacunarity were calculated to quantify the heterogeneity of the plot structure. Fractal dimension Dδ is useful for studying irregularity: it characterizes the way objects such as trees and forest canopies occupy space. The box-counting method has been used to estimate fractal dimension of objects (Pradal et al., 2009; Zeide and Pfeifer, 1991). The adaptation of this method to 3D scenes composed of voxels involved building a sequence of 3D arrays of voxels of decreasing size, and for each sequence to identify and count which voxels composed the object. The estimator of the fractal dimension of the object was defined as:

Once a scanned tree from the library was chosen, a scale factor was applied to each point of the TLS scan to obtain the same HGT of each tree found in the plot. Then, DBH, HGT, and total leaf area (LA) were given as inputs to L-Architect to create a reference tree. The LA was calculated using allometric relationships (≡f(DBH, HGT)) developed specifically for fir and spruce trees in Newfoundland (Lavigne et al., 1996; Newton, 2006). The parameterization of the variables used in the different algorithms and calculations followed the averaged values found in Côté et al. (2013; ref. Table 3) for the immature and mature tree samples used for validation. Each reference plot was created by positioning each reference tree in the reference plot according to the tree map location and rotating each tree randomly in azimuth to avoid spurious frequencies in the scene. These reference plots are the optimal structural representation of the plot based on the available in situ measurements.

1 = K

β x β − 1 −⎛ β α ⎞ ⎛ ⎞ e ⎝ ⎠ , x ≥ 0, α ⎝α⎠

Λ(δ ) =

Z (2) (δ ) σ 2 (δ ) =1+ , (1) 2 (Z (δ )) (Z (1) (δ ))2

(5)

where σ (δ) is the sample variance and Z (δ) is the mean of the probability function. The lacunarity index was taken as the normalized lacunarity statistic integrated across all scales (Frazer et al., 2005): 2

(2)

Λ=

th

where di and hi are the stem diameter and height of the i whorl (zascending ordered). Additional tree-level attributes were calculated by using the 3D array of voxels. The crown base height (CBH), the maximum crown projected area (CPA) and the height at CPA (HCPA) were calculated as in Seidel et al. (2011). The crown volume (VOL) was estimated by summing the projected area of each vertical layer of 10 cm in height starting from CBH. The crown asymmetry (ASYM) was estimated by calculating the horizontal distance between the center of the crown at the height of maximum crown projection area and the stem location at ground level (Seidel et al., 2011). The plot-level attributes were calculated from the 3D array of voxels using the same algorithm as used for trees (Sen and Srikanth, 2008). The fraction cover (FC) was defined as the ratio of 10 cm × 10 cm cells containing material over the total number of cells covering the projected area of the PSP. The leaf area index (LAI) and plant area index (PAI) were also calculated by summing the surface area of foliage or

1 Λ(1)

(1)

∑ Λ(δ )

(6)

Similar procedures were used to calculate structural attributes from modeled TLS plots in order to validate the structural consistency of the plot creation method. 3.3. Calculation of error statistics The use of L-Architect for modeling trees and plots was validated by comparing selected attributes of the reference plots with modeled plots: i.e., plots produced by an alternative modeling approach using simulated TLS data. The assessment considered the coefficient of determination (R2) and the normalized root mean square error (NRMSE) defined as: N

NRMSE =

103

∑n = 1 (x nmod − x nref )2 / N ref ref x max − x min

(7)

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

where N is the number of observations, and xnref, xmaxref, xminrefthe nth, maximal and minimal reference value of an attribute x; xnmod is the nth corresponding value calculated from the modeled trees or plots. Lower values of NRMSE indicate less residual variance. In addition to R2 and NRMSE, Taylor diagrams (Taylor, 2001) were used as concise graphical statistical summaries to compare tree- or plotlevel attributes from the reference plots (ref) and modeled plots (mod). The statistics selected were the correlation (r), the standard deviation (σ), ratio of dispersion σmod/σref and the centered RMSE (CRMSE) defined as:

trees, (2) create the validation plots using simulated TLS scans for selected trees from the reference plots, and (3) assess how the reference and validation plots compare. The 20 reference plots were selected randomly from the available pool of 77 plots from field work, making sure that 10 plots were fir dominated and the other 10 plots were spruce dominated. Within each of these 20 plots, five reference trees were randomly selected: 51 firs and 49 spruces. These 100 trees produced with L-Architect were then scanned by a TLS simulator. Terrestrial laser scans were simulated with the physically based ray-tracing renderer PBRT (Pharr and Humphreys, 2010) in order to create a library of 100 simulated TLS trees. The TLS instrument was modeled as a perspective camera with a 90° field of view operating at an angular resolution of 0.072° (equivalent to a 1250 × 1250 pixels image). The laser source was simulated as a spot light with a 90° field of view at the same location and pointing in the same direction as the camera. Two scans were generated for each tree from opposite directions (back and front of the tree). The distance between the tree and the instrument was equivalent to half the total height, which ensured capturing the full tree. The two scans were then coregistered. The tree-level structural attributes of the original 100 reference trees, selected for the library, and their simulated TLS scans were compared to validate the tree construction method. This tree-level comparison was considered as a test of internal tree model consistency, as was done in Côté et al. (2009). In addition to compare validation and reference plots, surrogate plots constructed with L-Architect and either in situ or ALS data were also compared with the reference plots at both tree and plot levels. As part of the validation process, a sensitivity analysis was performed on the number of trees used for modeling plots. As suggested in Côté et al. (2012), improvements to the accuracy of modeled attributes can be gained if the library of coregistered tree scans contains several scans for each species. The optimal number and distribution of tree scans in the library may vary according to the complexity of the environment: number of tree species, DBH, HGT distribution, etc. Here, the impact of having a reduced number of scans was assessed. The sensitivity analysis of the validation plot construction was performed by successively removing 25%, 50%, and 75% of the scans at random from the simulated TLS coregistered scan library. Error statistics were used to assess at what point a reduced library of scanned trees would suffice to maintain the high level of structural detail at both tree and plot levels.

N

CRMSE =

∑ [(xnmod − μmod ) − (xnref

− μref )]2 / N

n=1

(8)

where μref is the mean value of the attributes of the reference and, μmod is the corresponding values from the attributes of the plots modeled with the alternative approach. The Taylor diagram represents three different statistics simultaneously based on the similarity of the relationships between these quantities and the law of cosines. These statistics make it easy to determine how much of the overall difference is attributable to a difference in variance and how much is due to poor pattern correlation. 3.4. Validation of the plot creation method Options to validate plot reconstructions with L-Architect are limited, and they are far from trivial. Field destructive measurements are the only mean to obtain true validation of trees and plot attribute estimates. A field validation was accomplished at the tree level in Côté et al. (2013) with destructive sampling of six trees. Typically, an exhaustive sampling of individual tree structure takes one to three days per person per tree, depending on tree size and the level of detail required for the assessment (Côté et al., 2013; Landry et al., 1997). At the plot level, accurate measurement of geometric attributes repeated over hundreds of thousands structural units requires significant financial and time resources. Furthermore, these destructive measurements are simply not achievable for some attributes such as branch orientation, canopy gaps, or material distribution. In addition, measurement errors from field estimates are also a factor to take into account for a reliable evaluation. Another validation strategy was to generate point clouds from realistic 3D scenes of forest with a TLS simulator. In this way, we made use of the structural details specified by the reconstruction method to assess if the results are consistent with the original scene. However, Disney et al. (2018) do not recommend testing a canopy architecture model using simulated point clouds derived from an equivalent reconstruction method. This quantitative self-consistency test defines model-to-model deviations that can be circumvented. Another approach would be to compare a series of independent architectural models, which is a laborious effort that may be left for a model cook-off effort on its own; similar to the RAMI exercises for radiative transfer models (Widlowski et al., 2015). Instead, we adopted a strategy introduced by Côté et al. (2009) for tree level validation in a virtual environment, and adapted it for plot-level simulations. We simulated the measurement configuration of the TLS scans with a physically based ray-tracing algorithm on a subset of independent randomly chosen reference tree models. This resulted in realistic (imperfect) point clouds of individual trees to drive the plot creation method. Then we produced “second generation” tree models in the validation plots using the simulated TLS scans. In this way, the structural attributes of the validation plots could be compared in an unambiguous manner with their respective reference plots. Whilst not being a full comprehensive comparison (Calders et al., 2015; Gonzalez de Tanago et al., 2017; Stovall et al., 2017), this process provided a sound self-consistency check for the plot creation method. Validation of the plot creation method therefore involved three steps: (1) using L-Architect, produce a series of 20 reference plots using tree-level in situ measurements and a library of TLS scans of individual

3.5. Creation of surrogate plots Contrary to reference plots, surrogate plots provide forest plot data with a high level of detail where TLS scans are not available. Surrogate plots can be used with “wall-to-wall” remote-sensing data to develop models for mapping forest structural attributes for a complete survey area. In this study, the potential to produce surrogate plots with LArchitect was evaluated (i) using in situ measurements and (ii) using simulated ALS data and the software TIFFS v8.0 (Toolbox for LiDAR Data Filtering and Forest Studies) to extract the information needed for the model: i.e., tree positions, HGT and CPA (Chapman et al., 2010). For the surrogate plots, tree selection was made using the library of tree models. The best fit tree was selected from the available modeled trees based on HGT, CPA, and DBH, excluding those modeled trees from the reference plot being processed (Fig. 3). The following constraint was minimized:

|DBHin situ − DBHmod | |CPAin situ − CPAmod | |HGTin situ − HGTmod | + + DBHin situ CPAin situ HGTin situ (9) where DBHin situ, CPAin situ and HGTin situ were measured at the plots, and DBHmod, HGTmod, and CPAmod were for modeled trees produced with L-Architect. In this instance, where the constraint variables were measured in situ, the plots are referred to as surrogate ground plots. Errors statistics (Section 3.3) were used to compare the structural attributes of the surrogate ground plots with those of the reference plots. 104

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Fig. 3. Creation of surrogate ground and ALS plots with L-Architect. Attribute definitions are found in the list of acronyms; f(HGT,CPA) is a regression model.

a good compromise for individual tree crown detection between operational “coarser” density usually ≤4 points/m2 (White et al., 2013) and unmanned aerial vehicles with typical density ≥50 points/m2 (Wallace et al., 2014). The canopy surface is the outer portion of forests viewed by remotesensing instruments where numerous functional characteristics of forests are linked to developmental stage. Canopy surface models were produced to assess the simulated ALS data and other factors to account for the dissimilarity with the reference plots. Canopy height models (CHM) were created interpolating the maximum height from the first returns using a 0.5 m grid overlay. A natural neighbor interpolation algorithm was applied to the extracted heights to produce CHMs. Three attributes were calculated: canopy surface area (CSA), canopy surface volume (CSV), and the standard deviation of the CHM (σCHM). In addition, the lacunarity index (ΛCHM) was also calculated using the gliding-box algorithm described in Frazer et al. (2005). Following the simulation of ALS data, several tests were made to find the most adequate parameters in TIFFS to maximize detection of individual tree crowns. Tree detections were processed for each simulated ALS scene on a grid of 15 cm × 15 cm. This resolution was a compromise between the ability to depict the required details without an excessive amount of data. The use of a smoothing filter adjusted according to the type of forest cover (small vs. larger crowns) improved tree detection. Species was assigned to the dominant species measured in the corresponding plot. Surrogate plots were produced using the following constraint:

The performance was assessed with 1590 reference trees and corresponding modeled trees (comprised of 823 fir and 767 spruce trees). Alternatively, surrogate plots can be created using inputs that are derived using above canopy remote-sensing data such as ALS. This requires the detection of individual tree crowns and the estimation of several tree crown attributes (Fig. 3). In this study, we chose to simulate ALS data on the reference plots to gain better control over the ALS point density and provide the necessary inputs for modeling surrogate plots, hereafter referred to as surrogate ALS plots. Radiative transfer models have been used to simulate the return signal of ALS (Disney et al., 2010; Gastellu-Etchegorry et al., 2015; Hancock et al., 2012). They require the definition of input parameters to describe the structural and spectral properties of all forest components and the ground. In addition, the laser scanner and photo-detector characteristics as well as the positions and orientations of the sensor (during the flight) need to be known a priori. Due to the absence of configuration details for a specific sensor and the lack of spectro-directional and illumination-related properties of the modeled species, a simplified approach was adopted, similar to the one proposed by Morsdorf et al. (2009) but adapted to discrete-return systems. In this study, ALS data were simulated for the reference plots with PBRT. A camera with an orthographic projection that preserves the relative distance between objects (no foreshortening effect) was used. A distant light source was used to simulate a source of illumination from the same direction at every point in space. Both the camera and the light source were located at the same position over the scene. Applying this configuration, PBRT was executed to cast a predefined number of rays distributed onto the reference plot. Ray tracing is a point-sampling process: the rays used to assess light intensities are infinitely thin. However, each pixel of a rendered image has a finite width. 64 samples per pixel were used in order to apply the Mitchell's “best-candidate” sampler (Mitchell, 1991) implemented in PBRT. It allowed sampling the plot in a more natural way than using uniform sampling. For each pixel, the closest ray hits' (x,y,z)-position from the camera location was recorded. This approach served as a proxy for recording the first return, or echo, of a simulated laser pulse with a footprint proportional to the pixel size. We tested four different point densities against the capacity to detect the individual tree crowns of the reference plots, namely 4, 16, 64, and 256 points/m2. Using 4 points/m2 to approximately match the existing ALS coverage of Newfoundland, documented in Luther et al. (2014), proved to be too coarse to provide tree-level HGT and CPA for this forest environment. Preliminary analysis showed that nearly a third of the trees were not detected. From the results of this analysis, we chose to simulate the ALS data with a density of 16 points/m2, which is

|DBHALS − DBHmod | |CPAALS − CPAmod | |HGTALS − HGTmod | + + DBHALS CPAALS HGTALS (10) where HGTALS and CPAALS were derived from simulated ALS data using TIFFS. DBHALS of each tree was estimated using a relationship with tree height and crown area. A nonlinear regression model was developed for each species using MATLAB® Statistics Toolbox (MathWorks®) and ground measurements from over 1907 fir and 984 spruce trees in Newfoundland (fir: adjusted R2 and RMSE of 0.75 and 3.13 cm; spruce: adjusted R2 of 0.73 and RMSE of 2.59 cm; coefficients found in Table 1). The following model was used to estimate DBHALS (Filipescu et al., 2012):

DBHALS ≡ f (HGTALS , CPAALS) = β1 (HGTALS − 1.3) β2 ·β3 HGTALS − 1.3·CPAALS β4 (11) As in previous sections, error statistics were calculated by comparing the attributes in the surrogate ALS plots with those attributes 105

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

with attributes of the validation plots using the full library of 100 tree scans guided by three attributes estimated from in situ measurements, namely DBH, HGT, and CPA. Weak correlations were observed for Tcoef and ASYM, whereas correlations over 0.59 were observed for the other attributes. The ratio of dispersion σrec σref for VOL, CPA, Btot and Tcoef had the most discrepancies (> 33%) with the reference values, and the other tree-level attributes were within 20%. At plot level (Fig. 6, bottom), the mean correlation and NRMSE values were 0.91 [R2 = 0.83] and 0.16, respectively. Weaker correspondence were observed for Wβleaf with r = 0.69 [R2 = 0.48] (NRMSE = 0.57). Moderate correspondences were observed for Wβwood, VOL, ASYM and Tcoef with 0.74 ≤ r ≤ 0.89 [0.55 ≤ R2 ≤ 0.79] and 0.18 ≤ NRMSE ≤ 0.47. The ratio of dispersion σrec σref for Wαleaf, Wβleaf, Wβwood and Tcoef had the most discrepancies (> 1.5). In each case, the high agreement between the reference plots attributes and the validation plots attributes validated the modeling approach using a few key in situ measurements and the LArchitect model.

Table 1 Coefficient values for the nonlinear model of DBH for balsam fir and black spruce. Species

Coefficient

Initial values

Estimated values

bFa bF bF bF bSa bS bS bS

β1 β2 β3 β4 β1 β2 β3 β4

2.038 0.657 1.024 0.218 2.038 0.657 1.024 0.218

2.974 0.273 1.015 0.251 2.393 0.522 1.005 0.192

a

bF and bS stands for balsam fir and black spruce respectively.

calculated from the reference plots. Here, the error statistics for the attributes calculated from the surrogate ALS plots incorporate the errors associated with estimating the canopy inputs with simulated ALS data and TIFFS. The error assessment of the ALS surrogate plots was based only on plot-level attributes as it was not possible to match individual trees.

4.2. Sensitivity to the number of tree scans The sensitivity of the plot creation method to the number of tree scans indicated an average loss over the set of 24 plot-level structural attributes for R2 (resp. NRMSE) of 0.01 (0.005) using 75% of the available scans, 0.04 (0.008) using 50% of the available scans, and 0.08 (0.02) using only 25% of the available scans. The sensitivity analysis showed values of R2 higher than 0.92 and NRMSE below 0.10 for HGT (Table 2). Excluding results for the two model input variables DBH and LA, the analysis indicated sizeable differences on R2 (> 0.05) when using only 25 trees compared with 100 trees for 10 of the tested attributes, in particular the vertical profiles of material, branching (Btot, BR, Tcoef) and crown-related attributes (ASYM, CBH, HCPA, HLC). However, the performance was better with Sensitivity 25% for CW, VOL and Λ. No significant differences were observed for the other seven attributes with differences in R2 and NRMSE lower than 0.05 and 0.02 respectively.

4. Results 4.1. Validation of the modeling approach An example of a reference plot and corresponding validation plot modeled with simulated TLS scans is shown in Fig. 4. The reconstruction of 1590 reference trees, dispatched to ten independent threads, elapsed a total of 10.5 h on a Dell Precision T7500 workstation equipped with an Intel® Xeon® Processor X5680 (3.33 GHz) and 24 GB of DDR3 RAM. This corresponded to an average of 4 min per tree (most within 2.5 min), ranging between 32 s and 35 min depending on the tree size. Most of the attributes of the simulated TLS trees were strongly correlated (r > 0.75) with the reference tree attributes with the exception of ASYM and Tcoef. The Taylor diagram of Fig. 5 provides statistics comparing the tree-level attributes of the 100 original reference trees with their associated simulated TLS trees. The ratio of dispersion σrec σref was around 1.0 for all attributes except for Tcoef and Btot for which dispersion was higher. The tree-level comparison of 15 structural attributes indicated an overall agreement between the reference plots and the validation plots using simulated TLS trees with NRMSE values below 16%. Similarly, the plot-level comparison of over 24 structural attributes showed average NRMSE values below 17%. The Taylor diagram of Fig. 6 (top) provides statistics comparing tree-level attributes of the reference plots

4.3. Surrogate plots from in situ data The comparison between the reference plots and surrogate ground plots constructed with the modeled trees gave an average correlation (r) of 0.65 for the 15 tree-level structural attributes (Fig. 7). Weak correspondences were obtained with the reference trees (r < 0.37) for the Tcoef, ASYM, and HLC attributes. Moderate correspondences with the reference trees were observed for VOL and CBH with correlations of 0.46 and 0.51 respectively. Except for Tcoef, HLC, and CBH,

Fig. 4. Balsam fir-dominated permanent sample plot rendered with the physically based ray-tracing software (PBRT), with the reference plot [left] and the corresponding validation plot produced with simulated TLS scan [right]. 106

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Fig. 5. Taylor diagram comparing tree-level attributes of the 100 original reference trees and their associated TLS simulated scanned tree. The attribute definitions are found in the list of acronyms. Statistics are dimensionless, and symbol o represents the origin at (r, CRMSE/σref, σmod/σref) = (1, 0, 1).

attributes. Mean R2 of the CHM attributes was 0.84 and NRMSE was below 0.09. Overall, the surrogate ALS plots have a significant correspondence (R2 > 0.42, mean of 0.65) for a majority of structural attributes.

corresponding values of NRMSE were all below 15%. The ratio of the dispersion (σrec σref ) varied with the attributes but both options generally showed values less than one. At the plot level, significant correspondence with the reference plots was obtained for most of the variables with a mean and median value for R2 of 0.68 and 0.77 respectively (Table 3). The surrogate plots showed weak correspondence with the reference plots for the β-parameter of the Weibull fits, ASYM, HLC, and Λ. Moderate coefficients of variation were observed for Wαwood (R2 = 0.50), VOL (R2 = 0.52), and Tcoef (R2 = 0.60). Mean NRMSE value was 0.17.

5. Discussion Estimating structural attributes using remote-sensing data is a challenge considering the complexity and variety of tree and stand architectures. Regardless, a large number of methods exist to link variables available from remote sensing with forest structural attributes. The methods described in this study improve our ability to estimate detailed forest structural attributes by taking advantage of TLS data and an architectural model. Terrestrial laser scanning data contain relevant detailed 3D structural information, and L-Architect reproduces the detailed forest structure from TLS data. The combination of TLS data and L-Architect offers the ability to produce surrogate plots in boreal forests: detailed 3D representations of the stand structure of forest areas beyond what can be measured by field crews dealing with operational constraints. Several architectural models can be used to simulate trees from TLS data sets (Côté et al., 2011; Hackenberg et al., 2014; Raumonen et al., 2013). However, most of these models involve forward modeling only, i.e., use of TLS data to simulate individual trees. The simulation of all trees composing a plot is also possible (Côté et al., 2012) but requires extensive in situ measurements. The method presented here goes a step further using inverse modeling to create surrogate plots with simulated ALS data. Forward modeling produces a library of trees using TLS scans similar to how functional structural plant models simulate individual trees (Cieslak et al., 2011; Fourcaud et al., 2008). Then, inverse modeling with L-Architect segments all individual trees of an area and estimates at least three attributes for each tree: position, HGT, and CPA. In our test case, the estimates were derived from a simulated ALS data

4.4. Surrogate plots from simulated ALS data When comparing known values of the reference plots with the results from TIFFS, we obtained average per plot absolute differences of two detected trees (2%), 7 cm for the mean crown radius (17%), and 1.46 m for the mean tree height (18%). The Taylor diagram of Fig. 8 provides a comparison of 24 selected attributes from the reference plots with those of the surrogate ALS plots. For most attributes, the mean correlation values were 0.81 (r > 0.65). Weak correspondence were observed for Wβleaf, HLC, ASYM, CBH, Tcoef and Λ with 0.003 ≤ R2 ≤ 0.17 and 0.22 ≤ NRMSE ≤ 0.57 (Table 3). Moderate correspondences were observed for Wβwood, Wαwood, HCPA, DBH, Wαleaf, Ksurf and VOL with 0.42 ≤ R2 ≤ 0.60 and 0.16 ≤ NRMSE ≤ 0.27. The performance statistics for structural attributes also display much contrasted results between attributes: some attributes are very well predicted, whereas others have only weak correspondence. For instance, CPA, Btot, Wtot, CW, PAI, FC, HGT, LAI, Dδ, and LA showed the most favorable correspondence and error with 0.62 ≤ R2 ≤ 0.88 and 0.12 ≤ NRMSE ≤ 0.29 (Table 3). However, Wβleaf, HLC, ASYM and CBH showed much weaker correspondence (R2 < 0.05). Lastly, displaying the correlation for four attributes in Fig. 9 identifies how this correspondence works for these key CHM 107

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Fig. 6. Taylor diagrams comparing tree- [top] and plot- [bottom] level attributes of the reference plots with their validation plots created using TLS simulations of 100 tree scans. The attribute definitions are found in the list of acronyms. Statistics are dimensionless, and symbol o represents the origin at (r, CRMSE/σref, σmod/ σref) = (1, 0, 1).

108

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Table 2 Results of the sensitivity analysis to test the difference between the reference plots and the validation plots modeled with simulated TLS tree scans. Attributes

ASYM (m) CBH (m) CPA (m2) HCPA (m) HGT (m) CW (m) HLC (m) VOL (m3) Btot BR DBH (cm) Wtot Ksurf (m2) Tcoef Dδ FC Λ LAI LA (m2) PAI Wαleaf Wαwood Wβleaf Wβwood

Sensitivity 100% (100 tree scans)

Sensitivity 75% (75 tree scans)

Sensitivity 50% (50 tree scans)

Sensitivity 25% (25 tree scans)

R2

NRMSE

R2

NRMSE

R2

NRMSE

R2

NRMSE

0.667 0.722 0.714 0.88 0.97 0.803 0.73 0.545 0.835 0.857 1 0.871 0.993 0.688 0.912 0.854 0.757 0.999 1 0.996 0.889 0.877 0.48 0.791

0.181 0.256 0.197 0.11 0.072 0.16 0.245 0.205 0.251 0.142 0 0.229 0.04 0.466 0.116 0.142 0.105 0.011 0 0.03 0.116 0.195 0.383 0.266

0.576 0.716 0.709 0.916 0.96 0.767 0.723 0.536 0.828 0.878 1 0.858 0.992 0.659 0.89 0.819 0.775 0.999 1 0.996 0.907 0.891 0.439 0.745

0.207 0.251 0.209 0.108 0.083 0.172 0.236 0.212 0.261 0.149 0 0.24 0.046 0.483 0.12 0.163 0.1 0.009 0 0.031 0.114 0.196 0.374 0.283

0.674 0.639 0.731 0.864 0.922 0.698 0.511 0.799 0.833 0.838 1 0.864 0.992 0.639 0.88 0.802 0.75 0.999 1 0.996 0.784 0.805 0.262 0.693

0.183 0.252 0.212 0.103 0.101 0.204 0.22 0.148 0.252 0.168 0 0.229 0.044 0.487 0.137 0.184 0.111 0.009 0 0.033 0.142 0.199 0.387 0.294

0.203 0.445 0.666 0.717 0.959 0.89 0.552 0.643 0.779 0.766 1 0.847 0.991 0.567 0.919 0.862 0.856 0.999 1 0.995 0.722 0.778 0.186 0.671

0.497 0.27 0.203 0.142 0.076 0.115 0.301 0.172 0.249 0.162 0 0.221 0.04 0.469 0.103 0.133 0.083 0.011 0 0.031 0.153 0.17 0.411 0.367

R2 = 0.32) correlated with the references. Lacunarity (Λ) expresses the heterogeneity in the distribution of material volume and this statistic is more sensitive to minor differences in the structure than the fractal dimension (Dδ) explained by the gap distribution. Overall, the results demonstrate the capacity to construct detailed forest plot structure using TLS and L-Architect. Given the realistic representation of plot structure using in situ measurements or simulated ALS data, we believe that the surrogate plots provide an alternative to detailed in situ measurements for forest inventory. This study assessed the sensitivity of the modeling approach to the number of modeled trees required from the library. Generally, the greatest improvements were observed when the number of trees increased from 25 to 50. Beyond 50 trees, the improvements were modest. Moreover in some cases, the number of trees did not improve the results significantly (e.g., CW, VOL, and Λ; see Table 2). Fortunately, increasing the number of trees does not have a major impact on the computing requirements of the method. Therefore, we concluded that using 100 trees in the library was a good compromise to allow flexibility and efficient computation considering the relatively simple environment dominated by only two species (fir and spruce). Another level of analysis, not investigated in this study, consist of repeating the creation of surrogate plots by varying the configuration of LiDAR simulations on the reference trees (TLS) and plots (ALS). The independent sets of surrogate plots would allow estimating the accuracy and precision of the reconstruction method through a series of different LiDAR acquisition settings: e.g. viewing geometry, random noise level, point density, beam divergence. For this study, we used L-Architect because of its ability to reproduce fir and spruce in mature boreal forests but the proposed method could be adapted for other tree modeling approaches. For instance, in more complex environments, one option for generating realistic tree models could be using a Bayesian approach where stochastic tree growth models are optimised against TLS measurements to generate a large number of similar “clones” (Potapov et al., 2016, 2017). In the future, public libraries of tree models could be developed for different species from different environments. Such libraries could then be used to create virtual plots provided that a standard file format was used for storing the 3D structures. Several potential improvements could improve the results or

set. The ability to produce detailed 3D surrogate plots with only a few input variables available from above-canopy remote-sensing data sets represents a considerable step forward in our ability to estimate a wide range of structural attributes. Surrogate plots are modeled representations, however, they provide details that exceed what can be measured in situ within a reasonable amount of time. Some attributes are very complex to estimate accurately even for a field crew doing in situ measurements (e.g., 3D distribution of wood, foliage or gaps, tree topology). Therefore, an error assessment of the surrogate plots was not practical in a real environment. This motivated the use of a virtual, or controlled, environment (Côté et al., 2009; Leblanc and Fournier, 2014; Bremer et al., 2018) to evaluate the modeling approach and associated method of creating surrogate plots. Consequently, the errors of the structural variables were estimated based on known reference values of simulated scenes. The differences observed between the reference and the surrogate plots vary among the assessed variables. At the tree-level, the best modeled structural attributes of the surrogate ground plots, in decreasing order, included DBH, Ksurf, LA, and HGT, with R2 ranging from 0.75 to 0.92. The tree-level variables most challenging to model were VOL (R2 = 0.21), Tcoef (R2 = 0.13), HLC (R2 = 0.11), and ASYM (R2 = 0.03). These attributes correlated weakly with the reference essentially because we did not implement a procedure to adjust the crown shape and structure of the selected model according to individual tree measurements and its neighbours with TIFFS. The selection of a tree model without adjusting or reshaping its form and structure explains in part the low correlation values observed for these crown attributes and taper coefficient. In addition, not considering the tree neighbourhood configuration in the selection process (regardless of the availability of a tree reshaping procedure), was a significant factor influencing the resulting tree morphology, which contributed to increase the observed discrepancies. At the plot level, the vertical distribution of canopy components was correctly modeled (fitted) through the α-parameter of the Weibull function (R2 = 0.50–0.58). In contrast the β-parameter was notably sensitive to small variations in the spatial pattern (R2 = 0.23–0.24). However, the horizontal distribution of gaps (FC) showed strong agreement with the reference plots (R2 = 0.81). Lastly, the general spatial statistics were highly (Dδ, R2 = 0.81) and poorly (Λ, 109

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Fig. 7. Taylor diagram comparing tree- [top] and plot- [bottom] level attributes of the reference plots with the surrogate ground plots modeled using in situ measurements. The attribute definitions are found in the list of acronyms. Statistics are dimensionless, and symbol o represents the origin at (r, CRMSE/σref, σmod/ σref) = (1, 0, 1).

110

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

observed a loss of accuracy of around 7% for NRMSE (10% for R2) compared with the surrogate plots created with in situ measurements. This loss in accuracy is due in part to errors in the estimated HGT and CPA. Airborne laser scanning has been successfully used to estimate a wide range of tree descriptive attributes beyond position, HGT, and CPA (Korhonen and Morsdorf, 2014; Lin and West, 2016; Vauhkonen et al., 2014). For instance the light resource competition between trees is known to affect the preferred orientation of branches and crown plasticity (Pirk et al., 2012). Adding new variables could include competition indices and crown shape metrics. These variables could improve matching trees from the remote-sensing data with those from the library. A second way to improve the method would be to alter the selected trees in accordance with known differences from the input variables. Currently, selected trees are scaled following a simple procedure. However, it would be reasonable to reshape trees dynamically within their environmental conditions, e.g., according to light conditions and proximity to other trees (Pirk et al., 2012). Although this reshaping might deform the basic tree structure, we could also include realistic growth and procedural rules (e.g., Greenlab; Yan et al., 2004). A third way to improve the method is to add more species and ecosystems to the library of trees. L-Architect is currently limited to conifers or trees with a monopodial structure. Adding more conifer species would not require fundamental changes to the model. However, deciduous species exhibit sympodial structure, which makes them more complex to derive from TLS data. The use of fast skeleton extraction (Hackenberg et al., 2014; Raumonen et al., 2013; Ravaglia et al., 2017) would allow modeling the main branching structure and stem taper for a large number of either coniferous or deciduous tree species. Lastly, this study used ALS data for deriving input variables to the inverse modeling approach. The use of very high spatial resolution satellite images or their combination with ALS data is another possibility. The ability to create surrogate plots brings several new perspectives to the remote sensing of forest areas. The plot creation method using LArchitect and inverse modeling enables estimation of tree and plot attributes beyond those that can be measured easily on the ground. There

Table 3 Performance statistics comparing structural attributes of the reference plots with those of the surrogate plots. Attributes

ASYM (m) CBH (m) CPA (m2) HCPA (m) HGT (m) CW (m) HLC (m) VOL (m3) Btot BR DBH (cm) Wtot Ksurf (m2) Tcoef Dδ FC Λ LAI LA (m2) PAI Wαleaf Wαwood Wβleaf Wβwood

Surrogate ground plots

Surrogate ALS plots

R2

NRMSE

R2

NRMSE

0.163 0.639 0.680 0.764 0.862 0.767 0.436 0.515 0.862 0.828 0.986 0.853 0.964 0.596 0.812 0.812 0.317 0.970 0.969 0.974 0.576 0.503 0.242 0.226

0.296 0.197 0.206 0.118 0.129 0.168 0.234 0.208 0.247 0.136 0.040 0.228 0.080 0.416 0.121 0.153 0.187 0.049 0.051 0.053 0.188 0.179 0.290 0.222

0.018 0.034 0.622 0.487 0.815 0.726 0.014 0.547 0.692 0.603 0.527 0.706 0.541 0.148 0.832 0.797 0.170 0.818 0.875 0.754 0.532 0.477 0.003 0.423

0.309 0.412 0.285 0.201 0.131 0.245 0.422 0.249 0.288 0.196 0.176 0.256 0.172 0.567 0.117 0.155 0.215 0.131 0.128 0.151 0.267 0.221 0.342 0.161

enlarge the range of applications of the plot creation method. Adding more input variables has potential to improve the results without adding much procedural complexity. For instance, in our study we used CPA and HGT from the ALS data as the input variables to find the trees from the library. We showed that, using available software (i.e., TIFFS) and high resolution simulated ALS data, we obtained a detailed and reliable structural representation of forest canopy. However, we

Fig. 8. Taylor diagram comparing plot-level attributes of the reference plots with the surrogate ALS plots modeled using input derived using simulated ALS and TIFFS. The attribute definitions are found in the list of acronyms. Statistics are dimensionless, and symbol o represents the origin at (r, CRMSE/σref, σmod/σref) = (1, 0, 1).

111

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Fig. 9. Comparison between the reference plots and the surrogate ALS plots generated with simulated ALS data for: lacunarity, canopy surface area, canopy surface volume, and the standard deviation of the canopy height model (Stdev CHM).

Our results demonstrate the potential of L-Architect to reproduce the detailed forest structure at plot level with good correspondence (mean R2 = 0.83). The comparison between the surrogate plots and the reference plots showed that the use of a library of TLS coregistered tree scans was well adapted for natural forest. The surrogate plots created using the trees modeled by L-Architect reproduced the attributes of the reference plots with NRMSE (R2) values between 4 and 42% (0.16–0.97) with median of 18% (0.77) using in situ ground measurements. The same approach using inputs estimated from ALS gave NRMSE (R2) values between 12 and 57% (0.00–0.88) with median of 22% (0.54). The creation of surrogate plots provides a unique solution for linking remote-sensing observations to detailed physiological characteristics and ecological indicators of forest ecosystems.

are also possibilities to expand the attributes modeled by L-Architect. For instance, a 3D estimation of leaf area density, such as proposed by Pimont et al. (2018), could support ecophysiology studies. The availability of a large number of attributes brings possibilities for modeling new relationships including complex attributes like wood fiber quality or stem branchiness. For instance Luther et al. (2014) and Blanchette et al. (2015) explored how remote-sensing metrics can predict wood fiber quality. The use of fine-scale structural attributes of surrogate plots as predictor variables has the potential to improve on these relationships. The method and procedures described herein for creating surrogate plots and associated fine-scale structural attributes improves upon current LiDAR remote sensing capabilities by combining detailed information of TLS data and structural attributes available from ALS to produce detailed 3D representations of forest areas.

Acknowledgments 6. Conclusion The Newfoundland Fibre Project is a research initiative of the Canadian Wood Fibre Centre (CWFC), Canadian Forest Service of Natural Resources Canada, Corner Brook Pulp and Paper Limited (CBPPL), FPInnovations, Grenfell Campus of Memorial University of Newfoundland (MUN), NL Department of Natural Resources, Nova Scotia Community College and the University of Sherbrooke. We thank all participants of the project, notably, Dr. Wade W. Bowers of Grenfell Campus (MUN). We thank Danny Blanchette and Emilie Lessard of the University of Sherbrooke for their help during the data acquisition and processing of the tree map locations and TLS scans. We are also grateful to Janet Bourque, Matthew Glover, Heidi Kavanagh, and Michael Noonan of Grenfell Campus (MUN) for their help acquiring and

In this paper, we present a method for creating surrogate plots using in situ or ALS data making possible the selection of fine-scale structural attributes for developing empirical relationships with above-canopy remote-sensing. The method uses TLS data and the L-Architect model to create a library of trees (forward modeling) and then uses tree-level input variables, from ground survey or remote-sensing data, to simulate surrogate plots (inverse modeling) with detailed representation of the structural attributes. The method was validated in a virtual environment as the required measurements for the comparison with a reference data set was impossible to obtain at the correct level of detail and accuracy. 112

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

processing TLS scans. Thanks to Caroline Simpson and the anonymous reviewers for kindly commenting on, and incidentally improving the quality of, this manuscript. This work was supported by the Atlantic Canada Opportunities Agency (ICF 786-16766-195492 to W.W.B.), Centre for Forest Science and Innovation (221208 to W.W.B.), Natural Science and Engineering Council of Canada (CRDPJ-390394 to R.A.F. and J.E.L., RGPIN-2014-04508 to R.A.F.), NL Department of Innovation, Business and Rural Development (WES2010.06.02.045W to W.W.B.) and Research and Development Corporation (5404.1221.102 to W.W.B.). This research was also supported by the AWARE project (NSERC CRDPJ-462973-14, grantee N.C. Coops, UBC), in collaboration with the CWFC, FPInnovations and CBPPL.

Filipescu, C.N., Groot, A., MacIsaac, D.A., Cruickshank, M.G., Stewart, J.D., 2012. Prediction of diameter using height and crown attributes: a case study. West. J. Appl. For. 27, 30–35. Fourcaud, T., Zhang, X., Stokes, A., Lambers, H., Körner, C., 2008. Plant growth modelling and applications: the increasing importance of plant architecture in growth models. Ann. Bot. 101, 1053–1063. Fraser, D.A., 1966. Vegetative and reproductive growth of black spruce (Picea mariana (Mill.) B.S.P.) at Chalk River, Ontario, Canada. Can. J. Bot. 44, 567–580. Frazer, G.W., Wulder, M.A., Niemann, K.O., 2005. Simulation and quantification of the fine-scale spatial pattern and heterogeneity of forest canopy structure: a lacunaritybased method designed for analysis of continuous canopy heights. For. Ecol. Manag. 214, 65–90. Gastellu-Etchegorry, J.P., Demarez, V., Pinel, V., Zagolski, F., 1996. Modeling radiative transfer in heterogeneous 3-D vegetation canopies. Remote Sens. Environ. 58, 131–156. Gastellu-Etchegorry, J.-P., Yin, T., Lauret, N., Cajgfinger, T., Gregoire, T., Grau, E., Feret, J.-B., Lopes, M., Guilleux, J., Dedieu, G., Malenovský, Z., Cook, B., Morton, D., Rubio, J., Durrieu, S., Cazanave, G., Martin, E., Ristorcelli, T., 2015. Discrete anisotropic radiative transfer (DART 5) for modeling airborne and satellite spectroradiometer and LiDAR acquisitions of natural and urban landscapes. Remote Sens. 7, 1667–1701. Godin, C., 2000. Representing and encoding plant architecture: a review. Ann. For. Sci. 57, 413–438. Godin, C., Caraglio, Y., 1998. A multiscale model of plant topological structures. J. Theor. Biol. 191, 1–46. Gonzalez de Tanago, J., Lau, A., Bartholomeus, H., Herold, M., Avitabile, V., Raumonen, P., Martius, C., Goodman, R.C., Disney, M., Manuri, S., Burt, A., Calders, K., 2017. Estimation of above-ground biomass of large tropical trees with terrestrial LiDAR. Methods Ecol. Evol. 9, 223–234. Govaerts, Y.M., Verstraete, M.M., 1998. Raytran: a Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media. IEEE Trans. Geosci. Remote Sens. 36, 493–505. Hackenberg, J., Morhart, C., Sheppard, J., Spiecker, H., Disney, M., 2014. Highly accurate tree models derived from terrestrial laser scan data: a method description. Forests 5, 1069–1105. Hancock, S., Lewis, P., Foster, M., Disney, M., Muller, J.-P., 2012. Measuring forests with dual wavelength LiDAR: a simulation study over topography. Agric. For. Meteorol. 161, 123–133. Hopkinson, C., Chasmer, L., Young-Pow, C., Treitz, P., 2004. Assessing forest metrics with a ground-based scanning LiDAR. Can. J. For. Res. 34, 573–583. Kangas, A.S., 2010. Value of forest information. Eur. J. For. Res. 129, 863–874. Kankare, V., Vastaranta, M., Holopainen, M., Räty, M., Yu, X., Hyyppä, J., Hyyppä, H., Alho, P., Viitala, R., 2013. Retrieval of forest aboveground biomass and stem volume with airborne scanning LiDAR. Remote Sens. 5, 2257–2274. Kirkpatrick, L.A., Weishampel, J.F., 2005. Quantifying spatial structure of volumetric neutral models. Ecol. Model. 186, 312–325. Koch, B., 2010. Status and future of laser scanning, synthetic aperture radar and hyperspectral remote sensing data for forest biomass assessment. ISPRS J. Photogram. Remote Sens. 65, 581–590. Korhonen, L., Morsdorf, F., 2014. Estimation of canopy cover, gap fraction and leaf area index with airborne laser scanning. In: Maltamo, M., Næsset, E., Vauhkonen, J. (Eds.), Forestry Applications of Airborne Laser Scanning: Concepts and Case Studies. Springer, Dordrecht, The Netherlands, pp. 397–417. Landry, R., Fournier, R.A., Ahern, F.J., Lang, R.H., 1997. Tree vectorization: a methodology to characterize fine tree architecture in support of remote sensing models. Can. J. Remote. Sens. 23, 91–107. Lavigne, M.B., Luther, J.E., Franklin, S.E., Hunt Jr., E.R., 1996. Comparing branch biomass prediction equations for Abiesbaisamea. Can. J. For. Res. 26, 611–616. Leblanc, S.G., Fournier, R.A., 2014. Hemispherical photography simulations with an architectural model to assess retrieval of leaf area index. Agric. For. Meteorol. 194, 64–76. Lenz, P., Bernier-Cardou, M., MacKay, J., Beaulieu, J., 2012. Can wood properties be predicted from the morphological traits of a tree? A canonical correlation study of plantation-grown white spruce. Can. J. For. Res. 42, 1518–1529. Lessard, E., Fournier, R.A., Luther, J.E., Mazerolle, M.J., van Lier, O.R., 2014. Modeling wood fiber attributes using forest inventory and environmental data for Newfoundland's boreal forest. For. Ecol. Manag. 313, 307–318. Lin, Y., West, G., 2016. Retrieval of effective leaf area index (LAIe) and leaf area density (LAD) profile at individual tree level using high density multi-return airborne LiDAR. Int. J. Appl. Earth Obs. Geoinf. 50, 150–158. Livny, Y., Pirk, S., Cheng, Z., Yan, F., Deussen, O., Cohen-or, D., Chen, B., 2011. Texturelobes for tree modelling. ACM Trans. Graph. 30 (53), 1–10. Luther, J.E., Skinner, R., Fournier, R.A., van Lier, O.R., Bowers, W.W., Coté, J.F., Hopkinson, C., Moulton, T., 2014. Predicting wood quantity and quality attributes of balsam fir and black spruce using airborne laser scanner data. Forestry 87, 313–326. Mahoney, C., Hall, R., Hopkinson, C., Filiatrault, M., Beaudoin, A., Chen, Q., 2018. A forest attribute mapping framework: a pilot study in a Northern Boreal forest, Northwest Territories, Canada. Remote Sens. 10, 1338. Meades, W.J., Moores, L., 1994. Forest Site Classification Manual: A Field Guide to the Damman Forest Types of Newfoundland. Canada–Newfoundland Forest Resource Development; Forestry Canada and Newfoundland, Labrador. Department of Forestry & Agriculture, Western Newfoundland Model Forest, Inc., Corner Brook, NL, Canada. Měch, R., Prusinkiewicz, P., 1996. Visual models of plants interacting with their environment. In: Proceedings of SIGGRAPH 96 (New Orleans, Louisiana, August 4–9, 1996). Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, pp. 397–410. Mitchell, D.P., 1991. Spectrally optimal sampling for distribution ray tracing. SIGGRAPH

References Avery, T.E., Burkhart, H.E., 2002. Forest Measurements. McGraw-Hill, Boston, MA. Béland, M., Widlowski, J.L., Fournier, R.A., 2014. A model for deriving voxel-level tree leaf area density estimates from ground-based LiDAR. Environ. Model. Softw. 51, 184–189. Blanchette, D., Fournier, R.A., Luther, J.E., Côté, J.-F., 2015. Predicting wood fiber attributes using local-scale metrics from terrestrial LiDAR data: a case study of Newfoundland conifer species. For. Ecol. Manag. 347, 116–129. Boudon, F., Nouguier, C., Godin, C., Pradal, C., 2003. GEOM Module Manual I. User guide, CIRAD, Montpellier. Bouvier, M., Durrieu, S., Fournier, R.a., Renaud, J.-P., 2015. Generalizing predictive models of forest inventory attributes using an area-based approach with airborne LiDAR data. Remote Sens. Environ. 156, 322–334. Bremer, M., Wichmann, V., Rutzinger, M., 2018. Multi-temporal fine-scale modelling of Larix decidua forest plots using terrestrial LiDAR and hemispherical photographs. Remote Sens. Environ. 206, 189–204. Calders, K., Newnham, G., Burt, A., Murphy, S., Raumonen, P., Herold, M., Culvenor, D., Avitabile, V., Disney, M., Armston, J., Kaasalainen, M., 2015. Nondestructive estimates of above-ground biomass using terrestrial laser scanning. Methods Ecol. Evol. 6, 198–208. Calders, K., Origo, N., Burt, A., Disney, M., Nightingale, J., Raumonen, P., Åkerblom, M., Malhi, Y., Lewis, P., 2018. Realistic forest stand reconstruction from terrestrial LiDAR for radiative transfer modelling. Remote Sens. 2018 (10), 933. Ceccato, P., Gobron, N., Flasse, S., Pinty, B., Tarantola, S., 2002. Designing a spectral index to estimate vegetation water content from remote sensing data: part 1 theoretical approach. Remote Sens. Environ. 82, 188–197. Cescatti, A., 1997. Modelling the radiative transfer in discontinuous canopies of asymmetric crowns. II. Model testing and application in a Norway spruce stand. Ecol. Model. 101, 275–284. Champion, I., Porte, A., Bert, D., Loustau, D., Guedon, M., Jean-Courdier, F., Lagane, F., Lambrot, C., Lardit, A., Sartore, M., 2001. Tree architecture in remote sensing analytical models: the Bray experiment. Int. J. Remote Sens. 22, 1827–1843. Chapman, J., Hung, I.-K., Tippen, J., 2010. Evaluating TIFFS (toolbox for LiDAR data filtering and forest studies) in deriving forest measurements from LiDAR data. Math. Comp. For. Nat. Res. Sci. 2, 145–152. Chen, J.M., Leblanc, S.G., 1997. A four-scale bidirectional reflectance model based on canopy architecture. IEEE Trans. Geosci. Remote Sens. 35, 1316–1337. Cieslak, M., Seleznyova, A.N., Prusinkiewicz, P., Hanan, J., 2011. Towards aspect-oriented functional–structural plant modelling. Ann. Bot. 108, 1025–1041. Colombo, S.J., 1986. Second-year shoot development in black spruce Picea mariana (Mill.) B.S.P. container seedlings. Can. J. For. Res. 16, 68–73. Côté, J.-F., Widlowski, J.-L., Fournier, R.A., Verstraete, M.M., 2009. The structural and radiative consistency of three-dimensional tree reconstructions from terrestrial LiDAR. Remote Sens. Environ. 113, 1067–1081. Côté, J.-F., Fournier, R.A., Egli, R., 2011. An architectural model of trees to estimate forest structural attributes using terrestrial LiDAR. Environ. Model. Softw. 26, 761–777. Côté, J.-F., Fournier, R.A., Frazer, G.W., Olaf Niemann, K., 2012. A fine-scale architectural model of trees to enhance LiDAR-derived measurements of forest canopy structure. Agric. For. Meteor. 166–167, 72–85. Côté, J.-F., Fournier, R.A., Luther, J.E., 2013. Validation of L-Architect model for balsam fir and black spruce trees with structural measurements. Can. J. Remote. Sens. 39, S41–S59. Dassot, M., Constant, T., Fournier, M., 2011. The use of terrestrial LiDAR technology in forest science: application fields, benefits and challenges. Ann. For. Sci. 68, 959–974. Disney, M.I., Lewis, P., North, P.R.J., 2000. Monte Carlo ray tracing in optical canopy reflectance modelling. Remote Sens. Rev. 18, 163–196. Disney, M.I.I., Kalogirou, V., Lewis, P., Prieto-Blanco, A., Hancock, S., Pfeifer, M., 2010. Simulating the impact of discrete-return LiDAR system and survey characteristics over young conifer and broadleaf forests. Remote Sens. Environ. 114, 1546–1560. Disney, M.I., Boni Vicari, M., Burt, A., Calders, K., Lewis, S.L., Raumonen, P., Wilkes, P., 2018. Weighing trees with lasers: advances, challenges and opportunities. Interf. Focus 8, 20170048. Falkowski, M.J., Hudak, A.T., Crookston, N.L., Gessler, P.E., Uebler, E.H., Smith, A.M.S., 2010. Landscape-scale parameterization of a tree-level forest growth model: a knearest neighbor imputation approach incorporating LiDAR data. Can. J. For. Res. 40, 184–199.

113

Remote Sensing of Environment 219 (2018) 99–114

J.-F. Côté et al.

Stovall, A.E.L., Vorster, A.G., Anderson, R.S., Evangelista, P.H., Shugart, H.H., 2017. Nondestructive aboveground biomass estimation of coniferous trees using terrestrial LiDAR. Remote Sens. Environ. 200, 31–42. Taylor, K.E., 2001. Summarizing multiple aspects of model performance in a single diagram. J. Geophys. Res. Atmos. 106, 7183–7192. van der Zande, D., Hoet, W., Jonckheere, I., van Aardt, J., Coppin, P., 2006. Influence of measurement set-up of ground-based LiDAR for derivation of tree structure. Agric. For. Meteorol. 141, 147–160. van der Zande, D., Jonckheere, I., Stuckens, J., Verstraeten, W.W., Coppin, P., 2008. Sampling design of ground-based LiDAR measurements of forest canopy structure and its effect on shadowing. Can. J. Remote. Sens. 34, 526–538. van Leeuwen, M., Nieuwenhuis, M., 2010. Retrieval of forest structural parameters using LiDAR remote sensing. Eur. J. For. Res. 129, 749–770. Vauhkonen, J., Ørka, H.O., Holmgren, J., Dalponte, M., Heinzel, J., Koch, B., 2014. Tree species recognition based on airborne laser scanning and complementary data sources. In: Maltamo, M., Næsset, E., Vauhkonen, J. (Eds.), Forestry Applications of Airborne Laser Scanning: Concepts and Case Studies. Springer, Dordrecht, The Netherlands, pp. 135–156. Verstraete, M.M., Pinty, B., Myneni, R.B., 1996. Potential and limitations of information extraction on the terrestrial biosphere from satellite remote sensing. Remote Sens. Environ. 58, 201–214. Wallace, L., Lucieer, A., Watson, C.S., 2014. Evaluating tree detection and segmentation routines on very high resolution UAV LiDAR data. IEEE Trans. Geosci. Remote Sens. 52, 7619–7628. Watt, P.J., Donoghue, D.N.M., 2005. Measuring forest structure with terrestrial laser scanning. Int. J. Remote Sens. 26, 1437–1446. Weng, C., Jackson, S.T., 2000. Species differentiation of North American spruce (Picea) based on morphological and anatomical characteristics of needles. Can. J. Bot. 78, 1367–1383. White, J.C., Wulder, M.A., Varhola, A., Vastaranta, M., Coops, N.C., Cook, B.D., Pitt, D., Woods, M., 2013. A Best Practices Guide for Generating Forest Inventory Attributes From Airborne Laser Scanning Data Using an Area-based Approach. Natural Resources Canada, Canadian Forest Service - Canadian Wood Fibre Centre, Victoria, BC (Information Report FI-X-010). Widlowski, J.-L., Pinty, B., Lavergne, T., Verstraete, M.M., Gobron, N., 2006. Horizontal radiation transport in 3-D forest canopies at multiple spatial resolutions: simulated impact on canopy absorption. Remote Sens. Environ. 103, 379–397. Widlowski, J.-L., Côté, J.-F., Béland, M., 2014. Abstract tree crowns in 3D radiative transfer models: impact on simulated open-canopy reflectances. Remote Sens. Environ. 142, 155–175. Widlowski, J.-L., Mio, C., Disney, M., Adams, J., Andredakis, I., Atzberger, C., Brennan, J., Busetto, L., Chelle, M., Ceccherini, G., Colombo, R., Côté, J.-F., Eenmäe, A., Essery, R., Gastellu-Etchegorry, J.-P., Gobron, N., Grau, E., Haverd, V., Homolová, L., Huang, H., Hunt, L., Kobayashi, H., Koetz, B., Kuusk, A., Kuusk, J., Lang, M., Lewis, P.E., Lovell, J.L., Malenovský, Z., Meroni, M., Morsdorf, F., Mõttus, M., Ni-Meister, W., Pinty, B., Rautiainen, M., Schlerf, M., Somers, B., Stuckens, J., Verstraete, M.M., Yang, W., Zhao, F., Zenone, T., 2015. The fourth phase of the radiative transfer model intercomparison (RAMI) exercise: actual canopy scenarios and conformity testing. Remote Sens. Environ. 169, 418–437. Wulder, M.A., Hall, R.J., Coops, N.C., Franklin, S.E., 2004. High spatial resolution remotely sensed data for ecosystem characterization. Bioscience 54, 511–521. 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. Xu, H., Gossett, N., Chen, B., 2007. Knowledge and heuristic-based modeling of laserscanned trees. ACM Trans. Graph. 26 (19), 1–13. Yan, H.P., Kang, M.Z., De Reffye, P., Dingkuhn, M., 2004. A dynamic, architectural plant model simulating resource-dependent growth. Ann. Bot. 93, 591–602. Zeide, B., Pfeifer, P., 1991. A method for estimation of fractal dimension of tree crowns. For. Sci. 37, 1253–1265.

Comput. Graph. 25, 157–164. Morris, R.F., 1951. The effects of flowering on the foliage production and growth of balsam fir. For. Chron. 27, 40–57. Morsdorf, F., Nichol, C., Malthus, T., Woodhouse, I.H., 2009. Assessing forest structural and physiological information content of multi-spectral LiDAR waveforms by radiative transfer modelling. Remote Sens. Environ. 113, 2152–2163. Myneni, R.B., Hoffman, S., Knyazikhin, Y., Privette, J.L., Glassy, J., Tian, Y., Votava, P., Nemani, R.R., Running, S.W., 2002. Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sens. Environ. 83, 214–231. Newfoundland Forest Service, 2011. Permanent Sample Plot Program Procedures and Specifications (Revision). Newfoundland Forest Service, Forest Ecosystem Management Division, Corner Brook, NL. Newton, P.F., 2006. Intraspecific Competition Processes and Their Management Within Black Spruce (Picea mariana (Mill.) B.S.P.) Stands. University of British Columbia, Vancouver, BC, Canada. Palubicki, W., Horel, K., Longay, S., Runions, A., Lane, B., Měch, R.I., Prusinkiewicz, P., 2009. Self-organizing tree models for image synthesis. ACM Trans. Graph. 28 (58), 1–10. Pharr, M., Humphreys, G., 2010. Physically Based Rendering: From Theory to Implementation, 2nd ed. Morgan Kaufmann Publishers Inc., Burlington, MA. Piene, H., 1983. Nondestructive estimation of foliar biomass in balsam fir. Can. J. For. Res. 13, 672–677. Pimont, F., Allard, D., Soma, M., Dupuy, J.-L., 2018. Estimators and confidence intervals for plant area density at voxel scale with T-LiDAR. Remote Sens. Environ. 215, 343–370. Pinty, B., Lavergne, T., Dickinson, R.E., Widlowski, J.L., Gobron, N., Verstraete, M.M., 2006. Simplifying the interaction of land surfaces with radiation for relating remote sensing products to climate models. J. Geophys. Res. D Atmos. 111, 1–20. Pirk, S., Stava, O., Kratt, J., Said, M.A.M., Neubert, B., Měch, R., Benes, B., Deussen, O., 2012. Plastic trees: interactive self-adapting botanical tree models. ACM Trans. Graph. 31 (50), 1–10. Plotnick, R., Gardner, R., Hargrove, W., Prestegaard, K., Perlmutter, M., 1996. Lacunarity analysis: a general technique for the analysis of spatial patterns. Phys. Rev. E Stat. Phys. 53, 5461–5468. Potapov, I., Järvenpää, M., Åkerblom, M., Raumonen, P., Kaasalainen, M., 2016. Databased stochastic modeling of tree growth and structure formation. Silva Fenn. 50 (1), 1413. Potapov, I., Järvenpää, M., Åkerblom, M., Raumonen, P., Kaasalainen, M., 2017. Bayes Forest: a data-intensive generator of morphological tree clones. GigaScience 6, 1–13. Pradal, C., Boudon, F., Nouguier, C., Chopard, J., Godin, C., 2009. PlantGL: a Pythonbased geometric library for 3D plant modelling at different scales. Graph. Model. 71, 1–21. Prusinkiewicz, P., Lindenmayer, A., 1990. The Algorithmic Beauty of Plants. SpringerVerlag, New York, NY, USA. 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. Remote Sens. 5, 491–520. Ravaglia, J., Bac, A., Fournier, R.A., 2017. Extraction of tubular shapes from dense point clouds: application to stem taper estimation from laser scanned data. Comput. Graph. 66, 22–33. Rowe, J.S., 1972. Forest Regions of Canada. Natural Resources Canada. Canadian Forest Service (Publ. No. 1300). Runions, A., Lane, B., Prusinkiewicz, P., 2007. Modeling trees with a space colonization algorithm. In: Proceedings of the Third Eurographics conference on Natural Phenomena. Eurographics Association, Prague, Czech Republic, pp. 63–70. Seidel, D., Leuschner, C., Müller, A., Krause, B., 2011. Crown plasticity in mixed forests—quantifying asymmetry as a measure of competition using terrestrial laser scanning. For. Ecol. Manag. 261, 2123–2132. Sen, D., Srikanth, T.K., 2008. Efficient computation of volume fractions for multi-material cell complexes in a grid by slicing. Comput. Geosci. 34, 754–782.

114