Burned Area Mapping in the North American Boreal Forest Using Terra-MODIS LTDR (2001–2011): A Comparison with

An algorithm based on a Bayesian network classifier was adapted to produce 10-day burned area (BA) maps from the Long Term Data Record Version 3 (LTDR) at a spatial resolution of 0.05° (~5 km) for the North American boreal region from 2001 to 2011. The modified algorithm used the Brightness Temperature channel from the Moderate Resolution Imaging Spectroradiometer (MODIS) band 31 T31 (11.03 μm) instead of the Advanced Very High Resolution Radiometer (AVHRR) band T3 (3.75 μm). The accuracy of the BA-LTDR, the Collection 5.1 MODIS Burned Area (MCD45A1), the MODIS Collection 5.1 Direct Broadcast Monthly Burned Area (MCD64A1) and the Burned Area GEOLAND-2 (BA GEOLAND-2) products was assessed using reference data from the Alaska Fire Service (AFS) and the Canadian Forest Service National Fire Database (CFSNFD). The linear regression analysis of the burned area percentages of the MCD64A1 product using 40 km × 40 km grids versus the reference data for the years from 2001 to 2011 showed an agreement of R2 = 0.84 and a slope = 0.76, while the BA-LTDR showed an agreement of R2 = 0.75 and a slope = 0.69. These results represent an improvement over the MCD45A1 product, which showed an agreement of R2 = 0.67 and a slope = 0.42. The MCD64A1, BA-LTDR and MCD45A1 products underestimated the total burned area in the study region, whereas the BA GEOLAND-2 product overestimated it by approximately five-fold, with an agreement of R2 = 0.05. Despite MCD64A1 showing the best overall results, the BA-LTDR product proved to be an alternative for mapping burned areas in the North American boreal forest region compared with the other global BA products, even those with higher spatial/spectral resolution.


Introduction
Forests provide a range of beneficial services, such as habitats for organisms, climate regulation, biodiversity, soil conservation and watershed protection.Currently, forests cover approximately 30% of the total Earth's land surface [1].The boreal forest extends across the northern continents of North America, Europe and Asia.In certain areas, the boreal forest extends more than 2000 kilometers from North to South.This forest is the largest biome on our planet, and stores approximately 40% of the total carbon in terrestrial ecosystems [2,3].Boreal forests store carbon in surface vegetation, but in addition they have accumulated carbon for millennia in associated soils, permafrost deposits, wetlands and peat lands that has been sometimes underestimated [4,5].
Since the advent of Earth observation from space, remote sensing has become a valuable tool for the scientific community and natural resource managers, because it enables the periodic collection of data in different bands of the electromagnetic spectrum on very wide and inaccessible areas of the Earth [6].
Recently, the international community has made considerable progress in the development of initiatives for global monitoring because such information is critical to global issues such as climate change and land use change [7].The Moderate Resolution Imaging Spectroradiometer (MODIS) and Long Term Data Record (LTDR) programs from the National Aeronautics and Space Administration (NASA) and the European VEGETATION program are examples of such initiatives that provide daily imagery and higher level land and atmosphere products for global mapping.
 The MODIS program is part of the Earth Observation System (EOS) that provides daily data in 36 spectral bands (from 0.412 to 14.235 μm) with three native spatial resolutions (250 m, 500 m and 1,000 m).Daily observations are acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard both the Terra and Aqua satellite platforms and transferred using the Tracking and Data Relay Satellite System (TDRSS).Subsequently, the data from the observations are sent to the Earth Observation System (EOS) Data and Operations System (EDOS).The Level 1A, Level 1B, geolocation and cloud mask products and the higher-level MODIS land and atmosphere products are produced by the MODIS Adaptive Processing System (MODAPS) and are then parceled out among three Distributed Active Archive Centers (DAACs) for distribution [8]. The VEGETATION program is the result of collaboration among various European partners (Belgium, France, Italy, Sweden and the European Commission) and is conducted under the supervision of the French National Center for Space Studies (CNES).This program consists of two observation instruments in orbit, VEGETATION 1 and VEGETATION 2, and ground facilities.
