Agreement Index for Burned Area Mapping: Integration of Multiple Spectral Indices Using Sentinel-2 Satellite Images

Identifying fire-affected areas is of key importance to support post-fire management strategies and account for the environmental impact of fires. The availability of high spatial and temporal resolution optical satellite data enables the development of procedures for detailed and prompt post-fire mapping. This study proposes a novel approach for integrating multiple spectral indices to generate more accurate burned area maps by exploiting Sentinel-2 images. This approach aims to develop a procedure to combine multiple spectral indices using an adaptive thresholding method and proposes an agreement index to map the burned areas by optimizing omission and commission errors. The approach has been tested for the burned area classification of four study areas in Italy. The proposed agreement index combines multiple spectral indices to select the actual burned pixels, to balance the omission and commission errors, and to optimize the overall accuracy. The results showed the spectral indices singularly performed differently in the four study areas and that high levels of commission errors were achieved, especially for wildfires which occurred during the fall season (up to 0.93) Furthermore, the agreement index showed a good level of accuracy (minimum 0.65, maximum 0.96) for all the study areas, improving the performance compared to assessing the indices individually. This suggests the possibility of testing the methodology on a large set of wildfire cases in different environmental conditions to support the decision-making process. Exploiting the high resolution of optical satellite data, this work contributes to improving the production of detailed burned area maps, which could be integrated into operational services based on the use of Earth Observation products for burned area mapping to support the decision-making process.


Introduction
From the 1980s onwards, because of climate change conditions [1] and the increase in drought and land degradation phenomena [2], fires have become an important cause of land cover modifications around the world [3], affecting landscape patterns and functions [4][5][6], ecosystem processes [7,8], and air quality [9]. Identifying burnt areas provides key information on the extension of wildfires, the spatial location of the events, and the vegetation loss, which is useful to understand and quantify post-fire phenomena, such as air emissions, soil erosion, and vegetation recovery.
Earth Observation (EO) has proven to be a suitable technology for identifying areas affected by fire [10], from global [11,12] to local scales [13]. From this perspective, EO data have been used in the of medium-resolution optical satellite data analyzed with methodological approaches that make use of predefined sets of parameters and a constant number of SIs. The production of detailed burned area maps exploiting high-resolution optical satellite data, such as S2, and taking advantage of parameter tuning is yet to be accomplished. The operational service procedures could integrate the AIX approach to increase detailed wildfire mapping capacity for sites with different environmental conditions in all the fire seasons, balancing the errors, and optimizing the accuracy. This study aims to propose a procedure for burned area mapping through the combination of multiple SIs that i) identifies the best SI to be used to define SI thresholding based on the fire season and site conditions; ii) use an adaptive SI thresholding method which is flexible enough to be adapted to the local conditions, and iii) introduce a spatially explicit agreement index that can be consistently used for areas with different environmental conditions in all fire seasons.
After introducing the topic of burned area mapping through EO data and the use of SI, we demonstrate the effectiveness of the procedure in four case studies, consisting of distinct wildfires that occurred during different seasons in areas with different environmental characteristics. We describe the input data used and the method applied for analyzing S2 images and for computing the proposed index. The results section shows the SIs' performances and the final burned area map, and then we discuss the results of the analysis taking into account the advantages and shortcomings of the procedure. Finally, we summarize the derived conclusions.

