Postﬁre Tree Structure from High-Resolution LiDAR and RBR Sentinel 2A Fire Severity Metrics in a Pinus halepensis -Dominated Burned Stand

: Tree and plant structures remaining after ﬁres reﬂect well their degree of consumption, and are therefore good indicators of ﬁre severity. Satellite optical images are commonly used to estimate ﬁre severity. However, depending on the severity of a ﬁre, these sensors have a limited ability to penetrate the canopy down to the ground. Airborne light detection and ranging (LiDAR) can overcome this limitation. Assessing the di ﬀ erences between areas that have been burned in di ﬀ erent ﬁre severities based on satellite images of plant and tree structures remaining after ﬁres is important, given its widespread use to characterize ﬁres and ﬁre impacts (e.g., carbon emissions). Here, we measured the remaining tree structures after a ﬁre in a forest stand burned in SE Spain in the summer of 2017. We used high-resolution LiDAR data, acquired from an unmanned aerial vehicle (UAV) six months after the ﬁre. This information was crossed with ﬁre severity levels based on the relativized burnt ratio (RBR) derived from Sentinel 2A images acquired a few months before and after ﬁre. LiDAR tree structure data derived from vertical canopy proﬁles (VCPs) were classiﬁed into three clusters, using hierarchical principal component analysis (HPCA), followed by a random forest (RF) to select the most important variables in distinguishing the cluster groups. Among these, crown leaf area index (LAI), crown leaf area density (LAD), crown volume, tree height and tree height skewness, among others, were the most signiﬁcant variables, and reﬂected well the degree of combustion undergone by the trees based on the response of these variables to variations in ﬁre severity from RBR Sentinel 2A. LiDAR metrics were able to distinguish crown ﬁre from surface ﬁre through changes in the understory LAI and understory and midstory vegetation. The three tree structure clusters were well separated among each other and signiﬁcantly related with the RBR Sentinel 2A-derived ﬁre severity categories. Unburned and low-severity burned areas were more diverse in tree structures than moderate and high severity burned ones. The LiDAR metrics derived from VCPs demonstrated promising potential for characterizing ﬁne-grained post-ﬁre plant structures and ﬁre damage when crossed with satellite-based ﬁre severity metrics, turning into a promising approach for better characterizing ﬁre impacts at a resolution needed for many ecological processes.


Introduction
Wildfires burn heterogeneously through forested landscapes; while some patches may have the trees barely scorched, others will have severe damage in their canopies, with all leaves and small branches fully burned. The degree of consumption of the plant canopies has been used as a surrogate calculate fuel load. Nevertheless, the integration of optical remotely sensed imagery and LiDAR data provides more accurate information than either sensor type alone for classifying post-fire residual vegetation [11,12] or crown severity levels [13].
In this paper, we characterize the fine-grained tree structures, including their canopies, of trees in a large Pinus halepensis-dominated fire that occurred in south-eastern Spain (Yeste, Albacete) in the summer of 2017. We classified the damages to each individual tree by using LiDAR metrics derived from vertical canopy profiles (VCP). Our main objective was to characterize tree morphologies (from LiDAR) and verify the differences that exist in them after fire in relation to fire severity levels based on Sentinel 2A RBR metrics. More specifically, we evaluated the following: (i) which structural metrics (LiDAR-derived) are most important for separating individual trees into distinct groups, (ii) how these tree groups are distributed within fire severity categories, and also (iii) how these groups are related to RBR Sentinel 2A fire severity levels (how different were tree groups among them within each fire severity level, and how different were tree groups within them across different fire severity levels).
The main novelties of this work are as follows: (i) the use of different approaches (i.e., voxels, height bins and original points) to estimate vertical crown profile metrics at the tree level; (ii) using LAD profiles and LAI values to estimate crown properties; (iii) using the breakpoints method to calculate the canopy base height (CBH) for deriving the crown volume, and (iv) using multitemporal passive optical multispectral imagery to relate spectral fire severity indices with LiDAR data.

Study Sites
The study area was the Yeste fire (province of Albacete, SE Spain), that occurred in summer 2017. The fire started on 27th July, one day before an incursion of warm, tropical Africa air, and was controlled on August 1st, after burning 3217 ha (see [53] for further details). High-density LiDAR flights were carried out six months after the fire on three areas of 5-6 ha size with different levels of fire severity and in an unburned area adjacent to the fire ( Figure 1a). The unburned area was occupied by Pinus halepensis Mill. (60%) and shrublands. The low severity area was occupied by mixed Pinus pinaster Ait. and Pinus halepensis Mill. (74%) forests and open forests of Juniperus oxycedrus L. with Pinus halepensis Mill., and the moderate severity area was completely occupied by Pinus halepensis Mill. Forest. The highest fire severity area was entirely occupied by mixed Pinus pinaster Ait. and Pinus halepensis Mill. forest. The elevation ranges from 722 m to 1112 m and the slope from 14.6 • to 20 • across the four areas. Crosstabulation between LiDAR flights and RBR Sentinel 2A fire severity categories allowed us to broadly characterize each LiDAR flight regarding the distribution of RBR fire severity levels ( Figure 1b). Location of field parcels to validate tree delineation in LiDAR flights characterized by moderate and high severity. (c) Crosstabulation between LiDAR flights (columns) and RBR Sentinel 2 fire severity levels (rows) to broadly characterize each LiDAR flight by the RBR fire severity levels. LiDAR flight categorized as low severity had 91% of its area in RBR low severity level; LiDAR flight categorized as moderate severity had 13%, 49% and 38% of its area in low, moderate and high RBR severity levels, respectively; and LiDAR flight categorized as high severity had 51% and 48% of its area in moderate and high RBR severity levels, respectively.

