Use of LiDAR data during multi-annual periods for estimating forestry variables

Aim of study: To test the use of LiDAR data from a single acquisition in order to estimate volume over_bark variations in a 5-yr period of Pinus radiata D. Don. Area of study: Province of Bizkaia in the Autonomous Community of the Basque Country (Spain). Material and methods: Two field plot measurements were made in 2011 and 2015 and two wood volume models (one for each year) were fitted using the metric variables of the 2012 LiDAR points cloud. The models were applied to a 26.59 m raster covering the study area and the increase in volume at each pixel was calculated by subtraction. Main results: The increase in estimated wood volume, when added to the volume of timber extracted in the area during the 5-yr period under consideration, yielded an average increase of 13.74 m3 ha-1 yr-1, which corresponds to the average growth of the P. radiata in that area. The harvest area estimated using this procedure largely coincides with the actual harvest area in the same period. The value of R2 (85%) of the wood volume model for 2011 is similar to that obtained in other studies. However, as expected, the one obtained for the wood volume model for 2015 (80%) is significantly lower. Research highlights: The increase in wood volume can be estimated using a single LiDAR flight and field data from the 5-yr period provided that data from plots subjected to this kind of harvest is included in the models. Additional keywords: LiDAR forest inventory; wood volume increase; wood volume LiDAR model. Abbreviations used: CCF (Canopy Cover Fraction); CV (Coefficient of Variation); Dm (Mean Diameter); G (basal area, m2 ha-1); Hm (mean height, m); Ho (top height, m); IFN4 (4th National Forest Inventory); LiDAR (Light Detection and Ranging); MAE (Mean Absolute Error); ME (Mean Error); RMSE (Root Mean Square Error); Vob (stem volume over bark, m3 ha-1). Authors ́ contributions: MAV conceived and designed the experiments and contributed to writing, data production, and analysis. EM wrote the introduction and reviewed the work. FR revised the statistical calculations and fitted models. Citation: Valbuena, M. A.; Mateos, E.; Rodríguez, F. (2017). Use of LiDAR data during multi-annual periods for estimating forestry variables. Forest Systems, Volume 26, Issue 3, e019. https://doi.org/10.5424/fs/2017263-11468 Received: 28 Mar 2017. Accepted: 10 Jan 2018. Copyright © 2017 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution (CC-by) Spain 3.0 License. Funding: Gobierno Vasco (Projects SAI10/147SPE10UN90); Vicerrectorado de Investigación de la Universidad del País VascoEuskal Herriko Unibertsitatea (UPV/EHU) (Project NUPV14/11). Competing interests: The authors have declared that no competing interests exist. Correspondence should be addressed to Manuel A. Valbuena: mavalbuena66@gmail.com; valbuena@irakasle.eus Forest Systems 26 (3), e019, 13 pages (2017) eISSN: 2171-9845 https://doi.org/10.5424/fs/2017263-11468 Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria O. A., M. P. (INIA)


