Analysis of Mediterranean Vegetation Fuel Type Changes Using Multitemporal LiDAR

: Increasing ﬁre size and severity over the last few decades requires new techniques to accurately assess canopy fuel conditions and change over larger areas. This article presents an analysis on vegetation changes by mapping fuel types (FT) based on conditional rules according to the Prometheus classiﬁcation system, which typiﬁes the vertical proﬁle of vegetation cover for fuel management and ecological purposes. Using multi-temporal LiDAR from the open-access Spanish national surveying program, we selected a 400 ha area of interest, which was surveyed in 2010 and 2016 with scan densities of 0.5 and 2 pulses · m − 2 , respectively. FTs were determined from the distribution of LiDAR heights over an area, using grids with a cell size of 20 × 20 m. To validate the classiﬁcation method, we used a stratiﬁed random sampling without replacement of 15 cells per FT and made an independent visual assessment of FT. The overall accuracy obtained was 81.26% with a Kappa coefﬁcient of 0.73. In addition, the relationships among different stand structures and ecological factors such as topographic aspect and forest vegetation cover types were analyzed. Our classiﬁcation algorithm revealed that stands lacking understory vegetation usually appeared in shady slopes, which were mainly covered by beech stands, whereas sunny areas were preferentially covered by oak stands, where the understory reached greater height thanks to more light availability. Our analysis on FT changes during that 6 year time span revealed potentially hazardous transitions from cleared forests towards a vertical continuum of canopy fuels, where wildﬁre events would potentially reach tree crowns, especially in oak forests and southern slopes with higher sun exposure for lower fuel moistures and increased ﬂammability. Accurate methods to characterize forest canopy fuels and change over time can help direct forest management activities to priority areas with greater ﬁre hazard. Multi-date canopy fuel information indicated that while some forest types experienced a growth of the shrub layer, others presented an understory decrease. On the other hand, loss of understory was more frequently detected in beech stands; thus, those forests place lower risk of wildﬁre spread. Our approach was developed using low-density and publicly available datasets and was based on direct canopy fuel measurements from multi-return LiDAR data that can be accurately translated and mapped according to standard fuel type categories that are familiar to land managers.


Introduction
Longitudinal studies using multitemporal remote sensing series can be very helpful for the purpose of monitoring and analyzing the dynamics of ecosystem structural features [1]. Traditionally, most remote sensing-assisted studies on the vegetation dynamics were based on spectral sensors. For instance, vegetation phenology has been studied using Moderate Resolution Imaging Spectroradiometer (MODIS) [2], land cover with Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper (ETM) and Sentinel-2 [3], fire risks with Landsat or National Oceanic and Atmospheric Administration-Advanced Very High is defined as "an identifiable association of fuel elements of distinctive species, form, size arrangement, and continuity that will exhibit characteristics fire behavior under defined burning conditions" [49]. Particularly in Mediterranean regions, the Prometheus classification system [50] has been widely used, which is an adaptation of the Northern Forest Fire Laboratory system to the Mediterranean conditions [40]. According to the Prometheus classification system, there are seven different FTs (Table 1). Furthermore, this FT classification is an easy way of typifying the vertical profile of the vegetation cover. The capacity of remote sensing to reliably classify and map FTs has been tested using LiDAR [51], and also in combination with multispectral imagery, such as QuickBird [52], airborne imagery [53], Landsat-8 OLI [54,55], or ASTER [56]. However, no research has focused on change detection and analyzing FT changes using exclusively multi-date discrete-return LiDAR, even though these types of LiDAR data alone have been proven sufficiently accurate to characterize forest successional stages across large spatial extents [57]. The aims of this research are (1) to compute LiDAR metrics to measure and categorize forest and shrubland conditions into FT categories according to the Prometheus classification system; (2) to assess the FT classification accuracy; (3) to assess the relationships in the spatial distribution of the FTs with ecological factors; and (4) to determine FT change over time from the viewpoint of the forest stand vertical profile, characterized as FTs, using LiDAR data.

