Understanding Tree-to-Tree Variations in Stone Pine (Pinus pinea L.) Cone Production Using Terrestrial Laser Scanner

: Kernels found in stone pinecones are of great economic value, often surpassing timber income for most forest owners. Visually evaluating cone production on standing trees is challenging since the cones are located in the sun-exposed part of the crown, and covered by two vegetative shoots. Very few studies were carried out in evaluating how new remote sensing technologies such as terrestrial laser scanners (TLS) can be used in assessing cone production, or in trying to explain the tree-to-tree variability within a given stand. Using data from 129 trees in 26 plots located in the Spanish Northern Plateau, the gain observed by using TLS data when compared to traditional inventory data in predicting the presence, the number, and the average weight of the cones in an individual tree was evaluated. The models using TLS-derived metrics consistently showed better ﬁt statistics, when compared to models using traditional inventory data pertaining to site and tree levels. Crown dimensions such as projected crown area and crown volume, crown density, and crown asymmetry were the key TLS-derived drivers in understanding the variability in inter-tree cone production. These results underline the importance of crown characteristics in assessing cone production in stone pine. Moreover, as cone production (number of cones and average weight) is higher in crowns with lower density, the use of crown pruning, abandoned over 30 years ago, might be the key to increasing production in combination with stand density management.


Introduction
Three-dimensional arrangement of vegetation is of great importance in natural and in managed forest ecosystems. Canopy structure drives many functions such as light interception, habitat for fauna, and microenvironment regulation. Considering the inaccessibility and the complexity of large tree crown structures, it received little attention in the past. In the past decade, terrestrial laser scanning was increasingly used in many fields of forest science. These tools provide very accurate representation of the environment and allow unprecedented measurement of forest attributes at the stand and at the tree scale [1]. Traditional inventory metrics such as tree diameter and height can be extracted from TLS data, as well as precise positions of the trees within the stand [2,3]. A plethora of crown shape, crown volume, or crown deformation metrics can also be estimated [4,5], opening the way to accurately studying crown plasticity on large trees [6][7][8]. Moreover, material arrangement inside the crown attributes are rarely considered, given the difficulty in measuring them. When used, crown diameter (measured as the vertical limit of the crown in one or two directions) or volume, estimated using crown diameter, depth, and shape, was found to explain some of the tree-to-tree variability (Rodrigues et al. 2014). These measures remain approximate, and detailed crown characteristics could help in furthering our understanding of why some trees in a stand produce more cones when compared to others. We believe that TLS technology can help in modeling crown production and understanding its functional drivers.
The objective of the present study is to better understand the relationship between tree characteristics and cone production in Mediterranean stone pine. This was achieved by quantifying the gain in using detailed crown metrics in estimating cone production at the individual tree level (number of cones per tree and average cone weight). Models based on traditional variables (tree size and stand characteristics) were compared to models that relied on crown metrics extracted from TLS data. The resulting models should help owners and managers to better predict cone production.

Study Region
The study sites were located on the Spanish Northern Plateau, a 50,000 km 2 plain area defined by the Duero river basin, with altitudes ranging from 600-900 m. This territory is characterized by a continental Mediterranean climate, with an average temperature of 12.6 • C, annual rainfall over 440 mm, and severe summer drought, resulting in an aridity index [28] of 0.60 (dry sub-humid). Soils have a high content of sand (>90%), with low content of nutrients and low water retention capacity. Only in the limestone areas (less than 15% of total area) do the soils have a higher clay content.
The forests within the region are pure pinewoods of Pinus pinea and Pinus pinaster, and mixtures of both species with Quercus ilex, Quercus faginea, and Juniperus thurifera. Pure Pinus pinea forests were retained, which represent more than 50% of the total forest area of the region, occupying more than 70,000 ha. Management of these forests is orientated to promote joint timber and cone production; thus, typical interventions consist of heavy and early thinnings to favor horizontal crown expansion, promotion of even-aged stands by application of shelterwood systems, and rotation lengths over 120 years.