The VEGETATION 1 instrument, launched on 24 March 1998, is aboard the SPOT 4 satellite.VEGETATION 2 is aboard the SPOT 5 satellite, which was launched on 4 May 2002 [9].The aim of the VEGETATION instruments is to provide daily global images for the observation of long-term regional and global environmental changes at a spatial resolution of 1,165 m × 1,165 m using the following four spectral bands: 0.43-0.47µm (blue), 0.61-0.68µm (red), 0.79-0.89µm (near infrared (NIR)) and 1.58-1.75µm (middle infrared (MIR)).The VEGETATION images are processed, archived and distributed by the Flemish Institute for Technological Research (VITO) and are available at no cost on its website [10].Fire is an important disturbance regime in forests and is the primary process that organizes the physical and biological functions of the boreal biome [15,16].Due to the slow vegetation recovery process, an increase in boreal forest burning can substantially modify the carbon balance [17].
Burned area mapping provides information on fire seasonality, frequency of occurrence, location and quantification of the burned area, which is essential for developing environmental management policies.However, recent published studies on BA products derived from the daily coarse spatial resolution imagery [18][19][20][21][22][23] on a continental/global scale do not provide a quantitative comparison that spans the 2001-2011 decade at the final product level.Besides, as far as we know, there has not been any burned area products developed using the LTDR dataset generated from MODIS imagery.The BA-LTDR product (from AVHRR imagery) was only compared with other burned area products before the year 2000.Therefore, the goals of the present study are the following: (1) To generate a burned area time series for the North American boreal region (2001-2011) from the Terra-MODIS LTDR dataset Version 3 that extends the previous time series (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) from the NOAA-AVHRR LTDR dataset Version 3 using the same methodology [19].
(2) To evaluate the accuracy of the 2001-2011 time series using reference data from the Alaska Fire Service (AFS) and the Canadian Forest Service National Fire Database (CFSNFD).(3) To compare BA-LTDR product with the Collection 5.1 MODIS Burned Area (MCD45A1) product [24][25][26], MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1) [27] and the Burned Area GEOLAND-2 (BA GEOLAND-2) product derived from VEGETATION, which was an improvement of the Global Burned Area (GBA2000) product and the burned area product developed by the consortium of three academies under the supervision of the Joint Research Centre of the European Commission (L3JRC) [18].

Study Region
The study area selected for this analysis is shown in Figure 1.This site was selected because it was the original location of the BA-LTDR development, and there is a dataset of independent reference information for this region that allows for the validation of the BA products against field survey data.These reference data are provided by the Alaska Fire Service (AFS) [28] and the Canadian Forest Service National Fire Database (CFSNFD) [29], which collect forest fire data from various sources, including fire locations and perimeters.The boreal forest exists as a nearly continuous belt of coniferous trees across North America and Eurasia.This biome covers approximately 12 million square kilometers, two-thirds in Russia and Scandinavia and the remainder in Canada and Alaska.The North American boreal forest region stretches from coast to coast across Canada and into Alaska.This region represents approximately 36% of Earth's boreal forest; Canada contains 25%, and the remaining 11% is located in Alaska [30].
Frost resistant conifer species dominate in boreal forests, typically spruce, pine, larch and fir.The spruces are the most widespread and characteristic of boreal forests.The soil in these regions is typically poor because conifer needles acidify the soil and this has accumulated annual increments of carbon for millennia [4].Hence, the biodiversity of the boreal forest is quite low, although there are more species in the North American boreal forest than in other regions of the boreal belt [31].This region has a continental climate that is marked by considerable variations in temperature.Summers are short and warm, with a monthly mean temperature above 10 °C during at least one of the four months and a maximum value of approximately 25 °C.Over the long winters, however, the temperatures vary from −20 to −40 °C for at least six months [30].
Forests are constantly exposed to and modified by natural disturbances, such as fire, insects and diseases.Natural disturbances are an essential part of the forest renewal process.Wildfires play a critical role in maintaining the ecological integrity of boreal forests because they are part of the natural evolution of the forest.Fire not only stimulates the regeneration of many plant species but also recycles phosphorus and removes accumulated organic matter [32].
However, in the last few decades, fires have become more severe and the burned area across the boreal forest has increased because of a warmer climate [5,33,34].For example, in 2010, 7,319 forest fires were reported across Canada, which is approximately equal to the 10-year average (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009).Although the number of fires was the same, the area burned was much higher in 2010, with three million hectares burned.This increase in the extent of forest fires in boreal regions is producing a vast increase in CO 2 .Nevertheless, the impact of the boreal forest fires on climate change is not clear, and some authors suggest that it may not accelerate climate warning [35] or perhaps it could cool the high-latitude atmosphere [34].
According to the AFS and the CFSNFD, the burned area in the North American boreal forest varies substantially from year to year [36].For the last three decades, years with high fire activity, such as 1989, 1994, 1995, 1998, 2002, 2004, 2005 and 2010, registered up to four million ha of burned area in a single year, whereas other years, such as 1984, 1985, 1987, 1997, 2000 and 2001, registered as low as 1.20 million ha.
During the summer season, the months of June and July have the greatest fire incidence, accounting for two-thirds of the total burned area.Fires greater than 200 ha represent approximately 5% of the total number of fires and account for more than 95% of the total area burned.The number and area of local fires per year greatly depends on latitude and distance from the ocean [37].In Alaska, the fire frequency is highest for the closed-canopy boreal forests, followed by spruce-lichen woodland and finally tundra.This trend is similar in Canada, which shows a latitudinal decrease in fire frequency from South to North, following the progression from closed-canopy boreal forest to tundra, with spruce-lichen woodland as a transitional zone.Also, this latitudinal zoning follows latitudinal climate variability [38].
The primary driver of fire activity in the boreal forests of North America is the climate [39], but other factors, such as topography, vegetation and land use, also have a vital influence on spatial burn patterns [40].In interior boreal regions, where large fires are most frequent, the influence of vegetation on fire size appears to be linked with flammability [41].Topography also controls fire severity.As the topography becomes more complex, the heterogeneity of fire severity increases, and the influence of vegetation decreases [42].The natural burn/regeneration cycle is less over 75, reaching 100 years, depending on the species.For example, Spruce regeneration under pine on poor soil is about 90 years after burning [43].