Study Area
The procedure was developed at a local scale in Italy, and four burned areas were chosen from the archive of fire events made available by the Carabinieri Forestali (CUFAA, https://www.carabinieri.it/arma/oggi/organizzazione/organizzazione-per-la-tutela-forestaleambientale-e-agroalimentare): One fire was located in the municipality of Bussoleno, and one was in the municipality of Trivero, both of which were in the Piedmont region in northern Italy; one was in the municipality of Roccagorga in central Italy, and the last was in the municipality of Gravina in Puglia, located in southern Italy (Figure 1). The CUFAA wildfire database responds to Italian law (http://www.parlamento.it/parlam/leggi/00353l.htm) and reports information about the burned area surface extent and the burned vegetation. To test the methodology across various contexts, we selected study areas that differ in environmental conditions (climate, lithology, geomorphology, vegetation), land use, and fire season. In this study, the polygon representation of the four burned areas obtained from the fire events' archive was used in the procedure. The four polygons represented in Figure 1 are the reference maps used to detect the best SI, to identify the SI threshold, to compute the AIX score, and to validate the results.
Bussoleno and Trivero are situated in the Alpine mountain range. The fire events occurred in 2017 and 2015, respectively, both during the fall season and after an extended drought season. Roccagorga and Gravina in Puglia are situated in the Apennine mountain range, and the fire events occurred during a harsh summer dry season in 2017. The characteristics of the study areas and the information about the fire events are reported in Table 1.

Satellite Images
EO data used for the identification of burned areas in the study areas were acquired by Sentinel-2A and Sentinel-2B satellites. The Copernicus Sentinel-2 (S2) is a satellite mission carrying the Multi-Spectral Instrument (MSI) multispectral sensor with a high spatial resolution (10 m, 20 m and 60 m), high revisit capability (5 days with two satellites), an orbital swath width of 290 km and a moderately large band set (13 spectral bands) from the visible to shortwave infrared (https://sentinel.esa.int/web/sentinel/ missions/sentinel-2/instrument-payload/resolution-and-swath) (Table S1).    Pre-fire and post-fire S2 L2A products were downloaded from the Theia catalog (https://theia.cnes. fr/atdistrib/rocket/#/home) for each study area (Table 2). False-color maps generated from selected pre-fire and post-fire images are shown in Figure 2. The spectral bands available and the spatial resolution of 10 m for the visible-near-infrared range (VIS-NIR) and 20 m for red-edge-shortwave infrared range (red-edge-SWIR) enable detailed detection when applying SIs which are usually used for identifying burned areas (Table S1).

Methods
All the spectral bands of the images, distributed in the MUSCATE format as the bottom of the atmosphere (BOA) reflectances, which were orthorectified, terrain-flattened, and atmospherically corrected with MACCS-ATCOR Joint Algorithm (MAJA) [54,55], were processed for spatial resampling at 20 m and for the masking of topographic shadows. The images correspond to the bottom of the atmosphere reflectances, which were corrected for topographic effects to avoid the spectral variability originating from the local topography. The areas corresponding to topographic shadows (identified by product quality flags representing the topographic mask and low-sun mask, with the latter considering the pixel aspect, slope, and sun zenith angle) were not considered in the analysis and, therefore, also excluded from the reference maps. In this study, only the Trivero and Bussoleno areas showed topographic shadows and needed masking.
Resampled and masked images were used to compute SI before and after the fire event. Five SIs were used to recognize burned areas, which were chosen to cover the VIS/red-edge/NIR/SWIR spectral ranges (Table 3): Table 3. Spectral indices adopted, equations with Sentinel-2 Multispectral Instrument (MSI) bands, and bibliographic references.

Spectral Index Equation Reference
NBR The original formulae for the Normalized Burn Ratio (NBR) and Normalized Difference Vegetation Index (NDVI) have been modified to make the behavior of the indices consistent: The values of all the indices increased if the pixel was affected by fire. The NDVI formula uses B8A instead of B8 since it provides more precise spectral information.   Table 2.
Image differencing was performed [61] through the arithmetic difference between post-fire and pre-fire index values (∆ post-pre images) to detect abrupt changes and retrieve burned pixels. Regarding the difference SIs (i.e., dNBR, dNBR2, Mid-Infrared Bispectral Index (dMIRBI), Burned Area Index for Sentinel-2 (dBAIS2), and dNDVI), a set of thresholds was applied to obtain a first set of burned/not burned maps.
Index thresholds were defined, starting from a spectral sensitivity analysis to assess the efficacy of SI to separate burned and unburned land. The separability index (M) was applied as follows (Equation (1)) [33,62,63]: where µb and µu are the mean values of the difference SIs for burned and unburned pixels, and σb and σu are the standard deviations of the difference SIs. M values higher than 1 indicate a good separability between burned and unburned land [62]. To calculate the M index, the analysis was performed on an area of around 10 times the surface of the burned polygon defined in the reference map.
The threshold values were identified, starting from the polygons obtained from the difference SI with the highest M value. The median SI values of the pixels on the polygon contour and adjacent to polygon contour (inside and outside) of the reference map were calculated ( Figure 3). The empirical threshold value of the difference SI with the highest M value chosen by an expert is a value derived by the statistic distribution of the pixels inside, outside, and on the contour of the reference map.
Remote Sens. 2020, X FOR PEER REVIEW 8 to polygon contour (inside and outside) of the reference map were calculated ( Figure 3). The empirical threshold value of the difference SI with the highest M value chosen by an expert is a value derived by the statistic distribution of the pixels inside, outside, and on the contour of the reference map. Threshold values for the other SIs were calculated through an adaptive procedure by applying the regression equation of each difference SI vs. the difference SI with the highest M value, considered as the best SI. After that, the binary images classified as burned (pixels above the threshold value) and unburned (pixels below the threshold value) could be defined. The accuracy assessment for each difference SI was calculated using confusion matrices comparing the resulting maps against the reference maps. Omission and commission error estimation was performed.
To define the final burned area map, an Agreement IndeX (AIX) was proposed by combining the SIs as described in Equation (2): where P is a burned pixel, SI is the spectral index, and N is the number of SIs. The optimum combination of SIs to use to obtain the final map is the one that optimizes the error of each difference SI (both omission and commission) and the overall classification accuracy.
For each combination of n SIs corresponding to a different AIX value, the overall accuracy and omission, commission, and total errors were calculated from a comparison with the reference map. In addition, an Agreement Index Score (AIS) was computed according to Equation (3): where OE is the omission error, CE is the commission error, and OA is the overall classification accuracy. The AIS serves to identify the AIX value, which corresponds to a combination of n SIs that balances the errors and optimizes the accuracy. The final burned area map is the result of the pixels for which the combination of at least n SIs out of N identifies them as burned pixels. Figure 4 shows the workflow of the main phases of the procedure: (1) pre-processing and (2) the Threshold values for the other SIs were calculated through an adaptive procedure by applying the regression equation of each difference SI vs. the difference SI with the highest M value, considered as the best SI. After that, the binary images classified as burned (pixels above the threshold value) and unburned (pixels below the threshold value) could be defined. The accuracy assessment for each difference SI was calculated using confusion matrices comparing the resulting maps against the reference maps. Omission and commission error estimation was performed.
To define the final burned area map, an Agreement IndeX (AIX) was proposed by combining the SIs as described in Equation (2): where P is a burned pixel, SI is the spectral index, and N is the number of SIs. The optimum combination of SIs to use to obtain the final map is the one that optimizes the error of each difference SI (both omission and commission) and the overall classification accuracy. For each combination of n SIs corresponding to a different AIX value, the overall accuracy and omission, commission, and total errors were calculated from a comparison with the reference map. In addition, an Agreement Index Score (AIS) was computed according to Equation (3): where OE is the omission error, CE is the commission error, and OA is the overall classification accuracy. The AIS serves to identify the AIX value, which corresponds to a combination of n SIs that balances the errors and optimizes the accuracy. The final burned area map is the result of the pixels for which the combination of at least n SIs out of N identifies them as burned pixels. Figure 4 shows the workflow of the main phases of the procedure: (1) pre-processing and (2) the identification of burned areas. All the analyses were performed using R-CRAN and Quantum GIS 3.4 software.
Remote Sens. 2020, X FOR PEER REVIEW 9

Results
Regarding the separability index (M), dNBR2 was the index that showed the highest discriminant power for the two study areas of Bussoleno and Gravina in Puglia and dNBR was the index that showed the highest discriminant power for the Roccagorga study area, whereas the

Results
Regarding the separability index (M), dNBR2 was the index that showed the highest discriminant power for the two study areas of Bussoleno and Gravina in Puglia and dNBR was the index that showed the highest discriminant power for the Roccagorga study area, whereas the dMIRBI index showed the highest M value for the Trivero study area. The resulting separability index (M) values for all the study areas are reported in Table 4; the highest values for each study area are shown in bold. Having obtained the highest M value, dNBR and dNBR2 were selected as the best SIs to serve as reference for the summer dry season study cases (Roccagorga and Gravina in Puglia), while dMIRBI and dNBR2 were selected as reference SIs for the fall dry season (Trivero and Bussoleno). These SIs were used to identify the empirical threshold for the study areas, and the median values of the pixels on the boundary and adjacent to the boundary (inside and outside) of the reference maps were calculated (Table 5). SI thresholds, estimated using regression equations in the adaptive thresholding processing step, were computed starting from the thresholds defined for the best SI and reported in Table 5. Table 6 reports the SIs' performances by providing the estimated threshold values of all the difference SIs and the omission and commission errors and overall accuracies for each burned area map, with the highest values shown in bold. Figures S1-S4 show the ∆ post-pre and the classified burned/not burned images of each difference SI for the four study areas. Regarding the Bussoleno burned area maps, the highest overall accuracies of 0.57 were achieved by dBNR2 and dMIRBI, whereas for the Trivero study case, dMIRBI and dNBR2 achieved the highest overall accuracies (0.65 and 0.61, respectively). Both areas showed a high commission error for all the indices, especially Trivero. Regarding the Roccagorga burned areas maps, an overall accuracy of 0.96 was achieved by dNBR and dNDVI, while the dMIRBI results showed an overestimation of the burned surfaces, and dNBR2 showed the lowest omission error. The study case of Gravina in Puglia showed high overall accuracy for all the indices (dNBR, dNBR2, and dBAIS2 achieved 0.93) as well as the lowest commission errors; in particular, dNDVI achieved 0.17. Table 7 reports the classification performances combining a different number of SIs. For all the study areas, the increment of the AIX value corresponded to an increase in the omission error and a reduction in the commission error. In contrast, a higher overall accuracy value can result from the combination of different SIs. The AIS has the capability to combine both the total error (omission error and commission error) and the overall accuracy. Three study areas out of four resulted in a higher AIS with the combination of four SIs. The Gravina di Puglia study area reported a higher AIS (1.43) when only one SI was considered. The resulting number of SI combinations balanced the omission, commission, and total errors.   Figure 5 shows the spatial distribution of the AIX classes and the final maps of burned areas; that is, the pixels which the combination of n SIs identified as belonging to the "burned" class.
Remote Sens. 2020, X FOR PEER REVIEW 12 Figure 5. Spatial distribution of the Agreement IndeX (AIX) and maps of burned areas for the four study cases.

Discussion
In this research, we compared various case studies in four areas with diverse site conditions to evaluate the performances of different SIs for the burned area mapping of wildfires, which occurred in different seasons and under different environmental conditions in terms of climate, geomorphology, and vegetation cover. In fact, two wildfires occurred during the winter season in areas with vegetation cover characterized by grasslands and deciduous forest. Other considered wildfires occurred during the summer dry season: The Roccagorga study case mostly involved a natural area also covered by shrublands and grasslands, while Gravina in Puglia could be considered an ideal study case, where the wildfire burned a natural forest with both deciduous and evergreen tree species. The wildfires which occurred during the fall and winter seasons in the Alps also highlighted the need for a proper topographic shadows mask [17], with the capacity to remove dark pixels from further processing and reduce the commission error in the classification phase.
Research studies generally evaluate the performances of SIs on ideal case studies [44][45][46]-i.e., wildfires occurring during the summer dry season, when the vegetation is still photosynthetically active-but the results could differ from other case studies, especially for SIs that rely on vegetation greenness (e.g., NDVI), and that could not perform well if a fire event occurs at the end of the senescence phase. Our results demonstrate that SIs have different performances under different environmental conditions. In some circumstances, the same SI performs best in one study area and worst in another study case. Vegetation SIs, such as dNDVI and dBAIS2, have a reduced detection capability in deciduous tree species during the fall and winter dry seasons due to the poor photosynthetic activity. Indeed, dNDVI showed a high omission error, and both indices had large commission errors for Bussoleno and Trivero study cases since the burned area was mainly represented by deciduous forests at the end of the senescence phase. dBAIS2 generally has a discrete detection capability for grasslands, while it has a good detection capability for Mediterranean mixed vegetation (e.g., the study case of Gravina in Puglia). dNBR was revealed to be a good SI; in particular, it obtained the best performances for herbaceous environments. dMIRBI exhibited good performance during the fall dry season, whereas dNBR2 showed generally good performance in all the cases and was shown to be suitable for fall and winter seasons.
The integration of SIs was performed on a pixel base with the numerical support of the AIX. This is a spatially explicit methodology that could be adapted to different site conditions as appropriate. The two considered summer study cases reported a high overall accuracy, exceeding the value of 0.9 for all SI combinations. The results obtained for the Gravina in Puglia area-here considered an ideal study case-showed high AIS, high overall accuracy values, and low commission error values, even using a single SI. Since the commission error generated by each SI was very low, burned area identification for summer wildfires could benefit less from SI integration than a fire that occurs in other seasons. On the contrary, fall season study cases showed the need for an in-depth investigation to avoid a different phenology leading to erroneous evaluations. The support of the AIX allowed us to obtain higher overall accuracies and, at the same time, to lower the omission errors compared to considering the SIs separately. Another research study [33] attempted to integrate multiple SIs: After performing an SI separability analysis, SI values were reduced to a common domain, converted to positive and negative evidence [38], and integrated using fuzzy membership functions; finally, burned area pixels were identified, applying different operators, such as ordered weighted averaging [51]. Alternatively, research studies using multiple SIs adopted a set of empirical conditional statements, also defined as decision rules, to identify burned area maps [31] even with the support of a temporal harmonizer algorithm [64], or use distribution-modeling algorithms trained by hyperspectral indexes and hotspots to model and map burned areas [47]. Compared to other methodologies described in the literature, the approach presented in this study has the strengths of significantly reducing the use of pre-defined set of parameters, thanks to the adaptive thresholding, and rejecting the use of SIs whose mapping performances are not adequate for specific site conditions, thanks to the AIX.
The proposed procedure first identifies the most suitable SI to be used as a reference, then calculates thresholds for the other SIs based on regression models, and computes a single burned area mask from each SI. Later, all the burned area masks are combined to calculate the AIX, which is finally used to map the burned areas by optimizing omission and commission errors. The availability of reference maps supports the analysis by providing information for the identification of the best SI to be used for the index separability analysis, allowing us to determine the best SI threshold and calculate the AIS and classification accuracy metrics to support the selection of the number of SIs to be used for the final burned area map generation and for validation purposes. The procedure could be used in a training phase to evaluate the performances of SIs on a large set of wildfire cases for which a reference map is available. One may think that the pre-fire and post-fire images could be selected during the vegetation growing season, without the burn date; this could reduce, on the one hand, the commission error due to the slight variation in pre-fire and post-fire spectral signatures, but, on the other hand, this would consistently increase the omission error due to the vegetation restoration process [21]. It is, therefore, necessary to select the best-performing SIs, which may depend on the vegetative status, and to integrate all of them to obtain the most accurate mapping, balancing the omission and commission errors, and optimizing the overall accuracy. Compared to the methodologies available in the literature, the proposed approach allowed us to use a variable number of SIs, and it benefited from the best-performing SIs for the site-specific conditions, with the AIX selecting the number of SIs on a single pixel base.
In the procedure training phase, a reference map can be used to calculate the separability index and identify the SI that results in the most accurate burned area mapping and its threshold value. The identification of the best SI for specific site conditions allowed us to use it as the reference SI when using adaptive SI thresholding. Once the training phase is completed, the site conditions related to the time of the season and the phenology may be used to establish which SI and threshold value should be used, opening the way to a newer paradigm for the development of burned area mapping operational services. Here, the extensively adopted methodology for burned area mapping operational services was based on a defined set of SIs and the threshold to be used, independent of the site conditions [11,12]. Alternatively, operational services may use an automatic thresholding algorithm developed to exploit the statistical bimodal distribution analysis [49]. Further analysis is still required to create a dictionary that suggests the SI and thresholds to be used under the different site conditions, fire regimes, and ecosystems, to move the proposed approach from training mode to operational mode. The availability of high-resolution land cover maps as ancillary information to identify vegetation cover type prior to a wildfire taking place could help to identify the parameters from the dictionary (SI, threshold, and AIX) to be used for the burned area mapping. Alternatively, on-demand high-resolution land cover mapping operational services can be adopted. The information regarding the land cover (especially dominant tree species) and vegetation status in the pre-fire stage can also integrate the burned area maps with indispensable information for the estimation of atmospheric emissions during wildfires [65].

Conclusions
In recent years, the high spatial and temporal resolution of EO data, such as the S2 satellite, has enabled progress in the development of detailed and timely post-fire mapping. The availability of detailed and prompt mapping of post-fire areas is important to support post-hazard management strategies, including soil erosion estimation, air emission control, and environmental loss recovery.
SIs are widely used for burned area mapping, although the choice of the best SIs, the threshold values and their performances are often assessed without taking into account the fire season, the local geomorphology, the vegetation types and their phenological stage.
This study proposed an approach for improving the detection of burned land using satellite optical multispectral data, taking advantage of different spectral bands through the combination of multiple SIs. The results are encouraging considering the high level of accuracy of the final maps achieved, reducing the problems associated with under and overestimation. Moreover, the transferability of the procedure in different geographical contexts was verified by testing on different study cases, with satisfactory results, which highlight the potential for this method to be applied over larger areas.
The adaptive SI thresholding method and the proposed AIX could be applied to different site conditions and could represent a tool to link SIs, thresholds, and environmental conditions. The possibility to test the methodology on a large set of wildfire cases could allow a catalog to be created that encompasses the aspects mentioned above which are representative of the existing case studies and that could support operational services based on the use of EO products for burned area mapping in order to support the decision-making process.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/11/1862/s1, Table S1: Spectral characteristics of the Sentinel-2 MSI sensor, Figure S1: ∆post-pre images and burned/not burned images of each difference SI for the Bussoleno study area, Figure S2: ∆post-pre images and burned/not burned images of each difference SI for the Trivero study area, Figure S3: ∆post-pre images and burned/not burned images of each difference SI for the Roccagorga study area, Figure S4: ∆post-pre images and burned/not burned images of each difference SI for the Gravina in Puglia study area.