Study Area
The study area was a 2 × 2 km tile (southwest corner coordinates in UTM: 504,000; 4,660,000) located in La Rioja (Northern Spain, Figure 1). This area was selected due to the availability of data both from the first (2010) and second (2016) PNOA LiDAR coverages. In addition, the vegetation types comprised the full range of seven different FT from the Prometheus classification (Table 1). No forest fires occurred in the study area between 2010 and 2016. In terms of climate, the area selected is mainly continental-Mediterranean (mean annual rainfall: 853 mm and mean annual temperatures: 7.3 • C [58]). Given its topography and characteristic geographical situation, it hosts a great diversity of microclimates and ecosystems. The vegetation is mainly composed of meso-xerophilous oak, beech, and pine forests. The most representative species are Quercus ilex ssp. ballota, Quercus pyrenaica, Quercus petraea, Fagus sylvatica, and Pinus sylvestris [59]. The forest vegetation cover types and the orthophoto of the study area are shown in Figure 1.

LiDAR Data Acquisition and Processing
The LiDAR data corresponding to years 2010 and 2016 were downloaded from the Spanish Geographic Institute's website [60]. Scan densities calculated from flight parameters [61] were 0.5 pulses·m −2 for the 2010 survey and 2 pulses·m −2 for 2016. Both datasets were collected in August, hence under leaf-on conditions, with high photosynthetic productivity and most favorable conditions for the occurrence of wildfires. Both datasets had the same angle error specifications (≤0.005 • for Roll and Pitch and 0.008 • for Yaw) but different field of view (FOV): 35 • in the 2010 dataset and 26 • in the dataset from 2016. In addition, both LiDAR datasets could detect a maximum of 4 returns for each pulse. Data processing was carried out as a combination of a code developed in RStudio v.1.3.1093 [62] and FUSION software v4.10 (Forest Service U.S.D.A, Seattle, WA, USA) [63]. The original LiDAR data were filtered, and ground pulses were classified and interpolated into a digital terrain model (DTM) at 5 m × 5 m spatial resolution and then used to calculate heights above ground of individual LiDAR returns. For a 20 × 20 m spatial resolution grid covering the study area, we calculated LiDAR return relative densities at the following strata: 0-0.3 m; 0.3-0.6 m; 0.60-2.0 m; 2.0-4.0 m; and above 4.0 m, coinciding with the strata of the Prometheus classification system (Table 1). Then, for each cell we calculated the stratum with highest (Mode), second highest (2nd Mode) and third highest (3rd Mode) number of LiDAR returns; the maximum elevation of vegetation (Max.Elev); the forest canopy cover (Cover); and the proportion of ground returns (Ground). All returns were used for calculating these metrics, except for the case of Cover, which was estimated from the percentage of first returns above 4.0 m.

Fuel Type Classification and Mapping
Using the LiDAR height and the cover metrics, an FT was assigned to each cell based on conditional rules according to the Prometheus classification system, following the criteria in Table 2. Briefly, the general procedure was based on searching for the stratum with the highest frequency of LiDAR returns (Mode), and whenever a mode was found at the uppermost strata, then the next mode was searched for, up to the 3rd Mode in the case of forest FTs (FT5, FT6, or FT7). First of all, if Ground proportion was greater than 60%, the cell was directly assigned to pastures (FT1). Then, if Cover was lower than 50%, the cell was classified as either pastures or shrub FTs (FT1, FT2, FT3, or FT4), defining the exact FT according to the stratum that concentrated the highest number of returns (Mode). In case that the stratum with the greatest number of returns was the stratum above 4.0 m high, which would be due to the presence of few scattered trees, the second height interval with more returns (2nd Mode) determined the associated FT (Table 2). A similar procedure was employed when Cover was greater than or equal to 50%, a criterion that was associated to forest FTs (FT5, FT6, or FT7). A similar procedure was then followed, recursively looking at the 1st, 2nd and 3rd Mode, allocating to the suitable FT according to the abundance of vegetation material across strata ( Table 2). For forest areas with 1st or 2nd Mode at 0.6 m or above, an additional criterion was introduced to discriminate whether vertical continuity of plant material would allow ground fires to spread toward tree crowns. We empirically determined that having already established that there are both trees and understory in these areas, a Max.Elev above 12 m would ensure no vertical continuity because shrubs would not reach to crown base height. Thus, they were either forest with shrub layer (FT6) or forest with vertical continuity (FT7), depending on Max.Elev ( Table 2). The study area was divided into 20 × 20 m regular grid cells, and then we applied to each cell this set of recursive rules using RStudio (Boston, MA, USA) [62] so each cell was assigned to a FT. The same procedure was repeated for both the LiDAR files from 2010 and 2016, and then we studied the changes of FTs by comparing them.