LiDAR Data
Post-fire high-density LiDAR data were acquired using the TerraSystem-Lidarpod system formed by the Velodyne HDL-32e sensor with a range distance of 110 m. This system has 32 pairs of laser sensors (wavelength of 905 nm) with a rotation rate of 100 Hz, and dual-sensor first and last returns, collecting a total of 700,000 pulses per second. The average pulse spacing was between 0.01 and 0.04 m. The flights were made at an altitude of 40 m above the ground (footprint on the ground was 33 × 204 m), and in parallel lines 40 m apart with an overlap of 50%. The flight speed was between 5 and 7 m/s. At this flight height, the spatial distance (spacing) among points was 0.05 m, with a mean density of 300 points m −2 . The horizontal accuracy was <2 cm and the vertical accuracy <10 cm [54].

LiDAR Pre-Processing
The first step consisted in aligning forward and backward flight trajectories from the unprocessed Global Navigation Satellite System (GNSS) data and raw advanced navigation files from the Inertial Measurement Unit (IMU) (.anpp) using the Kinematica software [55]. Later, in the LasTools software [56], point clouds were filtered to eliminate duplicates and isolated points.
Next, ground returns were classified using the progressive Triangulated Irregular Network (TIN) (i.e., a Delaunay TIN) densification algorithm implemented in the lasground function of the LasTools software, using as parameters (i) the terrain type (-wilderness) with a step size of 6 m and (ii) the granularity ('-hyper_fine'). Moreover, we fine-tuned the algorithm by specifying (i) a threshold in meters at which spikes are removed from the ground (3 m), and (ii) the maximal offset (c) Crosstabulation between LiDAR flights (columns) and RBR Sentinel 2 fire severity levels (rows) to broadly characterize each LiDAR flight by the RBR fire severity levels. LiDAR flight categorized as low severity had 91% of its area in RBR low severity level; LiDAR flight categorized as moderate severity had 13%, 49% and 38% of its area in low, moderate and high RBR severity levels, respectively; and LiDAR flight categorized as high severity had 51% and 48% of its area in moderate and high RBR severity levels, respectively.

LiDAR Data
Post-fire high-density LiDAR data were acquired using the TerraSystem-Lidarpod system formed by the Velodyne HDL-32e sensor with a range distance of 110 m. This system has 32 pairs of laser sensors (wavelength of 905 nm) with a rotation rate of 100 Hz, and dual-sensor first and last returns, collecting a total of 700,000 pulses per second. The average pulse spacing was between 0.01 and 0.04 m. The flights were made at an altitude of 40 m above the ground (footprint on the ground was 33 × 204 m), and in parallel lines 40 m apart with an overlap of 50%. The flight speed was between 5 and 7 m/s. At this flight height, the spatial distance (spacing) among points was 0.05 m, with a mean density of 300 points m −2 . The horizontal accuracy was <2 cm and the vertical accuracy <10 cm [54].

LiDAR Pre-Processing
The first step consisted in aligning forward and backward flight trajectories from the unprocessed Global Navigation Satellite System (GNSS) data and raw advanced navigation files from the Inertial Measurement Unit (IMU) (.anpp) using the Kinematica software [55]. Later, in the LasTools software [56], point clouds were filtered to eliminate duplicates and isolated points.
Next, ground returns were classified using the progressive Triangulated Irregular Network (TIN) (i.e., a Delaunay TIN) densification algorithm implemented in the lasground function of the LasTools software, using as parameters (i) the terrain type (-wilderness) with a step size of 6 m and (ii) the granularity ('-hyper_fine'). Moreover, we fine-tuned the algorithm by specifying (i) a threshold in meters at which spikes are removed from the ground (3 m), and (ii) the maximal offset in meters (0.05 m) up to which points above the current ground estimate are included. Once the points were classified as "ground or not ground", the ellipsoidal heights of the point clouds were transformed into orthometric heights, using the geoid file from the National Geographic Institute (NGI) [57]. After that, the orthometric height was normalized so that the height of the vegetation that was above the ground line (0 m) could be obtained ( Figure 2). Finally, using the normalized LiDAR cloud, the canopy height model (CHM) was produced with a 0.1 m resolution.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 22 in meters (0.05 m) up to which points above the current ground estimate are included. Once the points were classified as "ground or not ground", the ellipsoidal heights of the point clouds were transformed into orthometric heights, using the geoid file from the National Geographic Institute (NGI) [57]. After that, the orthometric height was normalized so that the height of the vegetation that was above the ground line (0 m) could be obtained ( Figure 2). Finally, using the normalized LiDAR cloud, the canopy height model (CHM) was produced with a 0.1 m resolution.

