30m resolution Global Annual Burned Area Mapping based on Landsat images and Google Earth Engine

Heretofore, global burned area (BA) products are only 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 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 proposed an automated pipeline for generating 30-meter resolution global-scale annual burned area map from time-series of Landsat images, and a novel 30-meter resolution global annual burned area map of 2015 (GABAM 2015) is released. GABAM 2015 consists of spatial extent of fires that occurred during 2015 and not of fires that occurred in previous years. Cross-comparison with 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 error of GABAM 2015 are 13.17% and 30.13%, respectively.


Introduction
Accurate and complete data on 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 (Chuvieco et al., 2016). Remotely sensed satellite imagery has been widely used to generate burned area products. Burned area products at 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. Main global burned area products include GBS (Carmona-Moreno et al., 2005), Global Burned Area 2000 (GBA2000) (Tansey, 2004), GLOBSCAR (Simon, 2004), GlobCarbon (Plummer et al., 2006), L3JRC (Tansey et al., 2008), MCD45 (Roy et al., 2005), GFED (Giglio et al., 2010), MCD64 (Giglio et al., 2016), and Fire cci (Chuvieco et al., 2016).
The recently released Fire cci product is produced based on MODIS images and has the highest spatial resolution (250m) of all the existing global burned area products (Chuvieco et al., 2016;Pettinari and Chuvieco, 2018). However, the requisites of the climate modelling community are not yet met with the current global burned area products, as these products do not provide enough spatial detail (Bastarrika et al., 2011). Imagery collected by the family of Landsat sensors is useful and appropriate for monitoring the extent of area burned and provide spatial and temporal resolutions ideal for science and management applications. Landsat sensors can provide a longer temporal record (from 1970s until now) of burned area relative to existing global burned area products and potentially with increased accuracy and spatial detail in most areas on the earth (Stroppiana et al., 2012). Great importance has been attached to developing burned area products based on Landsat data in the past 10 years (Bastarrika et al., 2011;Stroppiana et al., 2012;Hawbaker et al., 2017a). Up to now, there is no Landsat based global burned area product, however, some regional Landsat burned area 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 time series change detection technique (Goodwin and Collett, 2014). Fire scars are automatically detected and mapped using  (Eidenshink et al., 2007). MTBS products are generated based on the difference of Normalized Burned Ratio (NBR) calculated from pre-fire and postfire images, in which the burned area 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), produced Landsat derived burned area products across the conterminous United States (CONUS) from 1984-2015, and its products have been released in April 2017 (Hawbaker et al., 2017a). The main differences between the MTBS and BAECV is the BAECV products are automatically generated based on all available Landsat images.
In summary, global burned area products are only available at coarse spatial resolution while 30-meter resolution burned area products are limited to specific regions. The majority of coarse spatial resolution algorithms developed to produce global burned area 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 MODIS burned area product (MCD45) is developed from the bi-directional reflectance model-based expectation change detection approach (Roy et al., 2005).
One of the difficulties to produce Landsat based burned area products is that the traditional approaches successfully applied to extract global burned area from MODIS, VEGETATION, etc. dont 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 . Another difficulty is that global 30-meter resolution annual burned areas mapping needs to utilize dense time-series Landsat images, and the required datasets can be hundreds of thousands of Landsat scenes, resulting impractical processing time. Although some researches have been addressed to detect burned area regionally from Landsat time series (Goodwin and Collett, 2014;Hawbaker et al., 2017b;Liu et al., 2018), results of 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 (Gorelick et al., 2017), it is now possible to perform global-scale geospatial analysis efficiently without caring about pre-processing of satellite images. In this study, we focus on an automated approach to generate global-scale high resolution burned area map using dense time-series of Landsat images on GEE, and a novel 30-meter resolution global annual burned area map of 2015 (GABAM 2015) is released.

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 global burned area map and also the completeness of quality assessment, a stratified random sampling method (Padilla et al., 2015;Boschetti et al., 2016;Padilla et al., 2017) 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 type of land cover.
Firstly, the Earth's Land Surface was partitioned based on the 14 land cover classes according to the MCD12C1 product (Friedl and Sulla-Menashe, 2015) of 2012 using University of Maryland (UMD) scheme.
These types were then merged into 8 categories based on their similarities , i.e. Broadleaved Evergreen, Broadleaved Deciduous, Coniferous, Mixed forest, Shrub, Rangeland, Agriculture and Others. Table 1 shows the reclassification rule from UMD land cover types to new classifications. As "Others" category consists of the biomes less prone to fire, only other 7 land cover categories are considered to create the geographic stratifications in this work.  (Giglio et al., 2013), the most widely used inventory in global biogeochemical and atmospheric modeling studies (Giglio et al., 2016). 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 • × 0.25 • GFED grid, and each pixel summed the total areas of BA (BA density, km 2 ) occurred in the grid cell during the whole year. The BA density of GFED4 2015 was then divided into 5 equal-frequency intervals  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 with larger sample sizes allocated to strata with higher BA extent (Padilla et al., 2014). According to the strategy of stratified sampling, 120 samples (24 for each BA density level) were randomly selected to generate training dataset, and 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 randomly stratified sampling, but trying to keep a distance (at least 200km) 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 2 shows the distribution of training and validation samples over the different land cover types.

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.1. 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 Google Earth Engine platform, whose ImageCollection IDs are "LAND-SAT/LC08/C01/T1 SR" and "LANDSAT/LC08/C01/T2 SR". These data have been atmospherically corrected using LaSRC (Vermote et al., 2016), and include a cloud, shadow, water and snow mask produced using Fmask (Zhu and Woodcock, 2014) and Green bands were composited in 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 (Koutsias and Karteris, 2000). 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 confusing pixels which were difficult to identify whether they were burned scars, further check was performed by examining the Landsat images on nearest date of the previous year or higher resolution images on nearest date in 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 (Bastarrika et al., 2011); and burned pixels located in burning flame or covered by smoke were also excluded to prevent potential contamination of burned samples. Land surface reflectance of the collected samples in BLUE, GREEN, RED, NIR, SWIR1, SWIR2 bands were extracted for further analysis.  Burned areas are characterized by deposits of charcoal, ash and fuel, and the reflectance of the burned pixels generally increases along with the wave length while the burned pixels have similar reflectance in 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 (Bastarrika et al., 2014), 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 (Boschetti et al., 2010). Consequently, in this work, we made use of most common spectral indices for Landsat image previously suggested in BA studies, and their formulas are summarized as Table 3. Some of these spectral indices were specifically developed for burn detection as they are sensitive to charcoal and ash deposition, such as normalized burned ratio (NBR (Key and Benson, 1999)), normalized burned ratio 2 (NBR2 (Lutes et al., 2006)), burned area index (BAI (Martín, 1998)), mid infrared burn index (MIRBI (S. Trigg, 2001)). In addition, other indices that are not burnspecific may also be useful to map burned areas when cooperating with burn-specific indices. For instance, although 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 (Stroppiana et al., 2009). The global environmental monitoring index (GEMI (Pinty and Verstraete, 1992)) is an improved vegetation index, specifically designed to minimize problems of contamination of the vegetation signal by extraneous factors, and it is considered very important for the remote sensing of dark surfaces, such as recently burned areas (Pereira, 1999). The soil adjusted vegetation index (SAVI (Huete, 1988)), which is originally designed for sparse vegetation and outperforms NDVI in environments with a single vegetation (Veraverbeke et al., 2012), is also helpful to improve separability of burns from soil and water (Stroppiana et al., 2012). The normalized difference moisture index (NDMI (Wilson and Sader, 2002)), which is sensitive to the moisture levels in vegetation, is also relative to fuel levels in fire-prone areas. We also evaluated the relative importance of 14 Landsat features (8 spectral indices in Table 3 and the surface reflectance in 6 bands of Landsat-8 image) when applied to classify burned areas using random forest algorithm (Pedregosa et al., 2011) (as shown in Figure 3). Figure 3 shows that spectral indices were generally more important than the surface reflectance in most bands, except for NIR band and SWIR2 band, which are sensitive to removal of vegetation cover and deposits of char and ash (Pleniou and Koutsias, 2013). We also found that NBR2, BAI, MIRBI and SAVI had the greatest relative importance, whereas SAVI was not initially developed for burned area detection. However,

