30 m Resolution Global Annual Burned Area Mapping Based on Landsat Images and Google Earth Engine

Abstract

Heretofore, global Burned Area (BA) products have only been available at coarse spatial resolution, since most of the current global BA products are produced with the help of active fire detection or dense time-series change analysis, which requires very high temporal resolution. In this study, however, we focus on an automated global burned area mapping approach based on Landsat images. By utilizing the huge catalog of satellite imagery, as well as the high-performance computing capacity of Google Earth Engine, we propose an automated pipeline for generating 30-m resolution global-scale annual burned area maps from time-series of Landsat images, and a novel 30-m resolution Global annual Burned Area Map of 2015 (GABAM 2015) was released. All the available Landsat-8 images during 2014–2015 and various spectral indices were utilized to calculate the burned probability of each pixel using random decision forests, which were globally trained with stratified (considering both fire frequency and type of land cover) samples, and a seed-growing approach was conducted to shape the final burned areas after several carefully-designed logical filters (NDVI filter, Normalized Burned Ratio (NBR) filter, and temporal filter). GABAM 2015 consists of spatial extent of fires that occurred during 2015 and not of fires that occurred in previous years. Cross-comparison with the recent Fire_cci Version 5.0 BA product found a similar spatial distribution and a strong correlation ($R^2=0.74$) between the burned areas from the two products, although differences were found in specific land cover categories (particularly in agriculture land). Preliminary global validation showed the commission and omission errors of GABAM 2015 to be 13.17% and 30.13%, respectively.

1. Introduction

Accurate and complete data of fire locations and Burned areas (BA) are important for a variety of applications including quantifying trends and patterns of fire occurrence and assessing the impacts of fires on a range of natural and social systems, e.g., simulating carbon emissions from biomass burning [1]. Remotely-sensed satellite imagery has been widely used to generate burned area products. Burned area products at the global scale using satellite images have been mostly based on coarse spatial resolution data such as Advanced Very High Resolution Radiometer (AVHRR), Geostationary Operational Environmental Satellite (GOES), VEGETATION, or Moderate Resolution Imaging Spectroradiometer (MODIS) images. The main global burned area products include GBS (8 km) [2], Global Burned Area 2000 (GBA2000, 1 km) [3], GLOBSCAR (1 km) [4], GlobCarbon (1 km) [5], L3JRC (1 km) [6], MCD45 (500 m) [7], Global Fire Emissions Database (GFED) (0.5 degree) [8], MCD64 (500 m) [9], and Fire_cci (250 m) [10].

The recently released Fire_cci product was produced based on MODIS images and has the highest spatial resolution (250 m) of all the existing global BA products [1,11], and it can greatly support global climate modeling at moderate resolution thanks to its promising results. However, owing to the growing availability of satellite images of higher resolution, the requirement of BA products derived from those images is also increasing as these products can provide more detailed perimeters of BA, as well as small burned patches. For instance, imagery collected by the family of Landsat sensors is useful and appropriate for monitoring the extent of BA and provides spatial and temporal resolutions ideal for science and management applications. Landsat sensors can provide a longer temporal record (from the 1970s until now) of burned area relative to existing global BA products and potentially with increased accuracy and spatial detail in most areas on the Earth [12]. Great importance has been attached to developing BA products based on Landsat data in the past 10 years [12,13,14]. Up to now, there is no Landsat-based global BA product; however, some regional Landsat BA products have been publicly released in recent years. Australia released its Fire Scars (AFS) products derived from all available Landsat 5, 7, and 8 images using the time-series change detection technique [15]. Fire scars are automatically detected and mapped using dense time-series of Landsat imagery acquired over the period 1987–2015, and the AFS product only covers the state of Queensland, Australia. The Monitoring Trends in Burn Severity (MTBS) project, sponsored by the Wildland Fire Leadership Council (WFLC), provides consistent, 30-m resolution burn severity data and fire perimeters across all lands of the United States from 1984–2015 (only fires larger than 200 ha in the eastern U.S. and 400 ha in the western U.S. are mapped) [16]. MTBS products are generated based on the difference of the Normalized Burned Ratio (NBR) calculated from pre-fire and post-fire images, in which the BA boundary is delineated by on-screen interpretation, and the process of developing a categorical burn severity product is subjective and dependent on analyst interpretation. The Burned Area Essential Climate Variable (BAECV), developed by the U.S. Geological Survey (USGS), produces Landsat-derived BA products across the conterminous United States (CONUS) from 1984–2015, and its products were released in April 2017 [14]. The main differences between MTBS and BAECV is that the BAECV products are automatically generated based on all available Landsat images.

In summary, global BA products are only available at coarse spatial resolution, while 30-m resolution burned area products are limited to specific regions. The majority of coarse spatial resolution algorithms developed to produce global BA products use a multi-temporal change detection technique, because such satellite data have very high temporal resolution and are capable of monitoring fire-affected land cover changes. For example, the algorithm of the MODIS BA product (MCD45) is developed by the bi-directional reflectance model-based expectation change detection approach [7]. One of the difficulties in producing Landsat-based BA products is that the traditional approaches successfully applied to extract global BA from MODIS, VEGETATION, etc., do not work well due to the limited temporal resolution of the Landsat sensors. Moreover, the analysis of post-fire reflectance may be easily contaminated by clouds or weakened by quick vegetation recovery, particularly in Tropical regions [17]. Another difficulty is that global 30-m resolution annual BA mapping needs to utilize dense time-series Landsat images, and the required datasets can be hundreds of thousands of Landsat scenes, resulting in impractical processing time. Although some research has addressed detecting BA regionally from Landsat time-series [15,18,19], the results at the global scale have not been reported. However, thanks to Google Earth Engine (GEE), a new generation of cloud computing platforms with access to a huge catalog of satellite imagery and global-scale analysis capabilities [20], it is now possible to perform global-scale geospatial analysis efficiently as the pre-processing of satellite images becomes more user-friendly.