Tree Detection and Vertical Canopy Profiles (VCP)
CHM, individual tree detection and crown metrics computation were performed using the LidR package [58] in the R software [59] in four steps ( Figure 3). First, we applied a "spike-free" algorithm [60] to build the CHM at 0.1 m pixel. This method resolves the problem of the empty pixels and socalled "pits" by the interpolation of first returns with a TIN, and then rasterizing it onto a grid. It consisted of several layers of triangulation at different height bins (at 0, 2, 5, 10, 15, 20 and 30 m). Moreover, the pit-free algorithm combined with the subcircling tweak, which replaces each return with a circle with a small radius (i.e., option '-subcircle 0.25′), worked adequately without empty pixels or pits. This algorithm is implemented with the 'grid_canopy' function in the LidR package [58] in R software [59]. Afterwards, the CHM was smoothed using a Gaussian filter (window size of 3 × 3 and a sigma value of 0.5). Second, individual trees were detected from the smoothed CHM, applying a tree detection algorithm based on a local maxima filter using a moving window of 3 × 3 m and 4 m as the minimum tree height. Third, crown segmentation was achieved based on local maxima plus Voronoi tessellation, using a threshold of 4 m, following Silva's method (see [61] for further details), and the computed maximum height (tree height) and crown area (CrArea, m 2 ) in the

Tree Detection and Vertical Canopy Profiles (VCP)
CHM, individual tree detection and crown metrics computation were performed using the LidR package [58] in the R software [59] in four steps ( Figure 3). First, we applied a "spike-free" algorithm [60] to build the CHM at 0.1 m pixel. This method resolves the problem of the empty pixels and so-called "pits" by the interpolation of first returns with a TIN, and then rasterizing it onto a grid. It consisted of several layers of triangulation at different height bins (at 0, 2, 5, 10, 15, 20 and 30 m). Moreover, the pit-free algorithm combined with the subcircling tweak, which replaces each return with a circle with a small radius (i.e., option '-subcircle 0.25 ), worked adequately without empty pixels or pits. This algorithm is implemented with the 'grid_canopy' function in the LidR package [58] in R software [59]. Afterwards, the CHM was smoothed using a Gaussian filter (window size of 3 × 3 and a sigma value of 0.5). Second, individual trees were detected from the smoothed CHM, applying a tree detection algorithm based on a local maxima filter using a moving window of 3 × 3 m and 4 m as the minimum tree height. Third, crown segmentation was achieved based on local maxima plus Voronoi tessellation, using a threshold of 4 m, following Silva's method (see [61] for further details), and the computed maximum height (tree height) and crown area (CrArea, m 2 ) in the same R package (Table 1). Finally, we delineated crown contours using the raster R package [62]. The fourth step consisted of clipping all the returns of the normalized heights within simplified crowns to obtain the following ( Table 1): (i) cloud metrics (standard height and intensity distribution metrics), including all height values and for two vegetation strata (<4 m and >4 m), and the density of different vegetation layers using raw points; (ii) diversity metrics derived from the density of height bins at 0.5 m from the ground up to 20 m (i.e., Simpson and Shannon diversity indices, and evenness dominance indices), using the Biodiversity R package ( [63]); and (iii) leaf area density (LAD, leaf area contained in unit crown volume) profiles based on voxels using the LeafR package [47] (Figure 3 and Table 1).
Tree delineation validation in the field was carried out on four plots of 30 m diameter, located randomly in the moderate fire severity LiDAR flights. We located LiDAR-derived tree crowns (≥4 m) in the field, using the X,Y position of the treetops of each crown segment. The locations of the treetops of selected trees were recorded using a Trimble Geo 7x (Sunnyvale, CA, USA), with post sampling corrections performed by the GPS provider. We measured the diameter at breast height (DBH) for each tree (n = 113). The height of the trees was not measured in the field. We identified 100% of the trees selected from LiDAR in the field. Moreover, we found a very good match between the tree crowns delineated by low-density pre-fire LiDAR data [53] and the post-fire crown segments ( Figure S1). Finally, we related field-estimated DBH with the LiDAR-estimated tree height and crown area to assess the coherence of the results ( Figure S2). In this study, we have tested various segmentation algorithms (i.e., watershed, Dalponte and Lee using LidR package in R software), and we have ascertained that Silva's tree segmentation [61] algorithm was better than the other ones because it did not merge so many adjacent trees together. Overall, the tree delineation was successful where the trees were more dispersed and thinner, such as in more severely burned areas [13,30,61,64].
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 22 to obtain the following (Table 1): (i) cloud metrics (standard height and intensity distribution metrics), including all height values and for two vegetation strata (<4 m and >4 m), and the density of different vegetation layers using raw points; (ii) diversity metrics derived from the density of height bins at 0.5 m from the ground up to 20 m (i.e., Simpson and Shannon diversity indices, and evenness dominance indices), using the Biodiversity R package ( [63]); and (iii) leaf area density (LAD, leaf area contained in unit crown volume) profiles based on voxels using the LeafR package [47] ( Figure 3 and Table 1).

LAI and LAD Profiles
Based on vertical profiles of the heights within each crown polygon, we calculated LAD profiles and LAI values, following the methodology proposed by Almeida [47] (Figure 4). The point cloud of each tree was binned into voxels (canopy volume units) with an XYZ resolution fixed at 1 × 1 × 1 m. For each voxel, we calculated the number of pulses that entered and the number of pulses that passed through that voxel. Next, the MacArthur-Horn equation was applied to each voxel to compute its LAD (LAD i ). From the LAD profiles, the total LAI (the sum of LAD from 1 m height to the top of the tree), the relative understory LAI (percentage of the total LAI from 0 to 3 m) and the crown LAI (total LAI-absolute understory LAI) of each tree were obtained ( Figure 4 and Table 1). Since we only used points with a maximum of 15 • off-nadir view, we could work under the assumption that each LiDAR pulse was vertically incident (see Almeida [47] for further details).
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 22 post sampling corrections performed by the GPS provider. We measured the diameter at breast height (DBH) for each tree (n = 113). The height of the trees was not measured in the field. We identified 100% of the trees selected from LiDAR in the field. Moreover, we found a very good match between the tree crowns delineated by low-density pre-fire LiDAR data [53] and the post-fire crown segments ( Figure S1). Finally, we related field-estimated DBH with the LiDAR-estimated tree height and crown area to assess the coherence of the results ( Figure S2). In this study, we have tested various segmentation algorithms (i.e., watershed, Dalponte and Lee using LidR package in R software), and we have ascertained that Silva's tree segmentation [61] algorithm was better than the other ones because it did not merge so many adjacent trees together. Overall, the tree delineation was successful where the trees were more dispersed and thinner, such as in more severely burned areas [13,30,61,64].

LAI and LAD Profiles
Based on vertical profiles of the heights within each crown polygon, we calculated LAD profiles and LAI values, following the methodology proposed by Almeida [47] (Figure 4). The point cloud of each tree was binned into voxels (canopy volume units) with an XYZ resolution fixed at 1 × 1 × 1 m. For each voxel, we calculated the number of pulses that entered and the number of pulses that passed through that voxel. Next, the MacArthur-Horn equation was applied to each voxel to compute its LAD (LADi). From the LAD profiles, the total LAI (the sum of LAD from 1 m height to the top of the tree), the relative understory LAI (percentage of the total LAI from 0 to 3 m) and the crown LAI (total LAI-absolute understory LAI) of each tree were obtained ( Figure 4 and Table 1). Since we only used points with a maximum of 15° off-nadir view, we could work under the assumption that each LiDAR pulse was vertically incident (see Almeida [47] for further details).
Based on the LAD profiles within each crown polygon, we detected inflection points (i.e., break points [BP]) using the breakpoints function from the strucchange package [65] in R, and derived crown base height (CBH, m) ( Figure 4). This step eliminates the influence of stems and understory vegetation from crown LAD profiles. We derived the minimum (out-min) and maximum (out-max) height values around each BP. The minimum height (min-out) of the first BP with LAD values > 0 was considered the minimum CBH (minCBH), whereas the minimum height (min-out) of the BP with the maximum LAD (>p90 LAD) was identified as the maximum CBH (maxCBH). With the height values of minCBH and maxCBH and the maximum height (tree height), we derived crown lengths, and using the crown area (CrArea) we calculated two crown volumes (min and max CrVol, m 3 ), although here we only show results from the max CrVol ( Figure 4).  Based on the LAD profiles within each crown polygon, we detected inflection points (i.e., break points [BP]) using the breakpoints function from the strucchange package [65] in R, and derived crown base height (CBH, m) ( Figure 4). This step eliminates the influence of stems and understory vegetation from crown LAD profiles. We derived the minimum (out-min) and maximum (out-max) height values around each BP. The minimum height (min-out) of the first BP with LAD values > 0 was considered the minimum CBH (minCBH), whereas the minimum height (min-out) of the BP with the maximum LAD (>p90 LAD) was identified as the maximum CBH (maxCBH). With the height values of minCBH and maxCBH and the maximum height (tree height), we derived crown lengths, and using the crown area (CrArea) we calculated two crown volumes (min and max CrVol, m 3 ), although here we only show results from the max CrVol ( Figure 4).

RBR Sentinel 2 Fire Severity Index
Two Sentinel 2 MSI images, captured on 18th July 2017 (pre-fire) and 14th August 2017 (immediately post-fire), were downloaded from the Copernicus open access hub server at the processing level L1C [66]. Correction of the images for surface reflectivity (level L2A) (hereafter 'Sentinel 2A') was performed using the Sen2Cor processing tool in SNAP 6.0 software [67]. All bands were resampled to 10 m pixel resolution, using the reference band B2, and a water and cloud mask, available for each scene in SNAP 6.0, were overlaid [53].We used the Normalized Burn Ratio (NBR) to calculate the Relativized Burnt Ratio (RBR) = (dNBR/(NBR prefire + 1.001)) ( [68]), which was the spectral fire severity index most related to field-based fire severity measures (see [53] for further details). Three fire severity levels were established based on RBR thresholds (i.e., <0.3, 0.3-0.5 and >0.5), following the results obtained by Viedma [53].

Statistical Analysis
To classify fire severity we used the hierarchical clustering on principal components (HCPC) approach, which allows one to combine three standard methods commonly used in multivariate data analyses: (i) principal component analysis (PCA), (ii) hierarchical clustering and (iii) partitioning clustering (i.e., k-means method). The algorithm of the HCPC method was implemented in the FactoMineR package [69] in R software. The first step consisted of PCA calculation, which is very useful for analyzing large data sets with multiple variables (n = 4711 and 49 LiDAR-derived explanatory variables, described in Table 1). PCA was calculated on standardized variables by z scores (mean 0 and standard deviation 1), and we checked that all PCA assumptions were not violated (i.e., continuous observations, linear relationship between all variables [See Figure S2], large enough sample sizes, and no significant outliers). Moreover, we checked two tests to assess the suitability of our data for structure detection (e.g., PCA), as follows: (i) Kaiser-Meyer-Olkin Measure of Sampling Adequacy (KMO) and (ii) Bartlett's test of sphericity (BTS). KMO high values (close to 1.0) and small values (less than 0.05) of the significance level in BTS generally indicate that factor analysis may be useful for the data, whereas a KMO value less than 0.50 and values greater than 0.05 of the significance level of BTS indicate the opposite.
Next, the hierarchical clustering was performed on the PCA results, using Ward's criterion. Later, an initial partitioning was made by cutting the hierarchical tree and selecting the number of clusters following the criterion of maximum statistical distance. Finally, K-means clustering was completed to improve the initial partition obtained from hierarchical clustering. This HCPC approach allowed us to quantitatively identify the variables that most describe each cluster (by means of the v test: a value of the v test greater than 1.96 corresponds to a p-value less than 0.05; the sign of the v test indicates if the mean of the cluster is lower or greater than the overall mean), as well as the principal dimensions that are most associated with the clusters [70]. Later, after identifying the most important variables, the statistical separability among tree clusters was estimated by the non-parametric post hoc test of Kruskal-Wallis.
Finally, to check the accuracy of the clustering, we split the entire database into training (70%, n = 3300) and validation (30%, n = 1411) data sets, and we ran random forest (RF) (using randomForest R package) in the classification mode [69] to predict the categorical tree classes, using the most important variables identified in HCPC analysis. To assess the precision of the RF classification, the Cohen's Kappa coefficient (i.e., level of agreement and the percentage of data that are reliable), sensitivity (i.e., the proportion of positive results out of the number of samples which were actually positive (True Positives, TP)), specificity (i.e., proportion of negative results out of the number of samples which were actually negative (True Negatives, TN)) and accuracy (i.e., how often the classifier is correct, as judged by the equation (TP + TN)/Total) values were calculated [13].

Which Structural Metrics (LiDAR-Derived) Are Most Important for Separating Individual Trees into Distinct Groups?
Before running the HCPC analysis, we checked that all the PCA assumptions were accomplished (see Figure S3). The KMO gave a value of 0.87, which was much above the limit (0.5), and the significance level of BTS was <0.001. Accordingly, we were prevented from continuing with the analysis. Based on the HCPC analysis, we identified three clusters of trees from five significant PCs (eigenvalues > 1) that explained 85% of the total variance ( Figure 5). HCPC classification selected only 14 variables as the most important ones according to the v test (Table S1). The PCA was redrawn using this shortlist of variables ( Figure 5). The performance of RF classification for explaining the separability of the three clusters using the final 14 variables was rather high, as follows: the accuracy was 0.97 ± 0.02; the Cohen's Kappa coefficient was 0.96; the sensitivity (i.e., True Positives) was 0.97 ± 0.02, and the specificity (i.e., True Negatives) was 0.99 ± 0.01 for the validation data set.

Which Structural Metrics (LiDAR-Derived) Are Most Important for Separating Individual Trees into Distinct Groups?
Before running the HCPC analysis, we checked that all the PCA assumptions were accomplished (see Figure S3). The KMO gave a value of 0.87, which was much above the limit (0.5), and the significance level of BTS was <0.001. Accordingly, we were prevented from continuing with the analysis. Based on the HCPC analysis, we identified three clusters of trees from five significant PCs (eigenvalues > 1) that explained 85% of the total variance ( Figure 5). HCPC classification selected only 14 variables as the most important ones according to the v test (Table S1). The PCA was redrawn using this shortlist of variables ( Figure 5). The performance of RF classification for explaining the separability of the three clusters using the final 14 variables was rather high, as follows: the accuracy was 0.97 ± 0.02; the Cohen's Kappa coefficient was 0.96; the sensitivity (i.e., True Positives) was 0.97 ± 0.02, and the specificity (i.e., True Negatives) was 0.99 ± 0.01 for the validation data set. According to the most important structural metrics (LiDAR-derived) for separating individual trees into distinct groups, we ascertained that Cluster 1 (C1, n = 1852) was significantly different from the other two groups, due to its high density of returns <1 m (63%) and the lowest crown LAI (CrLAI), crown volume (CrVol) and crown area (CrArea) (Figure 6 and Table S2). C1 had a high mean treetop height, of ca. 12 m, but a very low average height (ca. 3 m) due to the dominance of vegetation < 1 m. Moreover, C1 showed the highest positively skewed distribution of LiDAR returns (i.e., presence of According to the most important structural metrics (LiDAR-derived) for separating individual trees into distinct groups, we ascertained that Cluster 1 (C1, n = 1852) was significantly different from the other two groups, due to its high density of returns <1 m (63%) and the lowest crown LAI (CrLAI), crown volume (CrVol) and crown area (CrArea) (Figure 6 and Table S2). C1 had a high mean treetop height, of ca. 12 m, but a very low average height (ca. 3 m) due to the dominance of vegetation < 1 m. Moreover, C1 showed the highest positively skewed distribution of LiDAR returns (i.e., presence of tall, residual canopies in trees dominated by short vegetation), and the highest average and standard deviation of intensity LiDAR, mainly due to of the contrast between canopies and bare soils. This heterogeneous vegetation distribution gave place to the lowest evenness (Table S2). Cluster 2 (C2, n = 1422) had a low proportion of returns <1 m (26%), with the highest proportion of returns between 10 and 15 m (42%), the highest CrLAI, CrVol and CrArea, and the lowest understory LAI (<2.5 m) of all three clusters ( Figure 6 and Table S2). C2 was composed of the tallest trees (mean treetop height 13.5 m) with a mean tree height of ca. 8 m. In addition, C2 showed the lowest negative skewed distribution of LiDAR returns, mainly due to the dominance of tall vegetation, and a lower average and standard deviation of intensity than C1 (Table S2). Finally, C3 (n = 1436) was formed by tall shrublands and small tree formations, with the lowest proportion of returns <1 m (19%) and the highest percentage of returns between 2 and 6 m (66%), showing medium CrLAI, CrVol and CrArea and the highest understory LAI of all three clusters ( Figure 6 and Table S2). C3 had the lowest treetop height (c.a. 7 m) with a mean tree height of ca. 3 m. Moreover, C3 showed negative height skewness (i.e., even height distribution), the lowest average and standard deviation of intensity, and the highest evenness (Table S2).
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 22 13.5 m) with a mean tree height of ca. 8 m. In addition, C2 showed the lowest negative skewed distribution of LiDAR returns, mainly due to the dominance of tall vegetation, and a lower average and standard deviation of intensity than C1 (Table S2). Finally, C3 (n = 1436) was formed by tall shrublands and small tree formations, with the lowest proportion of returns <1 m (19%) and the highest percentage of returns between 2 and 6 m (66%), showing medium CrLAI, CrVol and CrArea and the highest understory LAI of all three clusters ( Figure 6 and Table S2). C3 had the lowest treetop height (c.a. 7 m) with a mean tree height of ca. 3 m. Moreover, C3 showed negative height skewness (i.e., even height distribution), the lowest average and standard deviation of intensity, and the highest evenness (Table S2). The behavior of height skewness in relation to tree crown metrics such as CrLAI was outstanding. Height skewness increased sharply in those trees with lower CrLAI, because of the combination of short vegetation and tall, sparse vegetation (i.e., residual crowns). On the contrary, the skewness lowered in healthier trees where the proportion of vegetation layers was more normally or evenly distributed in the vertical tree profiles (Figure 7a). It is likely that the density of returns <1 m was a good indicator of crown status, as it showed a robust negative relationship with CrLAI ( Figure 7b).  The behavior of height skewness in relation to tree crown metrics such as CrLAI was outstanding. Height skewness increased sharply in those trees with lower CrLAI, because of the combination of short vegetation and tall, sparse vegetation (i.e., residual crowns). On the contrary, the skewness lowered in healthier trees where the proportion of vegetation layers was more normally or evenly distributed in the vertical tree profiles (Figure 7a). It is likely that the density of returns <1 m was a good indicator of crown status, as it showed a robust negative relationship with CrLAI (Figure 7b).

How Tree Groups are Distributed within Fire Severity Categories
The tree structures of unburned areas and those that burned under different RBR fire severity levels were rather different (Figure 8a-d). The unburned areas were mainly occupied by C3 (tall shrublands and small trees) (66%) followed by C2 (tallest trees with large and dense crowns) (32%). Low severity areas were occupied by C2 (50%) (the tallest and greatest trees) and C1 (trees with tall and residual crowns dominated by vegetation <1 m) (31%). Moderate and high severity areas were mainly occupied by C1 (81.5% and 91%, respectively) (Figure 8e). C1 was the tree class with the highest RBR, whereas C3 showed the reverse behavior, C2 being in an intermediate position ( Figure  8f). Accordingly, there was a good fit between the RBR fire severity levels and Lidar tree cluster groups. However, as could be expected, there was greater variability in tree morphologies in unburned and low severity areas than in moderate to high ones.

How Tree Groups Are Distributed within Fire Severity Categories
The tree structures of unburned areas and those that burned under different RBR fire severity levels were rather different (Figure 8a-d). The unburned areas were mainly occupied by C3 (tall shrublands and small trees) (66%) followed by C2 (tallest trees with large and dense crowns) (32%). Low severity areas were occupied by C2 (50%) (the tallest and greatest trees) and C1 (trees with tall and residual crowns dominated by vegetation <1 m) (31%). Moderate and high severity areas were mainly occupied by C1 (81.5% and 91%, respectively) (Figure 8e). C1 was the tree class with the highest RBR, whereas C3 showed the reverse behavior, C2 being in an intermediate position (Figure 8f). Accordingly, there was a good fit between the RBR fire severity levels and Lidar tree cluster groups. However, as could be expected, there was greater variability in tree morphologies in unburned and low severity areas than in moderate to high ones.

How Tree Groups are Distributed within Fire Severity Categories
The tree structures of unburned areas and those that burned under different RBR fire severity levels were rather different (Figure 8a-d). The unburned areas were mainly occupied by C3 (tall shrublands and small trees) (66%) followed by C2 (tallest trees with large and dense crowns) (32%). Low severity areas were occupied by C2 (50%) (the tallest and greatest trees) and C1 (trees with tall and residual crowns dominated by vegetation <1 m) (31%). Moderate and high severity areas were mainly occupied by C1 (81.5% and 91%, respectively) (Figure 8e). C1 was the tree class with the highest RBR, whereas C3 showed the reverse behavior, C2 being in an intermediate position ( Figure  8f). Accordingly, there was a good fit between the RBR fire severity levels and Lidar tree cluster groups. However, as could be expected, there was greater variability in tree morphologies in unburned and low severity areas than in moderate to high ones.

How Tree Groups Related to RBR Sentinel 2A Fire Severity Levels
When assessing the difference among the tree groups within each RBR Sentinel 2A fire severity level, we found that all tree clusters showed similar RBR values among them on the unburned and each of the RBR fire severity levels. On the contrary, all clusters showed significant differences among them regarding tree structures in each RBR level ( Figure S4), although the differences tended to be less acute as fire severity increased (Table S3 and Figure 9).  Moreover, when assessing the difference within each tree group across the RBR fire severity levels, we found that several metrics were able to significantly differentiate unburned from burned trees, as well as the different RBR severity levels within each tree group (Table S4 and Figure 10). For all clusters, CrArea, CrLAI, CrVol, understory LAI, as well as the proportion of understory/midstory vegetation (density of returns until 8 m), decreased as fire severity increased, whereas height skewness and density of returns <1 m increased. The loss of understory LAI was the tree property most significantly different between unburned and low severity trees in all clusters. Differences in CrArea, CrLAI and understory/midstory vegetation proportions were more significant for separating unburned trees from those trees less damaged by fire (C2 and C3) (Table S4 and Figure 10). Similar to those metrics, the LAD profiles were by themselves good proxies of fire severity (Figure 11). In this sense, we observed that as fire severity increased, the LAD values at low heights On the other hand, the main tree properties able to separate the different RBR fire severity levels were different for each of the tree clusters. For C1 they were CrLAI, average tree height, height skewness and midstory vegetation (2 m-8 m). For C2 they were CrArea and midstory vegetation (2 m-6 m), and for C3 they were CrArea, CrLAI, intensity (average and standard deviation) and evenness (Table S4).
Similar to those metrics, the LAD profiles were by themselves good proxies of fire severity (Figure 11). In this sense, we observed that as fire severity increased, the LAD values at low heights considerably decreased, due to the consumption and removal of understory and midstory vegetation by fire, allowing us to differentiate unburned from burned trees, as well as crown fires from surface fires.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 22 considerably decreased, due to the consumption and removal of understory and midstory vegetation by fire, allowing us to differentiate unburned from burned trees, as well as crown fires from surface fires. Figure 11. Lidar-derived leaf area density (LAD) profiles of the three structure-derived tree groups stratified by RBR Sentinel 2 fire severity levels with a confidence interval at 95% surrounding the lines. UB: unburned; LS: low severity; MS: moderate severity; and HS: high severity.

Discussion
Here, we have characterized tree structures (i.e., LAI values, crown properties, cloud metrics, vegetation cover distribution and height diversity) after a fire, which reflect well the degree of consumption undergone by the trees during a fire, and thus reflect the fire severity. Consequently, these metrics can affect post-fire regeneration and ecosystem recovery [64,71], allowing a much finer resolution as it descends to the level of a tree. The methodology presented, based on VCPs at the tree level, provides an alternative to traditional area-based approaches based on LiDAR data, or those based on satellite data alone, to identify the heterogeneity of burned trees within standard satellite fire severity levels. Moreover, the object-based nature of individual tree detection can also be advantageous when field data are not available, and can it make many structural measurements without the need for ground-truthing [37]. Recently, researchers have begun to use LiDAR as a primary data source to study forest canopy structures without reference to field data over large areas [11,12,72], demonstrating a potential use for ecological analysis.
The LiDAR metrics derived from VCPs demonstrated promising potential in characterizing finegrained, post-fire residual vegetation, thus providing useful information for post-fire impacts or restoration. The LiDAR system used has a very high pulse density, which allowed the production of vertical foliage profiles for each individual tree with high accuracy. According to Almeida [47], accurate LAD profiles are obtained when the grain size is less than 10 m and the pulse density is high (>15 pulses m 2 ). Here, we used a grain of 1 m and a mean pulse density of 300 pulses m 2 . Moreover, Figure 11. Lidar-derived leaf area density (LAD) profiles of the three structure-derived tree groups stratified by RBR Sentinel 2 fire severity levels with a confidence interval at 95% surrounding the lines. UB: unburned; LS: low severity; MS: moderate severity; and HS: high severity.

Discussion
Here, we have characterized tree structures (i.e., LAI values, crown properties, cloud metrics, vegetation cover distribution and height diversity) after a fire, which reflect well the degree of consumption undergone by the trees during a fire, and thus reflect the fire severity. Consequently, these metrics can affect post-fire regeneration and ecosystem recovery [64,71], allowing a much finer resolution as it descends to the level of a tree. The methodology presented, based on VCPs at the tree level, provides an alternative to traditional area-based approaches based on LiDAR data, or those based on satellite data alone, to identify the heterogeneity of burned trees within standard satellite fire severity levels. Moreover, the object-based nature of individual tree detection can also be advantageous when field data are not available, and can it make many structural measurements without the need for ground-truthing [37]. Recently, researchers have begun to use LiDAR as a primary data source to study forest canopy structures without reference to field data over large areas [11,12,72], demonstrating a potential use for ecological analysis.
The LiDAR metrics derived from VCPs demonstrated promising potential in characterizing fine-grained, post-fire residual vegetation, thus providing useful information for post-fire impacts or restoration. The LiDAR system used has a very high pulse density, which allowed the production of vertical foliage profiles for each individual tree with high accuracy. According to Almeida [47], accurate LAD profiles are obtained when the grain size is less than 10 m and the pulse density is high (>15 pulses m 2 ). Here, we used a grain of 1 m and a mean pulse density of 300 pulses m 2 . Moreover, using LAD profiles from VCPs, several crown metrics closely related to vegetation health (i.e., crown LAI, crown volume or understory LAI) could be obtained (e.g., [24][25][26]29,43]). The ability to calculate CBH and CrVol from breakpoints in VCPs is promising for its efficiency, and follows previous approaches that used inflection points from spline lines to estimate CBH [49], or the parameters of Weibull functions [38] to estimate canopy fuel load.

