Use of Multi-Temporal UAV-Derived Imagery for Estimating Individual Tree Growth in Pinus pinea Stands

: High spatial resolution imagery provided by unmanned aerial vehicles (UAVs) can yield accurate and efﬁcient estimation of tree dimensions and canopy structural variables at the local scale. We ﬂew a low-cost, lightweight UAV over an experimental Pinus pinea L. plantation (290 trees distributed over 16 ha with different fertirrigation treatments) to determine the tree positions and to estimate individual tree height ( h ), diameter ( d ), biomass ( wa ), as well as changes in these variables between 2015 and 2017. We used Structure from Motion (SfM) and 3D point cloud ﬁltering techniques to generate the canopy height model and object-based image analysis to delineate individual tree crowns (ITC). ITC results were validated using accurate ﬁeld measurements over a subsample of 50 trees. Comparison between SfM-derived and ﬁeld-measured h yielded an R 2 value of 0.96. Regressions using SfM-derived variables as explanatory variables described 79% and 86–87% of the variability in d and wa, respectively. The height and biomass growth estimates across the entire study area for the period 2015–2017 were 0.45 m ± 0.12 m and 198.7 ± 93.9 kg, respectively. Signiﬁcant differences ( t -test) in height and biomass were observed at the end of the study period. The ﬁndings indicate that the proposed method could be used to derive individual-tree variables and to detect spatio-temporal changes, highlighting the potential role of UAV-derived imagery as a forest management tool.