Permanent Plots for Monitoring Cone Production
In 1996, a network of 141 permanent plots was installed in the study region. The main objective was to analyze the effect of forest management on timber and cone production ( Figure 1). The complete ranges of site quality, stocking, and stand age were covered by the plots. The variable radius plots consisted of measuring the 20 closest trees to the plot center. At plot installation, diameter at breast height, total tree height, height to crown base, and crown diameter were measured on every tree, and tree coordinates were also recorded. The plots were remeasured on three occasions (2001, 2008 and 2016).
Cones were manually collected in all the plots every autumn from 1996 to 2005 by specialized climbers in the five trees closest to the center of the plot. The cones of each tree were counted and weighted. Since 2006, yearly cone collection is carried out on a subset of 30 plots. Only the cone production for 2016 and 2017 was used in the present work, as TLS acquisition was carried out in the fall of 2016. Furthermore, four of the 30 plots were not used as they were isolated from the other 26 plots, due to logistical constraints. Finally, a tree from one plot was dropped since the tree died after TLS data acquisition and cone collection, leading to 129 trees being used in the study (26 plots × 5 trees/plot-1 dead tree). Table 1 presents the summary of plot, tree, and TLS-derived variables. Cones were manually collected in all the plots every autumn from 1996 to 2005 by specialized climbers in the five trees closest to the center of the plot. The cones of each tree were counted and weighted. Since 2006, yearly cone collection is carried out on a subset of 30 plots. Only the cone production for 2016 and 2017 was used in the present work, as TLS acquisition was carried out in the fall of 2016. Furthermore, four of the 30 plots were not used as they were isolated from the other 26 plots, due to logistical constraints. Finally, a tree from one plot was dropped since the tree died after TLS data acquisition and cone collection, leading to 129 trees being used in the study (26 plots × 5 trees/plot-1 dead tree). Table 1 presents the summary of plot, tree, and TLS-derived variables.

Terrestrial Laser Scanner Data Acquisition and Metric Extraction
The scanner used was a Faro Focus 3D. Each scan was performed with a resolution angle of 0.036 • in both horizontal and vertical planes in combination with rotation angles of 360 • horizontally and of 150 • (from −60 • to 90 • ) vertically. A multiscan approach consisting of merging several scans to increase three-dimensional (3D) representation accuracy was used. The individual scans were registered using 6-8 30 cm spherical targets distributed throughout the plot, and co-registration between scans was carried out with Faro Scene 6.2.2.7 (Faro Technologies Inc., Lake Mary, FL, USA). Finally, the stray point and dark filters were performed to remove noise points. Plots were scanned from five positions in October 2016, with the Faro Focus 3D placed for one scan in the center of the plot. The other four positions were chosen to minimize occlusion of the target trees and ensure that the crown of each sample tree was at least seen from three sides.
The co-registered point clouds were then exported to Computree (Othmani et al. 2011) to segment ground returns from tree points. The digital terrain model (DTM) was also obtained from CompuTree as a raster with 30 cm cells. The points from each sample tree were then manually isolated in R (R Development Core Team 2011) and CloudCompare (CloudCompare (GPL software) 2015). The DTM and isolated trees were then used to extract crown and micro-topographical metrics.

Crown Metrics
Crown base of the sample trees was estimated using a segmented regression with the segmented package in R [29,30]. Convex hulls were fitted on the XY-coordinates of sub point clouds every 10 cm along the Z-axis, and the maximum distance between the center and the vertices of the convex hull was recorded [31]. Crown base is defined as the lowest breakpoint of theses distances following the Z-axis, which corresponds to the point where the distance between the center and the vertices of the convex hull start to increase because of the presence of branches ( Figure 2). The points above the crown base were then used to calculate the following previously published crown metrics [6]: Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 16

Terrestrial Laser Scanner Data Acquisition and Metric Extraction
The scanner used was a Faro Focus 3D. Each scan was performed with a resolution angle of 0.036° in both horizontal and vertical planes in combination with rotation angles of 360° horizontally and of 150° (from −60° to 90°) vertically. A multiscan approach consisting of merging several scans to increase three-dimensional (3D) representation accuracy was used. The individual scans were registered using 6-8 30 cm spherical targets distributed throughout the plot, and co-registration between scans was carried out with Faro Scene 6.2.2.7 (Faro Technologies Inc., Lake Mary, FL, USA). Finally, the stray point and dark filters were performed to remove noise points. Plots were scanned from five positions in October 2016, with the Faro Focus 3D placed for one scan in the center of the plot. The other four positions were chosen to minimize occlusion of the target trees and ensure that the crown of each sample tree was at least seen from three sides.
The co-registered point clouds were then exported to Computree (Othmani et al. 2011) to segment ground returns from tree points. The digital terrain model (DTM) was also obtained from CompuTree as a raster with 30 cm cells. The points from each sample tree were then manually isolated in R (R Development Core Team 2011) and CloudCompare (CloudCompare (GPL software) 2015). The DTM and isolated trees were then used to extract crown and micro-topographical metrics.