Which Structural Metrics (LiDAR-Derived) Are Most Important for Separating Individual Trees into Distinct Groups?
Here we show that a classification at the tree level using HCPC and random forests resulted in highly accurate clustering. Our approach distinguished three clusters characterized by different tree structures: the first was mainly composed of trees affected by crown fire; the second was formed by trees that were mainly burned by surface fire, and the third contained trees with tall shrubs and small tree formations that were mainly unburned. The most important variables differentiating the three clusters provided clues about how the ecologically meaningful LiDAR metrics were able to separate the trees by the level of damage (e.g., CrArea, CrVol and CrLAI) [10,32]. Other LiDAR metrics were good indirect indicators of tree damage (e.g., height skewness and density of returns <1 m), due to their negative relationship with crown LAI [18], whereas others were complementary to one another (e.g., average tree height, treetop height and density of vegetation layers), and using them alone could lead to erroneous results. For example, high-severity trees with tall treetops but low average tree heights, mainly due to the very low density of returns at high heights (>10 m) and the high density of returns at low heights (<1 m), would indicate important crown damage.

How Tree Groups Are Distributed within Fire Severity Categories
The three clusters based on tree characteristics were differentiated rather well in terms of their Sentinel 2A RBR metrics. However, the RBR severity levels were not homogeneous in terms of tree structure clustering, and each of these RBR levels had a mixture of tree structures [13]. While broad RBR thresholds were applied to detect different severity levels, a range of crown severities and morphologies was included in those thresholds, contributing to the variability in structural responses observed. We found that open structures with high damaged trees (cluster C1) were strongly dominant within the high fire severely RBR areas, whereas unburned and less severely burned ones showed more complexity in tree morphologies, suggesting that low and moderate severity fires increase heterogeneity [11]. Here, the large variation in tree structures in unburned and low to moderately severely burned areas highlighted the value of monitoring fire-induced damage with LiDAR, rather than optical satellite data [11,12,18,23], due to the spatial and spectral limitations that optical sensors have [20,21], for estimating the fine-grain heterogeneity in fire severity [11,22,23]. In unburned and less severely burned areas, optical satellite data have a limited ability to capture understory and midstory vegetation structures, particularly in high canopy cover forests, whereas LiDAR is more sensitive to these types of change [19].