Burned Area Products
Table 1 shows the primary characteristics of the four BA products analyzed in this study: MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR.Annual maps for the four evaluated remote sensing BA products and for the reference data were built at their respective native spatial resolutions for the study region for the period 2001-2011 and were re-projected to a Lambert Conformal Conic projection to minimize distortions and maintain equivalent square shapes [21].In the subsections below, we summarize the algorithm used in each product and the process used to build the annual maps.[19] and this study

MCD45A1
The MCD45A1 product is distributed monthly and is part of the MODIS 5 collection [25].This product uses daily data from both the Terra and Aqua platforms and is generated from a time series of daily land surface reflectance data.The BA algorithm is a bi-directional reflectance model-based change detection algorithm.The change method is applied to geo-located pixels over a long time series of reflectance observations.The algorithm begins by preprocessing the seven MODIS land surface reflectance bands, which are corrected for atmospheric effects.This step provides a maximum of one observation per pixel per day, although cloud, snow and water pixels may remain.The analysis of the MODIS reflectance bands shows burned/unburned discrimination.A daily reflectance is predicted based on the reflectance sensed for the previous 16 days.This information is used to model a bidirectional reflectance distribution function (BRDF) that allows for the management of the angular variations of reflectance.This model is used to predict changes in the reflectance by estimating subsequent observations in time.Next, a statistical measure (Z-score) is computed from MODIS bands 2 and 5 for each geo-located pixel to determine whether the difference between the observed and predicted reflectance indicates that the pixel has burned.This decision is based on a threshold (Zthre).Finally, the algorithm reduces commission errors by selecting from candidate pixels those that provide persistent evidence of fire occurrence.A detailed description of the algorithm applied to obtain the MODIS BA products can be found in the Roy et al. 2005 [44].
The monthly GeoTIFF subsets for the study period (2001-2011) that contain the burn date and the quality assessment pixel, with a 500 m × 500 m spatial resolution in a Lat-Long projection derived from MCD45A1 for the regions of Alaska and Canada, were downloaded from the MODIS website [45] and used to generate monthly BA maps.We considered all of the confidence levels (1-4) of the detection.These monthly maps were re-projected to a Lambert Conformal Conic projection with a pixel size of 500 m × 500 m and combined to generate the annual maps of burned/non-burned pixels.

MCD64A1
MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1) is an alternative burned area product based on an automated method using 500 m MODIS imagery and 1 km MODIS active fire observations [27].The algorithm applies a hybrid approach in several steps for mapping post-fire burned areas.First, it builds time series of valid daily band 1 (0.620-0.670 µm), 5 (1.230-1.250µm) and 7 (2.105-2.155µm) reflectances for each pixel.Then, a daily burn-sensitive vegetation index is computed (VI = (B5 − B7)/(B5 + B7)), so that the abrupt decrease in VI is used as the most important signal to detect burned areas.The algorithm analyzes daily VI data time series in two time spans (10 days before and 10 days after), in order to calculate descriptive statistics, and to define a measure of temporal separability, yielding composites images of these parameters.A mask with active fires is applied to filter pixels providing a mechanism to select a training set of burned/unburned pixels.Next, the algorithm computes conditional probabilities densities for each land cover class existing in the MODIS tile, and then calculates the posterior burned probability.The result of this step is an initial classification (burned/unburned) for each valid 500 m pixel.Finally, the algorithm uses neighbor conditions to refine the initial classification, analyzing the burned/non burned state of the eight adjacent neighbors.
The MCD64A1 monthly subset were downloaded from [46] and used to generate monthly BA maps for (2001-2011) for the study area.This product contains five layers in a sinusoidal projection: the burn date, the burn date uncertainty, the quality assurance, and first day and last day.The monthly maps were re-projected to a Lambert Conformal Conic projection with a pixel size of 500 m × 500 m and combined to generate the annual maps of burned/non-burned pixels.

