Individual Tree Attribute Estimation and Uniformity Assessment in Fast-Growing Eucalyptus spp. Forest Plantations Using Lidar and Linear Mixed-Effects Models

Fast-growing Eucalyptus spp. forest plantations and their resultant wood products are economically important and may provide a low-cost means to sequester carbon for greenhouse gas reduction. The development of advanced and optimized frameworks for estimating forest plantation attributes from lidar remote sensing data combined with statistical modeling approaches is a step towards forest inventory operationalization and might improve industry efficiency in monitoring and managing forest resources. In this study, we first developed and tested a framework for modeling individual tree attributes in fast-growing Eucalyptus forest plantation using airborne lidar data and linear mixed-effect models (LME) and assessed the gain in accuracy compared to a conventional linear fixed-effects model (LFE). Second, we evaluated the potential of using the tree-level estimates for determining tree attribute uniformity across different stand ages. In the field, tree measurements, such as tree geolocation, species, genotype, age, height (Ht), and diameter at breast height (dbh) were collected through conventional forest inventory practices, and tree-level aboveground carbon (AGC) was estimated using allometric equations. Individual trees were detected and delineated from lidar-derived canopy height models (CHM), and crown-level metrics (e.g., crown volume and crown projected area) were computed from the lidar 3-D point cloud. Field and lidar-derived crown metrics were combined for ht, dbh, and AGC modeling using an LME. We fitted a varying intercept and slope model, setting species, genotype, and stand (alone and nested) as random effects. For comparison, we also modeled the same attributes using a conventional LFE model. The tree attribute estimates derived from the best LME model were used for assessing forest uniformity at the tree level using the Lorenz curves and Gini coefficient (GC). We successfully detected 96.6% of the trees from the lidar-derived CHM. The best LME model for estimating the tree attributes was composed of the stand as a random effect variable, and canopy height, crown volume, and crown projected area as fixed effects. The %RMSE values for tree-level height, dbh, and AGC were 8.9%, 12.1%, and 23.7% for the LFE model and improved to 7.3%, 7.1%, and 13.6%, respectively, for the LME model. Tree attributes uniformity was assessed with the Lorenz curves and tree-level estimations, especially for the older stands. All stands showed a high level of tree uniformity with GC values approximately 0.2. This study demonstrates that accurate detection of individual trees and their associated crown metrics can be used to estimate Ht, dbh, and AGC stocks as well as forest uniformity in fast-growing Eucalyptus plantations forests using lidar data as inputs to LME models. This further underscores the high potential of our proposed approach to monitor standing stock and growth in Eucalyptus—and similar forest plantations for carbon dynamics and forest product planning.


