Vegetation Fuel Mapping at Regional Scale Using Sentinel-1, Sentinel-2, and DEM Derivatives—The Case of the Region of East Macedonia and Thrace, Greece

Abstract

The sustainability of Mediterranean ecosystems, even if previously shaped by fire, is threatened by the diverse changes observed in the wildfire regime, in addition to the threat to human security and infrastructure losses. During the two previous years, destructive, extreme wildfire events have taken place in southern Europe, raising once again the demand for effective fire management based on updated and reliable information. Fuel-type mapping is a critical input needed for fire behavior modeling and fire management. This work aims to employ and evaluate multi-source earth observation data for accurate fuel type mapping in a regional context in north-eastern Greece. Three random forest classification models were developed based on Sentinel-2 spectral indices, topographic variables, and Sentinel-1 backscattering information. The explicit contribution of each dataset for fuel type mapping was explored using variable importance measures. The synergistic use of passive and active Sentinel data, along with topographic variables, slightly increased the fuel type classification accuracy (OA = 92.76%) compared to the Sentinel-2 spectral (OA = 81.39%) and spectral-topographic (OA = 91.92%) models. The proposed data fusion approach is, therefore, an alternative that should be considered for fuel type classification in a regional context, especially over diverse and heterogeneous landscapes.

1. Introduction

1.1. Fuel Type Maps for Wildfire Risk Assessment

Wildfire regime is an inherent characteristic of the terrestrial biomes, having played a key role throughout Earth’s history in the evolution of the biota in these areas [1]. In recent decades, changes in the historically observed wildfire regimes, associated with extended duration and earlier start of fire season as well as with changes in the fire interval, intensity, and size [2], mask out the positive role of fire for biodiversity maintenance and raises concerns [3,4,5]. Indeed, there is evidence that disturbances in the distribution, composition, and structure of ecosystems, even in the ones previously adapted to the historically noted wildfire regimes, are likely to lead to a decline in the resilience of ecosystem functions to wildfires and their ability to provide ecosystem services to the society [6]. Changes in the wildfire regime have also led to significant wildfire disasters associated with human losses and major socioeconomic impacts on human assets, as in the case of Greece (2007 and 2018), Portugal (2003, 2005, and 2017), and Australia (2009 and 2019) [7].

The 2021 and 2022 fire seasons clearly demonstrated the impact of wildfire regime change in Southern Europe and in the Mediterranean biome in general. This also partly verifies earlier concerns [2] that the overall decrease observed in recent decades (1980–2010) in the annual number of fires and annual burned areas for several Euro Mediterranean countries might be only a short-term positive effect of the increased effort and resources committed in fire suppression and fire prevention efforts. Greece, in 2021, faced an unprecedented ecological disaster with 78,813 hectares burned—more than three times higher than the 2008–2020 annual average (22,809 ha), according to the European Forest Fire Information System [8]. Likewise, during the 2022 fire season, the burned area recorded in the western part of the Euro Mediterranean area (Portugal, Spain, and France) presented a three-fold increase compared to the mean burned area from 2006 to 2021, reaching up to 477,361 hectares [9].

There are several causal drivers shaping this regime modification and the occurrence of extreme wildfire events, including human-induced and natural climate change, land management, socioeconomic activities, fire management, fire exclusion, etc. While the underlying process leading to the generation of extreme wildfire events is complex and varies across space and time due to the interaction between these factors [10], the new reality urges for changes in the fire management paradigm, including, among others, more efforts put in up-to-date and reliable wildfire risk assessment and reduction strategies [7].

Wildfire risk assessment can support decision makers as well as disaster management to proactively respond to the challenges posed by the changing wildfire regime, even in the fire-prone countries of the Euro Mediterranean area [11]. Effective wildfire risk assessment and mitigation is also a prerequisite for achieving the EU forest strategy for the 2030 goal to improve the quantity and quality of EU multi-functional forests. EU’s biodiversity strategy for 2030 is also threatened by this changing regime, as indicated by the fact that almost half of the protected areas in southwest Europe were burned during the 2022 summer fire season [9].

Wildfire risk models must be ingested with up-to-date, spatial explicit representations of the real world, with biomass and the characteristics of live and dead vegetation (wildland or vegetation fuels) available for burning being the primary factor influencing fire behavior and risk assessment, and the only element of fire triangle that can be managed. Fuel types constitute fuel element complexes of species, structure, pattern, size, etc., that may cause a predictable fire behavior parameter (e.g., rate of fire spread, fireline intensity, and flame length) or difficulty of control under specified weather conditions [12], yet they can be used for standardizing the wide range of natural fuels’ physical characteristics and facilitating fire behavior modeling and risk assessment [13]. Fuel type maps are essential not only for drafting mitigation strategies and fire management policy tools [11] but especially at local to regional scales for fire behavior modeling, fire towers and water tanks construction, landscape level fuel management, firefighting resources allocation, smoke impact prediction, and monitoring of vegetation recovery after fire [13,14].

1.2. Fuel Type Mapping through Remote Sensing

Reliable but time- and resource-demanding field surveys were historically used for fuel type mapping at the local scale, gradually replaced by remote sensing methods, facilitated among others by the improvement in sensors characteristics and data availability [15,16]. Optical multispectral data have been the main source of information for fuel type mapping from local to continental scales [17]. Hyperspectral data have also been explored in a Mediterranean context, providing little, if any, benefits when compared to accuracy attained by conventional multispectral sensors [18]. Such a finding might be well linked (apart from the classification scheme) to the limited ability of passive sensors to penetrate the canopy in order to describe the vertical distribution of the fuels and understory characteristics [15,17,19].