BA GEOLAND-2 Version V1
There are two versions of the BA product from the SPOT-VEGETATION dataset.The first version, -V0‖, is an implementation of the L3JRC BA algorithm [18] for the region of continental Africa and was developed under the Global Burned Area (GBA) 2000 project [47].
The -V1‖ version (implemented in the GEOLAND-2 project) was modified with aggregation into a ten-day product with near-real-time dissemination for application on the global scale.This version enhances previous products by including data outside the primary fire season, shortening the preprocessing steps, improving the land-water mask and providing additional years than those available for the previous L3JRC product [48].
In the -V1‖ version, reflectance data are preprocessed by applying a cloud and snow mask based on a threshold at 0.45 µm and 1.66 µm wavelengths, respectively, and a fire smoke mask at 0.45 µm.Masks are combined to create a masked pixel product that is computed on a daily basis.The primary processing algorithm makes use of a temporal index, I, that allows for the classification of burned and unburned surfaces [18].A pixel is classified as -burned‖ if the pixel value in the 112 × 112 pixel array window, I, is lower than three mean values minus two standard deviations.Two further checks are made on reflectance values in the 0.83 μm layer, where the pixel reflectance must be less than 0.13 µm , and the 1.66 µm layers, where the pixel reflectance value must be greater than 0.125 μm.
After a 10-day period, a seasonal metric calculation is computed using four available successive 10-day periods of processed BA data.This method allows the moving average to be calculated.Furthermore, two consecutive moving average values must be available.Because of this requirement, significant gaps were observed in the product [49].
The 10-day BA maps derived from BA GEOLAND-2 version -V1‖ for the North American continent were obtained using the free Vegetation extraction software tool (VGTExtract) by VITO [50].Subsequently, these 10-day maps were re-projected to a Lambert Conformal Conic projection with pixel size 1,000 m × 1,000 m and combined to generate the annual maps of burned/non-burned pixels.
3.1.4.Generating the BA-LTDR Product from Terra LTDR LTDR Version 3 data was composited into 10-day maximum brightness temperature (MODIS band 31; 11.03 μm) images [19].These images were generated for the 2000-2012 study period over the entire study area.The BA-LTDR algorithm identifies pixels with high brightness temperature and low infrared values.Then, it calculates the differences between spectral variables within pre-fire and post-fire 10-day composite periods of potential fire dates.The Burned Boreal Forest Index (BBFI) was computed considering the MODIS Brightness Temperature band T 31 (11.03 μm) instead of the band T 20 (which is the equivalent to the AVHRR Brightness Temperature band T3 (3.75 μm)) because the band T 20 is unavailable in the published LTDR dataset Version 3 from MODIS.To detect the energy released by active fires, smoke, charcoal and fire scars, middle infrared sensing is the most likely approach.The mid infrared window (at temperature 500-1,000 K) has maximum spectral response for temperatures in this range.Nevertheless, the thermal infrared (TIR) (8-12 μm) band transmits energy for fires below 500 K [51].In order to evaluate if either band could be included to compute the BBFI and as LTDR dataset contains T3 and T4 bands before the MODIS era, the timing of the fire events for the year 1998 was computed.A mean delay was observed for the BBFI based on the T4 band vs. T3 of 1.6 10-day composites with a standard deviation of 1.5, which has hardly any impact on the algorithm.
A set of statistical variables was computed for the MODIS band 1 ρ 1 (0.645 μm) and band 2 ρ 2 (0.858 μm), the brightness temperature from band 31 T 31 (11.03 μm), the modified BBFI and the Global Environmental Monitoring Index (GEMI) spectral indices within each pixel for each 10-day composite of daily images for the fire event year (Y), the year before (Y − 1) and the year after (Y + 1).
The BA-LTDR algorithm searched the ten-day composite single date (T[x]) for the maximum BBFI value.From that single date, the five pre-fire 10-day composites (T[x − 5, x − 1]) and the five post-fire (T[x + 1,x + 5]) composites were also specified to calculate certain statistics for those variables used in the Bayesian network with the goal of computing the burn probability of a given pixel.
The pre-fire and post-fire statistics were computed for several indices, but the GEMI median was the best detector of differences between the two periods.The following thresholds for these variables were applied to avoid false detections: low values for the maximum BBFI value in (Y − 1), high values for the maximum of band 31 T 31 and the BBFI value in (Y) and (Y + 1), low minimum ρ2 values from the post-fire period in (Y) and a small difference in minimum ρ2 values from the post-fire period between (Y) and (Y − 1).The algorithm thresholds selected the candidate pixels to be evaluated.The Bayesian network classifier based on the training set, within the WEKA machine learning package [52], calculated the individual conditional probability of each selected pixel to belong to the burned and non burned classes [53].Following Bayes' theorem [54], the joint probability for a given class was computed as a product of individual conditional probabilities.The resulting probability value was computed by subtracting the normalized probability for the unburned class from the normalized probability for the burned class, where higher positive values indicate a higher probability to be classified as a burned pixel.
Finally, neighborhood conditions were applied for a spatial coherence to minimize commission and omission errors; isolated burned pixels were downgraded to the unburned class, and pixels with a low burned probability but surrounded by burned neighbors were upgraded to the burned class.
The 10-day probability maps with a 0.05° spatial resolution in a Lat-Long projection were generated by the Bayesian network algorithm from the Daily Global LTDR dataset [14].For these probability maps, a pixel was labeled as -burned‖ when its probability value was greater than 0. The resulting 10-day burn scar maps were re-projected to a Lambert Conformal Conic projection with a pixel size of 5 km × 5 km and combined to generate the annual maps of burned/non-burned pixels.