Sensitive features for burned surfaces
considering the potential contribution of those features with relative low importance in distinguishing burned scars, all of 14 features were selected as sensitive features to perform global burned area mapping in this study.

Name Formula
Normalized burned ratio is the surface reflectance in SWIR1 band, and ρ SW IR2 is the surface reflectance in SWIR2 band.

Burned area mapping via GEE
In this work, annual burned area map is defined as spatial extent of fires that occurs within a whole year and not of fires that occurred in previous years. Therefore, global 30-meter resolution annual burned areas mapping needs to utilize dense time-series Landsat images, and the pipeline of annual burned area mapping via GEE is described as Figure 4.
As shown in Figure 4, the pipeline mainly consists of three steps, model training, per-pixel processing and burned area shaping, and the following provides more details of each step.

Model Training
The random forest (RF) algorithm provided by GEE were 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.1 and section 2.2).
Random forest classifier with higher number of decision trees usually provides better results, but also causes higher cost in computation time. Since the input features of the algorithm includes the surface reflectance (SR) in 6 bands of Landsat-8 image as well as 8 spectral indices that have high sensitivity to burned surface, we limited the number of decision trees in the forest to 100 for 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.

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 a pixel, the occurrence of a single Landsat satellite could be more than 20 or 40 times (considering the overlap between adjacent paths) within a year, and it would double when contemporary satellites (e.g. Landsat-7 and Landsat-8) were utilized.
However, considering the failure of Scan Line Corrector (SLC) in the ETM+ instrument of Landsat-7 satellite, we only utilized USGS Landsat-8 Surface Reflectance collections ("LANDSAT/LC08/C01/T1 SR" and "LANDSAT/LC08/C01/T2 SR"). The quality assessment (QA) band of Landsat image, which was generated by FMask algorithm (Zhu and Woodcock, 2014), 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. At each pixel, the geometrically aligned dense timeseries 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 current year and that of the previous year.
For the reflectance stack of current year, 8 spectral indices were computed at each time period, and then the trained decision forest classifier in section 2.4.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 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 (N DV I 1 ) and NBR value (N BR 1 ) on that date. However, it is not usually possible to unambiguously separate in a single image the spectral signature of burned areas from those caused by unrelated phenomena and disturbances such as shadows, flooding, snow melt, or agricultural harvesting (Boschetti et al., 2015); the burned scars which occurred in previous years but not yet recovered (particularly in boreal forests) should also be excluded from the annual BA map of current year. In this sense, we also concerned the summary statistics of current year and previous year: N DV I 2 , the maximum NDVI value within the couple of years (current year and previous year); t 2 , the date of N DV I 2 ; and N BR 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.
1. N DV I 2 > T N DV I , the maximum NDVI value within the couple of years should be greater than a threshold T N DV I . 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 appear like burned but in fact were just lacking vegetation.
2. N DV I 2 − N DV I 1 > T dN DV I , the difference between the maximum NDVI and the NDVI when the pixel was most like burned scar should be greater than a threshold T dN DV I . This constraint ensures an evidence of vegetation decrease when burn happened.
3. N BR 2 − N BR 1 > T dN BR , the NBR value of a burned pixel should be less than the minimum NBR of the previous year, and the threshold T dN BR 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.
4. t 1 > t 2 or t 2 − t 1 > T DAY , the most flourishing date of vegetation should be earlier than the burning date, or the lagged days should be less than a threshold T DAY . For tree-covered surface, it usually takes a long time for the vegetation to recover more flourishing than the previous year, thus the burnlike 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 days is also acceptable.
We named the first two constraints as "NDVI filter", the third and fourth ones as "NBR filter" and "temporal filter", respectively. In this work, the thresholds in above constraints were chosen empirically, context, those pixels with annual burned probability greater than or equal to 0.95 ("probability filter") were selected as seeds for region growing.