Crown Metrics
Crown base of the sample trees was estimated using a segmented regression with the segmented package in R [29,30]. Convex hulls were fitted on the XY-coordinates of sub point clouds every 10 cm along the Z-axis, and the maximum distance between the center and the vertices of the convex hull was recorded [31]. Crown base is defined as the lowest breakpoint of theses distances following the Z-axis, which corresponds to the point where the distance between the center and the vertices of the convex hull start to increase because of the presence of branches ( Figure 2). The points above the crown base were then used to calculate the following previously published crown metrics [6]:  -Crown volume (CV) and surface area (SA): volume and area of a 3D alpha-shape fitted to the points belonging to the crown cloud points; -Light crown volume (LCV) and surface area (LCSA): the farthest distance from the crown center of each 10-cm slice of the crown was used to establish the maximum limits of the crown. A linear-log segmented regression was then used to find the transition between the shade and sun-exposed parts [32]. The volume and surface area of the crown above the transition zone were then calculated by fitting a 3D alpha-shape to the sun-exposed crown; -Height of the base of the sun-exposed crown (LCB): breakpoint of the linear-log segmented regression; -Crown length (CL): difference between total tree height and height to crown base; -Projected crown area (PCA): area of the projected crown on the XY-plane; -Maximum crown radius (maxCr): maximum of all the radii measured in each 10-cm slice of the crown; -Crown asymmetry (ASYM): ratio of the smallest and largest crown radius of the projected surface area; -Crown density: ratio of volume of non-empty 10 cm edged crown voxels to total crown volume, calculated for the complete crown (CDen) and the light crown only (LCDen).