How Tree Groups Related to RBR Sentinel 2A Fire Severity Levels
All clusters showed significant differences among them regarding tree properties in each RBR level, although as fire severity increased those differences tended to be diluted. On the other hand, within each cluster, increasing fire severity resulted in a significant loss in crown LAI, crown area and crown volume, as well as a large reduction in understory/midstory vegetation, in agreement with Karna et al. [32]. Furthermore, the reduction in understory LAI allowed for separating significantly unburned trees from low severity burned ones in the three tree clusters. However, the tree properties able to significantly separate the different RBR severity levels varied among tree groups, highlighting changes in midstory vegetation for clusters 1 and 2, and changes in crown area for clusters 2 and 3.
Furthermore, height evenness and skewness significantly decreased (higher uneven height distribution) as fire severity increased due to the presence of both residual tall stems and bare soils, as shown by Bolton et al. [18] and Bishop et al. [31], but contrary to the results of Karna et al. [32] who found that crown evenness significantly increased at moderate to high fire severity. In addition, we did not find significant differences in treetop height or mean tree height as fire severity increased [32]. On the other hand, it has been shown that higher LiDAR intensities have been associated with green foliage, while lower LiDAR intensities with non-photosynthetic woody material, including dead vegetation [10,13]. Here, we ascertained that intensity metrics played a secondary role in differentiating tree groups and in the assessment of differences in severity levels. However, we found that LiDAR intensity significantly decreased as fire severity increased in the case of small trees/tall shrubs (C3), although it was almost insensitive in the case of tall trees, as in clusters 1 and 2, maybe because they kept a great part of the top canopy undamaged.