Classification Accuracy
We validated rules sets and FT classification accuracy by randomly selecting 15 grid cells for each FT and year. There was only one exception in the case of low shrubs (FT2) in 2016, for which only 13 cells were available. The LiDAR returns at the location of those selected cells was extracted and randomly arranged to avoid any subjective biases. Then, the FUSION LDV (LiDAR Data Viewer; [63]) was used to assign FT classes to each cell by observing the relative position of vegetation features displayed on this 3D visualization environment. The validation was done by contrasting the rule-based classifications carried out in Table 2 against this visual classification, which was considered as the reference. Accuracy assessments of the classifications were performed using a confusion matrix, and accuracy measurements were summarized using overall accuracy, user's accuracy, producer's accuracy, and Kappa coefficient [64]. The results of the confusion matrix were weighted according to the proportion of area observed for each FT [65,66] (Equation (1)). The weighted proportions (p ij ) of the sample which were observed for FT class j and predicted as FT class i were where A j /A t is the ratio between the area (A j ) observed for each FT class j with respect to the total number of cells (A t = 10,000); n ij is the number of cells observed for class j and predicted to be class i, and n i· is the total number of plots validated for a class i.

Relationships to Ecological Factors
We analyzed the relationships between the presence of the different FTs and ecological factors, such as the topographic aspect and forest vegetation cover types ( Figure 1). We also considered a number of additional variables, meteorological data for instance, which are not detailed herein because no relevant relationships to FT changes were observed. Regarding topographic aspect as an ecological factor, we were interested in whether a given FT occurs more predominantly in sunny or shady areas; thus, this variable was converted into a categorical factor including either of those values. We classified the pixels of the DTM with aspect between 120 and 300 • as sunny, and else (<120 • or >300 • ) as shady (the north being at 0 • ), using the classification criterion of [67]. Unforested vegetation cover types in Figure 1 were irrelevant to our analysis; thus, we considered only the remaining as forest vegetation cover types: beech/pine/oak-dominated, or mixed deciduous forests. Then, we performed a three-factor log-linear analysis and a goodness of fit on whether the observed frequencies significantly differed from the expected frequencies. The three factors were FT as the dependent variable (7 FT of the Prometheus classification system), aspect (shady or sunny), and forest vegetation cover type (beech, pine, oak, mixed deciduous forest). Since we considered frequencies over combinations of categorical variables, statistical significance was tested using a Pearson chi-square analysis with a 0.05 significance level. We excluded the mixed deciduous category from the analysis to comply with the chi-square restriction that forbids more than 20% of the cells with frequencies lower than 5 [68]. The ratio of the log-linear parameter estimate to its standard error was used to obtain the frequency significance level (0.05), using SPSS Statistics v.26.

