BAMS: A Tool for Supervised Burned Area Mapping Using Landsat Data

A new supervised burned area mapping software named BAMS (Burned Area Mapping Software) is presented in this paper. The tool was built from standard ArcGISTM libraries. It computes several of the spectral indexes most commonly used in burned area detection and implements a two-phase supervised strategy to map areas burned between two Landsat multitemporal images. The only input required from the user is the visual delimitation of a few burned areas, from which burned perimeters are extracted. After the discrimination of burned patches, the user can visually assess the results, and iteratively select additional sampling burned areas to improve the extent of the burned patches. The final result of the BAMS program is a polygon vector layer containing three categories: (a) burned perimeters, (b) unburned areas, and (c) non-observed areas. The latter refer to clouds or sensor observation OPEN ACCESS Remote Sens. 2014, 6 12361 errors. Outputs of the BAMS code meet the requirements of file formats and structure of standard validation protocols. This paper presents the tool’s structure and technical basis. The program has been tested in six areas located in the United States, for various ecosystems and land covers, and then compared against the National Monitoring Trends in Burn Severity (MTBS) Burned Area Boundaries Dataset.


Introduction
Fires are one of the foremost factors of land surface disturbance in diverse ecosystems causing severe economic, ecologic, and atmospheric effects [1].Vegetation burning is a global scale environmental phenomenon that affects most global biomes like grasslands/savannah (tropical and subtropical), forest (Mediterranean, temperate, and boreal), and agricultural lands.In the last few decades satellite data have been used extensively to monitor biomass burning at the regional to global scale in order to overcome the problems related to large areas affected by fire, the very dynamic nature of the process, and the low accessibility of many key regions [2,3].
In the past few years several global burned area thematic products derived from low-to moderate-resolution satellite data based on various automatic publicly available algorithms have been released [4][5][6][7][8], but no global consensus exists yet on the validation procedures and error characterization methods [9][10][11].In spite of this, most papers dealing with validation of global burned area products recommend using medium-resolution images (Landsat TM or similar) as reference information [5,[10][11][12], since national statistics are usually either not available or not reliable, making field work unfeasible.In fact, the same recommendation is given by the standard protocol of the Committee on Earth Observation Satellites (CEOS) Land Product Validation (LPV) Subgroup for generating reference data [13].
The capability of Landsat TM and ETM+ data to provide information about burned scars has been widely recognized in scientific literature [14][15][16][17][18][19], due to its spatial coverage (185 × 185 km), medium spatial resolution (30 m for reflective bands), multispectral characteristics (covering the most important spectral areas for burned area mapping with one band in the near infrared and two in the shortwave infrared), good temporal range (TM sensor and equivalent data have been available since 1982), and temporal resolution (16 days) [18].Moreover, in 2008 the United States Geological Survey (USGS) made its newly acquired and archived Landsat data freely available.This trend-improving accessibility to medium-resolution images is being followed by other space agencies like the ESA or the Centre National d'Etudes Spatiales (CNES), which provide a window into the past and ease monitoring and modeling of global land cover and ecological change [20].
Mapping burned areas using Landsat data is a complex task: burned areas may be confused spectrally with other phenomena due to the spectral signature of burned vegetation, which varies as a function of factors like fire behavior, pre-fire surface properties, and time elapsed since the burning happened [21].Although automatic Landsat classification and thresholding techniques provide useful information on the spatial extent and characteristics of the burned patches [22][23][24][25], different authors have recognized the difficulty in developing a global algorithm that can cope with the great variety of global burning conditions: old and new burns, spatially fragmented small patches, low severity fires, etc. [12,26].The variety of characteristics of the burned areas and the confusion with other spectrally similar areas (flooded regions, dark soils, etc.) are behind the omission and commission errors of all burned area products [27].
The novel methodology presented in this paper combines the speed and consistency of the automatic detection algorithms (based on a two-step thresholding process of various burned area indexes) [22] with the power of pattern recognition techniques, which are more easily adaptable to local conditions.The resulting methodology consists of an interactive process where the user must assess the result of each phase and decide whether more supervising effort is necessary.This methodology has been implemented using the most popular GIS software (ArcGIS), as it uses the development framework behind it (ArcObjects) to develop both the visual interface and make image analysis computational tasks.The program is named Burned Area Mapping Software (BAMS) and it was mainly designed to generate burned area reference maps using Landsat TM, ETM+, and OLI-TIRS data, although it will easily be adjusted to other sensors.

BAMS Program Flow
The BAMS code has two processing options, one for a pair of Landsat scenes (one considered pre-fire and the other post-fire) and another option that previously creates temporal composites from a set of Landsat scenes (maximizing NDVI for pre-fire state and maximizing NBR for post-fire state) (Figure 1):

Generation of Reflectances
The program can readily access two data input sources: the USGS Landsat Terrain Correction (Level 1T) GeoTiff format [28] and the USGS Landsat Surface Reflectance Climate Data Record (SRCDR) [29] HDF format.In the case of the former, BAMS converts Raw Digital Numbers (DN) into at-sensor radiance and then into exoatmospheric Top of Atmosphere (TOA) reflectance using the simplified reflectance equation.For Landsat TM and ETM+, the following equation is used: where Lλ = Grescale Qcal + Brescale ρλ = Exoatmospheric Top of Atmosphere reflectance (TOA) Lλ = Spectral radiance at the sensors aperture d = Earth-Sun distance ESUNλ = Mean exoatmospheric solar irradiance θs = Solar zenith angle Grescale = Band-specific rescaling gain factor from the metadata Qcal = Quantized calibrated pixel value (DN) Brescale = Band-specific rescaling bias factor from the metadata Reflectance for OLI data is as follows: where: ρλ = TOA planetary reflectance, without correction for solar angle.Mρ = Band-specific multiplicative rescaling factor from the metadata Aρ = Band-specific additive rescaling factor from the metadata Qcal = Quantized and calibrated standard product pixel values (DN) θs = Solar zenith angle The Sun Zenith Angle (θs) and the band's specific gain and offset calibration parameters are extracted from the USGS metadata file (MTL.txt) and the Earth-Sun distance and the Solar Irradiance are as published in [30] (only required for TM and ETM+).The SRCLR is still in development; however, it will become more relevant in the future, when climate records move to higher spatial resolutions.For this format, no pre-processing is required as TOA reflectances are included in the product.
In order to exclude areas not observed, a simple cloud detection mask is implemented in the code (areas with reflectance values over 0.5 in the blue spectral band or brightness temperature below 10 °C).In the case of the Landsat 8 data, the Band Quality Assessment Band [31] is used instead, masking out areas labelled as cloud, cirrus, and snow/ice.Similarly, the program reads the Fmask variable [32] in the Landsat SRCDR format to mask out non-land areas.

Computation of Burned Area Spectral Indexes
Burned areas are characterized by deposits of char and ash, removal of vegetation cover, and fuel, as well as exposure of the underlying soil.However, the magnitude and direction of spectral changes caused by charcoal and ash deposition depend on the type and condition of the vegetation prior to burning and the degree of combustion [33].

Temporal Composites
BAMS makes it possible to produce temporal composites for users who require burned area information for more than two periods (to reconstruct fire history with multiple Landsat scenes, for example).Two temporal composites are created in the process: one minimizes the NBR index, which aims to identify the most affected burned areas observed at each time frame of the series.This criterion will create the post-fire image of the time series.The second criterion maximizes the NDVI index, and aims to identify the time framework when each pixel is less affected by fire.The final composite will be considered as the pre-fire image of the time series.

Burned Area Mapping Supervised Methodology
BAMS follows a two-phase burned area strategy, which has produced good results for mapping of burned areas with low and medium spatial resolution [22,23,26,39,40].This strategy aims at keeping a balance between commission and omission errors: In the first phase the goal is to reduce the commission errors by using strict criteria, so that only the more clearly burned pixels are retained (seed pixels), even at the cost of omitting some burned pixels within each burned patch.The second phase analyzes the vicinity of the seed pixels, applies a more flexible criterion, and accepts as burned those neighboring pixels with spectral characteristics similar to the seeds.This phase progressively increases the burned area until the whole burned patch is covered and aims at reducing the omission errors.
The two-phase strategy implemented by this tool is simple (see Figure 2).To start with, two Boolean rasters are generated, one for the seeds (Figure 2a, in red) and another for those pixels that fulfill the second-stage criteria (Figure 2b, in orange).From the second-stage criteria raster, groups are created from pixels connected by any of the eight neighboring sides (right, left, above, below, and diagonals).Therefore only those groups that intersected with the seeds are retained (Figure 2c,d).
A multi-index approach is taken for the first and second phase.It has shown to be effective in Landsat automatic methodologies [22,23].BAMS code makes use of 10 spectral variables (the post-fire BAIM, NBR, MIRBI, GEMI, and NDVI indexes, plus the temporal differences of those five indexes).For the spectral indexes NDVI, GEMI, and NBR, the fire decreases the pre-fire value, so the tool sets a maximum value (the lower the value, the higher the probability of being burned).For BAIM and MIRBI, where the fire increases their pre-fire value, a minimum value is set (a higher value indicates higher probability of being burned).Two independent sets of criteria are defined, one for the seed phase (strict criteria) and another for the second phase (more relaxed).To avoid the impact of isolated pixels, seeds with less than two pixels in the immediate eight neighborhoods are removed.The final burned area is obtained by keeping only burned patches generated by the second-phase raster-intersecting seeds selected in the first phase.Once this raster layer is obtained, it is transformed to an ArcGIS shape vector format.In doing so, an aggregation process may be applied to join the resulting perimeters within a neighborhood of 100 meters, although when using temporal composites for the input of the burned area algorithm, this option has to be applied carefully to avoid different fire scars of neighboring events, occurring on different dates within the year, being considered as one fire event.The tool also has the option to remove holes for those users who want to produce simpler perimeter polygons that contain the boundary of the fire only, without any internal island.The crucial point of the algorithm is the selection of threshold values for each of the input spectral indices.These values are extracted from the training polygons identified visually by the user from the post-and/or pre-fire color composites.In opposition to traditional supervised classification methodologies, where both burned and unburned training areas are required, BAMS only needs burned training areas, which makes the training process easier and faster than previous approaches, because the unburned category is always more heterogeneous due to the higher spectral diversity of unburned land covers.Two different burned training polygons have to be sampled to set the threshold values, one for the seeds and the other for the second phase.For each, the minimum of all the samples are retained for MIRBI and BAIM variables, whereas the maximum values are extracted for the rest (NDVI, GEMI, and NBR).It is very important to avoid false positives in this phase, since the burned area thresholds will be set with these values and the results will then be less accurate.
The thresholds extracted from the two user-defined training polygons are applied with the logical operation "AND".The result (after applying the two-phase algorithm described in Section 2.2) should be assessed before deciding whether more training areas are needed to complete the burned cartography, by extracting new thresholds from different iterations.

Batch Process
If the user needs to process a large number of images, the program could also be run in a batch mode.This option makes it possible to run the process automatically, using a set of thresholds previously defined in the supervised mode with different images.

BAMS Result
The output result of BAMS is a shapefile layer that follows the protocol for BA reference sites defined by the European Space Agency's (ESA) fire_cci project [41] which in turn agrees with the CEOS Cal-Val working group [42].This protocol defines a "Category" attribute that records the areas as burned (class 1), not observed (class 2), or unburned (class 3).These unobserved areas are commonly covered by clouds or by sensor missing data (as caused by the ETM+ SLC-off problems).As the effects of the fire may be observed for some time after fire occurrence, the "Pre-date" and the "Post-date" attributes (acquisition dates of the pre-fire and post-fire images used, respectively) are used to describe the time period over which the mapped burns occurred.Finally, the type of Landsat sensor used (Landsat 4, Landsat 5, Landsat 7, or Landsat 8) and the path and row are also recorded with the "PostImg" and "PreImg" attributes.

Methodology
The BAMS tool has been run in six test areas (see details in Table 1) to evaluate its performance for operational burned area mapping.All the scenes involved have been processed at the same time, while the burned area processing has been done separately for each path-row.These test areas are located in diverse ecosystems within the United States (Great Plains, Eastern Temperate Forest, Mediterranean California, Northwest Forested Mountains, and Taiga) (Figure 3).Land cover affected by the fire within these sites includes grassland, shrubland, deciduous, and evergreen forest.Burned path sizes are also very diverse (from approximately 30,000 ha to almost 250,000 ha).All the scenes have been downloaded from the Earth Explorer Web Page [43] and correspond to L1T level.In four of the six test areas, a Landsat 5 or Landsat 7 pair of scenes with a gap of 2-3 months between them were selected; the 72-15 scene (pre-fire image) was from a year before (no cloud/snow-free scenes were available closer in time).For the 30-35 scene, three post-fire images were used to obtain a temporal composite in order to resolve the missing data related to the Scan Line Corrector Failure (SLC-off) of the Landsat 7 ETM+ data.
In this test exercise, the BAMS program was run with all spectral variables (post-fire and temporal changes of the BAIM, NBR, MIRBI, GEMI, and NDVI variables), including the aggregation process and removal of islands.The results obtained with BAMS have been compared with the National MTBS Burned Area Boundaries Dataset [44] which was generated by visual on-screen digitization of multitemporal NBR of Landsat data.Even though this is not strictly a validation of our results (since the MTBS may also include errors), it provides some insights into the accuracy assessment, as those perimeters are routinely used for fire management purposes.Matrix indicators of coincidences between the two results are computed.The differences are considered as underestimation (percent of area not detected by BAMS with respect to the total area detected by MTBS) and overestimation (percent of burned area detected by BAMS and not by the MTBS with respect to the total area detected as burned by BAMS).We have also computed the Kappa values that measure the relationship between the beyond chance agreement and expected disagreement [45].
Also, the number of iterations required in each test site has been noted, equal to the number of times the user runs the algorithm after redefining the samples in case commission errors are detected (thresholds are too relaxed), or after adding new samples if omission errors are observed (thresholds are too strict) until a visual satisfactory result was reached.

Results
Less than four iterations were required for four out of the five test areas, with an average processing time of 10 min for each test area.A computer with an Intel Xeon 2.53 Ghz processor and 4GB of RAM was used.For the 30-35 test area, 11 iterations were necessary, due to a greater spectral heterogeneity of the burned signal, the temporal compositing, and the longer interval between the pre-and post-fire images (more than four months).For spectrally homogeneous burned areas the mapping was faster, as a lower number of training samples was required to detect the signal's variability.The main statistics of the tests carried out are shown in Table 2.In terms of the comparison of the BAMS results with the MTBS perimeters, the largest disagreements were found for the 18-34 scene with 18.6% of underestimated area and 13% overestimated area.These disagreements were clearly related to the topographic shadows.The fires took place in November and the low sun elevation angle (below 30°) resulted in shaded areas that obscured the burned signal, thus complicating the detection of the whole perimeter (see Figure 4).Moreover, the shadowed areas broke the connectivity of the burned signal, making it more difficult to obtain completed perimeters.
The 46-31 scene was also found to have high underestimated area (11.0%) due to the weaker burned signal and fragmentation in the left border of the largest fire that made it difficult to reach until the validation border (see Figure 5).It is remarkable that a high overestimation disagreement was found in the 30-35 scene (71%), which is related to agricultural areas with very dark soils, although we noticed that some of these regions were in fact burned (Figure 6) but not included in the MTBS dataset.This is a common problem when comparing burned area maps with official fire perimeters, as most administrations only take into account burned natural areas, but no crops.

Discussion
BAMS was created due to the need to quickly produce burned area reference maps using Landsat data in the context of the ESA's fire_cci global burned area product validation [41].In this project, 147.994 burned patches that affected 126.180 km 2 in several ecosystems were mapped after processing 242 Landsat pairs.The scientific team that is developing the USGS Landsat Burned Area Essential Climate Variable Product [46] is evaluating the software to complete the ESA's fire_cci validation dataset for assessing the accuracy of their product.In addition, more than 25 fire researchers have also used the tool for burned area mapping in diverse parts of the world so far, supporting the idea that a tool to speed up the way to obtain this cartography was needed.Some algorithms have been assessed for automatic or semi-automatic burned area mapping with Landsat data [23][24][25], but these are difficult to implement and the software has not been released in an operative way for the public.
One of the advantages of BAMS is the ease with which Landsat data is processed: the calculation of absolute reflectance values, the most used burned spectral indexes, and the multitemporal substractions directly from EarthExplorer L1T data format.BAMS can also read the more recent Climate Data Record product that is already processed and has the advantage of using the Moderate Resolution Imaging Spectroradiometer (MODIS) atmospheric correction routines.It has the computation of the Fmask masking grid integrated, something especially important when generating validation perimeters where the detection of clouds and missing data to identify nonvalid areas is crucial for a quantitative validation.The ability to do temporal compositing also provides some advantages.It makes it possible to work with LE7 SLC-off data that has gaps with no data (Figure 7), and also to obtain complete perimeters covered by clouds in individual scenes.However, the criterion employed for temporal compositing requires further research, because minimizing the NBR value may not be the best criterion for some types of vegetation, cloud cover, and illumination conditions.The supervising strategy of BAMS requires good knowledge of the methodology followed by the tool.The seeds have to be established in areas where there is a strong burned signal, to avoid spectral confusion with other land covers such as cloud, topographic shadows, dark soils, or wetlands.It is recommended that the seeds be located in small and quite homogeneous areas.In following iterations, more seeds may be set if the entire burned area is not detected.For the second phase, training polygons bigger than those used for the first stage are recommended to take the maximum spectral post and multitemporal variability, in areas that show low severity at the borders of the fires or within unburned islands.As the thresholds applied are the same for the whole scene, it is recommended to find the burned areas less affected by the fire in this phase and thus to obtain thresholds low or high enough to map the entire burned extent in fewer iterations.For both cases, it is crucial to ensure that the samples are not burned already in the pre-fire scene, or are not covered by clouds, because that will generate weaker thresholds that will result in high commission errors.The best performance of BAMS is achieved where the fire has severely burned areas (to establish the seeds) and spatially not fragmented in order to have connectivity between the seeds and the second-phase result.The joint use of the 10 variables, BAIM, NBR, MIRBI, GEMI, and NDVI (post-fire and multitemporal post-and pre-subtraction) avoids most of the confusion with unburned areas.Some typical problems have been observed in the test areas, however, including confusion with croplands (these areas have to be manually removed) [22,26] and effects related to the obscuration of the signal due to topographic shadows.If a sample is located within a burned shadow, all the shadows of the scene are detected; thus, it is more efficient to force a higher underestimation error and manually reshape the final perimeter.
The six test areas assessed in this paper have shown BAMS to be a useful and efficient tool for mapping burned areas in diverse ecosystems and land covers in a short period of time and for obtaining reasonable results in a semi-automatic way, if compared with the manually recorded MTBS burned database.For example, the complex test scene 30-35, with more than 10 iterations required, was processed in only 30 min.The more spectrally clear 70-15 scene, with almost 250,000 ha of burned areas, was mapped in less than 15 min with very good agreement; several hours would have been necessary to digitize those perimeters manually.The automatic focus of this tool provides another important advantage related to the systematic search for burned patches around the entire scene, making it more difficult to miss burned patches than with manual searches of the areas with burned signal.Although in the tests shown in this paper the internal islands have been removed, perimeters derived through BAMS may also help with mapping unburned islands within the fire perimeter, something that may cause a significant overestimation rate of the affected areas.This is especially important for users who use this data for damage assessment.

Conclusions
BAMS, an operative tool to obtain burned area maps from Landsat data, is described in this paper.It can be downloaded for free from [47].Scientific groups like ESA's Fire_cci project [41] and the USGS Landsat Burned Area Essential Climate Variable Product [46] have already used the tool to obtain burned patch perimeters for validating global burned area products.More than 25 fire researchers around the world have already used the tool successfully, thus inaugurating a new era of producing accurate burned area maps from free satellite data.Appendix: Description of the BAMS Software Tool BAMS was initially developed using Visual Basic 6 programming language and ArcObjects 9.x.After the expiration of the support period for this programming language in 2008, the code of ABAMS was rewritten in Visual Basic .NET and adapted to ArcObjects 10.x.BAMS requires ArcGIS 10.x Advance license level (formerly ArcInfo), the Spatial Analyst extension, and the NET Framework version 4 to work properly, with minimum hardware requirements of a Pentium 4 2.2 Ghz processor and 2 GB RAM.The installation is done by executing the setup program available from [47].Once installed, the BAMS Main menu allows one to choose one of two analysis options: scene processing or supervised burned area mapping (Figure A1).
The "Process scenes" button allows for batch processing of multiple scenes.The user should set the parent folder to where the unzipped scenes are kept (maintaining the original folder names).All the processable scenes are listed as in Figure A2, grouped by the path-row code.This folder may contain several path-rows, sensors, and formats among those supported (see Section 2.1 for more information).
Only the selected (ticked) scenes will be processed.If the "Create temporal composites" option is selected, the temporal composites for each selected path-row will be created (see Section 2.1.2).
In the "Supervised-BAMS" option, the user should select first the Landsat scenes or temporal composites to be used as the pre-and post-fire images.Once the scenes are selected, the program shows a post-fire temporal composite, computes the multitemporal variables required (temporal difference between the pre-and post-fire images) and displays the main interface (Figure A3).In this interface, the user can draw the polygons of burned samples that will be used to define the seeds and to extract the second-phase thresholds (see Section 2.2 for more information about the algorithm used).The interface is divided into five sections: with an Editor button that allows for editing of the training polygons (both for the seeds or the second phase), and a BAMS custom tool to change the default color composite and apply some basic enhancing options to the scenes.3. Map area (center section): The area where the scenes, the burned samples, and the iteration results will be loaded.4. Variables area (right section): In this section, the user can select the variables to define the first and second phases of the algorithm.By default, BAIM, NBR, MIRBI, NDVI, and GEMI post-and pre-variables are checked, but they can be unchecked if they don't show enough differences between the burned and unburned classes.These variables can also be loaded into the map area by selecting them from the list (left click) and clicking in the "Add to Map" contextual dialogue.5. Output perimeter options (bottom-right section): In this section, the aggregation polygons with a fixed distance to 100 meters and the option to remove the internal polygons to obtain only the boundary of the fires are available.
After the burned samples have been edited, the BAMS software automatically extracts the thresholds for both phases and applies the burned area mapping algorithm, adding the result into the Table of Content (TOC).At this moment, the user may assess the result and redefine samples in case commission errors are detected (thresholds are too relaxed) or add new ones if omission errors are observed (thresholds are too strict).
Finally, the "Batch BAMS" option allows the batch execution of the burned area mapping algorithm criteria obtained with the "Supervised BAMS" option to other sets of scenes (Figure A4).

Figure 2 .
Figure 2. Flow chart for the burned area mapping section (a) seeds, (b) second stage result, (c) seeds and second stage result superimposed, (d) result.Note that some of the burned areas fulfilling only the second stage criteria are not kept as they do not contain a burned seed.

Figure 4 .
Figure 4.The 18-34 Test site's detail, where the commission errors (in the left area due to very burned areas not being included in the validation dataset) and the omission errors (on the right side) related with topographic shadows.

Figure 5 .
Figure 5.The 46-31 Test site example where the omission errors related to the weaker burned signal and the fragmentation of the fire are shown.

Figure 6 .
Figure 6.The 30-35 Test site's detail shows two high overestimation disagreement sources.On the left side, real overestimation disagreement can be observed (due to confusion of mainly not burned croplands), and on the right side, false overestimation disagreement areas are shown (related to real burned areas not included in the MTBS burned database).

Figure 7 .
Figure 7.The 30-35 Test site's detail, where the temporal composite minimizing the NBR value (right column) of the four individual scenes (left column) allows avoiding the no-data gaps of individual Landsat 7 data.

1 .
Table of Contents(left section), where four layers are automatically displayed  Burned_seeds_path_row_count.shp:The shapefile where the training polygons for the seed phase will be saved (yellow lines)  Burned_second_path_row_count.shp:The shapefile where the training polygons for the second phase will be saved (red lines)  Post-fire: Color composite for the post-fire scene (SWIRL-NIR-RED composite band as default)  Pre-fire: Color composite for the pre-fire scene (same color composite as the previous) 2. Main toolbar (top section), where the usual ArcGIS map navigation tools are available, along

Table 1 .
Landsat scenes (LT5 is Landsat 5 TM and LE7 Landsat 7 ETM+) and characteristics of the test areas (the state name, ecoregion, main burned land cover(s), and extent).

Table 2 .
Results obtained in the six test areas.Iterations required, and statistics related with the matrix error as compared with the MTBS dataset (underestimation and overestimation rates and Kappa).