Accuracy Assessment
Burned area recorded by the Alaska Fire Service (AFS) and the Canadian Forest Service National Fire Database (CFSNFD) were used as ground truth.These databases represent the most complete information source of burned area in North America since the 1950s, and they have been often used as reference data to assess the accuracy of different burned area remote sensing products in this region [19][20][21][22][23].
The reference data were downloaded from the AFS website [28] and the CFSNFD website [29].The polygons of all of the fires registered in the CFSNFD and AFS databases that occurred from 2001 to 2011 were used to produce the annual vector layers of ground truth.Next, these vector layers were re-projected to a Lambert Conformal Conic projection with a pixel size of 500 m × 500 m to generate annual ground truth maps.In order to determine how the pixel is assigned a burned/non burned value, the method of the maximum area within the pixel was used.
All of the annual maps (BA products and reference) were clipped to the defined area (70°N, −168.5°E,45°N, 50°E) and masked with political boundaries for both the Canadian and Alaskan regions to precisely delimit the study region.The maps were re-sized to a Lambert Conformal Conic projection with a pixel size of 5 km × 5 km by pixel aggregation (an aggregated pixel represents the percentage of BA at the sub-pixel level) to compute the spatial and temporal accuracy.
The temporal accuracy of each product has been assessed considering the total annual burned area calculated in the process.Then, a timing distribution was represented on a chart.The next step was the computation of the burned area proportions of each product map vs. the reference map on a grid of cells to assess the spatial accuracy.The cells of 40 km × 40 km were uniformly distributed over the entire study region.A simple regression analysis of these proportions was carried out [55].Then, the error matrix (P ij ) [56] and its derived indices of interest (commission and omission errors) [57] were computed considering a 40 km pixel size [20] to prevent geo-reference-derived errors.The elements of the error matrix (P ij ) of a map of N pixels are calculated by the following equation: where (p ij ) k represents the percentage of pixel k classified as belonging to class C i but really belonging to class C j .In general, each pixel k will contain classified percentages of every class (c 1k , c 2k , …, c ck) and reference percentages of every class (r 1k , r 2k , …, r ck ).Moreover, the calculation of the elements of the error matrix of every pixel would imply solving a system of linear equations with more unknown quantities than equations (multiple solutions).A reasonable supposition to make such a system compatible (one single solution), would consist of comparing, for every class C i within each pixel k, the percentages classified (c ik ) and the reference percentages (r ik ) and assign the least value as the percentage of the class correctly classified (p ii ) k [58].The absolute value of the difference is distributed either in the percentages of class C i classified incorrectly (p ij ) k , or in the percentages of the other classes classified as class C i (p ji ) k .For the case of a dichotomized classification, as is the case of the evaluation of burned surface in this study, where only two classes of interest are involved, -burned‖ (C 1 ) and -non-burned‖ (C 2 ), the system of equations above is resolved for every pixel k by means of a simple algorithm [59]: Next, the Pareto boundaries were computed according to [60], in order to estimate if the accuracy of each burned area product could be limited by their spatial resolution.The original reference maps of burned areas (500 m spatial resolution) were resized to 1 km and 5 km of spatial resolution with pixel aggregate resampling.Each pixel of one of these resized maps contains a percentage p of burned area (p: [0.00-1.00]).Then, for each resized map, a set of N different dichotomic maps were built considering a pixel as burned in map i when p ≥ i/N (i = 1, …, N) and as non-burned otherwise.An error matrix was obtained for each dichotomic map and the commission/omission errors were computed.The set of the calculated commission/omission errors defined the Pareto boundary of the resized map (it represents the minimum omission and commission errors that could be attained at that resolution).Pareto boundaries with 1 km and 5 km of resolution were computed for each year and for the entire period (2001-2011), and were used to assess the performance of the respective burned area products with the same spatial resolution.The Pareto boundary of 500 m is represented by the Cartesian axes.
Finally, we analyzed the capabilities of each BA product to detect large burned areas.All of the burned areas greater than 100,000 ha were identified from the reference data.A set of sub-scenes (200 km × 200 km) that included at their center one (or several) of these large burned areas was created for all of the products and reference data.The total amount of BA was computed for each sub-scene.The detection of a minimum of 10% of burned pixels inside a sub-scene, with respect to the reference data, was the considered criterion to identify the presence of a large BA.

Annual Burned Area Distribution in the Period 2001-2011
The total BA computed in the 2001-2011 period with the reference data was approximately 28.56 million ha, with the year 2004 exhibiting the highest fire activity (5.62 million ha) and 2001 the lowest, with only 0.62 million ha burned.Along with the reference data, Figure 2 shows the annual distribution of the total BA in North America from 2001 to 2011 for the products analyzed: MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR.The total BA estimate over the entire study period for each of these analyzed products was 12.41, 19.68, 138.43 and 19.41 million ha, respectively.The MCD45A1, MCD64A1 and BA-LTDR products underestimated the total BA, whereas the BA GEOLAND-2 product significantly overestimated it.The determination coefficient between the annual burned area estimates for each product in relation to the reference data in the period considered was 0.90, 0.99, −0.08 and 0.97, respectively.