In this study, we focused on an automated approach to generate global-scale high resolution BA maps with dense time-series of Landsat images on GEE, in which all the available Landsat-8 images and various spectral indices were utilized to calculate the burned probability of each pixel using a machine learning model (random forest), and a seed-growing approach was conducted to shape the final burned areas after several carefully-designed logical filters. With learning-based global models, the classification becomes automatic without the need for regional adaptations after the training phase [21]. Due to the great variability of environments and burned conditions, several well-studied spectral indices for Landsat images were considered, including those specifically developed for burn detection, as they are sensitive to charcoal and ash deposition, such as NBR [22], NBR2 [23], the Burned Area Index (BAI) [24], the Mid-Infrared Burn Index (MIRBI) [25], and others that are not burn-specific, but useful to map burned areas when cooperating with burn-specific indices. For instance, although the Normalized Difference Vegetation Index (NDVI) is not the best index for burned area mapping, it is sensitive to vegetation greenness and therefore to the absence of vegetation in the case of burned areas [26]. The Global Environmental Monitoring Index (GEMI) [27] is an improved vegetation index, specifically designed to minimize problems of contamination of the vegetation signal by extraneous factors, which are considered very important for the remote sensing of dark surfaces such as recently-burned areas [28]. The Soil-Adjusted Vegetation Index (SAVI) [29], which was originally designed for sparse vegetation and outperforms NDVI in environments with a single vegetation type [30], is helpful to improve the separability of burns from soil and water [12]. The Normalized Difference Moisture Index (NDMI) [31], which is sensitive to the moisture levels in vegetation, is relative to fuel levels in fire-prone areas. By applying the proposed approach, a novel 30-m resolution Global Annual Burned Area Map of 2015 (GABAM 2015) was released, and the accuracy of this product was validated by using reference data derived from a stratified random sampling method and multiple data sources.

This work is related to the method of BAECV, but with many differences for global adaption, e.g., training data preparation, feature choosing, and burned seed generating. Particularly, some logical filters, i.e., NDVI filter, NBR filter, and temporal filter, are proposed to exclude the unreasonable confused surfaces, with the help of the MODIS Vegetation Continuous Fields (VCF) product. These filters are less related to the regional characteristics than the original BAECV algorithm, and thus more suitable for global BA mapping. Additionally, this was the first trial to produce a global BA map from a huge catalog of Landsat images in GEE, and it shows the potential to generate long time-series 30-m resolution global BA products automatically and efficiently (less than five days for an annual map), with the help of carefully-prepared training samples.

2. Methodology

2.1. Datasets

As many data were used in this work, we firstly give a brief description of all the involved data (as shown in

Table 1. Description of data sources. GFED, Global Fire Emissions Database; BA, Burned Area; MTBS, Monitoring Trends in Burn Severity.

Data	Usage	Source
MCD12C1 [32]	Stratified sampling for type of land cover	https://e4ftl01.cr.usgs.gov/MOTA/MCD12C1.006/
GFED4 [33]	Stratified sampling for fire frequency	https://www.globalfiredata.org/data.html
Landsat-8	BA mapping and validation	https://code.earthengine.google.com/dataset/LANDSAT/LC08/C01/T2_SR https://code.earthengine.google.com/dataset/LANDSAT/LC08/C01/T1_SR
MOD44B [34]	Adjustment constraint conditions for BA mapping	https://code.earthengine.google.com/dataset/MODIS/051/MOD44B
Fire_cci v5 [11]	Comparison	https://geogra.uah.es/fire_cci
CBERS-4 MUX	Validation	http://www.dgi.inpe.br/catalogo/
Gaofen-1 WFV	Validation	http://218.247.138.119:7777/DSSPlatform/productSearch.html
MTBS [16]	Validation	https://www.mtbs.gov/direct-download

), for the sake of clarity.

The following are more specific notes about the data selection for BA mapping and validation.

BA mapping: In this study, all the available Landsat images during 2014–2015 were the main data source used for global BA mapping. At a pixel, the occurrence of a single Landsat satellite could be about 23 or more times (considering the overlap between adjacent paths, particularly at higher latitudes) within a year, and it would double when contemporary satellites (e.g., Landsat-7 and Landsat-8) are utilized. However, considering the failure of the Scan Line Corrector (SLC) in the ETM+ instrument of the Landsat-7 satellite, we only utilized USGS Landsat-8 Surface Reflectance collections. The Quality Assessment (QA) band of Landsat image, which was generated by the FMaskalgorithm [35], was used to perform QA masking. Pixels flagged as being clouds, cloud shadows, water, snow, ice, or filled/dropped pixels were excluded from Landsat scenes, and only clear land pixels remained after QA masking.

Validation: Commonly, when satellite data are used as reference data, they should have higher spatial resolution than the data used to generate the BA product [36]. For the Landsat BA product, however, access to global higher resolution time-series satellite data is difficult, and [37] suggested a thorough validation scheme, in which high-resolution data were used to complement the independent Landsat-derived reference data. Consequently, in this study, some publicly-available satellite images of higher-resolution were utilized, while Landsat comprised the majority of the validation data source. Specifically, Landsat-8 (LC8) images were employed to generate reference data independently for most of the validation sites except those located in the United States (U.S.), South America, and China. In the U.S., the MTBS perimeters of 2015 were used as the supplemental reference data of LC8 images, and in South America and China, CBERS-4 MUX(CB4) and Gaofen-1 WFV(GF1) satellite images were used to create the perimeters of the burned area, respectively. The characteristics of CB4 and GF1 are illustrated in

Table 2. Characteristics of CBERS-4 MUX and Gaofen-1 WFV.

Sensors	Spatial Resolution at Nadir (m)	Swath Width at Nadir (km)	Spectral Bands ($\mathsf{\mu }$m)
			Blue	Green	Red	NIR
CBERS-4 MUX	20	120	0.45–0.52	0.52–0.59	0.63–0.69	0.77–0.89
Gaofen-1 WFV	16	192

.

2.2. Sampling Design

The spectral characteristics of burned areas vary in complex ways for different ecosystems, fire regimes, and climatic conditions. In terms of guaranteeing the accuracy of the global burned area map and also the completeness of quality assessment, a stratified random sampling method [38,39,40] was used to generate two sets of sites for classifier training and the validation of GABAM 2015, respectively. The training and validation sites were chosen randomly based on stratifications of both fire frequency and the type of land cover.

Firstly, the Earth’s land surface was partitioned based on the 14 land cover classes according to the MCD12C1 product [32] of 2012 using the University of Maryland (UMD) scheme. These types were then merged into 8 categories based on their similarities [41], i.e., broadleaved evergreen, broadleaved deciduous, coniferous, mixed forest, shrub, rangeland, agriculture, and others.

Table 3. Mapping between the original University of Maryland (UMD) land cover types and the new classifications for the geographic stratification.

New Classification	Original UMD Type
Broadleaved Evergreen	Evergreen Broadleaf Forest
Broadleaved Deciduous	Deciduous Broadleaf Forest
Coniferous	Evergreen Needleleaf Forest
	Deciduous Needleleaf Forest
Mixed Forest	Mixed Forest
Shrub	Closed Shrublands
	Open Shrublands