Burned Area Shaping
In this step, a region growing process was employed to shape the burned areas. Region growing has proved to be necessary for BA mapping in many studies (Bastarrika et al., 2011;Stroppiana et al., 2012;Hawbaker et al., 2017b), 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 (Stroppiana et al., 2012)), 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 confusing surfaces, e.g. shadows, 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, and by ignoring small fires with areas less than 1 ha (Laris, 2005), those fragmentary components (smaller than 11 pixels), which included most false seed pixels, were removed. Finally, an iterative procedure of region growing were 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 can 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.  uct, 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. 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 Landsat image (Figure 7b), yet the BAs in Figure 7c occupy 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 shows finer boundaries of BAs, compared with that of the Fire cci product.  Figure 9 shows the proportion of BA in 0.25 • × 0.25 • grids of different land cover categories in Table 1, for 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 after (July 26th, 2015) fire, respectively, displayed in false color composition (red: SWIR2 band, green: NIR band and blue: GREEN band); 7c shows the burned areas of annually composited Fire cci product, and 7d shows the burned areas generated by proposed method. and the coefficient of determination (R 2 ) for each land cover category (Figure 9a-9h) and for 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 to 10 or more) having the same proportion values.

Regression analysis
According to Figure 9, the intercept values of the estimated regression lines were close to 0 while the slopes were lower than 1, showing that GABAM 2015 generally underestimated  burned scars compared with Fire cci product. Moreover, the distribution and color of scatters in Figure 9i also show that a large number of grids were considered to have higher burned proportion by Fire cci product than by GABAM. The main reason for the inconsistence can be attributed to the difference in 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 • × 0.25 • grids were occupied by more than 90% BA in GABAM while grids with high proportion of BA were more common in Fire cci product. (R 2 = 0.62), and lowest strengths in agriculture land (R 2 = 0.31) and "Others" category (R 2 = 0.07).
In "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 owing to the uncertainty of both products, which will be further discussed in section 3.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.