Introduction
Forest plantations cover approximately 7% of the global forested area, including around 80 million ha in tropical and subtropical countries, and Eucalyptus spp. is one of the most widely cultivated species due to their fast growth rate [1][2][3]. In the context of sustainability, these plantations help to offset natural forest exploitation, restore some ecological services offered by natural forests, and mitigate climate change by storing carbon [4][5][6]. As economic endeavors, careful monitoring of forest plantation growth and productivity is critical for efficient and optimal management, with the objective of generating maximal yields while minimizing production costs and environmental disturbances. From a manager perspective who has to deal with large areas, having detailed information is necessary for this purpose and spatial modeling using remotely sensed data in fine resolutions could be useful in this regard.
Several remote sensing technologies have been used for forest monitoring in the past couple of decades. They provide comparable forest attributes estimations to traditional field survey-based methods while being less laborious and highly time-efficient, especially for areas that are difficult to access [7][8][9]. Light Detection and Ranging (lidar) has gained popularity among the remote sensing methodologies due to its ability to provide accurate and high-resolution 3D measurements of the forest. Airborne lidar data are well suited to calculate and/or predict forest structure attributes such as tree height [10][11][12], diameter at breast height (dbh) and basal area [13,14], stem volume [15,16] and aboveground biomass [17][18][19]. Many of the aforementioned applications have already been studied in the case of Eucalyptus plantations [20][21][22], although more attention needs to be given to the aboveground carbon (AGC) dynamics across different ages [8,23]. AGC information is important to understand the contribution of these forests in CO 2 sequestration, and its incorporation into management plans will help to mitigate climate change [24][25][26].
Modeling and mapping forest attributes using lidar are primarily based on either the area or individual tree-based approaches. In the former, forest attributes are estimated to an area level (e.g., mean canopy height and volume in a plot), whereas the individual tree approach estimates the attributes for each tree (e.g., tree height, dbh, and volume) [9,22,27]. The area-based approach has been more commonly used in forestry applications as it is comparatively more straightforward to integrate with conventional inventory measurements and works with less expensive small-footprint lidar data [28]. However, detailed forest information, such as tree counts, stem height, dbh, and tree biomass, is essential for effective management and quantitative analysis of forests, so the tree-level approach is highly appealing from an inventory standpoint [29][30][31][32]. Therefore, effectively modeling of the attributes at the tree level needs to be explored together with its auxiliary information (e.g., stem distribution and location) for understanding the real impact of this approach on improving management decisions [33].
Tree-level attribute modeling using lidar data are carried out by using allometric relationships between lidar-derived crown metrics (e.g., tree height, crown volume, crown area) and a field observed variable (e.g., tree volume). In Eucalyptus plantations, other critical factors may be considered, such as: (a) hierarchical structure (e.g., plots within forest stands, forest stands within regions, repeated measurements of trees) and (b) crossed grouping structures (e.g., species and genotype variation) [34][35][36]. These may occur due to the measurement of sample plots from multiple stands with different characteristics to generate a representative dataset for model fitting. In these cases, linear mixed-effects models (LME) are found to be more flexible, precise, and accurate compared to ordinary linear fixed-effects models (LFM) [37,38]. LME models account for "random effects" and therefore enable multiple models for different sites and genotypes. The LME assumption is that the individual-specific effects and the independent variables are uncorrelated [39,40]. LME models fitted with data from lidar-supported inventories have been applied to Eucalyptus stands before, with the primary focus being site index and volume estimation at the plot level [34,[41][42][43][44][45]. Nonetheless, combined individual tree level data and mixed-effect models have not been rigorously tested in fast-growing Eucalyptus plantations.
The tree-level estimations allow us to directly retrieve spatial distribution and uniformity of trees across the stand. In Eucalyptus plantations, this information is essential as stand heterogeneity may decrease productivity due to the impact of competition in tree yield [46][47][48]. Accurately identifying heterogeneous areas may also be used as a proxy for defining harvest and thinning goals [49]. One of the methods that have been used to assess tree uniformity is the Lorenz curves and Gini coefficient [50][51][52][53]. These methods have been mostly used with lidar to estimate plot-level forest structure in pine [53,54] and mixed-species [55] stands. However, it is still necessary to explore their application on lidar derived tree-level estimations in Eucalyptus-type stands, which are very different from pines.
In order to provide solutions for increasing forest plantation industry efficiency in monitoring and managing forest resources, the overall goal of this study was to develop and test a novel framework for modeling individual tree attributes, such as height, dbh, and AGC in fast-growing Eucalyptus stands using lidar-derived crown metrics and LME statistical modeling approach. We hypothesized a possible gain in accuracy over conventional LFM models given the known hierarchical structure of the field-measured data in these forest plantations. Herein, we also assessed forest uniformity across different stand ages using the Lorenz curves and Gini coefficient based on the tree-level estimated attributes from the best-fitted LME model. In addition, we discussed the implications of our framework for management decisions and precision forestry.

Study Area Description
The study area was located in southern Brazil within the Telêmaco Borba municipality in the state of Paraná (Figure 1), and it is owned by a pulp and paper company. There were 23 Eucalyptus stands where eight species and five genotypes were planted. Each stand was planted with a single genotype (i.e., stands were single-species and single-genotype). The stands managed to meet the company standards of high productivity and were free from significant damage from fire, disease, or insects. Trees were planted using a 3 × 3 m grid configuration, resulting in an average tree density of 1111 trees/ha. The site is characterized by an average altitude of 741 m with a relatively flat slope, annual average precipitation of approximately 1378 mm, and an annual average temperature of 18.4 • C [56].

Field Data Collection
A total of 25 rectangular (20 × 30 m; 600 m 2 ) plots were randomly established such that the sample plots represented the study area in an unbiased manner and captured forest variability among stands (summary of the plots is presented in Table S1). Individual trees were measured for diameter at 1.3 m above the ground (dbh), and a random subsample of 15% of the trees inside each plot was measured for total height (Ht). Tree heights that were not directly measured were estimated using hypsometric equations employing dbh as the independent variable [57], resulting Ht estimates with mean ± standard deviation (sd) bias of −0.03 ± 0.04 m (−0.23 ± 0.3%). Herein, plot locations were georeferenced using a decimeter accuracy differentially corrected geodetic GNSS (Trimble Pro-XR). Information regarding species, age, and genotype were also noted. AGC (kg/tree) was obtained through an allometric model, employing the logarithm of dbh and Ht as independent variables, based on an approach developed in [9] (see Supplementary Material). In summary, the stand ages ranged from 1.7 to 4.3 years, and the mean ± sd of the stand attributes Ht, dbh, AGC, and age were 16.09 (±3.68 m), 13.46 (±2.65 cm), 33.04 (±14.61 kg/tree), and 3.04 (±0.74 years), respectively.

Lidar Data Acquisition and Pre-Processing
Lidar data and high-resolution imagery were acquired at the same period of field data collection using a Harrier 68i system mounted on a CESSNA 206 aircraft in 2012. The lidar flight parameters are presented in Table 1. Digital Terrain Model (DTM), canopy height model (CHM), and lidar-derived crown metrics were the main products derived from lidar data processing and were created using the US Forest Service FUSION/LDV 4.0 software [58]. First, ground points were identified in the raw 3D point cloud data using Kraus and Pfeifer's algorithm [59], from which a 1 × 1 m resolution DTM was generated using the GridSurfaceCreate function. The height of each raw point was then computed using the Clipdata tool, and a 0.5 × 0.5 m CHM raster was generated using the CanopyModel. The points were then clipped in 3D to the sample plots using the ClipData tool. Lidar crown-level metrics were computed on this final height-normalized all-returns point cloud (item 2.4).

Individual Tree Detection and Crown Metrics Computation
The first step in the tree-based approach is to detect and delineate individual trees. The individual tree detection was performed by applying algorithms to the CHM rasters using the rLiDAR package [60] in R [61]. The entire process consisted of 5 core steps: (1) smoothing the CHM with a 3 × 3 mean filter to remove spurious local maxima caused by tree branches, (2) tree top detection from the smoothed CHM using the local maxima algorithm of the FindTreeCHM function (a 5 × 5 window was used in this case), (3) delineation of tree crowns using the Voronoi Tessellation algorithm based on the ForestCAS function, iv) extraction of all returns of the normalized 3D point cloud within the tree crowns, and v) computation of crown-level metrics-maximum height (HMAX, m), crown projected area (CPA, m 2 ), and crown volume (CV, m 3 ) for individual trees using the CrownMetrics function ( Figure 2). We assumed that these three metrics would be relevant for the tree attributes modeling based on observations made in previous studies [62,63]. In the accuracy assessment, we compared the number of lidar-detected trees per plot with the numbers of trees observed in the field and verified the tree location and crown boundary by overlaying it on the high-resolution RGB image.

