Fuel Type Classiﬁcation Using Airborne Laser Scanning and Sentinel 2 Data in Mediterranean Forest A ﬀ ected by Wildﬁres

: Mediterranean forests are recurrently a ﬀ ected by ﬁre. The recurrence of ﬁre in such environments and the number and severity of previous ﬁre events are directly related to ﬁre risk. Fuel type classiﬁcation is crucial for estimating ignition and ﬁre propagation for sustainable forest management of these wildﬁre prone environments. The aim of this study is to classify fuel types according to Prometheus classiﬁcation using low-density Airborne Laser Scanner (ALS) data, Sentinel 2 data, and 136 ﬁeld plots used as ground-truth. The study encompassed three di ﬀ erent Mediterranean forests dominated by pines ( Pinus halepensis , P. pinaster y P. nigra ), oaks ( Quercus ilex ) and quercus ( Q. faginea ) in areas a ﬀ ected by wildﬁres in 1994 and their surroundings. Two metric selection approaches and two non-parametric classiﬁcation methods with variants were compared to classify fuel types. The best-ﬁtted classiﬁcation model was obtained using Support Vector Machine method with radial kernel. The model includes three ALS and one Sentinel-2 metrics: the 25th percentile of returns height, the percentage of all returns above mean, rumple structural diversity index and NDVI. The overall accuracy of the model after validation was 59%. The combination of data from active and passive remote sensing sensors as well as the use of adapted structural diversity indices derived from ALS data improved accuracy classiﬁcation. This approach demonstrates its value for mapping fuel type spatial patterns at a regional scale under di ﬀ erent heterogeneous and topographically complex Mediterranean forests. ( Elev.min, Elev.max, Elev.mean, Elev.mode ( Elev. SQRT Elev. L Elev. L1, Elev. Elev. height variability metrics include deviation ( Elev. ( Elev.variance coe ﬃ cient of ( Elev.CV interquartile ( Elev.IQ ( ( Elev.kurtosis ret. returns ret.)


Introduction
Wildfires constitute a recurrent hazard in forested environments around the world, shattering many environmental and socio-economic resources each year. Fires shape the current landscape in the Mediterranean basin [1] modifying forest structure, composition and species richness [2]. An average of 45,000 fires is recorded yearly [3] in the Mediterranean region, with Spain being the second most affected country after Portugal [4]. For example, in 1994, 327,696 ha were burned due to an extreme dry summer and high temperatures, and was one of the most affected years reported by the Spanish fire database [5]. Forest fire risk mitigation is under the spotlight due to uncertain effects of climate and Sánchez et al. [42]. Some studies have explored the indirect fuel type mapping approach integrating ALS data with high-resolution multispectral images [8,29,43], with middle resolution passive sensors [44,45], its combination with SAR data [46] or ALS data alone [12], obtaining better results those that combine multi source datasets. The use of ALS data is particularly relevant for characterizing the vertical structure of fuel types as well as to discriminate fuel discontinuities [45]. Studies that have mapped fuel type have selected either a suite of statistical metrics commonly used for forestry applications, related to height distribution and canopy cover at different heights [44] or applying the height bin approach [8]. The height bin approach computes the proportion of ALS returns within different height ranges creating a 3D vegetation vertical profile. The use of structural metrics proposed to assess forest structure complexity such as rumple [47] or forest structural diversity as LiDAR height diversity index (LHDI), and LiDAR height evenness index (LHEI) [48], might help to create parsimonious models. Concerning the training sample in fuel type modelling using ALS data, two approaches can be found in previous studies; regional studies mostly used previous mapping to train the model [45,46] and more local ones used field plots as sample, usually from low to medium size, between 50 and 100 field plots [8,43,44].
Despite fuel mapping in Mediterranean ecosystems have been previously addressed [28,29], there is a lack of research on the integration of medium-resolution multispectral data with low point density ALS data for fuel types mapping at regional scale. To the best of our knowledge, this is the first attempt to classify fuel types at regional scale using low density ALS data, multispectral Sentinel 2 data, and field work that compares the performance of classification using both datasets together or separately. Sentinel 2 data have been used to characterize canopy fuel characteristics [23], fuel types [49] or mapping wildfire ignition [42], while the integration with ALS for fuel mapping at regional scale requires further analysis. Furthermore, the research analyzes whether differences in the presence of fuel types exists between areas previously burned in 1994 or unburned ones. Different authors have studied the regrowth of fire vegetation after fire [1,50], but the analysis of fuel type composition after fire and its comparison with unburned areas remains unknown in Mediterranean environments, characterized by a high topographic heterogeneity, vegetation species and climate. The study is developed within the SERGISAT project framework, whose aim is the analysis of the effects of fire severity and environmental variables on regeneration of areas affected by large forest fires in Spain, considering fraction cover, biomass and biodiversity.
In this context, the aim of this study was to classify and map fuel types according to Prometheus classification using low density ALS data, Sentinel-2A data and field work within three different Mediterranean forests dominated by pines (Pinus halepensis, P. pinaster y P. nigra), oaks (Quercus ilex) and quercus (Q. faginea). Furthermore, the following secondary objectives are addressed: (i) compare classification performance when combining ALS metrics with Sentinel-2 or when using these data sources separately; (ii) characterize fuel type spatial patterns under areas affected by a previous wildfire with high structural heterogeneity, topographical complexity, and different species representative of the Mediterranean region to create an integrated model for the Spanish Mediterranean arch; (iii) compare fuel type presence between burned and un-burned areas; (iv) analyze the use of synthetic ALS derived metrics to integrate structural complexity and diversity for obtaining parsimonious classifications, and depicting its importance in model performance; (v) compare two metric selection approaches and two non-parametric classification methods.

Study Area
The forests under study are located in areas affected by three wildfires as well as their surroundings.   Montmajor forested site includes 43,774.45 ha affected by fire, and 41,868.91 ha within the surrounding area. Forests are dominated by Quercus faginea, Quercus pubescens and Pinus nigra, accompanied by areas with presence of shrubs dominated by Rosmarinus officinalis, Rubus sp. and Genista scorpius. Pastures where characterized by the presence of Brachipodium sp. while areas with higher water availability where dominated by species such as Hedera helix. Furthermore, recurrence of fire burned 1000 ha in the central west part during 2005. In this area, the annual precipitation ranges from 430 to 890 mm in accordance with Digital Climatic Atlas of the Iberian Peninsula (DCAIP) (http://opengis.uab.es/wms/iberia/en_index.htm), and the average annual temperature is 13 • C. Forests are characterized by a hilly topography, with elevations ranging from 210 to 1240 m above sea level (a.s.l.) and slopes from 12% to higher than 45%. Cenozoic conglomerates, sandstones, clays and limestone dominate the lithology.
Requena forested site includes 38,750.15 ha affected by fire, and 36,884.14 ha within the surrounding area. Forests are mainly dominated by Pinus halepensis, while Pinus pinaster and Quercus ilex rotundifolia are present to a lesser extent. Furthermore, shrubland dominated by Rosmarinus officinalis, Quercus coccifera and Juniperus oxicedrus, and pasture areas characterized by the presence of Brachipodium sp., Ulex parviflorus and Rubia peregrina constitute a relevant land cover. Three fires affected recurrently the analyzed area after being burned in 1994. Two fires affected around 2900 ha in the eastern part in 2000 and 2003, and one small wildfire burned 100 ha in the western part in 2012. In this site, the annual precipitation ranges from 435 to 680 mm (DCAIP), and the average annual temperature is 13 • C. Forests are also characterized by a hilly topography, with elevations ranging from 180 to 1250 m a.s.l. and slope steepness ranging from 13% to higher than 45%. Sandstones and conglomerates from Mesozoic and Cenozoic Era dominate the lithology.
Yeste forested site includes 12,669.33 ha affected by fire, and 28,753.56 ha within the surrounding area. Forests are dominated by Pinus halepensis and, to a lesser extent, Pinus nigra and Pinus pinaster.
Shrublands are dominated by Rosmarinus officinalis, Juniperus oxicedrus and Quercus coccifera, while pastures are dominated by Brachipodim sp. and young trees from Quercus constitute the remaining vegetation. In this area, the annual precipitation ranges from 430 to 700 mm (DCAIP), and the average annual temperature is 14 • C. Forests are characterized by a hilly topography, with elevations ranging from 600 to 1670 m a.s.l., and slopes from 16.5% to higher than 50%. The lithology is dominated by dolomites, limestone and marls from Mesozoic Era.

Field Data
Field data was acquired in 136 field plots from three field campaigns performed during 2017 and 2018. A stratified random sampling technique was selected to establish the location of the field plots, considering a representative sample of terrain slope, aspect, fire severity, dominant pre-fire tree formation, height and cover of the forested areas. To perform this procedure, slope and aspect were derived using ALS data while fire severity was measured using passive remote sensing data from the Landsat archive. Slope and aspect where computed at 10 m resolution by resampling, using the bilinear interpolation in ArcGIS, the digital elevation model (DEM) generated at 1m resolution for ALS metric computation (see Section 2.3. ALS Data and Processing for further details). Fire severity was computed, within the SERGISAT project [51,52], using a simulation method based on the Composite Burn Index [53], called GeoCBI [54].
The Spanish Forest Map adapted by Rodrigues et al. [55] was used to determine pre-fire forest types and the height and cover of vegetation was derived from the ALS point cloud (see Section 2.3. ALS Data and Processing for information on data source and processing) to define homogeneous areas. Furthermore, road and path accessibility were checked to avoid inaccessible areas. A Leica VIVA ® GS15 CS10 GNSS real-time kinematic Global Positioning System was used to locate the centroid of the circular plots, which presented a variable radius of 5, 10 or 15 m. Plots with lower radii were mainly associated to homogeneous areas where really high tree or shrub covers, associated to post-fire regrowth processes, made difficult to move within the plot. Prometheus classification fuel type was in situ identified in each plot, determining one of the seven different fuel classes described in Table 1. An additional class was added to account for those areas without vegetation, named bare soil. Four photos in each of the cardinal points from the centroid of the plots were taken. Furthermore, the slope, elevation, flora inventory, tree and shrub characterization of height and cover, date and hour of acquisition were recorded. Table 2 include the number of field plot sampled by wildfire and fuel type.

ALS Data and Processing
The ALS data were acquired by the Spanish National Plan for Aerial Orthophotography (PNOA) between 2015 and 2016 using a Leica ALS80 discrete-return sensor operating at a wavelength of 1064 nm. Average flying altitude was 3150 m a.s.l. with a maximum scan angle of up to 25 • from nadir. The sensor used an average pulse repetition frequency of 230 KHz, average scanning frequency of 43.5 Hz and the average RMSE in z values obtained was 0.09 m. Point clouds were provided in 2 × 2 km tiles in LAS format in European Terrestrial Reference System (ETRS) 1989 Universal Transverse Mercator (UTM) with up to four return per pulse.
The processing of ALS point clouds started with the removal of noise points and overlapping returns, as the existence of vertical (z) or horizontal (x, y) displacements were found between flight lines by visual verification. The 699 tiles that cover the three areas affected by wildfires were filtered using the multiscale curvature classification algorithm [56] according to Montealegre et al. [57], implemented in MCC 2.1 command-line tool. The Point-Triangulated Irregular Network-Raster interpolation method [58] was selected according to Montealegre et al. [59] to generate a digital elevation model (DEM) with a 1-m grid resolution using ArcGIS 10.5 software. The point clouds were clipped to the spatial extent of the field plots and return heights were normalized by subtracting the DEM heights using FUSION LDV 3.60 open source software [60].
where p is the proportion of returns at regular intervals of 0.5 m or at defined Prometheus classification height intervals (h). The first step to compute LHDI and LHEI was the calculation of return proportion at regular intervals of 0.5 m according to Listopad et al. [48] ("LHDI_regular" and "LHEI_regular"). Furthermore, the indexes were computed using the Prometheus classification height ranges intervals at: 0.2 m, 0.6 m, 2.0 m and 4.0 m [64] ("LHDI_Prometheus" and "LHEI_Prometheus". Rumple index is the ratio of three-dimensional canopy surface model (CSM) to ground area. Rumple was computed as the ratio between the sum of the three-dimensional area of triangles from CSM grid points to the two-dimensional area of the grid cell surface. The CSMs were created for each plot using a 1.5 m pixel and a 3 × 3 smoothing algorithm, considering point cloud density of LiDAR-PNOA. The surface area of each CSM 1.5 m pixel is computed by creating triangles that fit the centroid of the pixel and those of the neighboring ones [47]. CSMs were created using the highest returns of each height range to account for canopy rugosity. Rumple was computed for Prometheus classification height ranges (0.2, 2 and 4) to characterize heterogeneity within each stratum, and for the forest canopy including the range between 0.2 up to 40 m.

Sentinel Data and Processing
Nineteen cloud-free images from the medium-high spatial resolution Sentinel-2A/B MSI sensor were used to cover the variability within the vegetative period for the three wildfires and their surrounding area under study ( Table 3). The images were downloaded from the European Space Agency's (ESA) Sentinel Scientific Data Hub (https://scihub.copernicus.eu/) with a processing Level-1C that includes geometric and radiometric corrections. Images were pre-processed using the Sen2cor plug-in implemented in SNAP software (https://step.esa.int/main/toolboxes/snap/) to convert Top-of-Atmosphere (TOA) reflectance values to Bottom-of-Atmosphere (BOA). The processing parameters for Sen2cor were chosen considering the environmental characteristics of the study area: Aerosol was set to rural, the summer period was chosen, a BRDF correction of 2 was applied, and a DEM at 90 m spatial resolution was considered to account for topographic correction. In this process, the bands were resampled to the original pixel resolution ranging from 10 up to 60 m. Table 3. Selected Sentinel-2 cloud-free images for the three wildfires and their surrounding areas used for generating the composite.

Wildfire
Dates Tile * stands for images that were captured in the same day within different tiles.
A yearly composite image was created for each wildfire site to better capture the vegetative period variability. In this sense, the mean reflectance value was computed for each band. Then, eight indices were calculated: the normalized difference vegetation index (NDVI) [65], the normalized difference infrared index (NDII) [66], the normalized difference water index (NDWI) [67], and the modified normalized difference water index (MNDWI) [68], the normalized burn ratio (NBR) [69], wetness, greenness, and brightness [70] (Equations (4)- (11)). Finally, the average value per plot was extracted for subsequent modelling steps.
where GREEN refers to band 3, RED is band 4, NIR refers to band 8, SWIR cirrus is band 10, SWIR 1 is band 11, and SWIR 2 is band 12 from Sentinel 2 A/B MSI.

Classification of Forest Fuels and Model Validation
The selection of ALS and Sentinel 2 metrics for forest fuel classification was performed using two approaches: (i) Spearman's rank correlation coefficient (p) that determines the strength and direction of the relationships between fuel types and remote sensing derived data, and (ii) all subset selection, which is an automatic selection procedure that selects the best subsets within a sample. Exhaustive, backward, forward and sequential replacement search approaches were tested following Domingo et al. [71]. The best subsets were subsequently used for model computation either when combining ALS and Sentinel 2 data or when modelling using each data source separately.
The performance of two frequently non parametric classification methods was tested to classify Prometheus fuel types: support vector machine (SVM), and random forest (RF). SVM was computed using linear (SVMl) and radial (SVMr) kernels. Cost and gamma were parametrized within an interval of 1-1000 and 0.01-1, respectively. RF classifier was parametrized by applying between 1 to 3000 trees to growth (ntrees) and between 1 to 3 metrics in each node (mtry) in accordance with Rodrigues et al. [72] and with our previous research Domingo et al. [71]. The model was computed using "randomForest" [73] and "caret" [74] packages in R environment and bias was corrected.
Models were validated using a stratified random sampling of 25% to include the different fuel types. The validation was executed 100 times to increase robustness in the results [75] and average performance values were computed. The classification overall accuracy, confusion matrices, user's and producer's accuracy were evaluated to compare and, subsequently, determine the best classification model [76]. Furthermore, metric importance for the models generated by combining ALS and Sentinel 2 data was analyzed by dropping one metric at a time and comparing overall accuracy between them. 2.6. Fuel type mapping and differences between burned and unburned areas Fuel type mapping was carried out using the most accurate model with a 10 m pixel resolution for the three study cases (Montmajor, Requena and Yeste). The differences between burned and unburned areas was performed by comparing the relative presence of a fuel type for each study case using a graphical assessment. The relative presence was computed using Equations (12) and (13).
RP fuel type surroundings = % area fuel type N in the surroundings % area fuel type N in the surroundings + % fuel type N in the wildfire (12) RP fuel type wildfire = % area fuel type N in the surroundings % area fuel type N in the surroundings + % fuel type N in the wildfire (13) where: RP stands for relative presence, N is fuel type. Table 4 shows the best classification model when combining ALS and Sentinel 2 data using the most explanatory metrics for fuel type classification determined by the Spearman rank correlation selection method. The model included three ALS metrics: the 25th percentile of return heights (P25), the percentage of all returns above mean (% all ret. above mean), rumple structural diversity index, and NDVI index from Sentinel 2 images, that showed correlation coefficients higher than 0.65 (Table A2 from Appendix A). The best classification method was SVM with radial kernel with an overall classification accuracy after validation of 0.59 (Table 4). The model was tuned with a cost value of 10 and a gamma value of 0.15. Lower accuracies where obtained with RF and SVM with linear kernel, that shows similar results. RF was parametrized with 2500 trees and 3 metrics in each node, while SVM with linear kernel was tuned with a cost value of 10 and a gamma value of 0.15. The analysis of metric importance (Appendix A, Table A3) determined that the dropping of % all ret. above mean decreased 0.14 in the overall accuracy, followed by rumple and P25, respectively. The dropping of NDVI slightly reduced overall accuracy.  Table 5 shows the performance of models computed using All subsect selection metrics and SVM with linear kernel. The selected ALS and Sentinel 2 metrics have overall accuracy values ranging from 0.51 when using the Exhaustive and forward approaches up to 0.57 when using the Backward selected metrics. As a result, none of these models computed using All subset selection metrics improved the performance obtained using the Spearman rank correlation selection method, and were disregarded for subsequent mapping. According to the results shown in Table 6, derived from the validation sample, there is confusion between fuel types 5 and 6. The structural similarities between both models as well as the low point density datasets may cause this confusion. Furthermore, we found confusion between fuel types 3 and 4, which may be associated to terrain complexity, and similarly confusion between fuel types 4 and 7. The high regrowth after fire frequently generate high density young forest formations with similar structural characteristics as Fuel type 7.  Table 7 shows the best classification model when using ALS data and the most explanatory metrics for fuel type classification, which were determined by the Spearman rank correlation. The model included three metrics: the 25th percentile of return heights (P25), the percentage of all returns above mean (% all ret. above mean), rumple structural diversity index. The best classification method was SVM with linear kernel with an overall classification accuracy after validation of 0.58, while SVM with radial kernel showed a similar performance (0.57). Both SVM models were tuned with a cost value of 10 and a gamma value of 0.01. Lower accuracies where obtained with RF with 0.54 overall accuracy, being parametrized with 500 trees and 2 metrics in each node. According to the results shown in Table 8, derived from the validation sample when classifying with ALS data, there greatest confusion is between fuel types 4, and 7 that generally represents dense canopies especially in areas affected by wildfires. Similarly, confusion is found between fuel types 1, 2 and bare soil associated with grasslands or low shrubs that are classified as bare soil or inversely.  Table 9 shows the best classification model when using Sentinel 2 data and the most explanatory metrics for fuel type classification, which were determined by the Spearman rank correlation. The model included four Sentinel 2 indexes and one band: NDVI, NBR, Wetness, Brightness and NIR. The best classification method was SVM with linear kernel with an overall classification accuracy after validation of 0.38, while SVM with radial kernel showed a similar performance (0.37). Both SVM models were tuned with a cost value of 10 and a gamma value of 0.01. Lower accuracies where obtained with RF was parametrized with 2500 trees and 3 metrics in each node. According to the results shown in Table 10, derived from the validation sample when classifying with Sentinel 2 data, there greatest confusion is between fuel types 4, 6 and 7 that generally represents dense canopies, but the optical data it is not able to differentiate the height. Similarly, confusion is found between fuel types 1, 2 and 3 associated with different shrub height.

Fuel Type Mapping and Differences between Burned and Unburned Areas
The models computed using data from only one data source, either ALS or Sentinel 2 data, were disregarded for subsequent mapping as presented lower overall accuracies. In this sense, the model that include the P25, % all ret. above mean, NDVI, and rumple was selected for further analysis and subsequent mapping for the three study sites (Montmajor, Requena, and Yeste). Figures 2 and 3 show the fuel type presence within each sub-study area and fuel type mapping, respectively. The predominance of fuel type 7 is confirmed in all the forested areas, being in Montmajor the one that reaches the highest value (62.68%). The Presence of fuel type 4 is significant in the study area, being more relevant in Requena (27.79%) and Yeste (18.08%), while reaching lower values in Montmajor (4.73%). Fuel type 5 have an intermediate presence in Montmajor and Yeste with values of 9.94% and 6.26%, respectively. The importance of fuel type 5 in Requena is residual. A residual presence of fuel type 2 and 3 was found all over the study area, while the presence of fuel type 0 and 1 varies from 6.3% in Requena, to 10.19% in Yeste, and up to 11% in Montmajor.

Fuel Type Mapping and Differences between Burned and Unburned Areas
The models computed using data from only one data source, either ALS or Sentinel 2 data, were disregarded for subsequent mapping as presented lower overall accuracies. In this sense, the model that include the P25, % all ret. above mean, NDVI, and rumple was selected for further analysis and subsequent mapping for the three study sites (Montmajor, Requena, and Yeste).  Generally higher presence of fuel type 4 is found within the areas affected by wildfire, and specially in Requena and Yeste, characterized by a drier climate and poorer soils than Montmajor. In contrast, higher presence of fuel type 5 and fuel type 6 are associated with the surrounding areas (e.g., northwest of Yeste and northeast of Montmajor), which denotes a higher predominance of the tree stratum that distinguish mature forests. On the other hand, fuel type 7 predominates in the areas affected by wildfires, which may be associated to a higher increase of shrub presence after-fire that    Generally higher presence of fuel type 4 is found within the areas affected by wildfire, and specially in Requena and Yeste, characterized by a drier climate and poorer soils than Montmajor. In contrast, higher presence of fuel type 5 and fuel type 6 are associated with the surrounding areas (e.g., northwest of Yeste and northeast of Montmajor), which denotes a higher predominance of the tree stratum that distinguish mature forests. On the other hand, fuel type 7 predominates in the areas affected by wildfires, which may be associated to a higher increase of shrub presence after-fire that remains within the forest structure. Aspect also plays an important role at local scales in vegetation regrowth. As can be observed in the southwest area of Requena a clear difference between north-east and south-west-oriented slopes exists. Thus, fuel type 4 is associated with south-west aspects while north-oriented areas have a predominance of fuel type 7. The presence of pastures is higher in Yeste and Montmajor, while in Requena its constrained to the central-east area. Furthermore, the presence of bare soil is low and associated to poor soils and rock presence (e.g., the south-east of Requena). Figure 4 shows that fuel type composition of forested areas affected by a wildfire 25 years ago still present differences with the non-burned area. Fuel type 7 has similar presence values in both, burned and non-burned areas, while the presence of fuels 5 and 6 is lower in burned areas with relative values around 0.25. Fuel type 4 has higher presence in burned areas than unburned ones with relative values around 0.6. On the other hand, fuel types 1, 2 and 3 show similar presence values which indicates a regrowth of these strata, except for case of Requena that still present lower values in the burned area. The presence of bare soil is higher in burned areas than in unburned ones in the three analyzed cases.

Discussion
Fuel type maps provide essential information to support preventive actions, fire management and fire modeling by forest managers [44]. This cartography is especially relevant in forested areas affected by wildfires, which can reach a higher structural diversity in an advanced stage of recovery [64]. Forest fires generate a partial or total modification of forest structure even after medium or advanced stage of recovery [55]. The fuel type classification performed in this research reveals the usefulness of low-density ALS data, integrated with Sentinel 2 images, to determine fuel types with moderate accuracy at regional scale in Mediterranean forested environments dominated by pines, oaks and quercus, characterized by high structural and topographical complexity.
The Spearman rho coefficients, considered as a good tool for determining the relationships between ALS and field metrics in accordance with Kristensen et al. [77], showed good results, agreeing with previous studies oriented to predict forest variables with ALS data in Mediterranean ecosystems [78]. According to the classification results, ALS-derived metrics have more importance in the models than multispectral Sentinel 2 data, which added a minor improvement to the classification than the ALS metrics (see Table A3 in Appendix A). Similar findings were reported previously for predicting fuel properties [35]. The percentage of all returns above mean showed the major importance on model performance, while 25th percentile of return heights and rumple index both showed similar and relevant importance. The inclusion of structural complexity metrics, such as rumple, in classification might be considered to generate parsimonious models, reducing the number of metrics to use in comparison to the height bin approach [8,43]. Although LiDAR height diversity index and LiDAR height evenness index diversity metrics were not included in the most accurate models, both showed a high correlation with Spearman values of 0.77 and 0.72, respectively in accordance with Listopad et al. [48] and Gelabert et al. [64].

Discussion
Fuel type maps provide essential information to support preventive actions, fire management and fire modeling by forest managers [44]. This cartography is especially relevant in forested areas affected by wildfires, which can reach a higher structural diversity in an advanced stage of recovery [64]. Forest fires generate a partial or total modification of forest structure even after medium or advanced stage of recovery [55]. The fuel type classification performed in this research reveals the usefulness of low-density ALS data, integrated with Sentinel 2 images, to determine fuel types with moderate accuracy at regional scale in Mediterranean forested environments dominated by pines, oaks and quercus, characterized by high structural and topographical complexity.
The Spearman rho coefficients, considered as a good tool for determining the relationships between ALS and field metrics in accordance with Kristensen et al. [77], showed good results, agreeing with previous studies oriented to predict forest variables with ALS data in Mediterranean ecosystems [78]. According to the classification results, ALS-derived metrics have more importance in the models than multispectral Sentinel 2 data, which added a minor improvement to the classification than the ALS metrics (see Table A3 in Appendix A). Similar findings were reported previously for predicting fuel properties [35]. The percentage of all returns above mean showed the major importance on model performance, while 25th percentile of return heights and rumple index both showed similar and relevant importance. The inclusion of structural complexity metrics, such as rumple, in classification might be considered to generate parsimonious models, reducing the number of metrics to use in comparison to the height bin approach [8,43]. Although LiDAR height diversity index and LiDAR height evenness index diversity metrics were not included in the most accurate models, both showed a high correlation with Spearman values of 0.77 and 0.72, respectively in accordance with Listopad et al. [48] and Gelabert et al. [64].
The comparison between classification methods shows that SVMr had the highest accuracy to classify Prometheus fuel types in accordance with García et al. [28]. RF showed and overestimation previous to the validation phase that has been previously reported when using low or medium sample sizes. Arellano-Pérez et al. [23] pointed out that RF overfitted the data for reduced sample plots (123 field plots) when modelling surface and canopy fuel characteristics with Sentinel-2A data. Hu et al. [79] described the same problem when predicting forest stock volume with Sentinel-2A images and 459 field plots. Our previous studies also showed that overfitting was produced when predicting different forest attributes (i.e., volume, biomass) [78] or residual biomass [80] using low density ALS data. Furthermore, specific studies about the overestimation of RF have been carried out as for example Janitza et al. [81], pointing out that few observations and larger number of predictor variables produce overestimations in RF algorithm. The performance of the classification when combining ALS and Sentinel 2 data with an accuracy value of 59% shows better results to the ones obtained by Huesca et al. [12], which used ALS-PNOA data and spectral mapping methods to classify Prometheus fuel types with a 44% overall agreement including five fuel types. Higher performance was obtained by García et al. [28], with a 88% overall accuracy agreement when predicting Prometheus fuel types. However, these authors applied decision rules to the output of a SVM classification, based on the mean height and the vertical distribution of LiDAR returns, and included ALS data of higher point densities (1.5 to 6 points m −2 ) and multispectral images of higher resolution (2 m grid). Alonso-Benito et al. [29] obtained also higher accuracies (from 84.27 to 85.43% of agreement) using also higher point density ALS data and images resolution. The comparison with studies that used a different classification system than Prometheus is not direct as these classifications are based on different parameters (i.e., different height and cover thresholds). In this sense, Alonso-Benito et al. [82] obtained different performance of classification using the same datasets, depending on the fuel type classification system applied, with a 10% difference in overall agreement between the NFFL and a specific fuel model developed for Canarias Island specific forest types. Other possible reasons to explain the higher performance in fuel type classifications in the aforementioned approach when comparing to our results may be associated to the higher structural complexity of vegetation in our study area, associated to the occurrence of wildfires [64]. Furthermore, our results were in accordance with those of Huesca et al. [12] when analyzing confusion between classes, which reported confusion between fuel types 5, 6 and 7 as well as between fuel types 3 and 4. García et al. [28] also reported confusion between fuel types 5, 6 and 7, but no problems for differentiating fuel types 3 and 4 were found. Riaño et al. [19] also found the same pattern of misclassification between Prometheus fuel types 3 and 4, as well as 5 and 6.
The mapping of fuel type within areas previously affected by wildfires and its comparison with control areas allow determining whether structural differences exists between both forested areas. A review published by Gómez et al. [83] determined, using passive remote sensing data in Spain, that pre-fire vegetation cover was reached after 7-20 years after fire. Rodrigues et al. [55] determined that vegetation recovery time after high-severity wildfires ranges from 21.5 years for high seeding trees, to 25.5 years for resprouter trees and up to 52.9 years for low seeding trees in Spanish Mediterranean forests. Our results are in accordance with Gelabert et al. [64], who pointed out that differences in height, canopy cover and structural diversity indexes were still present after 21 years of a fire in Pinus halepensis forests located within the central sector of the Ebro basin (northeast Spain). Furthermore, though taking into account the moderate accuracy of our models, the mapping results show concordance with our knowledge from field work campaigns. Fuel type 4 in burned areas was mainly related with high stem densities of thin trees in the case of Pinus forests, while Quercus forests where accompanied by bulky shrubs. The presence of fuel types 5 and 6 in burned areas was low, and those burned areas with dominance of tree strata were associated with high presence of shrubs with continuity to tree strata (fuel type 7). Furthermore, the lower presence of fuel types 1, 2 and 3 in Requena burned area may be associated with a high presence of stony areas and climate conditions. The present study shows the utility of integrating freely available ALS and multispectral data to classify Prometheus fuel type at a regional scale. More research should be done to increase discrimination between fuel type models by integrating other remote sensing datasets. The combination with high-resolution multispectral data or SAR data might improve classification accuracy [28,46]. In this sense, the use of orthophotos with NIR information acquired during recent years within the Spanish context might boost model performance. In addition, further analysis might focus on predicting fuel parameters within areas affected by wildfires as well as analyzing fuel type differences between burned and unburned areas.

Conclusions
This study assessed the usefulness of integrating low-point density ALS data and Sentinel 2 data to classify and map Prometheus fuel types in three different topographically complex Mediterranean forests, dominated by pines, oaks and quercus. These areas encompass forests affected by wildfires in 1994 and their surrounding areas. Spearman's rank coefficient has been the most powerful selection method to generate a representative and meaningful classification fuel types at regional scale. The SVM with radial kernel method produced the most accurate fuel type classification model, which included three ALS metrics: the 25th percentile of return heights, the percentage of all returns above mean, rumple structural diversity index, and NDVI index from Sentinel 2 images. The classification achieved an overall accuracy of 0.59 after validation. The use of ALS derived diversity structural indices boosts model performance as well as the fusion with multispectral data from Sentinel 2, allowing mapping at 10 m resolution. Forest fuel mapping provides valuable information to assist forest management. The structural characterization of forested areas affected by wildfires after two and a half decades constitutes a relevant input for analyzing post-fire forest treatments as well as to define actions to reduce fuel density at canopy level or lower strata.  Percentage of first returns above a height-break, above the mean or the mode e.g., % first ret. Above 0.20 Percentage of all returns above a height-break, above the mean or the mode e.g., % all ret. Above mean Percentage of all returns with a range of 0.5 m e.g., % all ret. between 1 and 1.5 m Canopy relief ratio CRR All returns above a height-break, above the mean or the mode x 100 e.g., (All ret. Above 0.20)/(total first ret.) by 100 Percentage of returns between 0.2 m and 95th percentile according to Naesset (2004)