Data sources
Accuracy assessment was carried out according to the 80 validation sites which were created in section 2.1, and the reference data were selected in these sites from multiple data sources, including fire perimeter datasets and satellite images. Commonly, when satellite data are used as reference data, they should have higher spatial resolution than the data used to generate the BA product . For Landsat BA product, however, access to global higher-resolution time-series satellite data is difficult, and a thorough validation of Landsat science products can be completed with independent Landsat-derived reference data while strengthened by the use of complementary sources of high-resolution data (Vanderhoof et al., 2017).
Consequently, in this study, some publicly available satellite images of higher-resolution were included in the validation scheme, while Landsat was the majority of validation data source. Specifically,  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 U.S., 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 perimeters of burned area, respectively. The characteristics of CB4 and GF1 are illustrated in Table 4. Note that the size of validation site varied by the type of data source, i.e. a WRS-II frame (about 185km × 185km) for Landsat images, a scene for CB4 images (about 120km × 120km) and a box of 100km × 100km for GF1 images. Using Landsat frames or image scenes as a unit of validation site is convenient for data downloading and processing; we chose a smaller box for GF1 to improve the data availability considering the extents of GF1 frames or scenes are not fixed due to the long orbital return period.

Reference data generation
In each validation site, all the available image scenes (LC8 1 , CB4 2 or GF1 3 ) acquired in 2015 were used.
LC8 images were ortho-rectified 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 following steps.