Fuel Type Mapping
Confusion matrices from mapped LiDAR FTs for 2010 and 2016 ( Figure 2, Tables 3  and 4) showed nearly identical accuracy results for both years, just slightly lower for the 2010 dataset, which had a lower scan density. We observed similar patterns at both years, with higher degrees of confusion between forest with shrub layer (FT6) and forest with vertical continuity (FT7), i.e., on whether there is vertical continuity of plant material in forests. Overall accuracies were 80.72% and 81.26%, and the Kappa coefficients were 0.73 and 0.73, for the years 2010 and 2016, respectively.
For the 2010 dataset (density 0.5 pulses·m −2 ) producer's accuracies were high, and Table 3 shows that most of the error was concentrated among the classes forest with shrubs (FT6), forest with vertical continuity (FT7), and low shrubs (FT2). The lowest values of user's accuracy were found in low shrubs (FT2) and forest with vertical continuity (FT7). For the remaining FTs, the values of weighted user's accuracy were high, reaching up to 0.87.
For the 2016 dataset (density 2 pulses·m −2 ) producer's accuracies were also high, except for forest with shrubs (FT6), which only reached 0.19. The highest levels of error were again found in the discrimination between forest with shrub layer (FT6) and forest with vertical continuity (FT7).

Fuel Type Changes
When analyzing the FT changes from 2010 to 2016, it was observed that in 74.0% of the total area the FT class remained unchanged between dates (Figure 3), being forest without understory (FT5) and pastures (FT1) the most stable FT categories with 90.7% and 78.3%, respectively (Table 5). Otherwise, most FTs experienced important changes, with interesting insights on their transition. A moderately high percentage of vegetation growth was observed, as 13.7% of the surface changed from a lower vegetation stratum to a higher vegetation stratum or even developed a tree stratum (denoted as 'Growth' in Figure 3a). Most changes toward higher vegetation strata occurred in oak forests, and also at the edge of other forested areas, suggesting ecotone changes and forest expansion toward pastures and shrublands (Figure 3). This vegetation growth was also observed as an area increase for medium shrubs (FT3; 198.3%) and forest with vertical continuity (FT7; 230.7%) ( Table 5). The increment in area of medium shrubs (FT3) showed shrub encroachment changes in areas that were formerly pasture (FT1), although we also observed transitions from low shrubs (FT2) toward pastures (FT1) ( Table 5). In addition, we observed a decrease from high shrubs (FT4) toward medium shrubs (FT3).   Regarding forested areas, we observed that increases in forest without understory (FT5; 110.4%) and forest with vertical continuity (FT7; 230.7%) were associated to decreases in forest with shrubs (FT6; 39.5%) and high shrubs (FT4; 63.9%), respectively. In addition, some forest without understory (FT5) changed to forest with shrubs (FT6) or to forest with vertical continuity (FT7), indicating that understory growth took place more frequently in the oak forests of the South West (Figure 3b; Table 5). On the contrary, a reduction in the understory height was observed in 11.79% of the study area (Figure 3b). This decrease principally occurred not only in transitions from forest with vertical continuity (FT7) to forest without vertical continuity (FT5-6), but also from forest with shrubs (FT6) toward forest without understory (FT5), which were detected over a wide area of beech forests, mainly situated in the Eastern part of the study area (shown in pale purple color in Figure 3b). Forest loss was recorded in a very small proportion, as only 0.48% of the FT changed from forest FTs (FT5-7) to non-forested FTs (FT1-4) (Figure 3a). Figure 4a shows the vertical profile of LiDAR returns within a given a cell, selected as it exemplifies this decrease in the understory presence in FTs corresponding to forests, which could be remarkably detected, especially given the increase in scan density.

