2.1. Study area
The research site is located in Redes Natural Park which occupies the eastern central area of the Principality of Asturias in the north of Spain (Figure 1). The annual rainfall ranges from 987 to 1351 mm, with an average annual temperature in the area of 10–11°C. The study site has a northern orientation with high slope (45-57%) and is between 600 and 700 metres above sea level. This study was conducted as a pilot forest management trial made up of three different treatments in sweet chestnut coppice stands: (1) Control, where there were no management operations (3300-3700 stems ha-1), (2) Treatment 1, which consisted of one thinning that left a living stock density of between 600 and 900 stems ha-1 and (3) Treatment 2, which involved a more intensive thinning that resulted in a living stock density of around 400 stems ha-1. Thinning treatments were carried out at the end of 2015, in winter when the sweet chestnut loses its leaves. In the Control and Treatment 1 plots Oak (Quercus Petrea) is also present.
2.2. Field data: Forest inventory and LAI measurements
Forest inventory data were collected from three established permanent plots (r=15 m) in the study area. Each tree in every plot was labelled and diameter (in cm) at breast height (d) (diameter at 1.3 m above the top of the stool) was measured with a tree calliper to the nearest 0.1 cm, as was total height (h), to the nearest 0.1 m, using a digital hypsometer. Using this information certain stand variables were calculated to characterize the study area. These plots were inventoried twice: first in winter 2015 (installation year) and then in summer 2019 (Figure 2).
LAI field measurements were collected in July 2019, using an LAI-2200C plant canopy analyzer, from 21 circular plots (r=10 m). The location of the plots was based on a grid system whereby where each plot centre was separated from those around it by 30 m. The locations of the plot centres were recorded using a GPS TrimbleExplorer XH™ (Trimble, Sunnyvale, CA, USA) with submetric accuracy and the total number of LAI plots per trial depended on the total surface area covered by each treatment.
The LAI-2200C is a portable instrument, in this case consisting of a control unit and an optical sensor, that directly provides an LAI value that may include vegetation which is not leaves and is thus referred to as Effective LAI (henceforth here, LAIe) (Chen and Black, 1992). This device has been widely used in LAIe studies and several extensive review papers have established it as an appropriate tool for LAIe field measurements (e.g. Bréda, 2008; Jonckheere et al., 2004; Morsdorf et al., 2006; Jensen et al., 2008; Solberg et al., 2009; Korhonen et al., 2011; Pearse et al., 2017; Fang et al., 2019). The LAI-2200C incorporates five concentric rings with central zenith angles of 7, 23, 38, 53 and 68 degrees in a “fish-eye” optical sensor, measuring in the blue band (320–490 nm). The measurements of LAIe were made below the canopy and also in a clearing, the latter to simulate the light falling directly on the crown. These two measurements were required in order to calculate the ratio between the two transmittances (above and below the canopy). This ratio was recalculated into LAIe as described in the instrument manual (LI-COR Inc., 2015). All measurements were performed under uniform overcast skies to reduce the effect of scattered light in the canopy and the sensor was equipped with a 90° view restrictor in accordance with the manufacturer’s instructions.
In the forest, LAIe field measurements were conducted using the plot layout employed in similar previous studies (Solberg et al., 2006, 2009) within the circular plots mentioned above. In each plot, one measurement was taken at the plot centre, and another four—one at each cardinal point—were made at 3 metres from the centre. Due to the differences in the canopy cover, to guarantee a low standard error of LAIe, supplementary measurements (n=18) were collected at random within the plot.
2.3. LiDAR data
LiDAR data were collected in leaf-on conditions in summer 2019, at the same time as the LAIe field data were measured, using a Velodyne VLP-16 LiDAR scanner mounted on an unmanned aerial vehicle (drone). The laser wavelength of the device was 905 nm with a field of range of ±15˚ Vertical FOV / 360˚ Horizontal FOV. The sensor recorded a maximum of 2 returns per pulse, and a minimum density of 25 points m−2 was achieved over the area. LiDAR data was recorded for the all the stands in the study area, which included all the LAIe plots.
To calculate the plot-level forest structural parameters, the following steps were followed in order to obtain the LiDAR metrics needed to model LAIe for sweet chestnut coppice (Figure 3).
The LiDAR point cloud was classified as either background or not background with FUSION software V4.00 (McGaughey, 2020), after which the point cloud was normalized against the ground surface height. The LiDAR standard elevation metrics were also computed with FUSION taking different fixed radii starting with the 10 m radius established for the LAIe plots in the field and increasing the radius by 2 m each time up to a limit of 30 m, following the methodology described by Pearse et al. (2017). This was done because the LiDAR and the LAI-2200c were not sampling the same canopy volume, LiDAR uses a vertical cylinder and the LAI-2220c an inverted cone. Moreover, since sweet chestnut resprouts from stump and this generates another canopy layer in addition to the principal one and because the volume of the cone depends on the height of the canopy, the LiDAR metrics were computed for different heights (2, 3, 4, 5 m). The optimum choice of the height threshold is crucial in this type of study in order to reliably identify ground hits as those that are below the height threshold and those above it as corresponding to the canopy (Zhao and Popescu, 2009). Following the testing of the various radii and height measurements, those most representative of the field measurements were selected to form the basis of the LiDAR standard elevation metrics.
Among the LiDAR standard elevation metrics extracted with FUSION, the KDE (kernel density estimation) metrics were considered essential to differentiate the layers within the canopy because it calculates the number of peak heights (minimum and maximum values and the range between them) using a Gaussian kernel to construct a probability density function for the peaks in each plot. In addition, the Profile Area, a new metric of FUSION V.4.00, was used, which corresponds to the area under the height percentile profile or curve (Hu et al., 2019). Both metrics are considered important in describing canopy vertical structure, so they were computed from ground to maximum height.
Another set of LiDAR metrics based on the canopy structure (LiDAR canopy metrics) were also calculated from ground to treetop following the methodology developed by Lefsky et al. (1999) and amended by Coops et al. (2007). Using the LiDAR point cloud data, for each circular plot the canopy was divided vertically from the ground to the maximum height into 1 metre intervals, resulting in a series of small cylinders, which were then classified into one of four canopy classes. First each cylinder was classified as empty or filled depending on the absence or presence of LiDAR points. After that, the empty cylinders were divided into “Open gap” if they were situated above the canopy and “Closed gap” if they were below it. In the case of the filled cylinders they were classified as “Euphotic” if they were located in the uppermost 65% of all filled cylinders and “Oligophotic” if below it (Figure 4).
Using this same point cloud cylinders, the Weibull density function was used to describe the Canopy Height Distributions (CHD) in the CVP. Weibull function is commonly applied to characterize foliage distribution for various species due to its flexibility in representing various distribution shapes (Lovell et al., 2003; Coops et al., 2007; Zhang et al., 2017; Melo et al.,2018). The two Weibull coefficients, i.e. scale (α1) and shape (β1) were fitted as explained in previous studies using LiDAR data (Lovell et al., 2003; Coops et al., 2007). Scale provides the position and the distribution of the movement in vertical scaling (the shape of the distribution density curve), while shape explains the capacity to increase or decrease the breadth of the distribution.
Density-based LiDAR metrics were also calculated as the proportion of points above the percentiles (1st, 5th, 10th, 20th, 25th, 30th, 40th, 50th, 60th, 70th, 80th, 90th, 95th, and 99th) as was proposed by Zhang et al. (2017). A summary of the LiDAR metrics with their corresponding descriptions are shown in Table 1.
Table 1. Description of LiDAR metrics for modelling LAIe.
LiDAR metrics
|
Description
|
LiDAR Standard Elevation Metrics
|
All returns
|
All returns from ground to maximum height of the LiDAR point cloud.
|
All returns above height threshold
|
All returns from height threshold to maximum height of the LiDAR point cloud.
|
Proportion return above height threshold
|
All returns from height threshold to maximum height of the LiDAR point cloud, expressed as a percentage.
|
First returns above height threshold
|
Only first returns from height threshold to maximum height of the LiDAR point cloud.
|
Proportion of first returns above height threshold
|
Only first returns from height threshold to maximum height of the LiDAR point cloud, expressed as a percentage.
|
Canopy relief ratio
|
Canopy cover generated by FUSION
|
Skewness and Kurtosis of heights
|
The skewness and kurtosis of the heights of first returns with respect to the total.
|
Percentile Heights
|
The percentiles of the canopy height distributions (1st, 5th, 10th, 20th, 25th, 30th, 40th, 50th, 60th, 70th, 80th, 90th, 95th, and 99th).
|
Mean, maximum and minimum height
|
Mean, maximum and minimum height above ground of all points.
|
KDE (kernal density estimator)
|
Number of peak heights with minimum and maximum values
|
Profile area
|
Area under the height percentile profile or curve
|
LiDAR Canopy Metrics
|
Return proportion
|
Proportion of returns in each 1 metre of height of the cylinders from the ground to the maximum height
|
Mean return proportion
|
Mean of the proportion of returns considering different ranges of interval, for example from the height threshold to the maximum height.
|
Density-based metrics
|
The proportion of points above the percentile heights compared to total number of points
|
Open and Closed gap zones
|
Proportion of empty cylinders located above and below the canopy, respectively
|
Euphotic and Oligophotic zones
|
Respectively, the proportion of cylinders located within the uppermost percentiles (above 65%) of all filled cylinders and those below the same point in the canopy vertical profile
|
α1 and β1 parameter of Weibull distribution
|
The scale parameter, α, and shape parameter, β, of the Weibull density distribution fitted to CHD.
|
2.4. LIDAR metrics selection and LAIe models
After extraction of the LiDAR metrics and using the selected radius and height threshold, the LiDAR data must be related to LAIe. Adding the LiDAR standard elevation metrics to those related to CVP computed from the LiDAR point cloud, generates a large dataset of variables that can be used for model development. To address this issue which is common to this type of study (e.g. Pearse et al., 2017; Zhao and Popescu, 2009), the LiDAR metrics to be used need to be decided prior to model development. To do this in this work, the metrics (Table 1) with low correlations with LAIe (r < 0.6) or highly correlated with other variables were excluded and then the candidate metrics were used in the regression analysis.
Stepwise linear and non-linear allometric modelling were chosen in contrast to non-parametric statistical learning approaches, because with small samples, it is usually wise to use simple methods (Pearse et al., 2016). However, the choice of the model for predicting LAIe from LiDAR data must be determined by the modeller, and in most cases it is dependent on the number of in-situ LAIe observations available (Zhao and Popescu, 2009).
As final model selection criteria, the coefficient of determination (R2) and the root mean square error (RMSE) were used. R2 indicates the proportion of variation explained by the model, while RMSE, which uses the same units as the dependent variable, gives an idea of the mean error when using the model.
Also, as validation is important in terms of evaluating the predictive performance of any models developed, the leave-one-out cross-validation was used. It considers the difference between the observed value and the predicted value whereby one piece of data in turn is omitted from the analysis to fit the model and the model is fitted to the remaining n-1 data. The cross-validation of each model was based on the analysis of RMSE and R2 statistics.
2.5. LiDAR-LAIe Map tool
The selected model was used to develop a tool to generate an automatic map of leaf area index at stand level. The processing of LiDAR data at stand level was programmed in a script, including the generation of a Digital Terrain Model (DTM) and CHM models with a pixel size of 10 m. Moreover, the calculation of all LiDAR metrics mentioned above and the selected LAIe model were incorporated into the tool (Figure 5). The process starts with the selection of the appropriate LiDAR data and the creation in the area of interest (AOI) followed by the calculation of LiDAR metrics and ends with the generation of different raster maps: mean height (Hm), using the 95th Percentile Height; coverage (FCC), using the percentage of first returns and LAIe raster applying the adjusted model.