Rangeland	Woody Savannas
	Savannas
	Grasslands
Agriculture	Croplands
Others	Water
	Urban and Built-up
	Barren or Sparsely Vegetated

shows the reclassification rule from UMD land cover types to new classifications. As the “others” category consists of the biomes less prone to fire, only the other 7 land cover categories were considered to create the geographic stratifications in this work.

Secondly, the globe was divided into 5 partitions based on the BA density in 2015 provided by the Global Fire Emissions Database (GFED) Version 4.0 [33], the most widely-used inventory in global biogeochemical and atmospheric modeling studies [9]. Specifically, GFED4 monthly products of 2015 were utilized to produce an annual composition (GFED4 2015), consisting of 720 rows and 1440 columns, which correspond to the global ${0.25}^{\circ }\times {0.25}^{\circ }$ GFED grid, and each pixel summed the total areas of BA (BA density, km${}^2$) occurring in the grid cell during the whole year. The BA density of GFED4 2015 was then divided into 5 equal-frequency intervals [41] with quantile classification.

By spatially intersecting the 7 land cover categories and 5 BA density levels, we obtained the final 35 strata with different fire frequencies and biomes. The samples were equally allocated to 5 BA density levels, but for different land cover categories, we also took into account the BA extent within each stratum: larger sample sizes were allocated to strata with higher BA extent [42]. According to the strategy of stratified sampling, 120 samples (24 for each BA density level) were randomly selected to generate the training dataset, and the spatial dimension of sampling units was based on Landsat World Reference System II (WRS-II). Similarly, 80 validation sites (16 for each BA density level) were also created by stratified random sampling, but trying to keep a distance (at least 200 km) from the training samples so as not to fall into the extent of training Landsat scenes. Figure 1 illustrates the distribution of 120 random Landsat image scenes and 80 validation sites over a map of BA density extracted from GFED4 2015, and

Table 4. Distribution of training and validation samples over the different land cover types.

Land Cover Type	Training Sample Count	Validation Sample Count
Broadleaved Evergreen	16	11
Broadleaved Deciduous	12	9
Coniferous	13	9
Mixed Forest	12	8
Shrub	18	12
Rangeland	25	15
Agriculture	24	16

shows the distribution of training and validation samples over the different land cover types.

2.3. Training Dataset

In terms of analyzing the characteristics of burned areas in Landsat images, 120 Landsat-8 image scenes were chosen according to the WRS-II frames generated by stratified random sampling in Section 2.2. All the Landsat-8 images used in this study were acquired from datasets of USGS Landsat-8 Surface Reflectance Tier 1 and Tier 2 in the Google Earth Engine platform. These data have been atmospherically corrected using LaSRC [43] and include a cloud, shadow, water, and snow mask produced using FMask [35], as well as a per-pixel saturation mask. For the purpose of burned area mapping, 6 bands of the Landsat-8 image were used, i.e., three visible bands (blue, 0.452–0.512 $\mathsf{\mu }$$\mathrm{m}$; green, 0.533–0.590 $\mathsf{\mu }$$\mathrm{m}$; red, 0.636–0.673 $\mathsf{\mu }$$\mathrm{m}$), Near Infrared band (NIR, 0.851–0.879 $\mathsf{\mu }$$\mathrm{m}$), and two Short Wave Infrared bands (SWIR1, 1.566–1.651 $\mathsf{\mu }$$\mathrm{m}$; SWIR2, 2.107–2.294 $\mathsf{\mu }$$\mathrm{m}$).

In this study, the burned area mapping algorithm was implemented on the GEE platform, and the maximum quantity of input samples was limited by GEE’s classifiers; thus, an average 90–100 sample points were collected by experienced experts from each Landsat-8 image, making the total quantity of sample points 12,881 (6735 burned samples and 6,146 unburned samples). Specifically, Shortwave Infrared (SWIR2), Near Infrared (NIR) and green bands were composited into a Red, Green, Blue (RGB) combination in order to visualize burned areas better, and burned samples, including fire scars of different burn severity and of various biomass types, were extracted from the pixels showing magenta color [44]. The unburned pixels were extracted randomly over the non-fire-affected areas covering vegetation, built-up land, bare land, topographic shadows, borders of lakes, etc. For those confused pixels for which it was difficult to identify whether they were burned scars, a further check was performed by examining the Landsat images on the nearest date of the previous year or higher resolution images on the nearest date in the Google Earth software. To ensure only clearly-burned pixels were selected, the burned samples were collected carefully to avoid pixels near the boundaries of burned scar [13]; and burned pixels located in burning flame or covered by smoke were also excluded to prevent potential contamination of burned samples. Note that each sample point covered a few (5 on average) Landsat pixels; these pixels were clearly burned or not affected by fire since they were not located near the boundaries of burned scar. Land surface reflectance of the collected samples in the blue, green, red, NIR, SWIR1, and SWIR2 bands were extracted for further analysis.

2.4. Sensitive Features for Burned Surfaces

Figure 2 shows the statistical mean reflectance (with standard deviations) of burned and unburned samples in Landsat 8 bands.

Burned areas are characterized by deposits of charcoal, ash, and fuel, and the reflectance of the burned pixels generally increases along with the wavelength, while the burned pixels have similar reflectance in the SWIR1 and SWIR2 bands, which is greater than that in other bands. However, the spectral character of post-fire pixels varies greatly (standard deviations in Figure 2) according to the type and condition of the vegetation prior to burning and the degree of combustion [45], and none of existing spectral indices can be considered the best choice for identifying burned surfaces without misclassification with other targets in all environments or fire regimes [46]. Consequently, in this study, we made use of the most common spectral indices for Landsat images previously suggested in BA studies, and their formulas are summarized as

Table 5. The formulas of spectral indices that are sensitive to burned areas.

Name	Abbreviation	Reference	Formula
Normalized Burned Ratio	NBR	Key and Benson [22]	$NBR=\frac{{\rho }_{NIR}-{\rho }_{SWIR2}}{{\rho }_{NIR}+{\rho }_{SWIR2}}$
Normalized Burned Ratio 2	NBR2	Lutes et al. [23]	$NBR2=\frac{{\rho }_{SWIR1}-{\rho }_{SWIR2}}{{\rho }_{SWIR1}+{\rho }_{SWIR2}}$
Burned Area Index	BAI	Martín [24]	$BAI=\frac{1}{{({\rho }_{NIR}-0.06)}^2+{({\rho }_{RED}-0.1)}^2}$
Mid-Infrared Burn Index	MIRBI	Trigg and Flasse [25]	$MIRBI=10{\rho }_{SWIR2}-0.98{\rho }_{SWIR1}+2$
Normalized Difference Vegetation Index	NDVI	Stroppiana et al. [26]	$NDVI=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$
Global Environmental Monitoring Index	GEMI	Pinty and Verstraete [27]	$GEMI=\frac{\eta (1-0.25\eta )-({\rho }_{RED}-0.125)}{1-{\rho }_{RED}}$,
			$\eta =\frac{2({\rho }_{NIR}^2-{\rho }_{RED}^2)+1.5{\rho }_{NIR}+0.5{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}+0.5}$