Individual Tree Attribute Modeling
To model the tree attributes (i.e., Ht, dbh, and AGC), we spatially linked the field and lidar-detected individual trees within each plot. A total of 1520 reference trees were used, and a preliminary exploratory analysis was conducted using scatterplots and correlation coefficient to identify trends in the data. In all models, the independent fixed-effects variables were the lidar-derived HMAX, CPA, and CV. To accommodate the potential random nature of the variables that were visible during the graphical analysis stage, four LME models [37,38] were developed to estimate the three dependent attributes with varying intercept and slope. The first three models were built using the tree species (SP), genotype (GM), and stand (SID) separately as random-effects variables, and the fourth was built using them in a nested structure. The models were fitted using the lmer function in the "lme4" R package [64]. The general expression of the LME is presented in Equations (1) and (2).
where: Y is the vector of response variables; β is the vector of fixed effects; u is the vector of random effects with E(u) = 0 and u~N(0, σ 2 u D); ε is the vector of residuals with E(ε) = 0 and ε~N(0, σ 2 ε I); X and Z are design matrices relating the observations y to β and u, respectively.
where: Y ij is the response variable for the i-th tree; β 0 and u 0ij are the fixed and random intercepts, respectively; β 1 , β 2 , β 3 are the fixed coefficients associated to HMAX, CPA, and CV, respectively; u 1ij , u 2ij , u 3ij are the random coefficients for the i-th tree and random-effect variable j; ε ij is the random error with E(ε) = 0 and ε ij~N (0, σ 2 ε I). The results were compared with a conventional linear fixed-effect model (LFE) with the same fixed-effects (i.e., independent variables HMAX, CPA, and CV) used in the LME models (Equation (3)). The LFE provided a model based on only fixed effects to test whether including random effects in the LME model significantly improved the results. The R lm function [65,66] was implemented for this purpose. The models were tested for normality and heteroscedasticity of the residuals with the Shapiro-Wilk [67] and Breusch-Pagan tests [68], respectively, for a stand-based estimate.
where: Y i is the response variable of the i-th tree; β 0 , β 1 , β 2 , β 3 are the intercept and slope parameters associated with the lidar-derived variables HMAX, CPA, and CV, respectively; ε i is the random error with E(ε) = 0 and ε i~N (0, σ 2 ε I). The accuracy of each model was assessed according to the coefficient of determination (R 2 ), absolute and relative root mean square error (RMSE and %RMSE, respectively), and absolute and relative bias (%bias) (Equations (4)-(8), respectively). The Akaike information criterion (AIC) was also calculated for comparing and ranking the proposed models [69,70]. To validate the model performance, we computed the aforementioned statistics again using a leave-one-out cross-validation (LOOCV). The acceptable precision was considered to be a %RMSE below 15% and bias near to zero to satisfy conventional forestry inventory standards [23].
where:Ŷ i = estimated tree attribute for the i-th tree; Y i = observed tree attribute for the i-th tree; Y = observed mean tree attribute; n = number of observations.