Micro-Topography
The surrounding topography of each sample tree was quantified using the following metrics calculated with the 30-cm cell DTM raster (Figure 1 Slope: local slope calculated around the sample tree [36] using the slope method in the SDMTools package in R [37]. The micro-topographic metrics were calculated using 25 cell matrices (approximately 75-cm radius around the sample tree): the center cell of the matrix and the 24 neighboring cells of the DTM.

Data Analysis
The variables were divided into three groups: plot, tree, and TLS-derived metrics (Table 1). During preliminary analysis, high correlations were observed between several LiDAR-derived variables. In order to reduce the chances of observing (multi) collinearity during the modeling process, an initial principal component analysis (PCA) was carried out on the TLS-derived variables [38]. This was done to evaluate how correlated these were, and if variables were similar, to form groups of variables to be used in the modeling approach [39]. PCA was selected as it can easily detect patterns in a dataset by decomposing the variance-covariance matrix, and it is widely used to reduce the number of explanatory variables [39]. Furthermore, no a priori hypotheses were used to group the TLS-derived variables.
Models were then fitted for the following variable groups: (1) plot, (2) tree, (3) TLS, (4) plot and tree and (5) plot and TLS. For TLS variables, only one variable per group determined by the PCA results was included in a given model, i.e., if group 1 had variables W and X, and group 2 was composed of variables Y and Z, a model could have variable W, W, and Y or W and Z but not W and X. In the fifth case (plot and TLS), plot-level variables that could be obtained by TLS were kept (i.e., stand density (SD, st·ha −1 ), stand basal area (BA, m 2 ·ha −1 ), mean quadratic diameter (MQD, cm), and dominant stand height (domH, m). For each analysis, all of the possible combinations of potential variables were tested, and the best model was the one with the lowest Akaike's information criterion (AIC, for zero-inflated negative binomial regression) or Akaike's information criterion corrected for small samples (AICc, for linear mixed model). Thus, a total of 63 models were compared for the plot-level variables (group 1), 15 for the tree level models (group 2), 449 for the TLS-derived variables (group 3), 1023 for the combined plot and tree levels (group 4), and 6735 for the combined plot and TLS levels (group 5).
where g(y i ) is defined as The mean µ i is then related to the covariates as follows: where X and β are the covariate matrix and estimated fixed effect vector, respectively. θ i is given by the logistic function where where Z and γ are the covariate matrix and estimated fixed effect vector, respectively. Finally, plot-level random effects were included as µ i and v i , where µ i ∼ N 0, σ 2 1i and v i ∼ N 0, σ 2 2i . The model was calibrated using the GLMMadaptive library in R [41].
Variable selection was carried out for the negative binomial (Equation (3)) and logistic regression (Equation (4)) separately for each of the variable groups. These were calibrated using the lme4 library in R [42]. Model comparison was carried out using the AICcmodavg package in R [43]. The logistic and negative binomial parts were then combined in the ZINB model, and variables with a significance level larger than 0.10 (i.e., p > 0.10) were discarded. Some variables were standardized in order to allow for convergence. The standardized variables are indicated with an asterisk, i.e., TPI* is the standardized value of TPI.

Average Cone Weight
Average cone weight (w ij ) was predicted using a linear mixed effect model including a plot and year crossed random effect (Equation (6)).
where T and δ are the covariable matrix and fixed effect parameter vector, respectively. s i , t i and ijk are the plot and year random effects and residual error, where s i ∼ N 0, σ 2 3i , t i ∼ N 0, σ 2 4i and ijk ∼ N 0, σ 2 ijk . These models were also calibrated using the lme4 package in R.

Results
The PCA analysis indicated that the TLS-derived variables could be grouped into four categories ( Figure 3). The first group contained variables quantifying the local topography around a tree: TPI, TRI, slope, and roughness. The second was composed of variables that measure crown dimensions: crown volume and surface area, light crown volume and surface area, average crown and average crown projected area radius, maximum crown radius, height to light crown base, and crown length. The third was related to crown asymmetry: crown asymmetry and smallest projected crown radius. Finally, the fourth group estimated the density of the crown: total crown density and light crown density. The mean intra-and inter-group Euclidian distances indicated that variables from a given group were closer than the groups ( Table 2). In the analysis using TLS metrics, only one of the variables of a given group could be found in a given model. In other words, if TPI was included in a model, the other topographic variables (i.e., TRI, slope, or roughness) were not included.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 16 crown volume and surface area, light crown volume and surface area, average crown and average crown projected area radius, maximum crown radius, height to light crown base, and crown length. The third was related to crown asymmetry: crown asymmetry and smallest projected crown radius. Finally, the fourth group estimated the density of the crown: total crown density and light crown density. The mean intra-and inter-group Euclidian distances indicated that variables from a given group were closer than the groups ( Table 2). In the analysis using TLS metrics, only one of the variables of a given group could be found in a given model. In other words, if TPI was included in a model, the other topographic variables (i.e., TRI, slope, or roughness) were not included.  Based only on plot-level variables, the average number of cones a tree has was conditioned by the mean quadratic diameter, stand density, site index, and stand basal area (Equations (7) and (8)).
When uniquely considering the tree-level variables, the best model indicated that the presence of cones increased with crown diameter (−1.0944, p = 0.0008, Equation (9)), and the number of cones increased with both tree diameter (0.0432, p = 0.0054, Equation (10)) and crown length (0.2307, p = 0.0024).
= exp (5.5133 − 1.0944 CD) = exp (−1.4369 + 0.0432 DBH + 0.2307 CL) When plot-and tree-level variables were combined, plot-level variables were included only in the negative binomial (abundance) part of the model (Equations (11) and (12)).   Based only on plot-level variables, the average number of cones a tree has was conditioned by the mean quadratic diameter, stand density, site index, and stand basal area (Equations (7) and (8)).
The probability of observing a tree with cones increased with the mean quadratic diameter (−0.1698, p = 0.0647). The number of cones, when present, increased with both site index (0.1143, p = 0.0781) and basal area (0.0521, p = 0.0274) and decreased with stand density (−0.0087, p = 0.0008).
When uniquely considering the tree-level variables, the best model indicated that the presence of cones increased with crown diameter (−1.0944, p = 0.0008, Equation (9)), and the number of cones increased with both tree diameter (0.0432, p = 0.0054, Equation (10)) and crown length (0.2307, p = 0.0024).
When TLS-derived metrics were used to predict the number of cones, crown dimension and asymmetry variables entered the model for the presence of cones, and crown dimension and crown density group variables influenced the number of cones (Equations (13) and (14)).
The probability of observing cones increased with the projected crown area (−2.0151, p = 0.0089) and the minimum crown radius (−0.9438, p = 0.0152). When the cones were present, their number increased with crown volume (scaled, 0.5019, p < 0.0001) and decreased with crown density (−2.3075, p = 0.0630).
The models with TLS-derived variables performed the best, based on AIC (Table 3). Of the two models with TLS covariates, the plot and TLS model had the lowest AIC. The AIC weights for both models containing TLS-derived variables summed to 0.999 (Table 3). These weights can be interpreted as the summed conditional probability for both models [44]. The next models were the ones with tree-level variables, with the best one having both plot-and tree-level variables. The model showing the worst fit was the one with only plot-level information. Table 3. Akaike's information criterion (AIC) and variance estimates for the zero-inflated negative binomial regressions (Equations (3) and (5)). The average weight of the cones was proportional to the site index when only plot-level variables were used (6.646, p = 0.0292, Equation (17)). The average weight was found to be proportional to tree height in the case of the models with tree-level variables (5.359, p = 0.0473, Equation (18)). When both plot-and tree-level variables were used, the model with the lowest AIC was the same as the tree-level model (Equation (18)). The average weight was proportional to the surface area of the light crown (0.2543, p = 0.0423, Equation (19)) and inversely proportional to the light crown density (−103.3775, p = 0.0433) when TLS-derived variables were considered. Finally, roughness (88.1099, p = 0.1497, Equation (20)) and stand density (−0.0940, p = 0.3214) were also included with the light crown surface area (0.2517, p = 0.0447) and light crown density (−88.5663, p = 0.0853) in the plot and TLS model.

Number of Variables in
As with the zero-inflated negative binomial models, the average cone weight models that included TLS-derived variables had the lowest AIC (Table 4), with the plot and TLS model having the lowest AICc. Moreover, the AICc weights supported the models with TLS-derived variables, with both models containing TLS variables accounting for 0.98 of the weight ( Table 4). The models with plot-and/or tree-level variables had AICc values that were very close, with a spread of 1.23 points between the best and worst. These models, thus, could not be very well discriminated (Burnham and Anderson, 2002). Table 4. Akaike's information criterion (AIC) and variance estimates for the mean cone average weight (Equation (6)).

Model
Number

Discussion
Inter-site and inter-annual stone pinecone production, defined as either cone presence, number of cones in a tree, or average cone weight, is fairly well understood. Tree size [21,45], competition [46,47] (Freire, 2009, Piqué et al., 2015, stand stocking [27,48], stand maturity, and site conditions [49] all influence the presence and number of cones found in a tree. Certain crown dimensions, such as crown length, which positively enters the model for cone production, can be considered a good indicator of inter-tree competition for light. However, the differences in production between two trees in the same stand, of similar status and dimensions, are not well understood. It was, thus, intended to use more detailed crown and microsite information obtained by TLS to better understand the between-tree variations. In order to quantify the improvements by the addition of TLS variables in understanding tree-to-tree variability, models using inventory-based plot-and tree-level variables were compared to one using TLS-derived metrics and one using all potential variables (inventory and TLS). The probability of a tree bearing cones and the total number of cones were positively related with tree and crown size and negatively with stand density. Basal area, which reduces cone production, was statistically significant when stand density was also present, and it can be considered as a proxy of both stand maturity and site index. As for average cone weight, it was found to increase with tree height and site index. These results, pointing to higher cone production in older trees with larger crowns growing under low stocking levels, were already reported for stone pine [19,50], as well as for other species [51][52][53].
Water availability in the soil is important for female flower emergence, cone survival, and enlargement [26]. Water availability is intrinsically related to topography. At larger scales, inter-site water availability is explained by total rainfall and soil texture. Irregularity of the terrain influences the intra-site water availability. This was captured by the roughness index obtained from the digital terrain model surrounding the tree. Higher roughness indices indicated more irregular terrain, which could be due to depressions and microbasins. This would, thus, indicate that trees found on such microtopographical locations have access to more water, leading to cones with more biomass.
Very few studies looked at the links between crown dimensions and cone production. Rodrigues et al. [27] found that trees with larger crowns have more cone biomass. Crown width or length, measured during field surveys as measured by Rodrigues et al. [27], are imprecise when compared to metrics that can be obtained from remote sensing such as TLS. Of all the tested TLS covariates, crown projected area and minimum crown radius replaced inventory-based crown diameter when trying to segregate trees with and without cones in the present study. While crown diameter and crown projected area both measure the same dimensions (larger crowns have more probability of bearing cones), trees with a smaller minimum crown radius for a given projected area have an asymmetric crown. Asymmetric crowns are largely associated with lateral inter-tree competition which prevents crown expansion on all sides, reducing the probability of observing cones in the tree. It is interesting to point out that adjacent trees on unmanaged forests can have intermingled crowns such that they ultimately share a single umbrella shaped crown in order to maximize light acquisition [18]. Such fused crowns were, however, not observed in the dataset, since our study area was mainly composed of forests intensively managed to promote cone production, where early and heavy thinning is applied aiming at maximum lateral crown expansion but preventing crown interlacing.
The number of cones in a stone pine tree is determined by both the volume and the density of the crown. TLS-derived crown volume captures irregularities and asymmetries in the crown shape, going beyond the common simplification of circular-shape projected area and regular ellipsoid approximation of crown volume. In both cases (either from ground-estimated or TLS-measured crown volume), trees with larger crowns have more cones, as the number of shoots that can bear cones increases with crown size. The decrease in the number of cones with crown density could be explained by two hypotheses. Firstly, stone pine crown architecture varies with tree ontogeny. The spherical crown with foliage covering the whole external surface of young trees shifts toward a smoothly rounded umbrella-like shape with foliage only covering the upper, fully sun-exposed part of the crown, while the lower part of the crown has very little foliage [25,26,54]. As the measured crown density decreases with tree age, crown density may act as a proxy of tree maturity. Moreover, since mature trees produce more cones, a decrease in their number with increasing crown density is expected. On the other hand, a more physiological-based explanation can also be detailed. In denser crowns, the tree has to allocate resources (e.g., water, photosynthates, nutrients) to less vigorous and shaded branches in the inner and lower part of the crown, which commonly does not bear female conelets [18,26]. There are, thus, fewer resources left for reproductive organs, and the number of cones, thus, decreases.
The average weight of the cones is also inversely proportional to the density of the sun-exposed part of the crown. This result reflects the importance of light for reproductive shoots and branches bearing cones, where the resources assimilated are directly allocated to the enlargement of the cones. The same physiological argument can be put forward, whereby resources in denser crowns have to be allocated to the maintenance of non-reproductive shoots. Conversely, self-shading and light occlusion over the maturing cones, as they are located below two vegetative shoots, can result in losses in the photosynthetic activity of the cones in very dense crowns [55]. This would result in a reduction on the available photosynthates for cone development and enlargement.
These results shed light on the relevance of crown architecture on the whole process of female flowering and further cone enlargement in the species. The importance of crown dimensions and density has practical implications if cone production is to be maximized. Crown pruning (or poda de olivación in Spanish) consists of clipping the non-vigorous and dead branches of the inner part of the crown, in order to promote the emergence of female flowers in the following spring [50]. While widely practiced in Spain up to 30 years ago, this treatment is currently abandoned since no scientific evidence of significant increment in cone production justifies the high cost of the practice. Our results are in agreement with recent experiences [56] giving support to the idea of enhanced floral emergence in less dense crowns. The relationship between crown density, flowering, and cone maturation is nevertheless needed to establish the increase in cone production after pruning, after which the economic viability can be assessed.
Most of the TLS-derived metrics that are included in the models cannot easily be measured through traditional inventories. For example, crown projected area and minimum crown dimensions can be measured from the ground (presence/absence of cones). However, light and total crown volume/density (number of cones and average cone weight) cannot be estimated by survey crews. These variables improve model fit statistics and help better understand the important drivers that are not yet considered in understanding cone production.

Conclusions
When compared, the models based on TLS metrics clearly outperformed the inventory-based models, indicating that crown dimensions are an important driver in understanding cone production in stone pine. In this sense, our work points to the potential of using TLS as a tool for rapid, accurate, and unbiased measurements of crown parameters not easily accessible from the ground. Moreover, the links between crown dimensions and cone production should be further investigated to better understand the physiological processes involved, in order to both propose management practices to maximize production and better integrate the processes in prediction models. Nevertheless, managers should favor growing conditions that enable the trees to develop large, symmetrical crowns. Crown pruning should also be considered to increase cone production. Finally, the TLS attributes that were found to be important (e.g., crown volume, density, terrain roughness) can also be derived from other remote sensing sensors such as aerial or drone LiDAR. Quick, large-area estimates of cone production could ultimately be achieved and help managers identify areas of high cone productivity. Thus, such sensors might reduce the costs in estimating cone production.