Soil-Adjusted Vegetation Index	SAVI	Huete [29]	$SAVI=\frac{(1+L)({\rho }_{NIR}-{\rho }_{RED})}{{\rho }_{NIR}+{\rho }_{RED}+L}$, $L=0.5$
Normalized Difference Moisture Index	NDMI	Wilson and Sader [31]	$NDMI=\frac{{\rho }_{NIR}-{\rho }_{SWIR1}}{{\rho }_{NIR}+{\rho }_{SWIR1}}$

${\rho }_{RED}$ is the surface reflectance in the red; ${\rho }_{NIR}$ is the surface reflectance in NIR; ${\rho }_{SWIR1}$ is the surface reflectance in the SWIR1 band; and ${\rho }_{SWIR2}$ is the surface reflectance in the SWIR2 band.

. Together, 14 Landsat features (8 spectral indices in Table 5 and the surface reflectance in 6 bands of the Landsat-8 image) were selected as sensitive features to perform global burned area mapping.

2.5. Burned Area Mapping via GEE

In this study, the annual burned area map was defined as the spatial extent of fires that occurred within a whole year and not of fires that occurred in previous years. Therefore, global 30-m resolution annual burned areas’ mapping needed to utilize dense time-series Landsat images, and the pipeline of annual burned area mapping via GEE is described as Figure 3.

As shown in Figure 3, the pipeline mainly consisted of three steps, model training, per-pixel processing, and burned area shaping, and the following provides more details of each step.

2.5.1. Model Training

The Random Forest (RF) algorithm provided by GEE was applied to train a decision forest classifier, and the global training data consisted of 6735 burned and 6146 unburned samples, which were manually collected from 120 Landsat scenes generated by stratified random sampling (in Section 2.2 and Section 2.3). The random forest classifier with a higher number of decision trees usually provides better results, but also causes higher cost in computation time. Since the input features of the algorithm include the Surface Reflectance (SR) in 6 bands of the Landsat-8 image, as well as 8 spectral indices that have high sensitivity to the burned surface, we limited the number of decision trees in the forest to 100 for a trade-off between accuracy and efficiency. Additionally, we chose “probability” mode for GEE’s RF algorithm, in which the output is the probability that the classification is correct, and the probability would be further utilized to perform region growing in the step of burned area shaping.

2.5.2. Per-Pixel Processing

In this step, Landsat surface reflectance collections from GEE, which consist of all the available Landsat scenes, were employed for dense time-series processing. At each pixel, the geometrically-aligned dense time-series Landsat image scenes provided a reflectance stack of 6 bands, which was then split into two stacks by date filters, i.e., a stack of the current year and that of the previous year.

For the reflectance stack of the current year, 8 spectral indices were computed at each time period, and then, the trained decision forest classifier in Section 2.5.1 produced a stack of burned probability using the 8 spectral indices and the reflectance of 6 bands. The maximum value of a probability stack indicates the probability that the pixel had ever appeared like a burned scar during the whole year. Four quantities were noted for each pixel, i.e., the date on which the maximum probability was observed ($t_1$), as well as the burned probability ($p_{max}$), NDVI value ($NDV{I}_1$), and NBR value ($NB{R}_1$) on that date. However, a single image is usually not adequate to spectrally separate the burned areas from confused surfaces caused by shadows, flooding, snow melt, agricultural harvesting, etc. [47]; the burned scars that occurred in previous years, but had not yet recovered (particularly in boreal forests) should also be excluded from the annual BA map of the current year. In this sense, we also were concerned with the summary statistics of current year and previous year: $NDV{I}_2$, the maximum NDVI value within the couple of years (current year and previous year); $t_2$, the date of $NDV{I}_2$; and $NB{R}_2$, the minimum NBR value within the previous year. Then, most of the unreasonable tree-covered burned-like pixels would be excluded unless they met all the following constraints.

• $NDV{I}_2>T_{NDVI}$, the maximum NDVI value within the couple of years should be greater than a threshold $T_{NDVI}$. We choose NDVI as it has been found to be a good identifier of vigorous vegetation, and this constraint is used to exclude areas that appeared as burned, but in fact were just lacking vegetation.
• $NDV{I}_2-NDV{I}_1>T_{dNDVI}$, the difference between the maximum NDVI and the NDVI when the pixel was most like a burned scar should be greater than a threshold $T_{dNDVI}$. This constraint ensures evidence of vegetation decrease when the burn happened.
• $NB{R}_2-NB{R}_1>T_{dNBR}$, the NBR value of a burned pixel should be less than the minimum NBR of the previous year, and the threshold $T_{dNBR}$ is the minimum acceptable decline of NBR. This constraint is useful to exclude false detections with periodic variation of NBR and NDVI, such as mountain shadows, burned-like soil in deciduous season, snow melting, and flooding.
• $t_1>t_2$ or $t_2-t_1>T_{DAY}$, the date when the vegetation becomes greenest should be earlier than the burning date or the lagged days should be greater than a threshold $T_{DAY}$. For a tree-covered surface, it usually takes a long time for the vegetation to recover more flourishing than the previous year, thus the burn-like pixels with $t_1<=t_2$ are likely attributed to a false alarm. However, as the recovering of burned trees can be fast in tropic regions, high post-fire regrowth within a reasonable amount of days is also acceptable.

We named the first two constraints as “NDVI filter” and the third and fourth ones as “NBR filter” and “temporal filter”, respectively. In this work, the thresholds in the above constraints were chosen empirically, $T_{NDVI}=0.2$, $T_{dNDVI}=0.2$, $T_{DAY}=100$ (days), and $T_{dNBR}=0.1$. Determining a globally optimal NDVI threshold is not easy or even impossible for various types and conditions of the vegetation, and we chose a low threshold $T_{NDVI}=0.2$ [48], not expecting to exclude directly all confused surfaces never covered by vegetation. Actually, the second constraint would also help to exclude non-vegetation with high NDVI, because the decline of NDVI, in the absence of vegetation variation, commonly would not meet the constraint. The change of NBR in pre-fire and post-fire images, defined as delta NBR or dNBR, has proven to be a good indicator of burn severity and vegetation regrowth (the higher the severity, the greater the dNBR) [49,50]. It was suggested that a dNBR greater than 0.1 commonly indicates a burn of low severity [23]; thus, we chose $T_{dNBR}=0.1$. Lastly, in the temporal filter, a fixed time for all kinds of trees to recover spectrally is also not available, and we just approximately chose an average time, 100 days.