On the contrary, active sensors have the potential to overcome this limitation and record and characterize the vertical structure of fuel complexes [20]. The use of light detection and ranging (LiDAR) data was demonstrated in several cases for the mapping and characterization of fuel complexes [21,22]. Nevertheless, the LiDAR data have well-known constraints related to availability, cost, and update frequency, and until the introduction of future spaceborne LiDAR systems [16], LiDAR sensors use is unlikely to be used operationally for wall-to-wall regional or national applications (except for few countries sharing such financial and technical resources). On the other hand, synthetic aperture radar (SAR) data are a viable alternative to LiDAR sensors, also having the advantage of penetrating canopy layers due to the longer wavelengths sensed compared to optical sensors [23]. In this regard, SAR data have been employed for wildland fuel characteristics retrieval [24,25,26,27], vegetation mapping [28,29], and forest attributes estimation [30,31], which are essential for fuel type mapping at a regional and global scale [17].

Remote sensing approaches often consider topographical gradients such as elevation, aspect, and slope to generate fuel maps, assuming that these attributes are related to vegetation, which is correlated well with fuel dynamics and fuel models [14]. Numerous studies highlight the importance of topography variables on land cover and vegetation type classification [32,33], as it is well known that species composition of vegetation depends to a certain extent on ecological site factors [34].

Recently, in fuel type mapping, there is also a tendency for information extraction based on data fusion from different sensors as a promising approach to overcome individual sensors/platforms limitations [16]. Radar and LiDAR data mainly complement the optical imagery for the discrimination of vegetation or fuel types. Chirici et al. [35] predicted forest fuel types through airborne laser scanning and IRS LISS-III imagery in Sicily (Italy). García et al. [36] used multispectral and LiDAR data fusion for fuel type mapping using a support vector machine and decision rules in the Natural Park of the Alto Tajo in central Spain. Ruiz et al. [37] applied an object-based approach for mapping forest structural types based on low-density LiDAR and multispectral imagery in the Natural Park of Sierra de Espadán in Spain. Mendes et al. [23] applied optical and SAR synergism for mapping vegetation types in the Endangered Cerrado/Amazon Ecotone in Brazil. Dobrinić et al. [38] used Sentinel-1 and Sentinel-2 time series for vegetation mapping using random forest classification in Croatia. Li et al. [39] integrated ALOS PALSAR L-band SAR and Landsat ETM+ optical data to estimate forest fuel loads in northern Sweden. D’este et al. [27] estimated fine dead fuel load using optical LIDAR and SAR data and field inventory data, applying three different algorithms in south-eastern Italy.

Even though these studies have presented positive results, a challenge still remains to retrieve spatial information on fuel types using cost-effective innovations and easily updated data sources [17,40,41,42]. Furthermore, there is little information about the synergetic use of different data, such as SAR, optical, and topographic-related data, and their potential contribution to forest fuel types mapping in a regional context over Mediterranean areas, implying high compositional and structural heterogeneity.

In this context, the overall aim of this study was to develop spatially explicit wildland fuel-type maps within the Region of East Macedonia and Thrace, Greece (REMTH), based on the fusion of multiple types of remote sensing data. Fuel types refer to vegetation classes (composition, structure) exhibiting similar fire-propagation behavior [13,40,43]. Consequently, fuel type maps are often developed using vegetation maps by assigning a unique fuel bed to each vegetation class [44,45]. The present research is carried out as part of an operational project for the establishment of a regional Risk and Resilience Assessment Center (RiskAC) aiming to act in the interface between science and policy in REMTH. Fuel-type maps will be used for decision-making for short-term predictions of wildfire risk at an operational level, contributing to wildfire suppression and management. The specific objectives of this study are to (1) provide region-wide, easy-to-update fuel type maps based on a nomenclature facilitating fire modeling; (2) evaluate the potential of synergistic use of Sentinel-2, Sentinel-1, and topography variables; and (3) explore the relative importance of spectral indices, backscattering information, and topography in fuel type mapping over a regional context.

2. Materials and Methods

2.1. Study Area

The focus of our study was set on the 7919.17 square kilometers of forest and natural areas located within the Region of Eastern Macedonia and Thrace (REMTH) occupying the north-eastern part of Greece. Elevation ranges from sea level to 2200 m above sea level (a.s.l.), and along the slope, gradient forms a highly diverse terrain (Figure 1).

The average annual temperature recorded over the 1960–1990 period was 16 °C, and the average annual precipitation was 481 mm. REMTH predominantly features a hot-summer Mediterranean climate (Csa) near the coastline, but further into the mainland of the region, and in the Rhodopi mountain range across the Greek–Bulgaria borders, the climate becomes continental, winter is harsh, has no dry season, and has a warm summer (Dfb) [46].