Evaluating the Spatial Accuracy of the Burned Area Products
For each product, Table 2 and Figure 3 show the results of the accuracy of the BA estimate obtained from the scatter plots of the different products versus the reference data using 40 km × 40 km grids.The average determination coefficient (R 2 ) over the 2001-2011 period for the MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR products was 0.67, 0.84, 0.05 and 0.75, respectively, with a slope of 0.42, 0.76, 0.22 and 0.69, respectively.For the MCD45A1 product (excluding the year 2001), the value of R 2 was between 0.56 (2002) and 0.76 (2011) and the slope was between 0.34 (2005) and 0.54 (2006).For the MCD64A1 product (excluding the year 2001), the value of R 2 was between 0.73 (2007) and 0.91 (2004,2009) and the slope was between 0.57 (2007) and 0.84 (2009).For the BA GEOLAND-2 product, R 2 was less than 0.20 in all the years and, therefore, the slope was not significant.For the BA-LTDR product, R 2 was in the range of 0.33 (2006) to 0.86 (2004,2010,2011) and the slope was in the range of 0.31 (2006) to 0.92 (2010).
Table 3 shows the commission and omission errors derived from the error matrix for each year.The average commission errors for the MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR products were 0.19, 0.11, 0.94 and 0.18, respectively, whereas the omission errors were 0.65, 0.39, 0.70 and 0.44, respectively.For the MCD45A1 product (excluding the year 2001), the commission error range was between 0.07 (2011) and 0.38 (2010) and the omission error range was between 0.57 (2003) and 0.74 (2009).For the MCD64A1 product (excluding the year 2001), the commission error range was between 0.06 (2004) and 0. 25 (2006) and the omission error range was between 0.28 (2004) and 0.50 (2007).For the BA GEOLAND-2 product, all of the commission errors were above 0.86 and the omission errors ranged from 0.53 (2001) to 0.78 (2011).For the BA-LTDR product, the commission errors were in the range of 0.08 (2004) to 0. 44 (2001) and the omission errors were in the range of 0.32 (2010) to 0.64 (2006).Figure 4 shows the Pareto boundaries (1 km and 5 km of spatial resolution) of burned areas in the North America region for the entire period (2001-2011) computed from the reference maps with n = 10, and the commission/omission errors for each burned area product (data from Table 3).This figure shows than commission/omission errors of BA-LTDR are near from its Pareto boundary (5 km) while commission/omission errors of GEOLAND-2 are farther from its Pareto boundary (1 km).The Pareto boundary of the MODIS burned area products (MCD45A1 and MCD64A1) is equal to the Coordinate axes, and MCD64A1 are closer than MCD45A1 with respect of their Pareto boundary.

Detection of Large Fires (>100,000 ha)
Table 4 shows the results of the study of the different 200 km × 200 km sub-scenes containing large fires (>100,000 ha) that were recorded in the period from 2001 to 2011 in North America.The total BA encompassed by all of these sub-scenes was 8.6 million ha (approximately 30% of the total area burned across North America in that period).The average relative percentage of burned area estimated for the MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR products for all of these sub-scenes with respect to the reference data was 41.2%, 81.9%, 35.8% and 69.2%, respectively.
The BA recorded was less than 10% in none of the MCD64A1 sub-scenes, three of the MCD45A1 sub-scenes (2002_4, 2007_2 and 2011_2) and in two of the BA-LTDR sub-scenes (2006_1_2 and 2007_2).For the BA GEOLAND-2 product, this area was over 20% in every case, reaching values of 166% (2007_2) and 206% (2010_6).The spatial distribution of the burned pixels for the MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR products (in their native resolution) for each of the considered sub-scenes along with the perimeter of the reference data are shown in Figure 5.