Preprocessing
All the images utilized to generate BA reference data were spatially aligned with mean squared error 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 (Long et al., 2016) was applied to ortho-rectify the time-series GF1 images, taking the LC8 panchromatic images (spatial resolution is 15 meters) 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 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 software were 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.3 were utilized in SVM for each LC8 image pair; but for CB4 and GF1 images, features used for classification consisted of the digital number (DN) values in four bands of an image pair (totally 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 software.

Reviewing and manually revision
The result of 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 U.S., 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 of small fires.

Validation results
To assess the accuracy of GABAM 2015, a cross tabulation (Pontius and Millones, 2011) between the pixels assigned by in our BA product and in the reference data was computed to produce the confusion matrix for each validation site. Afterwards, the global cross tabulation (Table 5) was generated by averaging all the cross tabulations. Finally, three statistics, i.e. commission error, omission error and overall accuracy, can be derived from the confusion matrix: • 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 12 )/(X 11 + X 12 + X 21 + X 22 ), the ratio between the area classified correctly and the total area to evaluate.
According to 1. In the validation sites located in tropical zones, clear burned evidences were frequently missed by Landsat sensor due to the quick recovery of the vegetation surface. This point will be further discussed in section 3.4.
2. Some pixels located within a burned area, but not showing 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 in those validation sites using MTBS perimeters, e.g. E c and E o of the validation site in Figure A.16 were 1.45% and 67.97%.

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 Landsat satellite with 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 to analyze the time-series patterns of land surface. Without using the evidence of active fire, it is not easy to identify the burned scars at global-scale with high confidence due to the wide variety of vegetation types, phenological characters and burned-like landcovers, and spectral characteristics within a burned scar (char, scorched leaves or grass, or even green leaves when the fire is not very severe (Bastarrika et al., 2011)). In this work, 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 year and, moreover, two categories is 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 behaviour, 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-specified algorithms should be helpful to improve the accuracy of high-resolution global annual burned area mapping. Furthermore, despite of 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 strong burned appearance, were not included in GABAM. This situation can be considered as 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 to statistics the area of biomes actually burned, and therefore to improve the simulation of carbon emissions from biomass burning.
At present form, however, GABAM suffers limitations in the following aspects.

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: • A lot of croplands have comparable spectral characteristics to burned areas when harvested or ploughed.
• The temporal behaviour of harvest or burning of cropland is similar to that of grassland fire, e.g.
sudden decline and gradual recovery of NDVI, and 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 short duration, difficult to be captured by satellite sensors. In this sense, the traditional burned area detection algorithms, which are frequently used to generate BA products from data source of the middle 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-10j) and the burned scars mapped by Fire cci (Figure 10k) and GABAM (Figure 10l). Small fire spots, showing light orange color, can be visually observed from Figure 10a, 10b and 10i, but burned scars surrounding these fire spots were not included in 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.

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 weaken 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 11p-11s), a new burned scar, which occurred during August 21 to October 12, can be identified at the center of image patch.
However, all the Landsat-8 images (Figure 11a-11o) acquired between the date interval from September 21 to November 24 were contaminated by cloud, thus the region covering this burned scar in these images were masked by QA band during the process of BA detecting.

Validation
For satellite data product validation, a commonly used method is to employ higher spatial resolution satellite data. For example, in order to validate MODIS derived data product (1 km spatial resolution), Landsat satellite data is commonly used. In this study, however, Landsat images were used as the main reference source to validate 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 (Strahler et al., 2006). Accordingly, extensive validation of GABAM is expected to be further performed by professional users.

Conclusions
An automated pipeline for generating 30m resolution global-scale annual burned area map 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 year, was derived from all available Landsat-8 images, and its commission error and omission error are 13.17% and 30.13%, respectively, according to global validation. Comparison with 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 huge catalog of Landsat images, and our future effort will be concentrated to produce long time-series 30m resolution GABAMs.  Table A.7 summarizes the information of these validation sites, including the location, source of reference data, commission error, omission error and overall accuracy.

ID Location
Reference data E c (%) E o (%) A o (%) Figure   1 China GF1