Climatic variability is also reflected in the natural environment of the region. Three altitudinal vegetation zones were recognized. The first vegetation zone (Quercetalia ilicis) in the study area is located above sea level to an altitude of 200–300 m. In this area, the vegetation zone is represented by the sub-zone Quercion ilicis, in which various plant communities appear that are partly degraded and partly territorially dependent (permanent plant communities). The floristic composition of these communities species is involved species such as Erica manipuliflora, Erica arborea, Arbutus unedo, Calicotome villosa, Spartium junceum, Quercus ilex, Fraxinus ornus, Phillyrea latifolia, Quercus pubescens, etc. This subzone is optimum for the growth of Pinus halepensis ssp. halepensis, as well as its Pinus halepensis ssp. brutia. The second vegetation zone (Quercetalia pubescentis) in the study area is located above the previous one and up to an altitude of 1000 m. The vegetation zone has 3 sub-zones in Greece, with 2 in the study area. The first sub-zone is Ostryo-Carpinion, which includes species such as Quercus cocciferra, Carpinus orientalis, Fraxinus ornus, Pistacia terebinthus, Ligustrum vulgare, Rhus coriaria, Cotinus coggygria, Quercus pubescens, Quercus frainetto, Acer monspessulanum, Sorbus torminalis, etc. The second sub-zone that is above the previous sub-zone is Quercion confertae and includes the oak forests of the area, also with the presence of species Castanea sativa, Tilia tomentosa, Acer obtusatum, Ostrya carpinifolia, Carpinus betulus, etc. The third vegetation zone (Fagetalia) in the study area is located above the previous one and up to the altitude of 1800 m. This vegetation zone has two sub-zones, but only one is present in the study area (Fagion moesiacae). The vegetation zone includes the forest of Pinus nigra, Fagus sylvatica, and Abies borisii-regis. The fourth vegetation zone (Vaccinio-Picetalia) is represented by the sub-zone Vaccinio-Piceion, which is located across the forest lines and includes the forest of Picea abies. The last vegetation zone (Astragalo-Acantholimonetalia) is located above the forest lines and includes grasslands that degraded due to overgrazing because these lands serve as summer pastures for herds [47,48,49].

2.2. Methodological Analysis Outline

For the development of the fuel type maps, three main steps were followed. At the first step, fuel type scheme development relied on field reconnaissance and measurements, as well as manual interpretation of VHR base maps (Google Earth) and existing auxiliary vegetation maps. The second step involved Sentinel-2 and Sentinel-1 dataset preprocessing, spectral information enhancement, backscattering information retrieval, and topographic variables extraction. In the third step, three models were developed using different datasets. The final step of the workflow was the classification accuracy assessment along with the variable importance analysis rating. The overall process diagram of the methodology used in the current research is shown in Figure 2.

2.3. Fuel Type Sampling

The forest fuel type classification scheme includes nine fuel models (FM) (

Table 1. Fuel models of the study area.

FM	Fuel Model	Average Fuel Height (cm)	Total Fuel Load (t/ha)	Litter Depth (cm)
FM 1	Shrubs up to 2 m height (Pseudomaquis)	182	20.50	2.1
FM 2	Forest litter layer of Oaks	3.5	2.32	3.1
FM 3	Forest litter layer of black pine	8.2	6.43	7.2
FM 4	Mediterranean grasslands	30	4.85	1.5
FM 5	Forest litter layer beech–black pine	5.5	4.9	2.9
FM 6	Forest litter layer of Norway spruce	4.9	2.8	6.3
FM 7	Forest litter layer of beech	4.3	1.82	4.3
FM 8	Forest litter layer of Calabrian pine	7.1	2.35	6
FM 9	Forest litter layer of deciduous shrubs	3.5	1.12	2.8

). The nine-class scheme was specified considering fire behavior modeling requirements, vegetation diversity within REMTH, and the need to adopt a simple and reproducible scheme for future updates. The successional stages, time since disturbance, and stand-specific characteristics (i.e., forest age) were not considered during scheme development to avoid complexity and increase confidence in the mapping workflow. An additional class (NV) was added to account for those areas without vegetation or other land uses (water, burnt areas, urban areas, bare soil, crops, irrigated crops, and floodplains).

Fuel parameters were determined following the field sampling procedures described by Mallinis et al. [18]. In brief, the methodology developed by Brown et al. (1982) [50] for inventorying biomass was used in representative, homogeneous fuel plots. The sampling approach included both 30 m transects (for the 1 h, 10 h, 100 h, 1000 h, and total fuel loads) and 100 m² square plots (for the herbaceous (live and dead) and shrub vegetation loads, foliage load, litter load, and depth).

A total of 63,742 sample points (Figure 3) were collected using visual image interpretation, utilizing recent very high spatial resolution (up to 0.5 m) satellite imagery, available via Web Map Service (WMS) standards, and vegetation and land-use maps from Greek Ministry of Environment and Energy. Vegetation and land-use maps, as a key part of forest management studies at a scale of 1:20,000, provide information on the distribution of forest species within the forest areas of the study area. The number of samples per fuel type and their spatial distribution are presented in Figure 3.

2.4. Remote Sensing Data and Preprocessing