However, for herbaceous vegetation, we should use only the first two constraints, as grassland usually recovers very quickly and can be burned year after year. Accordingly, the annual MODIS Vegetation Continuous Fields (VCF) 250-m Collection 5.1 (MOD44B) product [34] of the current and previous year, which contains the tree-cover percent layer and non-tree vegetation layer, were utilized to determine whether the pixel was dominated by trees or by herbaceous vegetation. Passing the filters of NDVI, NBR, and temporal context, those pixels with an annual burned probability greater than or equal to 0.95 (“probability filter”) were selected as seeds for region growing.

In addition, the global training samples in Section 2.3 were used to test the sensitivity of thresholds in the four constraints. By applying the per-pixel processing in Section 2.5.2, we got the values of $NDV{I}_2$ and $NDV{I}_2-NDV{I}_1$ in all the burned and unburned samples and the values of $NB{R}_2-NB{R}_1$ and $t_2-t_1$ in the burned and unburned samples located in tree-dominated regions, and their histograms are illustrated in Figure 4a–d. Then, parameter tuning was performed to check the percentages of missing burned samples and confused unburned samples when applying various values of $T_{NDVI}$ (0–0.4), $T_{dNDVI}$ (0–0.4), $T_{dNBR}$ ($-0.1$–0.3), and $T_{DAY}$ (0–150), as shown in Figure 4e–h. Generally, more burned samples were missed and less unburned were included when $T_{NDVI}$, $T_{dNDVI}$, $T_{dNBR}$, and $T_{DAY}$ increased, and the selected thresholds should balance the omission and commission errors. According to Figure 4, less than 3% of burned samples were missed, but around 10%, 30%, 70%, and 10% of unburned samples can be excluded, respectively, by choosing $T_{NDVI}=0.2$, $T_{dNDVI}=0.2$, $T_{dNBR}=0.1$, and $T_{DAY}=100$.

2.5.3. Burned Area Shaping

In this step, a region-growing process was employed to shape the burned areas. Region growing has been proven to be necessary for BA mapping in many studies [12,13,18], because spectral-based methods sometimes give ambiguous evidence (i.e., spectral overlapping between burned areas and unrelated phenomena with similar spectral characteristics, such as cloud shadows, ephemeral water, or dark soils [12]), and accepting all positive evidence can lead to confusion errors. Although candidate seeds were chosen with high confidence, false seed pixels were still frequently included in confused surfaces, e.g., shadows and borders of lakes. Different from the candidate seeds in the actual burned scars, those falsely-introduced seed pixels always distributed sparsely. Consequently, we aggregated the seed pixels into connected components using a kernel of 8-connected neighbors. By ignoring small fires with area less than 1 ha [51], those fragmentary components (smaller than 11 pixels), which included most false seed pixels, were removed. Finally, an iterative procedure of region growing was performed around each seed pixel. For each iteration, the 8-connected neighbors of the seed pixels were aggregated as burned pixels (new seeds) if their burned probabilities were greater than or equal to 0.5, and the iteration stopped when no more pixels could be aggregated as burned pixels. Figure 5 shows an example of region growing. One can see that only some pixels showing strong magenta color in the burned scars were chosen as seeds, while those showing light magenta color were labeled as candidates for region growing, including some actual burned pixels, as well as some false detections (right-middle in Figure 5b). However, after the processes of small seeds’ removal and region growing, the false detections were excluded, while those candidates near the seeds were aggregated to the final BA map.

2.6. Comparison with the Fire_cci Product

As 30-m resolution global burned area products are currently not available, we made a comparison between GABAM 2015 and the Fire_cci Version 5.0 products (spatial resolution was approximately 250 m) [11], which are based on MODIS satellite imagery. The monthly Fire_cci pixel BA products of 2015 were composited as an annual pixel BA product by labeling the pixels as burned ones once their values in the Julian day (the date of the first detection) layer were valid (from 1–366) in any of the 12 monthly products. Additionally, in order to perform regression analysis between two products of different spatial resolution, we also produced an annual grid composition of BA within 2015 from the composited annual pixel BA product by computing the proportion of burned pixels in each ${0.25}^{\circ }\times {0.25}^{\circ }$ grid. Note that the monthly grid BA products of Fire_cci were not used to composite the annual grid product, because summing up the areas of BA for each grid in all monthly products might result in repetitive counting at those pixels burned more than once within the year.

2.7. Validation

2.7.1. Data Sources

Accuracy assessment was carried out according to the 80 validation sites that were created in Section 2.2, and the reference data were selected in these sites from multiple data sources, including fire perimeter datasets and satellite images. As described in Section 2.1, LC8 images were employed to generate reference data independently for most of the validation sites except those located in the United States (U.S.), South America, and China. In the U.S., the MTBS perimeters of 2015 were used as the supplemental reference data of LC8 images, and in South America and China, CB4 and GF1 satellite images were used to create perimeters of burned area, respectively. Note that the size of the validation site varied by the type of data source, i.e., a WRS-II frame (about 185 km × 185 km) for Landsat images, a scene for CB4 images (about 120 km × 120 km), and a box of 100 km × 100 km for GF1 images. Using Landsat frames or image scenes as a unit of the validation site is convenient for data downloading and processing; we chose a smaller box for GF1 to improve the data availability considering the extent of GF1 frames or scenes not being fixed due to the long orbital return period.

2.7.2. Reference Data Generation

In each validation site, all the available image scenes (LC8, CB4, or GF1) acquired in 2015 were used. LC8 images were orthorectified surface reflectance products; CB4 images were ortho products; and GF1 images were not geometrically rectified. The procedure of generating reference BA can be summarized as the following steps.

• Preprocessing
  All the images utilized to generate BA reference data were spatially aligned with a mean squared error of less than 1 pixel. The ortho-rectified LC8 and CB4 images met the requirement of geometric accuracy, yet the GF1 images did not. Accordingly, an automated method [52] was applied to orthorectify the time-series GF1 images, taking the LC8 panchromatic images (spatial resolution was 15 m) as geo-references.