Limitations and Future Work
The main limitations of this study are related to the lack of a larger field sampling of individual trees, and the lack of other measures such as stem and twig size or biomass indicative of the combustion degree. In future work, we propose to use multispectral and thermal unmanned aircraft systems (UAS) and LiDAR data, as well as high-density multi-temporal LiDAR data (pre-and post-fire), to quantifying fire-induced forest changes (e.g., [22,28,31]). This will permit the calculation of the degree of consumption of the plants or soil litter during a fire with very high precision. Whenever possible, during the progression of current forest fires, the UAV LiDAR (LidarPod) should be flown in coordination and under permit of the relevant authorities, during late evening or early morning, when the firefighting helicopters or planes have halted their operations, in order to avoid the risk of collision. This would facilitate the maintaining of unburned patches that will later burn as true controls for fire severity and biomass consumption calculations. Crown volume and biomass losses and post-fire tree morphologies, as reconstructed from LiDAR, will be used to provide comparable metrics of fire severity. Additionally, when the UAV LiDAR cannot be used during fire progression, other LiDAR data available, such as the National Flights available for mainland Spain (LIDAR-PNOA), could be used to characterize pre-fire vegetation structure. Nevertheless, our high-density LiDAR data provided by the device LidarPod are not directly comparable with LIDAR-PNOA, with a very low density of point clouds (0.5-1 point/m 2 ), thus operational challenges remain in obtaining both pre-and post-fire LiDAR metrics with comparable data quality [14]. As the pre-and post-fire LiDAR data had different point densities, the canopy metrics could be affected in complex ways, interacting with differences in slope and elevation due to large differences in the Digital Terrain Models (DTM) derived [31].
By combining measures of LiDAR-derived post-fire tree structure with satellite-derived fire severity levels, variabilities in early fire effects could be characterized with high precision over large forested areas. Our results demonstrate the value of fusing RBR fire severity levels with LiDAR data for contextualizing different tree structures in terms of fire severity, and giving a physical-ecological meaning to reflectance-based spectral indices. These results can assist forest managers in better identifying the areas most severely burned for conservation and restoration plans after wildfires, and for wildlife habitat assessments and other post-fire applications, such as fire's effects on soils.