Assessing Forest Attributes Uniformity
Field and best-fitted LME model derived tree attributes were used to assess forest uniformity at the plot level. Herein, we used the Lorenz curve and Gini coefficient (GC) for their interpretability in describing forest structural differences in terms of uniformity of tree sizes [50,52,53,71]. We calculated the Lorenz curve and GC for describing the uniformity of the three estimated attributes (i.e., Ht, dbh, and AGC). The indices were calculated using the estimated values, and field measurements and the GC values were compared using the two-tailed Student's t-test. The plots were also split into the age classes of 1.5-2.5, 2.5-3.5, and 3.5-4.5 years for observing if the uniformity patterns changed for different ages. In the calculation, the trees were ranked in decreasing order with the intention to depict the dominance of larger trees in the y-axis [46,52], and the Lorenz curves were built by plotting the cumulative percentage of trees against the cumulative percentage of each tree attribute. The GC was then calculated as the ratio of the area contained under the empirical Lorenz curve by the area above the diagonal line of perfect equality (when all trees are identical in size). The GC thus theoretically ranges from 0 (all trees of equal size) to 1 (highest theoretical dispersion where only one observation holds the sum of the attribute against many tree sizes which are negligible, e.g., by not reaching the breast height). In forest science, the Lorenz curves are compared against the Lorenz curve of maximum entropy (depicted by a grey area in Figure 3), which represents the most uneven-sized situation with trees ranging from all types of sizes [52]. Adnan et al. [72] showed that maximum entropy could be deducted for fixed threshold values of GC. The authors demonstrated that these thresholds are different depending on whether the tree attributes used are one-dimensional (e.g., tree height or dbh) or two-dimensional (e.g., basal area or crown area) measures and determined the thresholds GC = 0.33 (one-dimensional) and GC = 0.50 (two-dimensional) which correspond respectively to each of the aforementioned groups. Accordingly, in this research, we employed a threshold of GC = 0.33 to depict the maximum entropy (i.e., plots with GC values below the threshold were considered as single-cohort even-sized forests) [72]. Moreover, the inflection point of the Lorenz curve depicts the mean value of the tree attribute [51,71,73]. The possible theoretical scenarios for describing the stands are uneven-sized stands with balance among tree cohorts (i.e., symmetric; Figure 3a), uneven-sized stands with a dominance of small trees (i.e., negatively asymmetric; Figure 3b), or dominance of large trees (i.e., positively asymmetric; Figure 3c) and even-sized stands with all trees of roughly the same size (curves that fall below the curve of maximum entropy, within the grey-shaded area- Figure 3d). Figure 4 offers a synopsis of the complete methodological workflow.