• BA detection
  BA perimeters were generated from the time-series images via a semi-automatic approach. Firstly, image pairs (pre- and post-fire) were manually selected from the time-series image by checking whether any new burned scars appeared in the newer images. For LC8 images, SWIR2, NIR, and green bands were composited in a Red, Green, Blue (RGB) combination; for CB4 and GF1 images, red, NIR, and green bands were composited in an RGB combination. The identification of BA might be difficult for CB4 and GF1 images due to the lack of shortwave infrared bands; thus, the Fire_cci BA product was used to verify the BA identification. Secondly, burned and unburned samples were manually collected from each selected image pair. The burned samples included only the newly-burned scars, which appeared burned in the newer image, but unburned in the older image; the unburned samples consisted of unburned pixels, partially recovered BA pixels, and also pixels covered by cloud or cloud shadows in either images. Afterwards, the Support Vector Machines (SVM) classifier in ENVI^TM (provided by Harris Geospatial in Broomfield, CO, United States) software was used to classify each image pair into burned and unburned pixels, and the detected burned pixels in all the image pairs were integrated to create a composited annual BA map. Note that the sensitive features in Section 2.4 were utilized in SVM for each LC8 image pair; but for CB4 and GF1 images, the features used for classification consisted of the Digital Number (DN) values in four bands of an image pair (in total, 8 DN values), as most of the burned-sensitive spectral indices cannot be derived from the RGB-NIR bands. Finally, the BA perimeters of 2015 were generated from the annual BA composition using the vectorization tool in ArcGIS^TM (provided by Environmental Systems Research Institute in Redlands, CA, United States) software.
• Reviewing and manually revision
  The result of the supervised classifier (SVM) and automated vectorization algorithm might not be perfect; thus, BA perimeters were further edited visually by experienced experts, via overlapping the vector layer of BA perimeters with the satellite image layers.

Additionally, in the U.S., the MTBS perimeters of 2015 were directly used as the main reference data, supplemented by the interpreted results of LC8 time-series images, which could help to avoid missing small fires.

2.7.3. Assessment

In this work, a cross-tabulation (

Table 6. Cross-tabulation between Global annual Burned Area Map (GABAM) 2015 and the reference data.

		Reference Data (pixel)
		Burned	Unburned	Total
GABAM 2015 (pixel)	Burned	$X_{11}$	$X_{12}$	$X_{11}+X_{12}$
	Unburned	$X_{21}$	$X_{22}$	$X_{21}+X_{22}$
	Total	$X_{11}+X_{21}$	$X_{12}+X_{22}$	$X_{11}+X_{12}+X_{21}+X_{22}$

) between the pixels assigned by in our BA product and in the reference data was computed to produce the confusion matrix [53], which was then applied to derive three statistics, i.e., commission error, omission error, and overall accuracy, to assess the accuracy of our BA product.

• Commission error ($E_c$): $X_{12}/(X_{11}+X_{12})$, the ratio between the false BA positives (detected burned areas that were not in fact burned) and the total area classified as burned by GABAM 2015.
• Omission error ($E_o$): $X_{21}/(X_{11}+X_{21})$, the ratio between the false BA negatives (actual burned areas not detected) and the total area classified as burned by the reference data.
• Overall accuracy ($A_o$): $(X_{11}+X_{22})/(X_{11}+X_{12}+X_{21}+X_{22})$, the ratio between the area classified correctly and the total area to evaluate.

To assess the dispersion of accuracy among the validation sites, the standard errors of the above three statistics were used:

$$\begin{array}{cc}\hfill s_{E_c}& =\sqrt{\sum _{t=1}^n{({E_c}^t-\overline{E_c})}^2/(n-1)}\hfill \\ \hfill s_{E_o}& =\sqrt{\sum _{t=1}^n{({E_o}^t-\overline{E_o})}^2/(n-1)}\hfill \\ \hfill s_{A_o}& =\sqrt{\sum _{t=1}^n{({A_o}^t-\overline{A_o})}^2/(n-1)}\hfill \end{array}$$

(1)

where n is the number of validation sites to be evaluated; ${E_c}^t$, ${E_o}^t$, and ${A_o}^t$ are the $E_c$, $E_o$, and $A_c$ at validation site t; $\overline{E_c}$, $\overline{E_o}$, and $\overline{A_o}$ are the average $E_c$, $E_o$, and $A_c$ at the n sites.

Additionally, taking into account the stratified sampling design, a stratified combined ratio estimator [38,54] was also used to assess the global accuracy:

$$\widehat{R}=\frac{{\displaystyle \sum _{h=1}^H{S}_h{\overline{y}}_h}}{{\displaystyle \sum _{h=1}^H{S}_h{\overline{x}}_h}}$$

(2)

where H is the number of strata; $S_h$ is the area of stratum h; ${\overline{x}}_h$ and ${\overline{y}}_h$ are the sample means of $x_t$ and $y_t$ at stratum h; $x_t=X_{11}^t+X_{12}^t+X_{21}^t+X_{22}^t$, $y_t=X_{11}^t+X_{22}^t$, and $X_{11}^t$, $X_{12}^t$, $X_{21}^t$, $X_{22}^t$ are from the cross-tabulation at validation site t.

Note that the area of a stratum was calculated by summing up all the areas of pixels within the stratum at a specific resolution in the geographic coordinate system. The area of an image pixel, whose longitude and latitude are expressed in radians, $\lambda$ and $\varphi$, can be approximately calculated:

$$\begin{array}{cc}\hfill A_{\lambda ,\varphi }& =R_{earth}^2\times r_{rad}\times \left[sin\left(\varphi +r_{rad}\right)-sin\varphi \right]\hfill \end{array}$$

(3)

where $R_{earth}=$ 6,371,007.2 m is the Earth’s authalic radius [55] and $r_{rad}$ is the spatial resolution of the image expressed in radians. In this study, we used $r_{rad}=4.3633\times {10}^{-6}$ (approximately 30 m).

3. Results and Analysis

3.1. Product Description

Employing the proposed approach, we produced GABAM 2015, which was projected in a geographic (Lat/Long) projection at ${0.00025}^{\circ }$ (approximately 30 m) resolution, with the WGS84 horizontal datum and the EGM96 vertical datum. The result consisted of $10\times 10$ degree tiles spanning the range 180 W–180 E and 80 N–60 S and can be freely downloaded from https://vapd.gitlab.io/post/gabam2015/. To make the visualization of GABAM better, burned area density was used instead of directly drawing the burned pixels on a global map, and it was defined as the proportion of burned pixels in a ${0.25}^{\circ }\times {0.25}^{\circ }$ grid. An overview of global distribution of burned area density, derived from the one-arc-second resolution GABAM 2015, is shown in Figure 8a, together with that of the Fire_cci product in Section 3.2.