Ecological Factors
We found relationships in the spatial distribution of the FTs with topographic aspect and forest vegetation cover types (triple interaction statistically significant at p-value < 0.001). Table 6 displays how these factors had an influence on the presence of FTs as detected in the 2016 classification, whereas Table 7 shows how they affected the observed changes. Sunny and shady aspects covered 38.9% and 61.1% of the study area, respectively. A total of 58.4% of the study area was forested (FT5, FT6, and FT7) in 2010 whereas in 2016 the forested area grew to a total of 61.4%. As a sciophilous species, beech had a vertical structure of forest without understory (FT5) and preferentially appeared in shady slopes (58.0%; p-value < 0.05) over sunny areas (9.5%; p-value < 0.05) ( Table 6). This was also the case for pine forests, with most areas being of forest without understory (FT5) in shady areas (6.2%; p-value < 0.05). In contrast, oak forests did not seem to be linked with any particular FT, with most of them being forest without understory (FT5; 9.4%), followed by areas of forest with vertical continuity (FT7; 4.3%), and only few being forest with shrubs (FT6; 1.6%), with statistical significances clearly showing their predominance for sunny areas with higher insolation (Table 6). On the other hand, mixed deciduous stands were indifferent to aspect and preferentially had a forest without understory (FT5) structure (Table 6). In addition, it was observed that 57.7% of the non-forest FTs (FT1-4) was located in sunny areas, also displaying a heliophilous preference.  Table 6. Observed frequencies of the following factors: forest fuel type (FT; FT5, FT6, and FT7), forest vegetation cover type (beech, pine, oak and mixed deciduous) and aspect (sunny and shady) in 2016 in the study area. Statistically significant frequencies (p < 0.05) have been identified with an asterisk (*). The analysis of changes in vertical stand structure between 2010 and 2016 revealed processes of understory loss in beech forests, but also transitions from cleared forests towards a vertical continuum in plant fuel in sunny slopes that can create greater fire hazard ( Table 7). The main transition was observed as a loss of understory in beech stands, with a very large patch, which changed from forest with shrubs (FT6) to forest without understory (FT5) (p-value < 0.05; Table 7). This transition was clearly visible when observing the vertical profile of LiDAR returns in those areas (Figure 4a), which showed a loss in the understory despite the increase in the scan density. In contrast, oak stands developed a continuous vegetation profile, in a great proportion changing from being forest with no understory (FT5) toward forest with vertical continuity (FT7; Figure 4b), preferentially on sunny slopes (11.7%; p-value < 0.05). The transition of oak forests from forest without understory (FT5) toward forest with shrubs (FT6) in sunny slopes also was significant (5.4%; p-value < 0.05). On the other hand, pine forests covered a relatively small proportion of the study area, compared to beech and oak forests, and there were no significant changes. Mixed deciduous forest also experienced little or no FT change between dates (Table 7). Table 7. Changes among structural stand and fuel types (FT5: forest without understory; FT6: forest with shrubs; FT7: forest with vertical continuity), forest vegetation cover types (beech, oak, pine and mixed deciduous forest) and aspect (sunny and shady). Statistically significant frequencies (p < 0.05) have been identified with an asterisk (*).