Sentinel-2 surface reflectance products (L2A) were acquired on August 2021, and a monthly composite with average values was created within the Google Earth Engine (GEE) environment to fully cover the spatial extent of the region of REMTH. A series of spectral indices (

Table 2. Mathematic formula and motivation for considering each spectral index in the fuel type classification process.

Spectral Index	Equation and Use	Reference
Normalised Difference Vegetation Index—NDVI	$NDVI=\frac{NIR-RED}{NIR+RED}$	[54]
	One of the most widely used vegetation indicators that highlights the condition of vegetation.
Enhanced Vegetation Index—EVI	$EVI=2.5\frac{NIR-RED}{NIR+6\ast RED-7.5\ast BLUE+1}$	[55]
	Strongly related to forest cover.
Normalized Difference Water Index—NDWI	$NDWI=\frac{NI{R}_{n2}-SWI{R}_1}{NI{R}_{n2}+SWI{R}_1}$	[56]
	Monitor changes in the water content of leaves.
Soil Adjusted Vegetation Index—SAVI	$SAVI=1.5\frac{NIR-RED}{NIR+RED+0.5}$	[57]
	Takes into account the terrain and corrects the effect of soil brightness in areas with low vegetation cover.
Optimised Soil-Adjusted Vegetation Index—OSAVI	$OSAVI=1.16\frac{NIR-RED}{NIR+RED+0.16}$	[58]
	Provides greater soil variation than SAVI for low vegetation cover, while it demonstrates increased sensitivity to vegetation cover.
Green Normalized Difference Vegetation Index—GNDVI	$GNDVI=\frac{NIR-GREEN}{NIR+GREEN}$	[59,60]
	The use of the green band in areas with sparse vegetation reduces the sensitivity of the vegetation index to changes in the vegetation cover fraction.
Renormalized Difference Vegetation Index—RDVI	$RDVI=\frac{NIR-RED}{\sqrt{NIR+RED}}$	[61]
	Related to leaf biomass and canopy structure.
Normalised Difference Vegetation Index—NDVI2	$NDV{I}_2=\frac{NI{R}_{n2}-RED}{NI{R}_{n2}+RED}$	[62]
	Yields higher correlations with biomass than the standard NDVI.
Sentinel Improved Vegetation Index—SVI	$SVI=\frac{NIR-R{E}_2}{NIR+R{E}_2}$	[63]
	Additionally considers the red edge part of the spectrum.
Modified triangular Vegetation Index—MTVI2	$MTV{I}_2=1.5\frac{1.2\left(NIR-GREEN\right)-2.5\left(RED-GREEN\right)}{\sqrt{{\left(2NIR+1\right)}^2-\left(6NIR-5\sqrt{RED}\right)-0.5}}$	[64,65,66]
	Sensitive to leaf and canopy structure and leaf area index and weakens the effect of the soil background.

) were considered as potential predictors based on finding of previous fuel mapping [51] or vegetation mapping studies [52,53].

Sentinel-1 images were also retrieved and processed using the GEE platform [67]. Sentinel-1 C-band images of interferometric wide swath mode, VV (vertical transmit–vertical receive) and VH (vertical transmit–horizontal receive) polarizations, acquired in August 2021. The Sentinel-1 data were initially calibrated into σ⁰ backscatter coefficient values (Equations (1)) and despeckled using a 5 × 5 Lee filter.

$${\sigma }^{0 }\left(db\right)=10\times log10\times \left({\sigma }^{0 }\right)$$

(1)

By default, the Shuttle Radar Topography Mission (SRTM) digital elevation model, at a resolution of 1 arc-second (∼30 m) [68], available as a GEE asset, was used for terrain correction [69]. Finally, the Sentinel-1 mosaic with VH and VV polarization bands covering the study area was used for the ratio (VV/VH) calculation (Equation (2)).

$$RATIO=\frac{VV}{VH}$$

(2)

Information on the topography was derived from the digital elevation model (DEM) for Europe (EU-DEM) created for the Copernicus programme at a 25 m spatial resolution [70]. Elevation, slope, and aspect were incorporated into the analysis process.

2.5. Classification Models

The random forest (RF) classification algorithm [71], a non-parametric, tree-based classifier, was employed for fuel type classification. To investigate the optimal approach, three different RF models were developed using the following datasets (

Table 3. Variables of random forest models: spectral index model (SP); spectral index and topographic variables model (SPT); spectral indices and topographic variables and backscatter information model (SPTS).

Variables of Random Forest Models
	SP	SPT	SPTS
Spectral variables	NDVI	NDVI	NDVI
	EVI	EVI	EVI
	NDWI	NDWI	NDWI
	OSAVI	OSAVI	OSAVI
	SAVI	SAVI	SAVI
	GNDVI	GNDVI	GNDVI
	RDVI	RDVI	RDVI
	NDVI2	NDVI2	NDVI2
	SVI	SVI	SVI
	MTVI2	MTVI2	MTVI2
Topographic variables		DEM	DEM
		SLOPE	SLOPE
		ASPECT	ASPECT
Backscattering variables			VV
			VH
			RATIO VV_VH

):

• Spectral indices model (SP) using the 10 spectral indices;
• Spectral indices and topographic variables model (SPT) using 13 variables;
• Spectral indices, topographic, and backscatter information (SPTS) model using 16 variables.

The “randomForest” [72] and “caret” [73] packages were used for developing an end-to-end processing workflow (training, prediction, and mapping) within R environment software [74].

The reference data were split into training and test data samples with a ratio of 70:30 using random stratified sampling to include the different fuel types. In the modeling process, two-thirds of the train data were used while one-third (out-of-bag samples (OOB)) were used for model validation, calculating the out-of-bag (OOB) estimate of error rate [71]. RF was operated with default settings for the number of input variables at each split ((mtry) equal to the square root of the number of input variables at each node), the node size was set to the default value of 5, and the number of trees to grow (ntrees) was set to 1000.

2.6. Accuracy Assessment

For the validation process and the RF performance, the test data were used along with the reference data collected during the field survey. For the accuracy assessment, the confusion matrix for the actual and predicted fuel type identities’ overall accuracy (OA) [75], user’s accuracy (UA) [75], producer’s accuracy (PA) [75], kappa and weighted kappa statistic [76], and the out-of-bag (OOB) estimate of error rate were also reported. User’s accuracy (UA) represents the proportion of classified sample plots that correctly represents that category on the ground. Producer’s accuracy (PA) describes how well features on the ground are correctly classified on a map [29,75]. Overall accuracy (OA) is calculated by dividing the total correctly classified sample by the total number of observations [29,75]. The kappa statistic reflects the difference between observed agreement and agreement expected by random chance [77]. The weighted kappa takes into consideration the different levels of disagreement between categories [78]. The OOB error is the OOB estimate of the error rate [71]. Moreover, McNemar’s test [79] was performed to statistically compare error matrices of the three developed models, focused on the binary distinction between correct and incorrect class allocations (Equation (3)) [80].