Figure 6 illustrates an examples of GABAM 2015 in Canada, and the annually-composited Landsat reference images with minimum NBR values of 2015 and 2014 are also included. This region is located in high latitude zones, and the burned scars may not completely recover within a year. Consequently, when new burning occurs around the unrecovered burned scars, we must determine which burned scars come from this year. Owing to the temporal filters, GABAM succeeded in clearing up such confusion. From Figure 6b, one can see that the burned scars mainly consisted of two components, separated by the river. Figure 6a, however, shows that burned scars on the right side of the river can be observed in 2014; hence, the result of GABAM 2015 only remained the component on the left side.

3.2. Comparison with the Fire_cci Product

3.2.1. Visual Comparison

Figure 7 shows an example of the two annual pixel BA products, and it can be seen that both products correctly detected the BAs in the Landsat image (Figure 7b), yet the BAs in Figure 7c occupied more pixels than those in Figure 7d. Due to the limitation in spatial resolution of the input sensor of the Fire_cci BA product, some of the mixed pixels (consisting of burned and unburned pixels) may be classified as burned ones. On the other hand, the result of GABAM 2015 showed finer boundaries of BAs, compared with that of the Fire_cci product.

3.2.2. Global Grid Map

Figure 8 illustrates the GABAM and Fire_cci annual grid composition of BA, consisting of the percentage of burned pixels in each ${0.25}^{\circ }\times {0.25}^{\circ }$ grid. Figure 8a,b shows a similar global distributions of BA density.

3.2.3. Regression Analysis

Figure 9 shows the proportion of BA in ${0.25}^{\circ }\times {0.25}^{\circ }$ grids of different land cover categories in Table 3, for the Fire_cci product (x-axis) and GABAM 2015 (y-axis), and regression analysis was also performed between the two products, providing a regression line (expressed as the slope and the intercept coefficient estimates) and the coefficient of determination ($R^2$) for each land cover category (Figure 9a–h) and for the global scale (Figure 9i). Moreover, as many points overlapped in the scatter graphs, we also rendered the scatters with different colors according to the number of grid cells (1–10 or more) having the same proportion values.

According to Figure 9, the intercept values of the estimated regression lines were close to zero, while the slopes were lower than one, showing that GABAM burned area was less than the Fire_cci product [41]. Moreover, the distribution and color of scatters in Figure 9i also show that a large number of grids were considered to have a higher burned proportion by the Fire_cci product than by GABAM. The main reason for the inconsistency can be attributed to the difference in the spatial resolution of data sources, and less pixels were commonly classified as BA in Landsat images, e.g., Figure 7. Specifically, only a few ${0.25}^{\circ }\times {0.25}^{\circ }$ grids were occupied by more than 90% BA in GABAM, while grids with a high proportion of BA were more common in the Fire_cci product.

Considering the coefficients of determination of estimated regression lines, the two products showed the highest linear relationship strengths in coniferous forest ($R^2=0.82$), rangeland ($R^2=0.75$), and shrub ($R^2=0.62$) and the lowest strengths in agriculture land ($R^2=0.31$) and the “others” category ($R^2=0.07$). In the “others” category, which is considered to be not prone to fire, the two products only included a few grids containing BA (with low burned proportions); thus, they were not likely to be correlated; the low correlation in agriculture land is owed to the uncertainty of both products, which will be further discussed in Section 4.

The quantity and color of scatters in Figure 9 indicate that most of burned areas were located in rangeland, and the global relationship (Figure 9i) of the GABAM and Fire_cci product was mainly determined by that in rangeland (Figure 9e), i.e., woody savannas, savannas, and grasslands.

3.3. Validation

To assess the accuracy of GABAM 2015, confusion matrices were computed for each validation site according to Table 6, and the global confusion matrix (Formula (4)) was then generated by averaging all these confusion matrices.

$$\left[\begin{array}{cc}{X}_{11}& X_{12}\\ X_{21}& X_{22}\end{array}\right]=\left[\begin{array}{cc}5,473,720& 823,170\\ 2,360,096& 43,661,559\end{array}\right]$$

(4)

The averaging commission error, omission error, and overall accuracy were finally derived from the global confusion matrix as described in Section 2.7.3: $E_c$ and $E_o$ of GABAM 2015 were 13.17% and 30.13%, respectively, while $A_o$ was 93.92%; and the standard errors of the three statistics, $s_{E_c}$, $s_{E_o}$, and $s_{A_c}$, were 10.83%, 19.58%, and 7.77%, respectively. When taking the stratified sampling design into account, the global accuracy estimator $\widehat{R}$ was 92.03%.

Generally, GABAM 2015 was expected to have a lower $E_c$, but a higher $E_o$. High omission error might result from several reasons:

• In the validation sites located in tropical zones, clear burned evidence was frequently missed by the Landsat sensor due to the quick recovery of the vegetation surface. This point will be further discussed in Section 4.
• Some pixels located within a burned area, but not showing a strong burned appearance, might be excluded by GABAM 2015 (e.g., Figure 7d), while they were considered as a part of a complete burned scar in the reference data. Particularly, high $E_o$ was found at those validation sites using MTBS perimeters, e.g., the $E_c$ and $E_o$ of the validation site in Figure A5 were 1.45% and 67.97%. Furthermore, this high omission error might result from the high commission error associated with MTBS perimeters [56].

Table 7. Information of site validation examples.

Land Cover Type	${\mathit{E}}_{\mathit{c}}$ (%)	${\mathit{E}}_{\mathit{o}}$ (%)	${\mathit{A}}_{\mathit{o}}$ (%)	${\mathit{s}}_{{\mathit{E}}_{\mathit{c}}}$ (%)	${\mathit{s}}_{{\mathit{E}}_{\mathit{o}}}$ (%)	${\mathit{s}}_{{\mathit{A}}_{\mathit{o}}}$ (%)
Broadleaved Evergreen	8.64	10.95	90.99	9.14	16.52	6.78
Broadleaved Deciduous	23.59	34.85	99.03	12.33	19.16	7.22
Coniferous	7.41	18.27	99.77	11.47	16.15	6.00
Mixed Forest	8.73	34.33	98.36	9.30	24.19	8.41
Shrub	13.00	3.78	99.49	11.05	16.05	8.78
Rangeland	11.91	23.06	91.79	13.55	17.91	9.04
Agriculture	10.91	45.38	94.41	10.50	28.09	7.33