Discussion
The recently published LTDR dataset Version 3 (for the year 2000 and later) did not incorporate the MODIS band equivalent to the AVHRR band T 3 (3.75μm); therefore, to apply the same methodology, it was necessary to modify the original BA algorithm.The MODIS Brightness Temperature band 31 T 31 (11.03 μm), which is equivalent to the AVHRR band T 4 , was used instead.The results were similar in terms of estimation and accuracy to those obtained in the literature [19] for the 1984-1998 period.
We compared global BA products (MCD45A1, MCD64A1 and BA GEOLAND-2) and the BA-LTDR product from 2001 to 2011 (when all the products were available) for the North American boreal forest region with the reference data (AFS and CSFND).The temporal accuracy analysis showed that the annual BA distribution for the BA-LTDR, MCD64A1 and MCD45A1 products fit the reference data (with a high correlation) but with an average estimation of 68%, 69% and 43%, respectively.In contrast, the BA GEOLAND-2 product did not fit the temporal profile of the reference and greatly overestimated the reference data (485%), as shown in Figure 2.These results are comparable to those obtained by other authors in the same region.Chang et al. [22] compared the L3JRC (the antecessor of GEOLAND) and the MCD45A1 global burned area products from 2000 to 2007 and showed similar results in the North American region in which L3JRC overestimated and MCD45A1 underestimated the burned area with respect to the reference data.Nuñez-Casillas et al. [20] made a comparative analysis of different BA products in the Canadian boreal forest in the year 2000 that showed an underestimation of the large burned areas in the MCD45A1 and BA-LTDR products but with a good fit for timing distribution and a slight overestimation of the GBA-2000 product (the antecessor of L3JRC).Giglio et al. 2010 [23] made a global comparison of different burned area products in the period 2001-2006; they showed that L3JRC significantly overestimates burned area in North America and the annual totals are uncorrelated with the reference data, while MCD45A1 underestimates burned area in this region but annual totals are highly correlated with the reference data.
An analysis of the positional accuracy of the estimations shows that the results of the BA-LTDR, MCD64A1 and MCD45A1 products are comparable in terms of determination coefficient of a linear regression model of the burned area percentage vs. the reference data and in terms of the commission errors.However, the BA-LTDR and MCD64A1 products showed a significant improvement over the MCD45A1 product alone regarding the omission errors.The analysis of the BA GEOLAND-2 product exhibited very low values for the determination coefficient and very high values for the commission errors.In the validation report of the BA GEOLAND-2 product [61], the authors obtained the same significance, for both commission and omission errors, for other climatic regions.
The analysis of the Pareto boundaries (Figure 4) shows that the most important contribution to commission and omission errors for the product BA-LTDR is due to the low resolution bias, although a slight improvement in the algorithm is still possible.For the other products with higher spatial resolution, the main contribution to both errors is due to the performance of the algorithms, which is really poor in GEOLAND-2.
The MCD45A1 and MCD64A1 products produced highly inaccurate results in the year 2001.The inaccuracy was due to the lack of input data available during the month of June because during an 18-day period (15 June to 2 July), the Terra MODIS sensor did not produce data; this failure occurred prior to the launch of the Aqua satellite.Therefore, the data from the MCD45A1 and MCD64A1 products for June 2001 and for May and July (at the end and beginning of these months, respectively) are unavailable.
The year 2006 presented the most erroneous records (followed by the year 2005) in terms of both accuracy and precision for the BA-LTDR product, despite being one of the best years for the MCD45A1 and MCD64A1 products.We investigated this occurrence to determine the possible causes.We analyzed the radiometric behavior of the different variables of the model and found that a major cause of the omission errors stemmed from the higher limit of the maximum BBFI threshold in the year before the fire event.This threshold is used to limit commission errors by avoiding the re-detection of pixels that were burned in previous years.This fact is observable in sub-scene 2006_1_2 (Figure 5).The BA-LTDR product does not detect even a single pixel inside the fire polygons in that scene.
Because the BBFI variable is computed by (1/ρ 2 (0.858 μm) + T 31 (11.03 μm)/2), it is necessary to determine the factors that influence this behavior.An analysis of the pixels inside the fire perimeters of sub-scene 2006_1_2 that were not identified as -burned‖ reveals that although the time evolution of the reflectance channels along the fire event year is correct, those for the thermal channels are incorrect.The time evolution does not fit the behavior of a burned pixel because the values of both thermal channels in the year preceding the fire event were substantially higher than normal and were even higher than those values for the fire event year, thereby preventing those potential burned pixels from being considered by the Bayesian network.We calculated the difference between the mean values of the maximum T 31 variable in the fire event year and the preceding year for a selected region of interest in sub-scene 2006_1_2 that included the fire polygons, obtaining a value of approximately −2 K. Upon further analysis, we observed that the maximum temperature was calculated in the 15th composite in the preceding year, and there was no continuity in the before-and-after composites for that year, indicating that it was not an existing BA but, rather, an outlier in the thermal channels of the LTDR dataset.
In sub-scene 2007_2, the BA-LTDR product did not detect a single pixel within the fire polygon of the scene.A similar result was noted in the MCD45A1 product (in this case, due to adverse atmospheric conditions at high latitudes).The BA GEOLAND-2 product was better able to detect that fire, although with a significant number of false pixels burned.
For the remaining sub-scenes, the best performance corresponded to the MCD64A1 and BA-LTDR products in both the estimated accuracy and spatial precision.The highest overestimation of the BA-LTDR product (with respect to the reference) obtained a value of 114.3% for scene 2010_1, whereas for the remainder of the cases, the estimate was lower or approximately 100%.The average estimate for the set of sub-scenes for the MCD45A1, MCD64A1, BA GEOLAND-2 and BA-LTDR products was 41%, 82%, 36% and 69%, respectively.However, in the case of the BA GEOLAND-2 product, a high percentage of the BA in the sub-scenes was due to the commission errors (Figure 5).

Conclusions
A BA-LTDR time series (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) for North America from the Terra-MODIS LTDR dataset was built successfully as a continuation of the previous time series, NOAA-AVHRR BA-LTDR (1984-1998), using the same methodology.The combination of the LTDR Version 3 dataset from Terra-MODIS and the modified Bayesian network algorithm, which used the MODIS Brightness Temperature from band 31 T 31 (11.03 μm), obtained better BA results than those from MCD45A1 and BA GEOLAND-2 in terms of spatial and temporal accuracy.However, the MCD64A1 product obtained the BA overall best results.
The BA GEOLAND-2 product overestimated the BA by approximately a factor of five with a coefficient of determination (R 2 ) from the regression model of the scatter plots of nearly zero in every year, whereas the BA-LTDR, MCD64A1 and MCD45A1 products underestimated the burned area by 32%, 31% and 57% with R 2 values of 0.75, 0.84 and 0.67, respectively.
In the balance between spatial accuracy and applicability to different years in the North American boreal forest region, the BA-LTDR could be considered as an alternative to the other global BA products, even those with a higher spatial/spectral resolution, such as MCD64A1 and MCD45A1, which exhibited certain deficiencies, particularly in the year 2001.Moreover, the BA-LTDR product has the advantage of providing existing data from 1981, making it the longest continuous BA time series from a remote sensing data record.available.We are grateful to Mike Flannigan and John Little for providing the fire perimeter dataset for Canada.We thank the agencies and services that processed and distributed the NASA and NOAA satellite data from which we obtained the majority of the images used in this study (MCD45A1, MCD64A1 and LTDR).SPOT 4-VEGETATION images were made available as part of the GEOLAND project.Constructive comments from anonymous reviewers were greatly appreciated.