$$Z=\frac{\left|f12-f21\right|-1}{\sqrt{f12+f21}}$$

(3)

A difference in accuracy between the confusion matrices of the classifiers is statistically significant (p ≤ 0.05) if the Z value is more than 1.96 [80].

2.7. Variable Importance

Furthermore, the relative importance of each variable was examined by calculating the mean decrease in accuracy (MDA). The MDA measure implemented in the RF classifier is the most well-known measure to estimate variable importance [81,82]. MDA is a straightforward technique [83,84] derived from the idea that if the variable is not important, then rearranging its values should not reduce prediction accuracy [85]. MDA is calculated by averaging the difference in OOB error before and after one variable is randomly permutated from the analysis.

3. Results

3.1. Forest Fuel Type Classification

The results indicate the SPTS model as the most accurate model, with overall accuracy (OA) of 92.76%, followed by SPT (OA = 91.92%) and SP (OA = 81.39%). Accuracy metrics are presented in

Table 4. Fuel type classification accuracy for the spectral indices model (SP); spectral indices and topography model (SPT); and spectral indices, topography, and backscatter information model (SPTS).

	SP	SPT	SPTS
Overall accuracy (OA)	81.39%	91.92%	92.76%
Kappa	0.75	0.89	0.90
Weighted kappa	0.76	0.90	0.91
Out-of-bag (OOB) error rate	18.58%	8.42%	7.50%

, and the resulting maps generated by the three models are illustrated in Figure 4. McNemar’s test was applied to assess whether statistically significant differences between classifications exist. In particular, based on McNemar’s test, there was a significant difference (Z >1.96) at the 5% significance level between the confusion matrices of the SP, SPT, and SPTS.

3.1.1. Spectral Indices Model

Regarding the fuel-type map produced by the SP model, the results indicate overall accuracy OA (81.39%) (Table 4). In terms of the individual classes’ accuracy (

Table 5. Confusion matrix for the spectral indices (SP) model: values are given in percentage of validation points.

	Actual										UA (%)
Predicted	FM1	FM2	FM3	FM4	FM5	FM6	FM7	FM8	FM9	NV
FM1	69.88	1.10	4.92	1.42	3.62	4.22	0.28	7.94	4.88	0.92	66.22
FM2	0.79	69.13	0.97	0.57	9.19	0.69	6.29	0.22	7.09	0.01	72.1
FM3	1.58	0.18	58.11	0.28	2.79	6.71	0.40	9.96	0.38	0.13	69.2
FM4	1.49	0.18	0.39	44.23	0.00	0.09	0.00	0.00	2.41	1.15	74.4
FM5	0.00	0.73	0.29	0.00	46.80	2.15	1.64	0.30	0.00	0.00	70.89
FM6	2.45	0.27	7.05	0.09	10.58	77.54	0.17	3.30	0.06	0.01	82.43
FM7	0.26	15.43	0.97	0.47	23.12	1.55	84.54	0.30	8.11	0.05	77.88
FM8	4.55	0.00	17.95	0.19	0.84	5.34	0.00	71.76	0.06	0.16	74.96
FM9	6.04	11.14	1.16	3.59	1.67	0.34	4.42	0.52	69.22	0.64	73.65
NV	12.96	1.83	8.20	49.15	1.39	1.38	2.27	5.69	7.79	96.93	88.96
OA (%)	69.88	69.13	58.11	44.23	46.8	77.54	84.54	71.76	69.22	96.93

), certain issues were detected in the class of FM4 Mediterranean grassland, which was erroneously assigned in other classes (PA = 44.23%) and mainly in other uses (NV). This can be partially related to the fact that the NV class also includes agricultural land, which may be spectrally similar to grassland cover. Misclassifications were also noticed for type FM5 forest litter layer beech–black pine (PA = 46.80%), which is confused with the FM6 forest litter layer of Norway spruce and FM7 forest litter layer of beech, and this is obviously due to the mixing rate of type FM5 forest litter layer beech–black pine. Errors also occur between the FM3 forest litter layer of black pine and the FM8 forest litter layer of Calabrian pine. Moreover, the FM7 forest litter layer of beech shows the highest accuracies.

3.1.2. Spectral Indices and Topographic Variables Model