shows the average accuracy of GABAM 2015 in various land cover categories, and more details of the validation can be found in Appendix A, which includes five examples of validation sites from various regions, with different data sources as reference data.

4. Discussion

Different from the satellite images of coarse spatial resolution, the temporal resolution of Landsat images is not high enough to capture the short-term events on the Earth. Specifically, the general revisit period of Landsat image is more than 10 days; hence, active fire will be observed by the Landsat satellite with a probability less than 10% (considering the cloud coverage). In addition, the gaps between Landsat images of adjacent time phases and the occurrence of cloud also increase the uncertainty in analyzing the time-series patterns of the land surface. Without using the evidence of active fire, it is not easy to identify the burned scars at the global scale with high confidence due to the wide variety of vegetation types, phenological characteristics, burned-like land covers, and spectral characteristics within a burned scar (char, scorched leaves, or grass, or even green leaves when the fire is not very severe [13]). In this work, the MODIS Vegetation Continuous Fields (VCF) product was applied to discriminate tree-dominated and grass-dominated regions, but the VCF product is neither precise in spatial resolution nor available before 2000, and moreover, two categories are far from enough to separate different burning types. Actually, much prior knowledge can be utilized to improve the accuracy of GABAM, if the globe is carefully divided into intensive regions according to the fire behavior, land cover types, and climate. For instance, most biomass burning in the tropics is limited to a burning season; around 10% of the savanna biome burns every year; burning cropland after a harvest is extremely prevalent, and so on. Consequently, region-specific algorithms should be helpful to improve the accuracy of high-resolution global annual burned area mapping. Furthermore, despite the high correlation between GABAM and Fire_cci, the area of detected BA was generally smaller in GABAM than that in Fire_cci, since some pixels located within a burned area, but not showing a strong burned appearance, were not included in GABAM. This situation can be considered as an underestimation of BA or omission error if only taking into account the connectivity and completeness of burned patches; on the other hand, however, the detailed perimeter of BA from GABAM can be useful for the statistics of the area of biomes actually burned, and therefore to improve the simulation of carbon emissions from biomass burning. In its present form, however, GABAM suffers limitations in the following aspects.

4.1. BA in Agriculture Land

It is difficult to detect BA in cropland with high confidence (low commission error and low omission error) from satellite images:

• Many croplands have comparable spectral characteristics to burned areas when harvested or ploughed.
• The temporal behavior of harvest or burning of cropland is similar to that of grassland fire, e.g., sudden decline and gradual recovery of NDVI, as well as periodic variation of NBR values year after year.
• Different from the wildfires in rangeland and forest, most of the fires in croplands are human-intended stubble burning, and they are commonly small and of a short duration, being difficult to capture by satellite sensors. In this sense, the traditional burned area detection algorithms, which are frequently used to generate BA products from the data source of a medium resolution (e.g., MODIS, AVHRR, MERIS), are likely to have high omission error in croplands for small cropland fire.

Figure 10 shows an example of cropland in Mykolayiv, Ukraine, including the Landsat-8 time-series (Figure 10a–i) and the burned scars mapped by Fire_cci (Figure 10j), GABAM (Figure 10k), and reference data (Figure 10l). Small fire spots, showing a light orange color, can be visually observed from Figure 10a,b,h, but burned scars surrounding these fire spots were not included in the Fire_cci product. On the other hand, without fire evidence or field validation, it is also difficult to tell whether the burned-like surfaces detected by GABAM were false alarms.

Due to these difficulties, discriminating true-burned areas from croplands is not a trivial task, and cropland masks can be employed to remove potential confusions.

4.2. Omission of Observations

Using Landsat images as input data for GABAM, the number of valid observations is a limiting factor for detecting fires, since the active- or post-fire evidence may be omitted or weakened due to the temporal gaps caused by temporal resolution, as well as cloud contamination. Especially in Tropical regions, where vegetation recovery is quite quick after fire, temporal gaps usually result in high omission error. Figure 11 shows an example of omission error in South America. From the CBERS-4 images (Figure 11o–r), a new burned scar, which occurred during 21 August–12 October, can be identified at the center of image patch. However, all Landsat-8 images (Figure 11a–n) acquired between the date interval from 21 September–24November were contaminated by cloud; thus, the region covering this burned scar in these images was masked by the QA band during the process of BA detection.

4.3. Validation

For satellite data product validation, a commonly-used method is to employ higher spatial resolution satellite data. For example, in order to validate the MODIS-derived data product (1-km spatial resolution), Landsat satellite data are commonly used. In this study, however, Landsat images were used as the main reference source to validate the Landsat-derived burned area product. Although the validation process was conducted by independent experienced experts with great caution, relying on Landsat for both product generation and validation limits our ability to assess inaccuracies imposed by the satellite sensor itself, such as radiometric calibration accuracy, spectral band settings, geolocation, and mixed pixels [57]. Accordingly, extensive validation of GABAM is expected to be further performed by professional users.

5. Conclusions

An automated pipeline for generating 30-m resolution global-scale annual burned area maps utilizing Google Earth Engine was proposed in this study. Different from the previous coarse resolution global burned area products, GABAM 2015, a novel 30-m resolution global annual burned area map of 2015, was derived from all available Landsat-8 images, and its commission error and omission error were 13.17% and 30.13%, respectively, according to preliminary global validation. However, completing a thorough validation of the GABAM product is beyond the scope of this study, and users of the product are encouraged to make further validation since GABAM 2015 is publicly available. Comparison with the Fire_cci product showed a similar spatial distribution and strong correlation between the burned areas from the two products, particularly in coniferous forests. The automated pipeline makes it possible to efficiently generate GABAM from the huge catalog of Landsat images, and our future effort will be concentrated on producing long time-series 30-m resolution GABAMs.

One of the limitations of GABAM is the uncertainty in agriculture land, where burned areas can be confused with harvested or ploughed surfaces in spectral and temporal characteristics, without the evidence of fire. In this case, a field survey or very high resolution images are required to achieve promising results. Another limitation of GABAM is that it does not provide the date when the burned area was first detected, because GABAM was generated using the maximum burned probability over the dense time-series. Nevertheless, the approximate date when the burned area occurred can be estimated if the algorithm in this paper is adapted to produce monthly or daily products. The availability of adequate good quality observations remains another limiting factor for Landsat-based BA detection, since the temporal gaps caused by relatively lower temporal resolution and cloud contamination can result in omission errors, particularly in tropic regions. However, for more recent global BA products generation, a combination of Landsat and Sentinel-2 datasets can be a good choice to cope with the observation limit after fixing the inconsistency between these two datasets, as the temporal resolution of combined datasets has become three days since 2016.