Individual Tree Detection and Crown Metrics
From the 1520 reference trees, 1468 (96.6%) were correctly detected (missed trees = 20; falsely detected trees = 32). A 1:1 comparison of the number of field-observed and lidar-detected trees by plot shows %RMSE and %bias less than 5% and 0.8%, respectively (Figure 5a). When assessing the lidar-derived crown-level metrics (Table S2), we found HMAX to be highly correlated (r > 0.78) with the response variables, whereas CV and CPA showed very low correlations (r < 0.3). This is further highlighted when analyzing the scatterplots between HMAX and field attributes across the random-effect variables ( Figure 6). Figure 6. Scatterplots of the lidar-measured maximum height (HMAX) against the field-measured height (Ht), diameter at 1.3 m above the ground (dbh), and aboveground carbon (AGC) (a, b and c, respectively). Colors represent the levels of species (SP) (a1,b1,c1), Genotype (GM) (a2,b2,c2), and Stand (SID) (a3,b3,c3) used as random-effect variables in the linear mixed-effect models. Plotted-lines were fitted to each group to better demonstrate the different patterns for each level.

Predictive Models
Overall, the LME performed better than the LFE models ( Table 2). The most accurate LME was built with the SID as a random-effect variable. The fitted equations using SID as a random-effect variable for Ht, dbh, and AGC explained 90%, 87%, and 91% of the attribute's variation, respectively, while the LFE models only did 84%, 62%, and 71% for the same tree attributes. The latter models also satisfied residual normality and homoscedasticity assumptions for the stand estimations (W > 0.83 and p-value > 0.09 in the Shapiro-Wilk test, and BP > 1.6 and p-value > 0.13 in the Breusch-Pagan test). Even though the LME with nested SID, SP, and GM had very similar results to using only SID, it lacked simplicity, as shown by the AIC values-9131.3 for the former and 9119.3 for the latter. Therefore, the LME model with SID as a random-effect was selected for the subsequent analyses. In all the best-fitted LME, the random-effects were used in the intercept and HMAX-associated parameters only, as doing so for the other variables did not improve performance (Table S3). The HMAX was the most important variable in all models, as commonly seen in allometric models using lidar metrics [10], while CV did not contribute significantly to the models (p < 0.05). The CPA was not significant only for the Ht estimation with the LFE model. The difference among LFE and the LME model with SID as a random-effect variable was also observed in the LOOCV assessment (Figure 7). The %RMSE decreased by approximately 1.5, 5, and 10 percent points for Ht, dbh, and AGC, respectively, when using the LME with SID as a random-effect variable rather than LFE. The differences in %bias were not as significant where the decrease was 0, 0.01, and 0.04 percent points with all models having a value lower than 0.3.
The Lorenz curves calculated from the field measurements and the best LME model estimations were overall similar (Figure 8). The main differences occurred in younger ages (1.5-2.5 years), especially for the AGC estimation. As expected, all the curves fell within the maximum entropy range. The mean attribute values (inflexion points checked with crossed symbol) of the younger stands were evenly distributed around the axis of symmetry (Figure 8a1,b1,c1). For stands of other ages, the values were mostly above this axis, representing a dominance of single cohorts of larger trees. The GC values all ranged around 0.1. Values slightly over 0.2 were found only in the coefficients calculated for the field measured AGC in the age classes 1.5-2.5 and 2.5-3.5 years. Significant differences between the GC using the field measurements and estimates with the LME with SID occurred mostly at the stands of 1.5-2.5 years, with the p-values being less than 0.002.

Figure 7.
Field measured (observed) and leave-one-out cross-validation (LOOCV) estimated tree attributes using linear fixed-and mixed-effects models (a,b, respectively) with the SID as a random-effect variable. LOOCV estimations for height (Ht), diameter at 1.3 m above the ground (dbh), and aboveground carbon (AGC) (a1,b1,a2,b2,a3,b3, respectively). The red dashed line presents the fit between observed and estimated values, and the blue line is the equality line (intercept = 0, slope = 1).

Figure 8.
Lorenz curves and Gini coefficient (GC) of the tree attributes as calculated from the field-measured observations (red) and estimated values (blue) using the linear mixed-effects model (LME) with the stand (SID) as a random-effect variable. The tree attributes were aboveground carbon (AGC), diameter at 1.3 m above the ground (dbh), and height (Ht) (a-c, respectively) at ages of 1.5-2.5, 2.5-3.5, and 3.5-4.5 (subscript 1, 2, and 3, respectively). Curves and coefficients were calculated for each sample plot. The p-value (p) for the Student's t-test of mean differences between the field and LME derived GC values.