Compared to the SP model, the SPT model showed higher overall accuracy (OA = 91.92%) (Table 4). According to the confusion matrix (

Table 6. Confusion matrix for the spectral index and topographic variables (SPT) model: values are given in percentage of validation points.

	Actual										UA (%)
Predicted	FM1	FM2	FM3	FM4	FM5	FM6	FM7	FM8	FM9	NV
FM1	86.69	1.64	2.99	0.47	2.23	0.00	0.11	4.19	3.80	0.45	81.89
FM2	1.23	86.67	1.64	0.19	0.00	0.00	1.87	1.05	3.99	0.00	86.9
FM3	0.53	0.37	79.05	0.19	2.79	0.17	0.40	9.51	0.13	0.05	83.32
FM4	0.00	0.00	0.10	82.99	0.28	0.60	0.45	0.00	0.57	0.06	96.59
FM5	0.00	0.00	0.00	0.00	59.33	2.07	0.40	0.00	0.00	0.00	87.30
FM6	0.00	0.00	0.10	0.19	8.36	96.64	0.34	0.00	0.00	0.00	96.64
FM7	0.09	4.29	1.06	0.57	24.79	0.52	93.66	0.00	3.23	0.00	88.69
FM8	5.34	0.18	12.16	0.19	0.00	0.00	0.00	82.25	0.25	0.17	83.94
FM9	2.10	6.76	0.58	1.51	1.95	0.00	2.77	0.15	85.12	0.21	87.27
NV	4.03	0.09	2.32	13.71	0.28	0.00	0.00	2.85	2.91	99.06	96.58
OA (%)	86.69	86.67	79.05	82.99	59.33	96.64	93.66	82.25	85.12	99.06

), the greatest confusion is between the fuel types FM5 forest litter layer beech–black pine and FM7 forest litter layer of beech, as well as FM3 forest litter layer of black pine and FM8 forest litter layer of Calabrian pine. Even though the underestimation of FM5 forest litter layer beech–black pine still remained (PA 59.33%, UA 87.30%), the accuracy results have improved compared to those of the SP model, with a significant increase in accuracy in the classification of FM4 Mediterranean grasslands type (PA = 82.99% and UA = 96.59%) and FM1 pseudomaquis (PA = 86.69% and UA = 81.89%), demonstrating the importance of topographic variables in the fuel type classification model.

3.1.3. Spectral Indices, Topography, and Backscattering Information Model (SPTS)

The SPTS model accuracies show a clear improvement over the SP model and a slight increase over the SPT model (OA = 92.76%) (Table 4). According to the confusion matrix (

Table 7. Confusion matrix for spectral indices, topographic variables, and backscattering information (SPTS) model. Values are given in percentage of validation points.

	Actual										UA (%)
Predicted	FM1	FM2	FM3	FM4	FM5	FM6	FM7	FM8	FM9	NV
FM1	88.88	1.83	2.12	0.19	2.51	0.00	0.11	3.00	4.24	0.29	84.44
FM2	1.05	87.67	1.54	0.09	0.28	0.00	1.81	0.45	3.55	0.00	88.56
FM3	0.35	0.37	80.79	0.38	2.51	0.09	0.28	8.76	0.19	0.03	84.8
FM4	0.00	0.00	0.10	85.16	0.00	0.52	0.34	0.00	0.51	0.07	97.09
FM5	0.00	0.00	0.10	0.00	63.23	1.98	0.28	0.00	0.00	0.00	88.67
FM6	0.00	0.00	0.10	0.09	8.36	96.99	0.40	0.00	0.00	0.00	96.66
FM7	0.00	3.84	1.06	0.57	21.17	0.43	94.05	0.00	3.29	0.00	89.64
FM8	4.47	0.18	11.97	0.09	0.00	0.00	0.00	85.77	0.06	0.16	85.58
FM9	1.84	6.03	0.39	1.04	1.67	0.00	2.72	0.07	84.99	0.22	88.41
NV	3.42	0.09	1.83	12.38	0.28	0.00	0.00	1.95	3.17	99.22	96.96
OA (%)	88.88	87.67	80.79	85.16	63.23	96.99	94.05	85.77	84.99	99.22

) and comparing to SPT models (Table 6), Sentinel-1 predictors improved the separability of the FM5 forest litter layer beech–black pine from the F7 forest litter layer of beech and of FM4 Mediterranean grasslands from NV non-forest areas. Moreover, the inclusion of backscattering information indicates an improvement in discrimination of the FM8 forest litter layer of Calabrian pine from FM1 pseudomaquis and NV non-forest-areas, of FM1 pseudomaquis from FM8 forest litter layer of Calabrian pine, and of FM3 forest litter layer of black pine from FM1 Pseudomaquis. However, confusion still remained in the FM5 forest litter layer beech–black pine, which was misclassified to other fuel types (PA = 59.33%) and mainly as the FM7 forest litter layer of beech. Overall, with the exception of FM5 Forest litter layer beech–black pine, the remaining categories were classified with PA and UA accuracies exceeding 80%.

In general, FM1 pseudomaquis, FM7 forest litter layer of beech, NV non-forest areas class presented overestimations (Figure 5), FM3 forest litter layer of black pine, FM4 Mediterranean grasslands, and FM5 forest litter layer beech–black pine were underestimated, albeit presenting high PAs and UAs, respectively.

3.2. Visual Inspection

For visual inspection of the maps, a representative area with different topography and tree species composition was selected. Figure 6 illustrates that the largest differences were obtained between SP and SPT or SPTS model. In this case of SP, the RF predictions also seemed to be noisier, underestimating the FM4 Mediterranean grasslands and overestimating FM1 pseudomaquis. The spatial patterns of fuel types are quite similar when using SPT or SPTS model; however, the separability of FM5 forest litter layer beech–black pine from the F7 forest litter layer of beech is improved on the SPTS model.

3.3. Variable Importance

The SP and SPT models were disregarded for subsequent analysis as they presented lower overall accuracies. In this sense, the SPTS model included 16 variables and was selected for further analysis of variable importance. The variable importance of SP and SPT models is given in Appendix A (Figure A1).

The importance of each variable was assessed using a mean decrease in accuracy (MDA) measure provided by the RF model. Figure 7 shows the variable importance ranking of the SPTS model based on MDA statistics. Relating to the MDA statistics (Figure 7), the most significant variables are the DEM and aspect, SVI and NDWI, and backscattering information, which were included in the seven most significant variables.