Introduction
Tree growth is an intermittent process characterized by temporal changes in stem form and dimension and is used to predict future forest conditions [1]. Changes in tree height and crown area are key to determining aboveground biomass growth for single trees, while at the stand level they also provide information about tree competition, productivity and forest health, amongst other variables. In this respect, forest managers require flexible, non-destructive methods for monitoring forest growth and estimating biomass.
Traditional field forest inventories are used for multi-temporal analysis based on tree height and crown growth measurements at different times. However, tree height growth can be demanding to estimate from the ground, because of the difficulty in determining the top of the crown [2,3], especially in species that do not display strong apical dominance. Crown growth is also difficult to determine when crowns overlap, as the edges of the crowns are often hidden by the branches of other trees [4]. Crown area is also difficult to estimate from field data; two cross crown diameters are usually measured and crowns are assumed to be ellipses, to simplify the calculations. Although field-measured data are commonly assumed to be ground truth values for remote sensing data, the associated measurement errors tend to be large [2,[4][5][6]. In addition, forest mensuration with traditional field techniques is expensive and time consuming, and measurements may take several years.
The suitability of remote sensing-based technologies for forest mensuration is widely acknowledged [2,7]. Since the 1990s, digital aerial photography (DAP), airbone laser scanning (ALS) (e.g., [8][9][10] or a combination of both [11][12][13][14][15] have been increasingly used to support forest inventories at different scales. Interest in using UAVs to acquire both ALS and DAP data has increased greatly in recent years [16][17][18]. Indeed, ALS and red, green, blue (RGB) sensors mounted on UAV platforms are becoming cost-effective tools for monitoring forest structure because the high spatial (achieved by the low flight height) and temporal resolution (allowed by the relatively low cost of the flight surveys) meet the requirements of forest managers [17,19]. Multi-temporal UAV-acquired data can be used to help assess tree growth quickly, accurately and economically, providing up-to-date information to support decision-making in forest management [20][21][22]. Due to the low cost and operational flexibility, light UAVs equipped with inexpensive consumer grade cameras have recently emerged as a feasible option for monitoring three dimensional (3D) forest structures [19]. In this case, an innovative computer vision technique-Structure from Motion (SfM)-enables extraction of 3D information from UAV flights, providing 3D photogrammetric point clouds based on feature matches within overlapping images [4,11,12,14,23,24].
Two main strategies have been adopted for DAP and ALS-based analysis: the area-based approach (ABA, a distribution-based technique, which typically provides data at stand level) and individual tree crown (ITC) delineation, in which individual tree crowns, heights and positions are the basic units of assessment). ABA has been used with ALS and DAP to estimate forest attributes over a wide range of forest types including Temperate (e.g., [25]), Boreal (e.g., [9,10,15,[26][27][28][29]), Atlantic (e.g., [30,31]), Tropical (e.g., [32]), Alpine (e.g., [33,34]) and Mediterranean forests (e.g., [35][36][37]. At the stand level, results from recent research on small-to medium-sized boreal and tropical forests have demonstrated the potential use of UAV-derived data for estimating forest biomass [19,38]. On the other hand, ITC has also been applied to DAP point clouds [39] and to ALS clouds ( [40][41][42]; see [43] pp. 166-169, for a review of the ITC approach to estimate aboveground biomass using ALS points clouds). ITC presents several advantages over ABA for estimation of above-ground biomass because it can be used to derive biomass when an allometric model is available at individual tree level [43]. Moreover, fewer reference data are required in the ITC approach, because a greater number of physical parameters are measured directly, while the final result of the ABA depends on good calibration with extensive, accurate, representative and expensive field-acquired data [44]. Finally, in many forests, stand-level results are usually not appropriate for forest management planning, especially for harvest management, pruning and forest fruit production, or in valuable stands in which precision forestry, which usually requires information about individual trees, is applied.
However, the ITC approach is traditionally very expensive and not widely used at the operational level, because it requires spatially dense ALS data or very high resolution DAP imagery. This problem has been partly solved by advances in SfM techniques that have enabled generation of high-density point clouds derived from top-of canopy by UAV-based DAP. Compared to ALS, UAV-based DAP is inexpensive and the point cloud can easily match ALS densities. However, the larger point density does not necessarily yield greater vertical accuracy, because DAP cannot penetrate vegetation. Therefore, for UAV-based DAP, the requirement for a precise, high-spatial-resolution Digital Elevation Model (DEM) to produce an accurate normalized canopy height model (CHM) must also be considered, especially in densely vegetated areas.
Despite the wealth of literature demonstrating that ALS (e.g., [45][46][47][48]), DAP (e.g., [28,49]) and UAV-SfM [50][51][52] can be used for ITC as an effective means of detecting change at fine scales, to the best of our knowledge, few studies have evaluated the usefulness of CHMs, created with SfM, to assess forest growth biomass from multiple UAV datasets. In this regard, recent advances in SfM may provide a means of constructing multi-temporal CHMs at a very high spatio-temporal scale, for accurate estimation of individual-tree growth.
The PINEA project (www.pineaproject.com) is focused on modelling growth for Pinus pinea under changing environmental conditions. The project established irrigation and fertilization experiments in Portugal in 2013. This study explores the performance of high spatial resolution UAV-derived imagery to monitor forest growth and biomass under different fertirrigation treatments. It uses an ITC approach, based on object-based image analysis (OBIA) techniques, to examine the temporal changes in individual-tree variables in P. pinea plantations.
The objectives of this study were as follows: (i) to investigate the combined use of SfM-derived individual-tree measurements (height and crown area) and nonlinear regression models to estimate individual tree diameter and aboveground biomass for the period 2015-2017 and (ii) to assess the use of SfM-derived individual-tree measurements to estimate spatio-temporal changes in biomass growth in relation to different fertirrigation treatments.

Study Area
This study was conducted in Esteveira, a privately owned forest close to Alcochete in Central Portugal (latitude 38.729789 • , longitude −8.834107 • ). The study site is an experimental umbrella pine (P. pinea) forest plantation covering an area of 16 ha. Trees were planted in rows of 8 × 10 m in 1992-1993. In order to form an open canopy of scattered trees and to homogenize the initial site conditions, systematic thinning was conducted in 2013, yielding a final density of approximately 63 stems ha −1 distributed in a 10 × 16 m grid. The study site is characterized by fairly flat terrain (slopes from 0 to 6%, elevations from 26 to 49 m above WGS84 ellipsoid) and no understorey.
During the spring of 2014, six adjacent rectangular plots were established in the field and a complete randomized block design including 2 blocks and 3 treatments was implemented: control (T0) and two different levels of fertirrigation (T1 and T2) (see Figure 1). A total of 290 trees were included within the study area (96 for T0, 94 for T1 and 100 for T2).
The fertirrigation treatments were carried out between mid-May 2015 and mid-September 2016. A fertilizer solution was injected directly into a drip irrigation system with the aid of a hydraulic driven injector in early May, July and August. After analysis of soil and foliar nutrient conditions, two levels of fertirrigation were selected: doses of 76 kg N ha −1 year −1 , 15 kg P ha −1 year -1 and 15 kg K ha −1 year −1 for treatment T1 and double these amounts for treatment T2 (152 kg N ha −1 year −1 , 30 kg P ha −1 year −1 and 30 kg K ha −1 year −1 ). The mean annual rainfall at the study site was 406.7 mm in 2015 and 580.5 mm in 2016 (water year started on 1 October and ended on 30 September). In total, 139-142 mm of additional water was supplied in T1 and 278-284 mm in T2 during the 2015-2016 and 2016-2017 growing seasons. The water was supplied at a constant daily rate of 0.92-0.94 and 1.85-1.90 mm in T1 and T2, respectively, to ensure that tree growth was not constrained by soil water deficit.

Field Measurements and Field Biomass Estimation
Individual-tree aboveground biomass (wa) was estimated using field measurements of diameter at breast height (d) and tree height (h) from a subsample of 50 trees from the centre line of one plot per treatment in 2015 and 2017 (Table 1). In 2015, h was measured to the nearest decimeter using a Vertex III hypsometer (Haglöf, www.haglof.se). In 2017, a telescopic pole was used to improve the height measurement accuracy, and h was recorded to the nearest centimetre. The change in method was justified after noting that the hypsometer tended to overestimate height [53], which was attributed to the lack of strong apical dominance and the umbrella-shaped crown of this species. In 2015, d was measured using tree callipers: two measurements were made at right angles to each other and averaged. At the initial measuring time, a mark was painted on the stem of each selected tree at breast height (1.3 m above ground level). In 2017, d was measured using a Pressler increment borer (Djos 400, diameter: 6 mm, height: 33 cm). One core sample was obtained at the paint mark on each tree. Cores were air-dried, stuck with glue onto wood guides and sanded using progressively finer grain papers, until the rings were clearly visible. Finally, tree ring widths were measured with a stereomicroscope (LEICA M80) and a LINTABTM table (Rinntech ® ) connected to a TSAP-WinTM 4.64 tree-ring analysis system (Rinntech ® ).
In both 2015 and 2017, wa was estimated as the sum of the biomass of the tree components (stem, bark, needles and branches) (Equation (1)).
where Y is the dependent variable (wa, kg), X 1 to X 2 are the independent variables (tree stem circumference-c (cm)-and h (m), respectively), β 0 to β 3 are the multiplicative parameters, λ 1 to λ 5 , are the exponential parameters and ε is the additive error term. For this study, a new biomass equation was developed for P. pinea in Portugal, improving the fit of the existing equation described in [54], by adding data from the destructive sampling of 35 newly harvested trees (outside of the study area), for a total of 75 trees (for more details on procedures and the field data used, see [54]). Model parameters were estimated using non-linear seemingly unrelated regression (SUR). Finally, for the 50 trees sampled per treatment, we predicted wa in by using h and c as inputs in Equation (1). We decided to use the height from the image as a ground truth value to predict wa in 2015 because of the aforementioned field errors in h measurements in 2015 and following the methodology of recently studies [50]. h is the field-measured individual tree height (m), d is the field-measured diameter at breast height (cm), wa is the field-derived individual above-ground biomass (kg). # In 2015, field-measured tree height (h) was replaced by SfM-derived tree height (h SfM ).

Remote Sensing Data Acquisition and Use
The airborne surveys were conducted by Terradrone Co. and 202 images (year 2017) were used to generate orthomosaics and digital surface models (DSMs) by the SfM image reconstruction process. The aerial data acquisition required greater control as it is important that the flights are always carried out under conditions of equal luminosity, a high degree of overlap and an covered area much larger than the selected area to avoid edge distortion [23,55].

3D Model Generation
The absolute orientation of the aerial photos was determined using aerotriangulation techniques performed by pix4D 3.1.22 (pix4D ® , Ecublens, Switzerland). Five (2015) and ten (2017) ground control points (GCPs) were measured in the field with topographic methods in order to georeference the SfM mosaic to a projected coordinate system. A LEICA GS15 GPS (Leica Geosystems, San Galo, Switzerland) (dual frequency, real time kinematic receiver with a planimetric precision of ±8 mm + 0.5 ppm and an altimetric precision of ±15 mm + 0.5 ppm) was used to capture the ground control photogrammetric target (marked with 60 × 60 cm matte white finish plastic boards detectable in RGB photographs). For reliable accuracy of GPS measurement, all GCPs were located in open areas with no canopy cover. At each point, GPS signals were logged in real time kinematic (RTK)-global navigation satellite system (GNSS) mode. The recordings were processed with real-time correction data retrieved from the fixed base station in Palmela (Lisbon) (latitude 38.5714631 • , longitude −8.9034504 • , and ellipsoidal elevation 245.989 m above the WGS84 reference ellipsoid). In 2015, GCPs were marked with surveying wooden stakes in order to relocate the centre of each point.
The photogrammetric point clouds were generated using the stereomatching algorithm implemented in pix4D 3.1.22. The matching parameters were set as follows: multiscale, image scale = 1/2 (half image size) and point density = 'optimal'. We also set the minimum number of matched images to 3. Postflight Terra 3D project from 2015 was migrated to pix4D 3.1.22 in order to improve the accuracy of the digital elevation models (DEMs). Non-ground points were removed using the automatic object classification implemented in Pix4D. DEMs were generated from the ground points by using a natural neighbour interpolation technique implemented in Pix4D (additional details of the algorithms are proprietary and were not disclosed by Pix4D). Finally, one CHM was constructed for 2015 and another for 2017 by subtracting the DEMs from the DSMs point cloud using the FUSION LiDAR Toolkit [56].

ITC Process and SfM-Derived Variables
Individual tree position (X and Y coordinates), height (h SfM ), and crown area (ca SfM ) were retrieved from the CHM (Figure 2). Resampling of the CHM to 20 cm resolution and subsequent smoothing with mean (5 × 5 window) and median filters (5 × 5 window) were conducted using the FUSION LiDAR Toolkit [56]. Crown delineation followed the procedure detailed in [42], but because of the forest stand and tree species characteristics, the initial maximum search domain in the iterative process (see Figure 3 in [42]) was changed from 5 to 25, and only one interaction was applied. Crown delineation results ( Figure 3) and tree top positions were then exported in ESRI™ shapefiles as vector polygons and points, respectively, for subsequent analysis.

Individual Tree Biomass Estimation
Most biomass equations, including Equation (1), require measurement of tree diameter or circumference, which is not available from UAV imagery. We therefore tested two approaches for estimating SfM-derived individual-tree biomass (wa SfM ) ( Figure 4). In the first approach, the multiplicative (power function) model in Equation (2), was fitted by using d as the dependent variable and h SfM and ca SfM as explanatory variables. Finally, the predicted diameter and h SfM were included as explanatory variables in Equation (1) to predict wa SfM for the subset of 50 trees. In the second approach, the multiplicative (power function) model in Equation (2) was also fitted to predict wa SfM for the 50 trees, but wa (calculated using the field-measured d and h in Equation (1)) was used as a dependent variable, and h SfM and ca SfM were the explanatory variables.
where Y is the dependent variable, X 1 to X 2 are independent variables β 0 is the intercept, λ 1 and λ 2 , are the exponential parameters to be estimated by non-linear regression analysis and ε is the additive term of the error. The models were fitted using the nls function implemented in the BASE package of R software (R Core Team, Vienna, Austria, 2016).

Accuraccy Analysis of the ITC and Individual Tree Variables
The ITC computations were validated in order to detect possible miscalculations in the image processing phase. We searched for false positives and/or false negatives using the method described by [57] and modified by [42] to link the SfM-detected trees to their corresponding field-measured trees. In order to assess the accuracy of SfM-detected tree heights, we plotted h vs h SfM by using the correctly detected and linked trees and then visually inspected the graphs. We also compared d with d SfM and wa with wa SfM in the subsample of 50 trees for both 2015 and 2017. Finally, we computed the coefficient of determination (R 2 ), the overall root mean square error (RMSE) and the relative root mean square error (rRMSE) to determine the accuracy of SfM for estimating tree heights, diameter and aboveground biomass in both 2015 and 2017.

Growth Analysis
A paired t-test was also conducted to compare SfM-predicted forest growth variables (tree height growth, ∆h SfM , and crown area growth, ∆ca SfM ) for the entire sample tree data set (N = 290) over one growth period (2015-2017) to determine whether the interval was sufficient to enable detection of statistically significant growth measurements in individual trees. The normality of the residuals of SfM-predicted forest growth variables was tested using the Shapiro-Wilk Test [58].
One-way ANOVA was used to analyze the effect of the fertirrigation treatments on the mean ∆h SfM and ∆ca SfM and to assess the validity of the UAV-based DAP method to detect differences in forest growth related to the treatments. We included a random plot effect, to account for the lack of independence of trees measured within a plot, and a random block effect. The experimental data were analysed using SAS 9.4 (SAS Institute, Cary, NC, USA, 2012), with the proc mixed procedure and restricted maximum likelihood and the Kenward-Roger correction for degrees of freedom. We also used Tukey's test for multiple pairwise comparisons to identify any significant differences between the treatments. Table 2 Table 2. In 2017, 5 independent GCPs were also used as check points to validate the accuracy of the geo-referenced point clouds (check points).

Field and SfM Biomass Estimation
The parameter estimates and goodness-of-fit statistics of the model used to predict wa (kg) as the sum of the biomass components (stem, bark, needles and branches), with c (cm) and h (m) as explanatory variables (Equation (1)) are shown in Table 3. w w is the individual tree stem wood biomass (kg), w b is the individual tree stem bark biomass (kg), w l is the individual tree needle biomass (kg), w br is the individual tree branch biomass (kg), w a is the individual tree aboveground biomass (kg), d is the diameter at breast height (cm), c is the stem circumference at breast height (π*d/100) (m), h is the individual tree height (m), Mef is the model efficiency statistic and RMSE is the root mean square error.
The model represents an improvement in terms of Mef relative to previous models used to predict individual tree aboveground biomass in P. pinea plantations in Portugal [54] and it is also more representative, as the sample size is larger, including 35 additional data points. The RMSE was 163.7 kg, with a model efficiency (Mef) of 0.96 (Equation (1), Table 3). The improved model was applied to estimate wa for the subsample of 50 trees. Table 4 shows the parameter estimates and goodness-of-fit statistics for the models used to predict d (cm) and wa SfM (kg) in 2015 and 2017 by approach 1 and the parameter estimates and goodness-of-fit statistics of the model used to predict wa SfM (kg) in 2015 and 2017 by approach 2 (see the description of both approaches in Section 2.6. Individual tree biomass estimation).

Accuracy Analysis of the ITC and Individual Tree Variables
We tested the accuracy of the ITC delineation by comparing and plotting h and h SfM . The linear regressions between these two variables yielded an R 2 value of 0.96 and an RMSE value of 0.19 cm with an rRMSE of 1.89% for trees with field-measured height ranging between 7.00 and 12.72 m (Figure 5a). We also checked the overall detection rate, which was 100%, with no false positives or negatives. UAV-based DAP methods tended to underestimate tree height, relative to field measurements (vertex and telescopic pole), for both 2015 and 2017. Although d is not directly measured in CHMs derived from UAV, SfM-derived variables, especially ca SfM , were strongly correlated with d ( Figure 5b). The diameter predicted by Equation (2) was regressed against d, explaining 79% of the variability for both years (Figure 4b). There was no appreciable bias throughout the observed diameter range. The linear regression model that predicted wa SfM by following approach 1 explained 87% (2015) and 85% (2017) of the variability (Figure 5c). However, this approach seemed to underestimate field above-ground biomass for 2017, especially for the largest trees (Figure 5c). Finally, the wa SfM predicted by approach 2 explained 86% and 85% of the variability in 2015 and 2017, respectively, and did not show any appreciable bias (Figure 5d).

Growth Analysis
There was a statistically significant difference between SfM-predicted forest growth variables (∆h SfM , ∆ca SfM and ∆wa SfM ) for the individual tree sample (N = 290) after two growing seasons (2015-2017). The mean annual ∆h SfM was 0.45 m, with a 95% CI between 0.43 and 0.47 m. The mean annual ∆ca SfM was 17.29 m 2 , with a 95% CI between 16.17 m 2 and 18.40 m 2 (Table 5). Finally, the average aboveground biomass production for the 2-year period at the individual tree level was 189.7 kg, with a 95% CI between 187.57 kg and 209.26 kg (estimate is derived from the second approach, Equation (2), Table 4).   Height growth was detected for all 290 trees, but with a reduction in crown area in 5 trees (Table 5, Figure 6a,b). We observed small differences in the mean ∆h SfM between the treatments, with an average difference of 7 cm between the control and the combined fertirrigation treatments. For ∆h SfM , one extreme positive outlier, corresponding to a small tree for treatment T1 in plot 2 (probably due to an error in obtaining the single tree height in 2015) was removed from the analysis. The differences in mean ∆wa SfM and particularly in ∆ca SfM were more evident (Figure 6b,c, respectively). The negative values in ∆ca SfM were checked in the field and were explained by a variety of causes, including a dead tree, strong competition with neighbouring trees, and trees with large fallen branches. In terms of stand level biomass, the differences amounted to 3.86 Mg ha −1 year −1 , 4.45 Mg ha −1 year −1 and 5.62 Mg ha −1 year −1 for T0, T1 and T2, respectively.  There was no evidence of any effect of the fertirrigation treatments on ∆h SfM and ∆wa SfM (p-values of 0.27 and 0.21, respectively, both based on an F-test on 2 and 2 degrees of freedom). Although tree growth was greatest in the fertilized and irrigated plots, the small sample size (only 6 plots and 2 degrees of freedom for the error) yielded a low power and a lack of ability to detect differences between treatments. However, we found a significant difference in ∆ca SfM (p = 0.04). Tukey-Kramer's test for multiple pairwise comparisons revealed that T1 and T2 were statistically significantly different from the control (T0), but not from each other.

Discussion
The spatial accuracy of the 3D cloud data is important to ensure a good description of tree crowns and was therefore a key issue in this study. Given the possible use of UAVs for monitoring vegetation dynamics, accurate direct georeferencing is an important processing step in forest areas [18,59,60].
Few studies have examined the effect of the sample size of GCPs on the accuracy of photogrammetric projections from UAVs [60,61]. We observed very small horizontal and vertical differences in accuracy with 5 or 10 GCPs (Table 1), as reported in other studies in flat terrain [61]. The mean RMSE X , RMSE Y , and RMSE Z values for 10 CGPs were 0.026, 0.007 and 0.020 m, respectively. These values are slightly better than those obtained using a rotatory-wing UAV with 10 CGPs and a flight altitude of 120 m [61], but slightly worse, in terms of horizontal and vertical accuracy (0.007 m, 0.005 cm, and 0.014 cm, respectively), than those reported for a fixed-wing UAV in a boreal conifer forest, with 13 GCPs, a flight altitude of 120 m, and an average slope of 10 • [19]. Our study findings confirm that horizontal and vertical accuracy increases with the number of GCPs. Furthermore, the RMSEz values increased with flight altitude, in agreement with the vertical RMSEz trend observed in a previous study [61].
The findings of the present study showed that individual tree variables (h SfM , d SfM and wa SfM ) could be estimated automatically from high resolution imagery, suggesting promising research lines for fast, reliable and cost-effective annual monitoring of P. pinea plantations that could help assess the efficiency of silvicultural treatments. The h SfM estimation was similar to or better than that of previous studies uding SfM point clouds in palm plantations [62], olive orchards [4,24], boreal forests [32] and Eucalyptus forest [18]. The proposed methodology could be helpful in P. pinea forest management, as estimation of annual tree height growth from field data (∆h) in this species is extremely difficult. Indeed, we were not able to validate the ∆h, because of the large measurement errors observed using hypsometers in 2015, reaching approximately 0.5 m [53]. In 2017, although we significantly improved the method of measuring h in the field in order to validate h SfM , the CHM still slightly underestimated h. These observations are consistent with those of other studies working with DAP [63][64][65], DAP-UAV [20,66,67] and ALS [46,47,68,69] point cloud data. The following factors may partly explain the differences between h and h SfM : (i) inaccuracies or bias in field-measured h caused by the difficulty in determining the top of the P. pinea trees due to their umbrella-shaped crown and lack of clear apical dominance; (ii) The developed-DEM from unsupervised ground-classified points filtering implemented in Pix4D 3.1.22 performed well in this flat area; however, the poor performance of the SfM photogrammetric technique for capturing the terrain surface below closed canopy cover, possibly leading to overestimation of the ground elevation [20,66,67] and underestimation of h SfM ; and (iii) smoothing of the image processing during the ITC process (necessary for noise reduction), possibly resulting in slight underestimation of h SfM .
Field measurements are traditionally considered ground-truth values and are used as reference values for comparisons with remote-sensing estimates. However, some field measurements, such as h and ca, have unavoidable associated error and large variation [2,[4][5][6] and thus may not be suitable for validation of the SfM-derived tree growth measurements, especially in dense stands or for trees which do not have well-defined top [70]. In our study, h SfM was less variable than field-based measurements in both 2015 and 2017. Field h was measured to the nearest 0.1 m, while h SfM was measured from tree apexes detected to the nearest 0.01 m. Field-measured h standard deviations were 1.05 m and 1.02 m in 2015 and 2017, respectively, while h SfM standard deviations ranged from 0.92 m to 0.98 m, in 2015 and 2017, respectively. Future research should consider the combined use of UAV and LiDAR to improve the DEM/CHM generation, as well as the combination of GPS and total station measurements to obtain ground-truth values of the forest variables and ground point under the canopy for improved quantification of the errors. A recent study has already demonstrated improvements of up to 40 cm in tree height measurement using LiDAR/UAV instead of UAV imagery in a dry sclerophyll eucalypt forest within areas of relatively low canopy closure [18]. Another solution might be to consider 1 or 2 points per square meter under canopy using a lightweight, cost effective rangefinder with the ability to log and co-register image-based point clouds.
The wa SfM was accurately estimated from UAV-based DAP using both approaches (approaches 1 and 2, described in Section 2.6. Individual tree biomass estimation). Individual tree variables derived from the automated processing of UAV imagery and ITC delineation, such as h SfM and ca SfM , were found to be significant explanatory variables for predicting d and wa. Studies conducted in Picea abeis (L.) H. Karst. and Pinus sylvestris L. stands in Sweden and in Pinus taeda L. stands in the SE United States found that ALS-derived h and crown diameter (cd) explained up to 87% and 91% of the variance associated with the estimation of d, with an RMSE of 3.8 and 4.9 cm, respectively [41,57]. Zhao et al. [71] reported an adj. R 2 value of 0.87 and an RMSE value of 5.2 cm for ALS-derived tree dimension variables including h, cd and crown base height in P. taeda stands. In this study, the diameter equation based on SfM-derived variables performed well, although the values of Mef are slightly lower than some of those reported for other species [41,57]. The performance of the wa SfM estimates for predicting tree biomass directly from SfM-derived variables (approach 2, adj. R 2 = 0.86-0.85, RMSE = 87.46-117.8 kg) was similar to that obtained by [41] (R 2 = 0.88, RMSE = 169 kg) and slightly better than those for P. taeda (R 2 = 0.80, RMSE = 237 kg) [71]. For wa SfM , the potential sources of error and variation may be similar to those already discussed for h SfM and ca SfM . Specifically, the tendency of SfM to underestimate h may be the main reason for underestimation of wa SfM with approach 1 in 2017. The allometric relationship between wa SfM and SfM crown-derived variables could be refined, either by harvesting more trees or through improvements in UAV imagery processing.
The results of the growth analysis demonstrated that the resolution of the photogrammetric point clouds obtained by stereomatching aerial photographs would be sufficient for estimating ∆h SfM , ∆ca SfM and ∆wa SfM in P. pinea plantations, because low-altitude flights can be carried out with UAVs. The relatively short period between data collection (2015-2017), and the relatively slow tree height growth rate of this species at this age, may also have led to some inaccuracy as measurement and precision errors may be within the physical limits of tree growth. The mean annual ∆h SfM computed using 96 trees within T0 treatment was 0.20 m year −1 , lower than the published site index 21 growth rate of 0.28 m year −1 derived from the forest growth model for this species [72]. The mean annual ∆wa SfM from high resolution images showed that, during the experimental period, fertirrigation increased the biomass growth rate of P. pinea plantations, relative to the unfertilized treatment (T0), but the small sample size led to a low level of statistical power for detecting differences after 2 years of treatment ( Figure 5). It was possible to estimate the annual biomass increment at tree level because the h SfM was considered as a ground-truth value in 2015, following the methodology of recent published studies [50]. Goodbody et al. [50] estimated the mean gross tree volume increments considering h and ca directly from CHMs derived from light UAV-based DAP and LiDAR, respectively. Use of the PINEA.pt growth model [73] implemented in the forest growth simulator (StandSIM module) [74,75] with an input tree file for the existing stand structure (Stem density = 54 tree ha −1 , Stand age = 24 and Site index = 21) in the study area showed that the mean annual aboveground biomass production (3.9 Mg ha −1 year −1 ) was similar to the estimated biomass growth in the control (3.86 Mg ha −1 year −1 ). Estimates were within the limits of forest growth models [74,75] and local independent field measured sample plot. UAV-DAP point clouds were therefore effective for updating individual tree growth and biomass between 2015 and 2017. The results also show a hierarchical growth habit in this species due to a lower level of apical dominance of the leader shoot than in most other conifers [76].
It was also possible to detect differences in biomass yield under different treatments. Our findings are consistent with those of previous studies, which investigated the effects of fertilization on the biomass allocation of hardwood species [77][78][79]. Albaugh et al. [79] used data from nine sites in which nutrient and water optimization studies were carried out with a 2 × 2 factorial design to determine maximum biomass production in response to a simple set of treatments. These authors found that above-ground, stem and total biomass production ranged from 3 to 30, 2 to 28, and 5 to 32 Mg ha −1 year −1 , respectively. In the present study, above-ground production (∆wa SfM = 5.62 Mg ha −1 year −1 ) was within the range reported for fertirrigation treatment in other species, although at different ages and considering only the effect of two years of treatment. For example, the values of wa SfM were lower than for 4-year-old P. taeda receiving fertilization and water treatments [80] (6.9 Mg ha −1 year −1 ) after three growing seasons in the southeastern United States and those reported by [81] (10 Mg ha −1 year −1 ) and up to 18 Mg ha −1 year −1 [82] in a 2-year-old stand and 8-year-old stand of P. taeda after the 6th and 4th year of treatment, respectively. The fertirrigation treatment (T2) increased above-ground biomass increment, with growth responses reaching up to 51% of that the control. This value is similar to that reported for a P. pinaster plantation in south-western France (5-year-old stand) after five years of treatment [83]. These differences may reflect differences in initial site fertility and climate, as well as the responsiveness of the species to nutrient treatments.
Finally, inferences from multi-temporal images could also be improved by carrying out the flights on almost the same date, under similarly good atmospheric conditions, and using similar image-acquisition parameters. In this study, this methodology was applied on flat terrain, in open canopies with no understory. In natural and managed closed-canopy forests and steeper terrain, poorer performance could be expected, as aerial photography cannot penetrate vegetation, and additional ALS data may therefore be needed. In order to ameliorate such errors in future studies, an approach using broad scale ALS acquisition at 10 year intervals, with regular 4-6 year UAS-DAP point cloud updates for mature stands, could be applied [50,84].

Conclusions
This study shows that multi-temporal UAV datasets and SfM techniques are promising tools for monitoring individual tree variables, forest and biomass growth and, ultimately, supporting decisions regarding the management of P. pinea forests. Future research with multi-temporal UAV datasets could benefit from the knowledge provided in this study, which reports the monitoring of spatial and temporal processes, analysis of individual tree-level competition, growth and effects of the treatments and forest management practices in P. pinea stands. The proposed methodology could help in the management of P. pinea forest and woodland-pasture ecosystems, as estimation of tree height growth from field data in this species is extremely difficult due to the lack of strong apical dominance and an umbrella-shaped crown. The study findings represent progress in the use of high resolution UAV imagery for estimate forest tree variables and forest growth, improving the temporal resolution of the forest inventories, reducing their cost and in some cases improving the accuracy of the forest variable estimates. In summary, this study has shown the potential of multi-temporal UAV data to characterize forest growth and to estimate biomass in Mediterranean P. pinea plantations.