Discussion
We developed a methodology to map FTs according to the Prometheus classification system using publicly available LiDAR data. It could be deduced from the 2010 confusion matrix ( Table 3) that a greater error occurred in categories with shrub stratum that may be harder to distinguish, i.e., forest with shrubs (FT6) or forest with vertical continuity (FT7). For the 2016 dataset (Table 4), still forest with understory (FT6, FT7) were among the FT most difficult to distinguish due to their structural complexity. Low shrubs (FT2) and forest with vertical continuity (FT7) were the FT with the lowest producer's accuracy in the 2010 dataset because they were rare in the study area given that only 0.93% and 1.66% of the study was classified as low shrubs (FT2) and forest with vertical continuity (FT7), respectively. The errors between low shrubs (FT2) and medium shrubs (FT3) may be the result of an error in the visual classification because in case of doubt, the more restrictive FT was chosen. In the case of the errors between low shrubs (FT2) and forest with shrubs (FT6), the source of error might be a Cover underestimation. Other sources or error, apart from the ones coming from a visual classification error, might be the terrain slope of our mountainous study area, and the different FOV for each LiDAR flight, which was greater in the 2010 dataset and, therefore, allowed a better shrub detection than for the 2016 dataset. That may explain why an increase in the scan density did not result in an increase in the overall accuracy since the higher scan angle may have compensated the lower density in the 2010 dataset.
In change analyses carried out post-classification, error matrices such as those in Tables 3 and 4 are a common endpoint of accuracy assessment [65]. However, Olofsson et al. [65] showed that these single date error matrices do not provide information relevant to assess the accuracy in change identification. Thus, the matrix in Table 5 is a summary of the vegetation structural changes using LiDAR between 2010 and 2016 in an area that contains the FTs typically found in Mediterranean ecosystems. Our analysis, however, does not directly evaluate the capacity of this method to detect changes. We nonetheless have an indication of potential to detect vegetation changes in a relatively short period (6 years). Analyzing vegetation changes using Prometheus FTs revealed that, while most of the area persisted as stable forest types, especially forest without understory (FT5) and pastures (FT1), many other areas showed significant changes between 2010 and 2016, hence the relevance of analyzing the changes of FTs. We observed changes of ecosystem development due to vegetation growth from lower to higher vegetation strata (13.7%; Figure 3; Table 5), sometimes even reaching the tree stratum with processes of forest expansion occurring at edges of forested vegetation types. These changes could be experienced not only because of a vertical growth, but also because of a growth in width so the canopy cover increases, and it surpasses the threshold of 50% of canopy cover to be considered as a forested FT. Other reasons could be that the increase in scan density allowed us to make a more precise estimation of the canopy cover, or it can also be just an error in the classification. Understory growth was frequently observed in the interior of oak forests (Figure 4; Table 7), and also at the border of pine and beech forests as ecotone changes. Special attention should be paid to these areas, where the fuel hazard has considerably increased over the time span between the two datasets, due to the growth of the shrub stratum that may lead to vertical continuity and crown fires [69]. Only few of the changes observed were due to forest loss (0.48%), for instance changes from forested FTs to pastures (FT1) or shrubs (FT2-4), which were mainly observed on pine forests (Figure 3a). These cells classified as forest loss may be interpreted as an error in the classification given that there are not forming patchy structures and because of the few number of cells classified as forest loss. The 11.8% decrease in understory proportion (grey areas in Figure 3a), which was observed in some beech areas, may be a result of the closure of the tree canopies inhibiting the growth of the understory (Table 7). Regarding the transitions from shrub FTs (FT4, FT3, and FT2) to pastures (FT1), they might be the result of grazing, droughts that could affect shrub density, or survival or large-sized herbaceous vegetation, which the classification algorithm interprets as small-sized shrubs.
We observed a clear preference of beech and pine for shady slopes (61.0% and 6.7%, resp.) as opposed to sunny areas (11.2% and 2.4%, resp.). On the contrary, oak forests were more abundant in sun-exposed areas (11.9%) as compared to shady slopes (3.4%), given their heliophilous nature. Regarding the vertical structure, forests with no understory (FT5) preferentially appeared on shady slopes (68.7%) over well-illuminated areas (19.8%; Table 6). Forests with low understory (FT6) appeared in similar proportions in either shady or sunny areas, while forest with vertical continuity (FT7) were more frequent in sunny slopes. Both FT6 and FT7 (11.4% together) were mainly represented by oak stands (3.9%). We also observed that changes toward FTs with risk of wildfires reaching tree crowns (FT7) occurred with higher frequency at these southern slopes with more sun exposure. Thus, competition for the light resource of these tree species and accompanying understory drives forest change toward the development of a more complex forest structure, especially in oak forests (Table 7). We thus recommend fire prevention measures to prioritize oak forests and southern slopes over shady areas, as our results showed FT changes to be clearly determined by terrain aspect.
Our study combines the following characteristics, which have not been present at the same time in previous studies. Firstly, it is important to highlight the simplicity of the method, which is based on conditional rules that may be tailored to specific vegetation types, as opposed to statistical adjustment. Indeed, LiDAR provides a direct measurement of vertical and horizontal fuel structure to which classification rules can be applied. Multireturn information provides data on the distribution of fuels from the ground up, which is not available from optical sensors. Secondly, the availability and characteristics of the data, as they were publicly available datasets, saves economic resources and yields good-quality results. Lastly, there is high overall accuracy for all the different FT of the Prometheus system, including the FTs which are more difficult to detect using remote sensors, which are forests with understory (FT6 and FT7). In summary, the major contribution of our work is that our FT classification algorithm efficiently detects vegetation changes, especially in the understory, over time. Authors in [51] employed an endmember approach applied also to a LiDAR dataset from the PNOA project (scan density of 0.5 pulses m −2 ) to map seven and five different FTs in a National Park area of 41,000 ha, obtaining a 39% and 55% of overall accuracy, respectively. Our method based on conditional rules applied to LiDAR metrics describing canopy and sub-canopy vegetation structure provided a 79% overall accuracy, although in a smaller study area (400 ha). Studies discriminating lesser number of classes are also essentially prone to show better overall results. Other studies have more typically employed error-minimizing statistical adjustments to training data, subsequently obtaining higher accuracies than studies like this contribution. Studies discriminating lesser number of classes are also essentially prone to show better overall results [51,70]. Authors in [71] also used small-footprint LiDAR with a nominal pulse spacing of 2.0 m to distinguish between only two different structure classes (single-story and multi-story) at regional scales, obtaining an overall classification accuracy of 97%. The results provided by [57] accomplished an overall accuracy of 90% distinguishing between four different stages of forest stand development using exclusively LiDAR data with a scan density of 10 pulses m −2 . Moreover, authors in [71] used only LiDAR data with scan density of 0.26 pulses m −2 to classify six and seven different forest successional stages with an overall accuracy of 95.54% and 90.12%, respectively, by means of a random forest algorithm. Authors in [39] integrated LiDAR and hyperspectral imagery to classify ten different structurally based vegetation type classes using a combination of image segmentation, principal components analysis, and unsupervised classification, obtaining an overall accuracy of 94%.
On the other hand, some authors approached the classification of FT by means of spectral sensors or combinations of data from both active and passive sensors, such as [72], who used a Quickbird image and ancillary data, with 75% overall accuracy and 0.69 Kappa coefficient, but with the particularity that no forest with shrubs (FT6) were present in their study area. In addition, [11,40] obtained slightly higher accuracies (81.82% and 82.8%, respectively), but it should be noted that they used far more complex methodologies than ours. Authors in [40] used a higher-scan density (between 1.5 and 6 pulses m −2 ), plus an Airborne Thematic Mapper multi-spectral sensor with 11 different bands, which they used to calculate Normalized Difference Vegetation Index and other spectral indices. Their results were similar to those obtained by [11], who used Landsat TM images and also involved ancillary data. Additionally, authors in [73] obtained similar results with a 3-dimensional discrete anisotropic radiative transfer simulation and multiple endmember spectral mixture analysis and spectral angle mapper applied to LiDAR-based FT classification. Our methodology produced similar results with simpler conditional rules, but further research could explore how these could be improved by tailoring the rules to specific vegetation types, using information from spectral sensors similar to those employed in these studies.

Conclusions
This cost-effective methodology has proven that LiDAR provides direct measurements which, in most cases, showed good agreement between calculated and visually assessed FT categories. As the demand for higher-quality LiDAR increases, so will opportunities for enhancing wildland fuel characterization and classification. Secondly, it has shown that multi-date LiDAR can provide insight and information with an adequate accuracy of changes in forest canopy fuels, particularly where vertical fuel structures develop and may increase fire hazard and the likelihood of a crown fire. Furthermore, the classification of vertical vegetation structures as FT and mapping methodology presented here has some advantages: it is a simple, inexpensive, and reliable application of the Prometheus classification system for wall-to-wall FT mapping of wildlands in countries with public LIDAR coverage. Moreover, this methodology could be applied in areas where the accessibility does not allow to collect field data and can also help to improve fuel management efforts where they are most needed, particularly where lives, critical wildlife habitat, or human infrastructure are at risk.

Data Availability Statement:
The data presented in this study are openly available in http:// centrodedescargas.cnig.es/CentroDescargas/index.jsp (accessed on 30 May 2020).