Boxes chart for backscatter (Figure 8) and spectral (Figure 9) profile of fuel type samples were developed to highlight the similarity of the fuel type classes.

4. Discussion

4.1. Fuel Type Classification Model’s Comparison

Wildfire is one of the main processes affecting the dynamics of Mediterranean forests [35]. An accurate and regularly updated description of the spatial distribution of wildland fuels is crucial for forest fuel management, fire management, firefighting, and mitigation of environmental and economic threats. Fuel mapping is essential for predicting fire behavior and managing vegetation, while wildfire behavior is partly determined by stand structure, composition, and vegetation characteristics. Currently, categorical fuel types corresponding to major forested vegetated/forested types are being used globally for fire research and management applications [86,87,88,89].

The strength of data fusion for forest type or fuel type classification has been successfully evaluated by many studies, but very few specifically addressed wildland fuel type mapping with the synergetic use of Sentinel-2, Sentinel-1 data, and topographic information. Most of these studies relied on the combined use of LiDAR with optical data [20,90]. The present work mapped vegetation/forest types implying categorical fuel types in a heterogeneous landscape of Northern Greece, combining three different datasets: Sentinel-2 spectral indices, topography, and Sentinel-1 backscattering information.

Optical data alone achieved less satisfactory results. The higher omission error noted in the case of FM4 Mediterranean grasslands, FM3 forest litter layer of black pine, and FM5 forest litter layer beech–black-pine is mainly related to the spectral similarity of these types with the NV non-forest areas, FM8 forest litter layer of Calabrian pine and FM6 forest litter layer of Norway spruce, and FM7 forest litter layer of beech, respectively.

Topographic variables included in the SPT model enhanced classification accuracy for all fuel types. Previous research related to vegetation patterns [34,91,92] and cover estimation [93] have also identified the influence of topography. In our case, the most noticeable improvement was observed for the FM3 forest litter layer of black pine, FM4 Mediterranean grasslands, and FM6 forest litter layer of Norway spruce. In Mediterranean areas, black pine generally colonizes areas at higher altitudes [94]; thus, the elevation predictor improved the FM3 forest litter layer of black pine separability from other fuel types on lower altitudes such as FM1 pseudomaquis, FM9 forest litter layer of deciduous shrubs, and NV non-forest areas. Topographic variables also supported the separability of FM4 Mediterranean grasslands from NV non-forest areas, as agricultural land occupies land on lower elevations. Similarly, the FM6 forest litter layer of Norway spruce gains classification accuracy due to topography-related variables since Norway spruce occupies higher altitude areas of the Rhodope Mountain range. However, misclassifications were noted again, especially in FM5 forest litter layer beech–black pine, due to spectral similarity with other classes.

The microwave energy of the C-band is mainly affected by the structural properties of a forest, such as the size, orientation, and spatial pattern of trees as a whole, and it scattered back to the sensor from the upper tree crown, branches, and leaves, and less from beneath [95]; thus, C-band backscatter signals are sensitive to forest biomass forest structure and vegetation water content [96,97]. Consequently, SAR predictors support the discrimination of forest/non-forest [98] areas as well different forest types of vegetation types with different water content and structure [29]. The synergetic use of optical and topographic variables with backscattering information slightly improved classification accuracy mostly of the FM5 forest litter layer beech–black pine, FM8 forest litter layer of Calabrian pine, FM1 pseudomaquis, FM4 Mediterranean grasslands, and FM3 forest litter layer of black pine. In the present study, FM5 forest litter layer beech–black pine was the most difficult class to map because, depending on the forest mix rate at the pixel level, the class shares spectral similarities with the FM7 forest litter layer of beech and FM3 forest litter layer of black pine [99]. However, Sentinel-1 variables contributed, even slightly, to separate these fuel types due to different canopy roughness and stand height between mixed and single species stands.

4.2. Variable Importance

With respect to variable importance ranking, the results demonstrate that all three datasets provide valuable information on fuel type properties; the MDA statistic indicates that not all variables contributed equally to the tree species classification.

Topographic metrics, mainly the aspect of the slope and altitude, have a considerable role in vegetation composition [100] and can distinguish different tree species [101]. The highest rank observed for elevation has been related in numerous studies to the altitudinal zonation of vegetation groups [32]. The high contribution of DEM was also confirmed by previous research on land cover classification [33,102] and tree species separability [101]. Aspect is presented as the second most important variable, as it is related to vegetation differences mainly due to differential exposure to the sun. Moreover, herb biomass significantly differed between aspects [100].

Concerning the slope, the south-facing slopes are drier than the north-facing slopes [103]. The wet and cool north-facing slopes are indirectly exposed to the sun and present low-intensity north-east rains during the winter season, thus holding more water and providing suitable conditions for the moisture-loving tree species [103]. Slope is also related to solar radiation across the surface, runoff, and gravity flow, all of which contribute to the presence of individual tree species as well as their number [104]; for instance, steep slopes are generally associated with broadleaves due to their developed root system [32], agricultural land and urban areas are usually located in relatively flat areas, and some tree species are likely to appear in sunny or shady-slope areas [105].