Figure 1 .
Figure 1.The North American boreal forest study region (70°N, −168.5°E;45°N, −50°E), where boreal forest is represented by green color and all the other land covers by brown color.

Figure 2 .
Figure 2. Annual distribution of the burned area estimate (million ha) in the study region for the analyzed products and the reference data.

Figure 3 .
Figure 3. Graphical comparison of the burned area estimate accuracy by year (slope and R 2 ) for the North American boreal region for the Moderate Resolution Imaging Spectroradiometer (MODIS) burned area product (MCD45A1), MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1), burned area GEOLAND-2 product (BA GEOLAND-2), and burned area product from Long Term Data Record (BA-LTDR).The center represents a value of 0.0, and the external circle represents a value of 1.0.

Figure 4 .
Figure 4. Pareto boundaries (5 km and 1 km) of burned areas in the North America region, and commission and omission errors of the different burned area products for the period 2001-2011.

Figure 5 .
Figure 5.The 200 km × 200 km sub-scenes with large fires for the Moderate Resolution Imaging Spectroradiometer (MODIS) burned area product (MCD45A1), MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1), burned area GEOLAND-2 product (BA GEOLAND-2), and burned area product from Long Term Data Record (BA-LTDR).The red line represents the perimeter of all of the fires registered by the reference data on the sub-scene.Id MCD45A1 MCD64A1 GEOLAND-2 BA-LTDR The LTDR program is funded as part of the NASA Earth Science Research, Education, and Applications Solutions Network (REASoN).The goal of this program is to produce a consistent long-term dataset at a spatial resolution of 0.05° from the Advanced Very High Resolution Radiometer (AVHRR) and the MODIS instruments for use in global change and climate studies.For the period 1981-1999, the LTDR was based on Advanced Very High Resolution Radiometer [11][12][13]al Area Coverage (GAC) data.The LTDR from AVHRR records contains land surface reflectance from AVHRR bands 1 (0.6 μm) and 2 (0.85 μm); Top of Atmosphere (TOA) Brightness Temperature bands from AVHRR bands 3 (3.7 μm), 4 (10.8 μm) and 5 (12 μm); Solar Zenith, View Zenith and Relative Azimuth angle bands; and a Quality Assessment Field band, which includes improved atmospheric correction and inter-calibration between sensors[11][12][13].The continuation data records of the LTDR series for the years 2000 and later were created by processing the daily data acquired by MODIS onboard NASA's Terra and Aqua satellites.These data contain land surface reflectance from MODIS bands 1 (0.645 μm), 2 (0.858 μm) and 7 (2.13 μm); Brightness Temperature bands from MODIS bands 31 (11.03 μm) and 32 (12.02 μm); and Normalized Difference Vegetation Index (NDVI) and quality flag bands.The LTDR series was created by applying a Bidirectional Reflectance Distribution Function (BRDF) correction to the spectral subset of bands in the MODIS Surface Reflectance Climate Modeling Grid (CMG) level 3 product [14].

Table 1 .
Moderate Resolution Imaging Spectroradiometer (MODIS) burned area product (MCD45A1), MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1), burned area GEOLAND-2 product (BA GEOLAND-2) derived from VEGETATION sensor on board of SPOT satellites, and burned area product from Long Term Data Record (LTDR) derived from the Advanced very High Resolution Radiometer (AVHRR) on board of National Oceanic and Atmospheric Administration (NOAA) satellites and the MODIS sensor on board of Terra satellite (BA-LTDR).

Table 2 .
Burned area estimate accuracy by year (slope, R 2 and intercept) for the North American boreal region for the Moderate Resolution Imaging Spectroradiometer (MODIS) burned area product (MCD45A1), MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1), burned area GEOLAND-2 product (BA GEOLAND-2), and burned area product from Long Term Data Record (BA-LTDR).

Table 3 .
Commission and omission errors by year for the North American boreal region for the Moderate Resolution Imaging Spectroradiometer (MODIS) burned area product (MCD45A1), MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1), burned area GEOLAND-2 product (BA GEOLAND-2), and burned area product from Long Term Data Record (BA-LTDR).

Table 4 .
Detection of large areas of the North American boreal region for the Moderate Resolution Imaging Spectroradiometer (MODIS) burned area product (MCD45A1), MODIS Collection 5.1 Direct Broadcast Monthly Burned Area Product (MCD64A1), burned area GEOLAND-2 product (BA GEOLAND-2), and burned area product from Long Term Data Record (BA-LTDR).