Introduction
Over recent years, there has been increasing interest in the evaluation of the forestry variables, such as biomass and timber volume, especially within the context of the Kioto protocol (Lindner & Karjalainen, 2007;Zhang et al., 2007). The quantification of these forestry parameters reveals the structure, operation and dynamics of forestry ecosystems and their assessment as an energy resource. For this reason, it is essential to determine these parameters in forestry management (Torres et al., 2010;Pan et al., 2011). In most European countries, the increase in these variables is greater than the amount harvested on an annual basis, converting forests into major carbon and biodiversity reserves (Mateos et al., 2016). One of the key indicators for assessing forestry conditions is the estimation of forestry parameters, such as timber volume, but the usual methods for estimating this parameter require prolonged fieldwork (Song et al., 2016). A large part of the forestry bibliography relating to this subject focuses on adjustment of timber volume equations by means of regression for geographical regions and specific tree species (Nilsson, 1996;Hyyppä et al., 2001;Popescu et al., 2003;Holmgren, 2004;Zianis et al., 2005;Dalponte et al., 2009;Latifi et al., 2010;Magnussen et al., 2010). The high cost of sampling, especially when this is done in inaccessible areas, limit in many cases the establishment of sufficient plots to determine real variability (Popescu, 2007).
The link between the parameters obtained in forestry inventories and LiDAR (Light Detection And Ranging) data through regression models offers an attractive approach to estimating forestry variables on a regional, continental or global scale (Li & Weiskittel, 2010;Naesset, 2011;Babcock et al., 2015). Over recent decades, LiDAR technology has become an extremely powerful tool for estimating the biophysical properties of the vegetation. This is an active, non-destructive remote sensing method that offers a better alternative for the study of the structural variables of the canopy due to the penetrating capacity of the laser pulse through gaps in the vegetation (Leiterer et al., 2015;Giannico et al., 2016;Koenig & Höfle, 2016). The recording of multiple returns allows us to obtain sufficient three-dimensional, morphological information at a broad range of space and time scales of forestry variables such as the height of trees (Dean et al., 2006;Popescu & Zhao, 2008), the basal area and the diameter of the crown (Giannico et al., 2016), volume and aboveground biomass (Popescu, 2007;Jochem et al., 2010;Zolkos et al., 2013), and the density and classification of groups of forest species (Brandtberg, 2007;Zhang et al., 2016). Moreover, the use of LiDAR data offers the possibility of reducing the number of required sampling plots (Maltamo et al., 2007;Montaghi et al., 2013) and the accuracy obtained with LiDAR is greater than that reached with traditional inventory techniques (Condes et al., 2013;Watt et al., 2013;Maltamo et al., 2014).
Now that LiDAR technology has demonstrated its efficiency and accuracy in estimating most forestry variables, today the challenge is to try to take advantage of this tool in order to predict the annual growth of the plant structure at a large scale. The ease with which data can be updated allows us to readjust the estimations and quickly assess changes or disturbances in forest stands by estimating their development over time, or establish comparisons between different areas with a simple update of the forestry cartography (Riaño et al., 2004). The bibliography contains numerous studies that used bi-temporal data from two LiDAR flights to estimate plant growth (Yu et al., 2004(Yu et al., , 2006Naesset & Gobakken, 2005;Hopkinson et al., 2008;Vepakomma et al., 2011;Fekety et al., 2014;Cao et al., 2016;Song et al., 2016). Some researchers have estimated the growth of trees in forests on the basis of the differences in height derived from LiDAR data in flights separated by a 2-yr interval (Naesset & Gobakken, 2005;Yu et al., 2006;Vepakomma et al., 2011;Cao et al., 2016). These studies allowed the growth in different forestry parameters to be calculated, such as the basal area or volume (Naesset & Gobakken, 2005) as well as estimations of the vertical and lateral growth of trees (Vepakomma et al., 2011). Other researchers have used multi-temporal LiDAR datasets to estimate the growth of forests (Yu et al., 2006;Hopkinson et al., 2008;Fekety et al., 2014;Cao et al., 2016;Song et al., 2016). However, we have not found any previous studies that use the data from a single LiDAR flight to estimate the increase in timber volume over a period of several years by fitting models that relate the timber volume at different moments in time with the metrics of the heights of a single LiDAR points cloud obtained at the beginning of the period. In view of the high cost of the flight, the optimisation of the long-term use of LiDAR data gathered during a campaign may represent a considerable economic advantage.

Material and methods
This study assesses the application of data from the last LiDAR flight made over the Autonomous Community of the Basque Country in 2012 to evaluate the growth of volume over bark (Vob) during 5 years on stands of Pinus radiata D. Don (PR) through the comparison of estimated volumes by applying fitted models to relate separately the field data obtained with 5 years of difference with the metrics of the LiDAR points cloud at the beginning of the period.
The area based approach (ABA) for forest inventory applications using LiDAR has been utilised in this study, following the methodology applied by other researchers (Hopkinson et al., 2008;Jochem et al., 2010;Yu et al., 2010). In order to get a better understanding of how LiDAR and the set of field data can be used to map Vob, we used two sets of field data gathered in 2011 and 2015 respectively. This involves analysing and validating the proposed Vob models, fitting the corresponding models for both sets of data in accordance with the metric variables of the LiDAR points cloud obtained in the 2012 LiDAR flight. The results will allow us to assess whether it is possible to estimate -with an acceptable margin of error-the growth in the timber volume based on field data obtained over a 5-yr period and the metrics of a single LiDAR flight made during this time could optimise the available resources and would represent a large decrease in the cost compared to traditional estimation methods.

Study area
The study area lies in the Province of Bizkaia in the Autonomous Community of the Basque Country (ACBC) in the north of Spain. Its center is the municipality of Muxika and data from the surrounding municipalities between the coordinates 2º 30' -2º 54' west and 43º 06' -43º 24' north are also used. The average height is 230 m. The mean annual rainfall is 1277 mm and the mean annual temperature is 13.8 ºC.
Frosts are infrequent and no physiological drought is apparent. With regard to the silviculture applied in the area, the rotation cycle is, on average, 35 yrs, varying between 30 and 40 according to the site quality and the price of wood at the time when harvest takes place. Clear cutting are applied followed by repopulation with initial densities of 1300-1600 trees ha -1 . Between the 7th and 9th year, non-commercial thinning takes place, eliminating 500 trees ha -1 and pruning up to 2.5 m on the trees that have not been eliminated. Between the 12th and 15th year, the first commercial thinning is carried out, leaving between 500 and 700 trees ha -1 . Between the 17th and 20th year, the second commercial thinning is applied, leaving between 400 and 500 trees ha -1 . The third thinning is carried out between the 22nd and 25th year, leaving between 250 and 300 trees for the clear cutting. The total production is on average 550 m 3 ha -1 in one rotation, and therefore the average growth is 14.85 m 3 ha -1 yr -1 . Volumes and growths of 3153 ha occupied by P. radiata in the municipality of Muxika were estimated. The area in which data was gathered and results were obtained can be seen in Fig. 1.

Field data
Two sets of data were used according to the year in which these were gathered. In order to fit the wood volume model from 2011, data was used from the Fourth National Forestry Inventory (IFN4) of the Province of Bizkaia (http://www.euskadi.eus/ inventario-forestal-2011/web01-a3estbin/es/). Initially, all the plots included in the IFN4 and located within the boundaries of the municipality of Muxika and neighboring municipalities were selected. From these plots, 111 located within the area occupied by P. radiata were selected, according to the 2010 Forestry Map of the Basque Country (http://www.euskadi.eus/mapaforestal-del-pais-vasco-ano-2010/web01-a2hestat/es/). From these 111 plots, 13 were eliminated as they had no data or did not have P. radiata as their main species. From the remaining 98, those that showed a clear incompatibility between the field data and the vegetation visible in the aerial photography were eliminated. This occurred because the coordinates of the centers of the plots were taken with GPS navigators that did not allow differential correction and had a potential error of more than 12 m. For the same reason, those plots in which the difference between the dominant height of the field trees and the 95th percentile of the height of the points of its LiDAR points cut was greater than 12 m were eliminated. Furthermore, those plots with a basal area (G) greater than 70 m 2 ha -1 were also eliminated as these were not considered possible for this species in this area. Following the application of these filters, only 8 additional plots were discarded. Therefore, 90 plots were used finally in this study. Their location can be seen in Fig. 1.
IFN4 plots are of variable radius, so that a tree will or will not be measured in accordance with its diameter and the distance to the center of the plot. For each selected Figure 1. The area in which the study was carried out. The area in which field data were gathered is shown in green. The area in which the study was conducted and the results were obtained is shown in red. The points correspond to the location of the IFN4 plots selected for the study. plot, the number of trees per hectare (N, trees ha -1 ), the dominant height (Ho, m) defined as the mean height of the 100 thickest trees per hectare, the mean height (Hm, m), the mean diameter (Dm, cm), the basal area (G, m 2 ha -1 ) and the wood volume over bark (Vob, m 3 ha -1 ) were calculated. For the last of these calculations, the Excel add-in CUBIFOR (Rodríguez et al., 2008) was used. With the aim of reducing data-gathering costs for subsequent years, a Vob tariff was adjusted in accordance with G and Ho. In this way, in order to calculate volumes in measured plots it is sufficient to measure both variables quickly and economically. The Vob estimated with this tariff will be considered as the field Vob for the entire study.

Field data
In order to be able to assess the possibility of using LiDAR data from a single flight when estimating the forestry resources in subsequent years, 55 field plots were measured during the month of December, 2015. So as to prevent any lack of accuracy in locating the center of the plot from affecting the estimation of the final results, those areas with a sufficient level of uniformity were selected so that a displacement of up to 10 m in the center of the plot would have only a slight effect on the final result.
With the aim of locating uniform areas, two influential variables were selected, the first being Canopy Cover Fraction (CCF), as a percentage of the surface covered by the canopy of trees, and the second being the height of the vegetation. The data gathering area was divided with a 10 m-sided grid and in each 100 m 2 square and the CCF was estimated as a percentage of first returns, higher than 3 m over the total of first returns of the LiDAR points contained in it. For each grid, the typical deviation and the average of its CCF value and those of the 8 surrounding squares were calculated and those squares in which the coefficient of variation (CV) was not greater than 30% were selected if the CCF was lower than 70% and those in which the CV was greater than 10% if the CCF was greater than 70%. To estimate the height of the vegetation, the data gathering area was divided with a 10 m-sided grid and in each grid cell, the height of the 95th percentile of the heights of the LiDAR points contained in it was calculated. For each grid cell, the mean and typical deviation of its value and that of the 8 surrounding cells was calculated. Those grid cells in which the CV was not greater than 20% were selected if the 95th percentile was less than 8 m; those grid cells in which the CV was not greater than 10% were selected if the 95th percentile was between 8 and 12 m, and those in which the CV was not greater than 5% were selected if the 95th percentile was greater than 12 m. Only those grid cells that met both conditions were candidates for containing a plot center for data gathering. In accordance with these conditions, 55 points were found and located in the field with a GPS navigator. From each point, the dominant height was measured as the average of the thickest three trees within a radius of 9.77 m from the point that corresponds to a plot with a surface area of 300 m 2 . The basal area was also measured by relascope inventory from the center of the plot. The measurement of these two variables allows the field work to be done cheaply and quickly. Later, the Vob was calculated for each plot applying the Vob tariff in accordance with Ho and G calculated with the plots from the IFN4. The values and statistics of the variables used in the measurements of the plots from the IFN4 (2011) and the new plots (2015) are shown in Table 1.

LiDAR Data
The LiDAR data used came from the flight over the Autonomous Community of the Basque Country in July, 2012, with a mean density of 1 pulse m -2 and up to 4 returns per pulse.
For the entire data-gathering area, those las format files that contained a selected plot from the IFN4 or some of the plots measured in the field in 2015 and all the sheets contained partially or completely in the euskadi.net/lidar/). Once the points clouds had been normalized, the corresponding cylinders were clipped to the projection of each one of the plots, the center of which was that of the IFN4 plots and radius 15 m. In this way, 90 sets of LiDAR points coinciding with the 90 IFN4 plots were obtained. For example, the cut of plot 559 can be seen in Fig. 2.
In order to make this clip, the PolyClipData tool of the FUSION program (http://forsys.sefs.uw.edu/fusion/ fusionlatest.html) was used. Subsequently, the statistical parameters that describe each of the points clouds corresponding to plots were calculated. The results of the LiDAR points cloud statistics were compared to the results of the field measurements and volumes for each plot. In this way a table was obtained in which each row contains the field data and Vob of each IFN4 plot and the statistics corresponding to its LiDAR points cut. The selected plots and their data were separated in two groups of 60 and 30 plots respectively. The first of these was used to fit the model that relates the Vob with the statistics of the LiDAR points cloud and the second was used to validate this model. In order to separate the two groups, a random number between 0 and 100 was generated for each plot. Those plots with assigned numbers of under 65 were used to fit the model; those plots with a number over 65 were kept in order to validate the model. The same classification was used for adjusting the Vob tariff in accordance with G and Ho. The location of the plots of the model and of the validation group can be seen in Fig. 3.
The same procedure was followed using the LiDAR data for 2011 and the plots measured in the field in 2015. Fifty-five cuts of LiDAR points were obtained with their statistics and 55 Vob data calculated for each plot measured in 2015. Forty plots were selected to fit the model and the validation group was formed by the 15 plots. It should be pointed out that none of the plots measured in 2015 was located in areas in which clear cutting had been performed during the 5-yr period 2011-2015. This meant that the fitted model did not consider the possibility of volume 0. Some plots were indeed located in areas that had been subjected to thinning during the 5-yr period.

Clear cutting and thinning data in Muxika during the period 2011-2015
Clear cutting and thinning data in the municipality of Muxika between 2011 and 2015 was gathered in order to be able to use this to estimate the growth in wood volume during the period 2011-2015. Data on harvest provided by the Provincial Council of Bizkaia (the administrative office that manages the forest stands of the Province of Bizkaia), can be seen in Table 2.
In order to adjust the Vob tariff per hectare in accordance with the stand variables G and Ho, 90 plots of the IFN4 were used. Of these, 60 were used to fit the model and 30 to validate this. Plots were assigned to the model group or the validation group according to the random number assigned to them. In order to build the tariff, the STATGRAPHICS (http://www.statgraphics. com) was used. Using the regression model selection option, Vob was used as a dependent variable and  Ho, Ho 2 , G, G 2 and Ho*G were used as independent variables (Table 3). The models were examined to determine whether all the terms should be retained in the final regression equations. Once the stepwise regression had been applied, the Mallows' Cp selection method was used as the best selection technique (Mallows, 1973). Only models whose parameters were significant at a given level (in our case, 5%) were taken into account. The best model was chosen with the following indicators: adjusted determination coefficient (R 2 ), root mean square error (RMSE) and mean absolute error (MAE). The heterocedasticity was examined visually, by plotting residuals as a function of predicted values and using the Breusch-Pagan test (Breusch & Pagan, 1979). Moreover, the non-autocorrelation condition was checked using the Durbin-Watson test (Durbin & Watson, 1971). The estimations of the different fitted models were analyzed by means of numerical and graphical analysis. Once the model had been fitted, the Vob in each plot from the IFN4 was calculated in accordance with Ho and G. This value was the one considered for each plot as a Vob value measured in the field and was the dependent variable in the stem volume models in accordance with LiDAR statistics for 2011 and 2015.

Fitting the stem volume model based on 2011 LiDAR data (Vob_2011)
Linear models were used to establish empirical relationships between field measurements and LiDAR variables. In order to fit the Vob for 2011, the Vob calculated in each IFN4 plot based on the model fitted with the independent variables G and Ho was used as the dependent variable. In the STATGRAPHICS regression model selection procedure, all statistics relating to the height of the points obtained with the FUSION CloudMetrics order were included as independent variables. In order to fit the model, 60 plots from the IFN4 were used. The other 30 plots were used to validate the model. The highest model efficiency in accordance with the LiDAR metrics was obtained by ordering the proposed combinations according to the best Mallows Cp. The same goodness-of-fit test utilized for the models in the previous section was used. Likewise, for the selected models, a residuals analysis was made for the selected models, with the aim of identifying evidence of normality deviations, lack of adjustment and/or heterocedasticity.

Fitting the 2015 stem volume model (Vob_2015)
In order to fit the Vob for 2015, the Vob was used as dependent variable, calculated in each one of the 55 plots measured in the field in 2015 with the previously adjusted tariff of G and Ho. In order to maintain the uniformity of the results obtained, the model was fitted with the same independent variable as the Vob_2011 model, in other words, the 70 percentile of the height of the points cloud was used. In order to fit the model, 40 plots were used. The other 15 plots were used to validate the model.

Application of LiDAR data to Vob models
In order to apply the wood volume models obtained, the area occupied by P. radiata in Muxika (3,153 ha) was separated and a grid made up of 26.59-m grid cells was superimposed. Each grid cell had a surface area equal to a circular plot with a 15-m radius. Therefore, each grid cell was of 706.86 m 2 . The cells that intersected the polygons corresponding to the surface areas of P. radiata according to the forestry map of the Basque Country, were singled out from the vector layer. In this way, a 26.59-m grid mask was obtained with the location of P. radiata forests in the municipality of Muxika. Parallel to this, the GridMetrics procedure of FUSION was applied to all the LiDAR files that intersected the municipality of Muxika. In order to apply GridMetrics, a grid cell size equal to that of the 26.59-m grid was used. The work process diagram can be seen in Fig. 4.

Vob tariff per hectare
The lowest Mallows Cp value was 0.51. This corresponded to the model (3) that had a single independent variable: Ho*G (Table 3). This is a model with a high R 2 (0.99), although it presents heteroscedasticity (p<0.05) for the Breusch-Pagan test. This occurs because the higher the volume value, the greater the absolute error.

Vob_2011
The lowest Cp value was 3.28 and this corresponded to the model that had the value of the height of the percentile 70 (h 70 ) of the points cloud (Eq. [4], Table 4) as single independent variable. This model had a good R 2 (0.85) but was heterocedastic with a p< 0.05 in the Breusch Pagan test. Once again, this occurs because the higher the volume value, the greater the absolute error.

Vob_2015
The stem volume model for 2015 based on LiDAR data is shown in Table 5.
The models obtained allowed us to estimate the total volume over bark with a different degree of accuracy, although in all cases it is quite high. The best regression uniformity of the waste of the three fitted models, no evidence was found of any non-compliance with the assumptions.
The mean values of volume with bark of the three models are shown in Table 8.

Results of applying the Vob models to the entire surface area of Muxika
By applying the volume per hectare of wood with bark models to the entire surface area of Muxika for 2011 and 2015, the results shown in Table 9 were obtained.
The results of the distribution of the Vob for 2011 and 2015 and the difference between both can be seen in Fig. 6.

Annual increase in volume over bark
The difference between the Vob estimated for 2015 and that estimated for 2011 is 81,603 m 3 for the entire surface area of P. radiata in the municipality of model corresponds to the Vob model based on stand variables (H 0 , G). This model presents the highest value of the adjusted coefficient of determination (R 2 adj = 99), followed by the volume estimation model in accordance with LiDAR metrics (Vob_2011, R 2 adj = 85) and lastly, the model based on the field data obtained in 2015 (Vob_2015, R 2 adj = 80) (Table 6). Likewise, the lowest values of RMSE, MAE and ME corresponded to the best model (Vob) obtaining in this the RMSE values (13.61-12.70), MAE (10.04-9.86) and ME = (0-3.49) obtained in the adjustment plots and the validation, respectively (Table 6). Models Vob_2011 and Vob_2015 exhibited very similar values in RMSE (~ 60 m 3 ha -1 ) and MAE (46 m 3 ha -1 ). Residual plots for Vob prediction indicated that the proposed model was appropriate, as can be seem from the predicted versus the observed values for Vob in Fig. 5.
Each of the adjusted equations showed good fit statistics. All the parameters were significant at a confidence level of 95% (Table 7). In the graphical evaluation of the assumptions of normality and      Muxika. Therefore, the annual mean growth of the volume over bark (ΔVob) obtained during the 5-yr period was 5.14 m 3 ha -1 yr -1 (Table 9). The forestry harvests data shows that the wood volume extracted during the period 2011-2015 was 135,624 m 3 (Table 2) and therefore, on average, 8.60 m 3 ha -1 were extracted. The estimated mean growth for the period 2011 -2015 is 13.74 m 3 ha -1 yr -1 , this being the value obtained by adding the estimated mean growth in the aforementioned period (ΔVob = 5.14 m 3 ha -1 yr -1 ) and the mean volume extracted 8.6 m 3 ha -1 .
As can be seen in Fig. 7, the fitted models predict the greatest increases for Vob of between 0 and 150 m 3 ha -1 . The increases diminish gradually for Vob values close to 300 m 3 ha -1 and become negative above that value.

Discussion
Both the R 2 and the RMSE of the model fitted for the Vob of 2011 are similar to those published by other authors under similar conditions (Packalen et al., 2011;Stone et al., 2011;Tonolli et al., 2011;Watt et al., 2013). The values of R 2 and RMSE of the Vob model for 2015 are more unfavorable than those obtained under similar working conditions. This is due to the fact that in the preparation of the Vob_2011 model, the field data and the LiDAR data were almost contemporary, i.e., they were taken just a few months apart. However, for the fit of the Vob_2015 model, the field data is more than four years older than the LiDAR data. This time difference causes a reduction in the ability to predict the fitted model, even more so if we take into consideration that P. radiata is a fast-growing species and is normally subjected to intensive forestry regimes that bring about considerable changes within a short time span. With regard to adapting the results obtained, the estimated mean growth being 13.74 m 3 ha -1 yr -1 , the volume will have increased 68.7 m 3 ha -1 on average during the period 2011-2015. Therefore, in those places where harvest was more than 68.7 m 3 ha -1 , the growth has failed to offset the volume extracted and as a result will have less volume in 2015 than in 2011. These are the grid cells with a negative growth value. On the other hand, the mean volume extracted in a thinning operation according to Table 2 is 39.2 m 3 ha -1 . It is for this reason that the areas that were subjected to thinning between 2011 and 2015 must correspond to cells with a volume difference of less than 29.5 m 3 ha -1 as the 68.7 m 3 of growth in the 5 yrs, minus the 39.2 m 3 extracted in thinning operations will have allowed this amount to accumulate on average in these areas. Those grid cells with a volume difference of under that amount add up to 1,552 ha. This is the surface area in which, according to the fitted models, some kind of harvest took place during 2011 and 2015 and coincides with those contributed as real registered harvest: 1068.1 ha. As shown in Table 1, the minimum Vob of the plots measured in the field data taken in 2015 was 90.84 m 3 ha -1 and, therefore, no plot in which clear cutting had taken place was measured. Consequently, the 2015 model cannot detect the clear cutting that took place during this period and assimilates these to intensive thinning, and therefore the surface area of harvest estimated by the models includes both clear cutting and thinning. This circumstance must be borne in mind for subsequent applications of the method and those plots subjected to clear cutting during the period of use of LiDAR data. It can be seen, moreover, that the areas in which the models predict greater harvest intensity correspond to those that had a higher volume in 2011. These areas are the ones in which third thinnings and clear cuttings (being the harvest of greater weight) are normally carried out. Both the assimilation of clear cutting to intensive thinning and the fact that the areas that were actually thinned or clearcut after July 2012 are captured in the raster models, which are based on a Lidar dataset from July 2012, can be explained when analyzing the diagram resulting from a comparison of Vob_2011 and the growth between 2011 and 2015 (Fig.  7).
One of the aspects that must be clarified in subsequent studies is the development in growth in annual periods between both ends of the 5-yr interval. In this study it has not been possible to deal with this question due to a lack of field data in this intermediate period. Moreover, due to a lack of funding, locating the IFN plots was not possible for the second gathering of field data and for this reason the samples for 2011 and 2015 were different. Even though this circumstance might reduce the validity of the results obtained, it was decided to begin the study to verify whether this methodology could be used in subsequent studies with higher funding levels. Once sufficient funds are available, field data must be taken on an annual basis and in the same plot sample.
With regard to the multi-annual use and estimation of growth using a single LiDAR flight, no published references have been found concerning similar methods. This is the main contribution of this paper: the use of LiDAR data in successive years through the updating of the mass conditions with fast, reliable and cost-effective annual field data gathering. In this way, the investment involved in LiDAR data gathering is justifiable as this can be amortized in several years.
Fitted models that relate the descriptive variables of the points cloud of a LiDAR flight with the wood volume field data obtained in field plots with a 5-yr interval, provides a good approximation to the growth values produced in masses of P. radiata in this region (13.74 m 3 ha -1 yr -1 ), as well as an estimation of the surface area subjected to harvest in this period. This allows the multi-annual use of LiDAR data gathered in a campaign with the consequent decrease in costs. To ensure the highest levels of accuracy in the estimation of the growth and volume of wood during the years following the gathering of LiDAR data, it is necessary to include in the sample of plots to be measured those in which both thinning and clear cutting operations have been carried out.