Conclusions
Damage to canopy structures can provide a sensitive indicator of the degree of fuel consumption by fire, and thus of fire severity. Here, we made a comprehensive use of the vertical canopy profiles from high-density LiDAR data to obtain several metrics related to biophysical single-tree characteristics, such as crown properties, LAI values, vegetation layer distribution and diversity. Moreover, using inflection points to derive the canopy base height, crown volume could be estimated. All these metrics at the tree-level allowed the distinguishing of three clusters of tree morphologies, which were separated rather well among each other. These tree structure-derived clusters were unevenly distributed in the different RBR Sentinel 2 fire severity levels. Unburned and low severity burned areas were more diverse than moderate and high severity burned ones. Furthermore, LiDAR metrics were able to distinguish crown fire from surface fire through changes in the understory LAI, and the understory and midstory vegetation. Moreover, the LiDAR metrics indicated that vegetation layers' distribution within a tree, and crown properties, were important for estimating the impact of fire on single trees, a feature that was not captured by the RBR Sentinel 2 index. Although RBR was not able to distinguish tree morphologies, especially in unburned and low severity areas, due to their inability to penetrate the forest canopy, we could ascertain that the RBR Sentinel 2 index is an adequate proxy of fire severity. The three tree structure clusters were rather well differentiated spatially and spectrally by the RBR fire severity levels, and significant intra-cluster differences in terms of LiDAR metrics were observed among RBR fire severity levels. When high-resolution LiDAR data, as we used here, is available, it provides a high spatial resolution, at the level of a single tree, which is needed for ecological impact studies. Tree-based fire severity maps at this level of detail could thus be used to develop more precise, site-specific fire impact studies and post-fire management plans.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/21/3554/s1, Figure S1: Snapshots of tree delineation validation in the field, Figure S2: Scatterplots relating tree top height (a) and crown area (b) derived from LiDAR data with diameter at breast height (DBH) of individual trees measured at ground-field, Figure S3: Scatterplots relating pairs of variables included in the Principal Component Analysis (PCA) and correlation matrix of the same set of variables, Figure S4: Normalized height point clouds of individual trees classified by using LiDAR derived metrics (C1-C3) and stratified by RBR Sentinel 2 fire severity levels, Table S1: Scores given by the v test derived from the PCA classification of individual trees and descriptive statistics. Here we show the most important variables (v test value > 10), Table S2: Mean values of the most significant variables used to classify individual trees (C1-C3), Table S3: Mean values of the most significant variables used to classify individual trees (C1-C3) grouped within each RBR fire severity levels to assess differences among tree groups (inter-cluster variability), Table S4: Mean values of the most significant variables used to classify individual trees (C1-C3) stratified by RBR fire severity levels (UB: unburned; LS: low severity; MS: moderate severity and HS: high severity) to assess differences within groups (intra-cluster variability).