Discussion
Our study explored tree-level information from lidar for enhancing the forest management decision-making processes and underscores the importance of incorporating random effects while generating predictive models for individual tree attributes in forest plantations. Herein, linear mixed-effects models were found to outperform ordinary least square regression models, which only consider fixed effects. By demonstrating a forest uniformity estimation procedure, we hope to translate our research approach into future forest inventory practices and subsequently bolster ongoing data collection and production endeavors.
The application of the Silva et al. [60] algorithm for tree segmentation was able to detect more than 95% of the trees in our work. This constitutes an outstanding result when compared to other studies involving the tree individualization approach [74] and serves as an example of improvements to lidar-based forest inventory. This relatively new paradigm is therefore perfect for tree-level attribute estimation in Eucalyptus plantations where such in-depth detail is necessary for accommodating small-scale variations that are untraceable via plot-based estimation procedures [21,[75][76][77][78][79]. Herein, the lidar-derived crown-metrics were found to provide a satisfactory estimation of Ht, dbh, and AGC, and the best results were obtained while employing a LME with SID as the random-effect variable. This variable is related to the site and would have likely enabled us to account for the existing grouping characteristics and growth variances of the tree attributes, improving the model performance. A similar trend was observed in several previous studies as well, where authors highlighted the better performance of LME models over LFE models [35,45,80,81]. One of the reasons might be related to the spatial correlation between the samples, which makes them not completely independent, and this introduces correlation into the model residuals [45,82,83]. This characteristic is not desired in LFE models [45,82,83]. On comparing our findings with similar studies in Eucalyptus stands, the major takeaway is that we were able to estimate forest attributes at the tree level and assess forest uniformity as opposed to only estimating forest attributes either at the tree or plot level. For instance, Packalen et al. [34] focused on using LME for estimating site index at the plot level, while Vauhkonen et al. [42] focused on evaluating and comparing area and tree-based lidar approaches for estimating stand volume with LME models in Eucalyptus plantations by applying genotype and stand as random-effect variables.
The lidar-derived crown-level metrics had different contributions in the models. The CPA and CV were not as important as the HMAX but were still found significant to most of the models. Previous studies provide some evidence that CPA and CV metrics might be important to biomass and other forest attribute estimations [62,63,[84][85][86]. Crown-level field measurements are not commonly collected in Eucalyptus plantations but can be considered in further studies to infer their real impact for lidar-based allometric modeling in these forests. There is a cost (time and money) related to measuring CPA and CV variables in the field, and a more straightforward way as described by Maltamo et al. [86] could be obtaining only the canopy base height as a reference. The crown information could be an indicator of growth potential and tree health and give insights on stem quality and tree competition affecting tree growth. Therefore, crown-related metrics, although less important in our study, need to be further understood to be incorporated into tree attribute models in Eucalyptus plantations.
In Eucalyptus plantations, tree heterogeneity might not be desirable as unbalanced competition could affect tree yield especially when it is present at an early stage [32] and decrease stand productivity [46][47][48]87]. Identifying heterogeneous areas can aid forest managers in silvicultural decisions and growth modeling but is limited due to constraints associated with time, labor and financial needs; assessing it from remotely sensed detected trees becomes an alternative in such cases. Hentz et al. [88] reported that the errors incurred due to tree height calculation discrepancies do not significantly affect the uniformity index results given these errors are distributed evenly throughout the study area. Therefore, the concerns associated with individual tree detection are only the major ones to be noted while choosing between lidar and field-based studies. For plantation studies, such as ours, lidar-derived crown-level metrics were effective for calculating uniformity. However, the same might not be recommended in forests with more complex structures [73].
We used the Lorenz curves and GC to assess forest uniformity at tree level and found that the estimations fairly depicted the tree attributes variation, mainly in stands over 2.5 years. Trees in the given Eucalyptus stand were found to be quite similar, which was evident from the low values of GC (<0.33) and the close proximity of Lorenz curves to the perfect equality line in Figure 8. The differences found in younger stands may be related to errors associated with tree detection and crown delineation at this age, which is commonly the main source of errors in tree-level approaches [79,89]. Further research may clarify the relationships between the theoretical threshold GC < 0.33 [72] and productivity loss in Eucalyptus plantations. This will allow us to verify if the heterogeneity affects tree growth rates. Another important point to consider is the scale-factor associated with GC calculation. As described by Adnan et al. [90], GC tends to zero when assessed at a finer scale (i.e., using a considerably smaller plot), while at coarser scales, GC can increase asymptotically toward a maximum value that best describes the size heterogeneity of uneven-sized forests.
Despite the advantages of using lidar data for forest inventory, there is a cost-benefit trade-off associated with obtaining and processing data at more precise and fine-scale localized levels, which comes from lidar flight cost and hardware limitations (e.g., memory and processing); this should be circumventable in the long run given the continuing progress of science and technology. Even then, the scalability, reusability, and transferability of the proposed approach is expected to outweigh the capital cost limitations. For instance, the same LME models presented here could be easily recalibrated to become applicable to different SID levels coming from new field measurements [91][92][93] (example in the Supplementary Material). In addition, the lidar data bring the benefit of directly extrapolating the estimations for the entire area. Our results also inform that LME models are viable alternatives for scenarios where numerous levels need to be considered (e.g., different sites and genotypes), and stratifying the data may not be deemed feasible as it may lead to considerably fewer samples per level. Additionally, the same methods can be extended to other types of forest plantations as well for determining stand variability and/or for improving management decisions.
The novel framework presented in this study can be perceived as a step forward within the precision forestry realm if implemented [94]. Tree-level lidar studies are currently not undertaken at large-scale plantation areas. Yet, fine-scale attribute estimations and their uniformity can be delivered through the operationalization of the proposed approach (Figure 9). Such information allows localized silviculture treatments associated with irrigation, disease detection, growth imbalances, and pest control. This also enhances forest management efficiency as we would have better estimates of the attributes and variability within the stands. Presently, while performing forest inventories, several forest industrial companies use only a small number of randomly established plots (~1 everỹ 15-20 ha), and they can be placed in an unrepresentative region with respect to stand variability, and this can introduce bias into their forest attribute estimations. Thus, sampling plot allocation and concomitant costs can be optimized using lidar data [95][96][97][98]. For harvest planning, managers can use detailed information to estimate potential yield, harvest time, and volume stock for specific regions in the stand or define selective logging objectives [35]. Our strategy could also facilitate the application of individual-tree level growth and yield models, calculation of competition indices, identification of areas of low productivity (e.g., due to low soil fertility, pests, or disease), and identification of dominant trees for site index classification. It is worth mentioning that the resolution of the point cloud data used to achieve the results in this study was not as high as the currently available state-of-the-art lidar technologies. That might make it feasible to apply these methods to large areas at a lower cost. Given all the benefits in relation to large-scale spatialization, enhanced individual tree attributes estimation, and forest uniformity assessment, we expect our workflow to be implemented extensively in the coming years not only in Eucalyptus plantations but also in diverse species varieties globally.

Conclusions
In this study, we demonstrated that Ht, dbh, and AGC individual tree attributes could be accurately estimated in a fast-growing Eucalyptus forest plantation using lidar-derived crown-level metrics combined with LME models. Our framework also allowed assessing forest uniformity at the tree level using the LME model-derived estimates. The LME model using only the stand as a random-effect variable, was able to effectively account for the grouping variability characteristics of the data and resulted in the best model. This further highlights the importance of considering random effects into the modeling approach for estimating individual tree attributes in forest plantations, as the observations usually used are not completely independent due to, e.g., spatial correlation. Understanding these effects is important for supporting the development of robust models by incorporating sources of variation and reducing error in tree attribute estimations. By presenting a transferrable and scalable state-of-the-art approach that emphasizes the need to integrate spatial data with management decision-making processes, we hope to boost outreach capabilities of precision forestry-related endeavors worldwide.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/21/3599/s1, Table S1: Descriptive statistics (mean ± standard deviation) of field measurements in the area of study. Ht = height; Dbh = diameter at 1.3 m above the ground; AGC = aboveground carbon, Table S2: Summary statistics of lidar derived crown-metrics maximum height (HMAX), crown projected area (CPA), and crown volume (CV), Table S3: Statistics for the linear mixed-effect model (LME) using the stand (SID) as random-effect variable considering intercept (u 0i j ) and slope (u 1i j , u 2i j and u 3i j ) combinations. RMSE = root mean square error, R 2 = Coefficient of determination; AIC = Akaike information criterion, Figure S1: Scatterplot comparing recalibrated linear mixed-effect model and fixed-effect models (LME and LFE, respectively) for estimating aboveground carbon (AGC). R 2 = coefficient of determination; %RMSE = relative root mean square error.