In regard to the highest-ranking spectral indices SVI and NDWI have a special interest in the classification of fuel types proving the importance of NIR_n2, RE₂, and SWIR spectral bands. Normalized difference water index (NDWI) as a common wetness index is sensitive to changes in the vegetation moisture content [106] and is often used in fire-related [107,108] and forest parameters studies [109,110]. Previous research proves that NDWI presents good spectral proxies for variations in vegetation density [111], and there is a correlation between NDWI and live fuel moisture content [25], facts that are associated with forest fuel type. Previous studies also use SVI for fuel type mapping [51]. The importance score noted for SVI supports the major role of the NIR (B8) and the red edge 2 (B6) bands of Sentinel imagery for forest fuel classification. SVI is a slight alteration to the traditional NDVI and considers a narrower waveband at the edge of the chlorophyll absorption feature (e.g., 740 nm), where the spectrum is characterized by lower absorption by chlorophyll, reducing the saturation effect. However, the reflectance still remains sensitive to chlorophyll absorption at its moderate-to-high values [112,113], and the sensitivity of the VI is enhanced to moderate–high vegetation densities [114], which is associated with forest fuel load and type.

Regarding backscattering information, Sentinel-1 variables improve the separability in classes that vary in terms of height (e.g., FM8 forest litter layer of Calabrian pine and FM3 forest litter layer of black pine from FM1 pseudomaquis and NV non-forest areas), or density (e.g., FM5 forest litter layer beech–black pine from FM7 forest litter layer of beech, and of FM4 Mediterranean grasslands from NV non-forest areas). Polarisation VV and VH present similar importance in terms of MDA; however, VV was rated slightly more important than VH in fuel type classification. VV backscatter is usually linked both to water content and surface roughness [115] and is dominated by direct contributions from the ground and the canopy [116]. VV polarization also performs better than VH polarization in land cover classification [117]. In contrast, in other studies, VH predictors have been proven to be more sensitive to tree/forest structure than VV predictors [95,118]; however, leaf-on SAR images used in the present study may influence the importance ranking of backscattering polarization because the leaves-on season both deciduous and evergreen contributes to crown volumetric scattering [115]. Ratio VV/VH was less important than the two polarisations since VV/VH reduces the dependence of the radar signal on soil moisture [119].

4.3. Limitations and Further Research Perspectives

The specific aim of this study was to investigate the synergistic use of Sentinel-2, Sentinel-1, and topography in classifying vegetation/forest types, as well as to explore the relative importance of spectral indices, backscattering, and topography information for an easy-to-update categorical fuel type mapping.

Our study considered static, structural, and single-month spectral characteristics of vegetation. Unfortunately, such variables restrict estimates of seasonal variations in phenologies. Vegetation phenology could provide important information on the successful discrimination of the fuel types [120] since phenological characteristics such as moisture content, leaf biomass, and chemical composition affect plant flammability. Future research should aim to investigate the potential of multi-temporal Sentinel-1 and Sentinel-2 datasets to develop a phenological fuel map considering the seasonality and inter-annual variability of forest fuel.

This study aimed to achieve a simple and robust model by using spectral indices, which have been used successfully for fuel mapping [51] and vegetation mapping [52,53] in previous studies. However, spectral data analysis is an important tool for classification problems [121], yet the choice of a suitable spectral index needs more consideration since it is contingent upon the specific application and environmental conditions [122]. Thorough research on the spectral profile of forest fuel types, as well as on spectral band importance in forest fuel mapping, could provide additional information to fuel type mapping, improve understanding of spectral indexes, and develop new effective spectral indexes.

A further point to be considered in future studies is the transferability of the model. The spatial variability of fuel types, along with uncertainties related to satellite imagery preprocessing, constrain model transferability between different areas. Overall, the developed regional model is currently specific to REMTH, but the general methodology is transferable to any other area with similar environmental characteristics. For improving the spatial transferability of fuel type classification models, future research might also evaluate the effectiveness of optical and SAR texture variables or features such as folded aspect [123] or vegetation management indicators that are independent of phenological stages. Finally, the development and delineation of a more detailed fuel type classification scheme, also considering successional stages and time since disturbance could facilitate fire behavior modeling, albeit it will be more challenging at extended scales requiring additional complementary information from LiDAR sensors [124].

The originality of this study methodology resides in the utility of Sentinel-1 data for the spatial distribution maps of wildland fuel in a Mediterranean landscape. The combination of information arising from backscattering data, optical data, and topographic metrics appears to be promising for accurate forest fuel mapping. The results obtained can also be used to guide the selection of suitable variables to develop an accurate fuel type classification model in a regional setting approach. Moreover, the developed models can be easily updated with newly available data, supporting decision-makers and disaster management with the necessary spatially distributed information on fuel types over the region.

5. Conclusions

Up-to-date mapping of fuel types is vital for wildland fire prevention. This study aimed to evaluate the synergetic use of three different datasets to develop accurate vegetation fuel maps and investigated the relative importance of spectral indices, backscattering information, and topography in fuel type mapping in a Mediterranean environment in Southern Greece.

In summary, the conclusions of this study are:

• The SPTS model yielded the highest overall accuracy (OA) of 92.76%, followed by SPT (OA = 91.92%) and SP (OA = 81.39%);
• The Sentinel-2 spectral indices produce accurate maps of forest fuel types;
• Topographic variables enhance model performance;
• Sentinel-1 data have a positive influence on the accuracy of forest fuel type classification;
• Variable importance measures highlighted the importance of the digital elevation model.

Overall, the results of this study provide evidence that encourages the synergistic use of Sentinel-2, Sentinel-1, and topography variables in fuel type mapping. Future work will involve the application of the generated fuel map as a key aspect of fire risk assessment across the region of East Macedonia and Thrace, providing automatic production of reliable results on fuel mapping and fire behavior at an operational level.