MODIS Sensor Capability to Burned Area Mapping—Assessment of Performance and Improvements Provided by the Latest Standard Products in Boreal Regions

This paper presents an accuracy assessment of the main global scale Burned Area (BA) products, derived from daily images of the Moderate-Resolution Imaging Spectroradiometer (MODIS) Fire_CCI 5.1 and MCD64A1 C6, as well as the previous versions of both products (Fire_CCI 4.1 and MCD45A1 C5). The exercise was conducted on the boreal region of Alaska during the period 2000–2017. All the BA polygons registered by the Alaska Fire Service were used as reference data. Both new versions doubled the annual BA estimate compared to the previous versions (66% for Fire_CCI 5.1 versus 35% for v4.1, and 63% for MCD64A1 C6 versus 28% for C5), reducing the omission error (OE) by almost one half (39% versus 67% for Fire_CCI and 48% versus 74% for MCD) and slightly increasing the commission error (CE) (7.5% versus 7% for Fire_CCI and 18% versus 7% for MCD). The Fire_CCI 5.1 product (CE = 7.5%, OE = 39%) presented the best results in terms of positional accuracy with respect to MCD64A1 C6 (CE = 18%, OE = 48%). These results suggest that Fire_CCI 5.1 could be suitable for those users who employ BA standard products in geoinformatics analysis techniques for wildfire management, especially in Boreal regions. The Pareto boundary analysis, performed on an annual basis, showed that there is still a potential theoretical capacity to improve the MODIS sensor-based BA algorithms.


Introduction
Wildfires cause deforestation and habitat loss, and they are responsible for releasing a huge amount of aerosol particles and greenhouse gases into the atmosphere. These emissions vary depending on the Burned Area (BA) extension and on the type of biomass present in the region where the fire occurs. For example, Equatorial Asia, which is responsible for only 0.6% of the Global Burned Area (GBA), generates CO 2 and CH 4 emissions of 8% and 23%, respectively. Meanwhile, boreal forests, responsible for 2.5% of GBA, emit 9% of global CO 2 and 15% of CH 4 emissions [1]. Annually, between 5 and 15 million ha are burned in boreal forests, mainly in Siberia, Canada and Alaska, and the projections of different climate models estimate from decreases to increases in BA, which, as suggested by Kitzberger et al. [2], generates uncertainty in boreal regions, where global warming may create contrary effects.

•
In relation to spatial accuracy, to estimate the main metrics derived from the confusion matrix (commission and omission errors) and determine the Pareto Boundary (PB) for the native spatial resolutions of each product from the reference data to separate the errors of each product from the intrinsic errors associated with its spatial resolution.

•
To intercompare the spatiotemporal performance of the latest versions of the Fire_CCI 5.1 and MCD64A1 C6 products and to analyze any possible improvements over previous versions (i.e., Fire_CCI 4.1 and MCD45A1 C5.1).

•
To quantify the contribution of burned area fragmentation to the classified map errors, linking the area under the annual Pareto boundary curve with the total annual errors of each product to its spatial resolution.

Study Region
The study area spans a large section of Alaska, extending 10 • in latitude (60 • N-70 • N) and 27.5 • in longitude (Figure 1). This area is dominated by boreal forest, a complex set of plant communities modulated mainly by fire, soil type and drainage. The boreal forest forms a mosaic of hardwood-conifer mixed stands with closed canopy in well-drained areas, while in those with permafrost, open spruce stands predominate. Boreal forest makes up 90% of Alaska's forests, an area of approximately 42 million ha [29].
Sensors 2020, 20, x FOR PEER REVIEW 3 of 24 • To intercompare the spatiotemporal performance of the latest versions of the Fire_CCI 5.1 and MCD64A1 C6 products and to analyze any possible improvements over previous versions (i.e., Fire_CCI 4.1 and MCD45A1 C5.1).

•
To quantify the contribution of burned area fragmentation to the classified map errors, linking the area under the annual Pareto boundary curve with the total annual errors of each product to its spatial resolution.

Study Region
The study area spans a large section of Alaska, extending 10° in latitude (60° N-70° N) and 27.5° in longitude (Figure 1). This area is dominated by boreal forest, a complex set of plant communities modulated mainly by fire, soil type and drainage. The boreal forest forms a mosaic of hardwoodconifer mixed stands with closed canopy in well-drained areas, while in those with permafrost, open spruce stands predominate. Boreal forest makes up 90% of Alaska's forests, an area of approximately 42 million ha [29].

Reference Data
Polygons delimiting the area burned by wildfires in Alaska are available from the Alaska Fire Service (AFS, Fort Wainwright, AK, USA). AFS has compiled a very complete and accurate database since 1940. In addition to the geographic coordinates of the fire site and perimeter, AFS contains information such as the name of the fire; start and extinction dates; estimated BA; cause (naturally (e.g., lightning), human negligence or maliciously) or municipality of origin. Fire perimeters have always been delineated using the best available data source, from traditional hand-drawing on topographic maps in the early decades, to the interpretation of recent fine-scale satellite images with spatial resolutions less than 30 m [30]. This information was used as the reference data (the ground truth) for the accuracy assessment of the MODIS-derived BA products.

Reference Data
Polygons delimiting the area burned by wildfires in Alaska are available from the Alaska Fire Service (AFS, Fort Wainwright, AK, USA). AFS has compiled a very complete and accurate database since 1940. In addition to the geographic coordinates of the fire site and perimeter, AFS contains information such as the name of the fire; start and extinction dates; estimated BA; cause (naturally (e.g., lightning), human negligence or maliciously) or municipality of origin. Fire perimeters have always been delineated using the best available data source, from traditional hand-drawing on Sensors 2020, 20, 5423 4 of 23 topographic maps in the early decades, to the interpretation of recent fine-scale satellite images with spatial resolutions less than 30 m [30]. This information was used as the reference data (the ground truth) for the accuracy assessment of the MODIS-derived BA products.
For the study period, from 2000 to 2017, AFS recorded 1868 fires [31]. The total BA exceeded 11.6 million ha, with an annual average of 0.65 million ha, although strong year-on-year fluctuations were found (see Figure 2): in 12 of the 18 years analyzed, a total BA of 0.5 million ha was not exceeded, with 2001 and 2008 recording the lowest levels, 0.09 and 0.04 million ha, respectively. On the other hand, in 2004 and 2015, the total BA exceeded 2 million ha, with values of 2.71 and 2.08 million ha, respectively.
Sensors 2020, 20, x FOR PEER REVIEW  4 of 24 For the study period, from 2000 to 2017, AFS recorded 1868 fires [31]. The total BA exceeded 11.6 million ha, with an annual average of 0.65 million ha, although strong year-on-year fluctuations were found (see Figure 2): in 12 of the 18 years analyzed, a total BA of 0.5 million ha was not exceeded, with 2001 and 2008 recording the lowest levels, 0.09 and 0.04 million ha, respectively. On the other hand, in 2004 and 2015, the total BA exceeded 2 million ha, with values of 2.71 and 2.08 million ha, respectively. . Five categories were considered for the sizes (BA extension in ha) of the fires: very small (<100 ha), small (≥100 ha and <1000 ha), medium (≥1000 ha and <10,000 ha), large (≥10,000 ha and <100,000 ha) and very large (≥100,000 ha).
Throughout the period considered, small and very small fires account for an average of 60.50% of the total registered fires, but they only represented 1.90% of the total burned area ( Figure 3). In contrast, large and very large fires accounted for 13.60% of the total fires and 82.04% of the total burned area. Between them, the 16 fires of more than 100,000 ha burned 20.07% of the total burned area in the study period ( Figure 3). . Five categories were considered for the sizes (BA extension in ha) of the fires: very small (<100 ha), small (≥100 ha and <1000 ha), medium (≥1000 ha and <10,000 ha), large (≥10,000 ha and <100,000 ha) and very large (≥100,000 ha).
Throughout the period considered, small and very small fires account for an average of 60.50% of the total registered fires, but they only represented 1.90% of the total burned area ( Figure 3). In contrast, large and very large fires accounted for 13.60% of the total fires and 82.04% of the total burned area. Between them, the 16 fires of more than 100,000 ha burned 20.07% of the total burned area in the study period ( Figure 3).
To construct the annual reference data maps, all the perimeters of fires occurring during the MODIS era (2000-2017) were downloaded from AFS. The annual vector layers were projected to the Albers Conical Equal Area projection using the maximum area method to assign a burned/non-burned class label [32]. The final size of each pixel in the reference maps was 50 m × 50 m.  To construct the annual reference data maps, all the perimeters of fires occurring during the MODIS era (2000-2017) were downloaded from AFS. The annual vector layers were projected to the Albers Conical Equal Area projection using the maximum area method to assign a burned/nonburned class label [32]. The final size of each pixel in the reference maps was 50 m × 50 m.

Burned Area Products
This study compared the main global burned area products using MODIS sensor data: Fire_CCI 5.1 from the ESA project of the same name, led by the University of Alcalá de Henares, and the official NASA MODIS Direct Broadcast Monthly Burned Area Product, developed by the University of Maryland. Both of their most recent versions (Fire_CCI v. 5.1 and MCD64A1 C6), along with the previous versions (Fire_CCI v. 4.1 and MCD45A1 C5.1), were considered to analyze possible performance improvements. Table 1 shows their main characteristics. The main concerns of ESA and NASA when updating their products are to extend the time series and improve the algorithms to obtain the best validation results. Besides, NASA as owner and developer of MODIS sensors periodically reprocesses the entire data archive to incorporate better calibration and improved upstream data into all MODIS products. To construct the annual BA maps for Alaska, the respective monthly composites of the four products were downloaded. As with the reference maps, all these maps were re-projected to the Albers Conical Equal Area, with a pixel size of 50 m × 50 m, and finally combined on an annual basis. Figure 4 shows the annual maps of each of the products generated for the year 2009.

Burned Area Products
This study compared the main global burned area products using MODIS sensor data:  Table 1 shows their main characteristics. The main concerns of ESA and NASA when updating their products are to extend the time series and improve the algorithms to obtain the best validation results. Besides, NASA as owner and developer of MODIS sensors periodically reprocesses the entire data archive to incorporate better calibration and improved upstream data into all MODIS products. To construct the annual BA maps for Alaska, the respective monthly composites of the four products were downloaded. As with the reference maps, all these maps were re-projected to the Albers Conical Equal Area, with a pixel size of 50 m × 50 m, and finally combined on an annual basis. Figure 4 shows the annual maps of each of the products generated for the year 2009. Three of the four products base their algorithmic strategy on a hybrid approach. They use the MODIS active fires product (hotspots) in combination with changes in daily surface reflectance to identify burned pixels. In contrast, the fourth product (MCD45A1) only uses daily surface reflectance imagery. The following is a brief description of the algorithms used in each of the products.  Three of the four products base their algorithmic strategy on a hybrid approach. They use the MODIS active fires product (hotspots) in combination with changes in daily surface reflectance to identify burned pixels. In contrast, the fourth product (MCD45A1) only uses daily surface reflectance imagery. The following is a brief description of the algorithms used in each of the products.

MCD45A1 Collection 5.1
This BA product uses the MCD45 algorithm to identify burned pixels at a 500 m spatial resolution. It detects changes in the time series of daily bi-directional reflectance in bands 2 (0.841-0.876 µm) and 5 (1.23-1.25 µm) of MODIS sensor [37]. MCD45A1, which is currently deprecated, was the official product until the release of MCD64A1 C6.
The algorithm compares the observed daily reflectance values for each pixel with the predicted values using a bi-directional reflectance model in a 16-day-minimum time window that selects the candidate dates for burning in the forward and backward directions. If both dates match, the pixel is considered burned. That date is then used as a seed to identify, through a contextual iterative process, whether neighboring pixels can be classified as burned or not. In a final step, unselected candidate pixels are considered burned if they have at least three neighbors burned, with their burn date being the average of their neighbors. The algorithm finally excludes those pixels already burned in previous seasons and years [34][35][36]. This BA product uses the MCD45 algorithm to identify burned pixels at a 500 m spatial resolution. It detects changes in the time series of daily bi-directional reflectance in bands 2 (0.841-0.876 µm) and 5 (1.23-1.25 µm) of MODIS sensor [37]. MCD45A1, which is currently deprecated, was the official product until the release of MCD64A1 C6.
The algorithm compares the observed daily reflectance values for each pixel with the predicted values using a bi-directional reflectance model in a 16-day-minimum time window that selects the candidate dates for burning in the forward and backward directions. If both dates match, the pixel is considered burned. That date is then used as a seed to identify, through a contextual iterative process, whether neighboring pixels can be classified as burned or not. In a final step, unselected candidate pixels are considered burned if they have at least three neighbors burned, with their burn date being the average of their neighbors. The algorithm finally excludes those pixels already burned in previous seasons and years [34][35][36].

MCD64A1 Collection 6
This product applies the MCD64 algorithm to identify burned pixels using 500 m daily surface reflectance products for bands 5 (ρ5: 1.23-1.25 µm) and 7 (ρ7: 2.105-2.155 µm), along with daily active fire products (hotspots), both derived from the imaging products of the MODIS sensor on board the Terra and Aqua satellites [38].
The algorithm initially calculates the maximum daily changes in the time series of a burn-sensitive vegetation index VI = (ρ5 − ρ7)/(ρ5 + ρ7) for two pre-and post-date temporal windows. From these dates, it assigns the burn/unburn label to the pixel using the MODIS daily active fire product. It then extracts a set of training samples and performs an initial supervised classification of all pixels based on Sensors 2020, 20, 5423 7 of 23 the normalized distance measurement of each pixel to the nearest pixel in the training set. The final classification is obtained by using contextual information (nearest neighbors) [13].

Fire_CCI 4.1
The Fire_CCI 4.1 product, which has now been discontinued, covers only the period from 2005 to 2011. It identifies burned pixels at a 300 m spatial resolution using time series of daily surface reflectances from the MERIS sensor on board the ENVISAT satellite, and the active fires product derived from MODIS sensor images. The algorithm initially constructs monthly composites of surface reflectances by selecting the candidate pixels to be burned using the MODIS hotspot dates as criteria. It calculates cumulative distribution functions to discriminate the most clearly burned pixels by means of near-infrared reflectance thresholds. Then, it selects seed pixels in a 5 × 5 pixel window centered on the hotspot, to grow the burned regions by contextual analysis of the neighboring pixels. A final filter removes isolated pixels, both burned and non-burned [33].

Fire_CCI 5.1
This is the latest version of the burned area product from the Fire_CCI project. The algorithm, similar to its predecessor, is based on a two-stage hybrid approach. In the first stage, the seed pixels are selected, guided by the daily active fires (thermal anomalies) derived from the MODIS sensors on board the Terra and Aqua satellites. In the second stage, the growth and delimitation of the burned area is performed using the daily surface reflectance products from MODIS sensor bands 1 (0.62-0.67 µm) and 2 (0.841-0.876 µm) at a 250 m spatial resolution [12].

Accuracy Assessment
To validate the burned area products, an accuracy assessment was performed against the AFS reference data for those periods when each of the products was available (Table 1). The AFS perimeters used in this study are derived primarily from satellite images with spatial resolutions of less than 30 m, including the Landsat and Sentinel-2 missions, the latter since 2016. Figure 5 shows a description of the workflow of this assessment exercise carried out on an annual basis. All annual burned area map time series were delimited to the study region with 50 m × 50 m pixels in the same projection as the reference data. Each pixel therefore represents an area of 0.25 ha. Pixels were labeled 0 (Not Burned) or 1 (Burned). This simplified the calculation of the total annual burned area (in hectares) of each product, which was obtained by multiplying the number of pixels labeled as 1 (Burned) by 0.25.
To assess the temporal accuracy, the percentage of annual burned area for each product (P Year (%)) was calculated with respect to the reference set (Equation (1)). (1) where BAP Year (i) is the value of the pixel i (0: Not Burned; 1: Burned) for the indicated year of BA product; BAR Year (i) is the value of the same pixel i for the reference data; and n is the total number of pixels in each burned area map. For each annual burned area percentage distribution, a centrality measure (total value for the entire period), calculated as the weighted average of the annual percentages of burned area (P BAP (%)) (Equation (2)), was obtained. Y 0 and Y E are the first and last years of the BA product, respectively. Scatterplots of the annual percentages for each product were constructed against the reference data to analyze the linear dispersion around the central value and to identify the extreme cases (years with estimate percentages well below or above the mean).

of 23
Sensors 2020, 20, x FOR PEER REVIEW 8 of 24 To assess the temporal accuracy, the percentage of annual burned area for each product (PYear(%)) was calculated with respect to the reference set (Equation (1)).
where BAPYear(i) is the value of the pixel i (0: Not Burned; 1: Burned) for the indicated year of BA product; BARYear(i) is the value of the same pixel i for the reference data; and n is the total number of pixels in each burned area map. For each annual burned area percentage distribution, a centrality measure (total value for the entire period), calculated as the weighted average of the annual percentages of burned area (P (%)) (Equation (2)), was obtained. Y0 and YE are the first and last years of the BA product, respectively. Scatterplots of the annual percentages for each product were constructed against the reference data to analyze the linear dispersion around the central value and to identify the extreme cases (years with estimate percentages well below or above the mean).
Subsequently, a plot was constructed of the annual burned area temporal distribution of each product and that of the reference data. Correlation analysis of each time series was performed with the reference data, calculating the coefficient of determination R 2 (the square of Pearson's correlation coefficient).
The spatial accuracy assessment for each annual BA map employed the confusion matrix method, which is commonly used to validate thematic maps [39]. Table 2 shows the confusion matrix for a pixel-level thematic classification with two classes (Burned and Non-Burned). The independent reference information (AFS) is located in the matrix columns, and the burned area map data for each Subsequently, a plot was constructed of the annual burned area temporal distribution of each product and that of the reference data. Correlation analysis of each time series was performed with the reference data, calculating the coefficient of determination R 2 (the square of Pearson's correlation coefficient).
The spatial accuracy assessment for each annual BA map employed the confusion matrix method, which is commonly used to validate thematic maps [39]. Table 2 shows the confusion matrix for a pixel-level thematic classification with two classes (Burned and Non-Burned). The independent reference information (AFS) is located in the matrix columns, and the burned area map data for each product in the rows. The diagonal elements are the correctly classified data (true Burned and true Non-Burned). The other cells indicate commission errors (CE), i.e., pixels classified as Burned that are not actually burned (false burned), or omission errors (OE), pixels that are actually burned but that have been classified as Non-Burned (false non-burned) [40]. Other commonly used metrics that can be derived from the confusion matrix are Overall Accuracy (OA) (the percentage of correctly classified pixels), Sensibility (S) (or producer's accuracy) which calculates the rate of true burned pixels (the proportion of burned pixels that were correctly identified) and Specificity (Sp), or the rate of true Non-Burned pixels (the proportion of properly identified Non-Burned pixels) ( Table 2). The OA and Sp are metrics that can create a false sense of correctness in classified maps, especially in this study, due to the asymmetry between the two classes considered (the Non-Burned class is the majority and its success rate would be very high). On the other hand, S is related to the omission error (S = 1 − OE), so only commission and omission errors of the Burned class were considered. The two errors are not comparable since, although they both represent percentages of the pixels labeled as Burned, in one case they are Burned concerning the classified map (CE) and in the other case for the reference map (OE). Therefore, for the total error (TE) calculation, the weighted sum of both errors was considered, taking into account the percentage (P) of BA identified by the classified map, according to Equation (3).
The annual distribution of OE and CE for each burned area product was calculated, as well as their average values (the totals of all years considered). To identify extreme values, CE scatterplots were constructed against the annual OE of each product, and the annual deviations from the average values were analyzed. Likewise, to analyze the possible relationship between the spatial accuracy and the annual burned area, scatterplots of commission errors and omission errors were constructed for each BA product, against the annual burned reference area.
To determine whether fragmentation of the burned areas limits the spatial accuracy of the BA products, PB were constructed at 250, 300 and 500 m spatial resolutions [41]. To obtain each PB, the AFS reference map was used, with a high spatial resolution (50 m × 50 m) per pixel. Using a pixel spatial aggregation process, the reference data was resized to maps of 250, 300 and 500 m, where the value of each pixel contained the percentage of burned area. If one of these mixed pixels is classified as Burned, in a strict binary classification, a commission error equal to a 1-pixel value is being made. However, if it is classified as Non-Burned, an omission error equal to the pixel value is produced. The final decision on how it is classified depends on a p parameter selected in the range [0, 1], which sets the minimum threshold for assigning a mixed pixel to the Burned class. Then, for each degraded map, the CE and OE pairs, obtained by varying this parameter between 0 and 1 in equidistant steps, were calculated. The set of pairs {(CE i , OE i )} resulting from this process is the PB (Figure 6). The PB represents the lowest possible errors, obtained in a strict binary classification. These errors are attributable exclusively to the fragmentation of referenced burned areas that occur in mixed pixels when the spatial resolution of the data decreases. All points on the PB represent ideal classifications for a given spatial resolution. The PB further demonstrates that minimizing the commission and omission errors simultaneously is contradictory: a classification with minimum commission errors would mean greater failure of omission and, conversely, minimizing the OE would increase the CE.
represents the lowest possible errors, obtained in a strict binary classification. These errors are attributable exclusively to the fragmentation of referenced burned areas that occur in mixed pixels when the spatial resolution of the data decreases. All points on the PB represent ideal classifications for a given spatial resolution. The PB further demonstrates that minimizing the commission and omission errors simultaneously is contradictory: a classification with minimum commission errors would mean greater failure of omission and, conversely, minimizing the OE would increase the CE. To assess the impact of the PB on the accuracy of a BA classified map, the area enclosed by it (the Area Under the Pareto Boundary, AUPB) was calculated, by analogy with the area of the ROC (Receiver Operating Characteristic) curve used in other burned area studies using machine learning or data mining [42][43][44][45][46][47][48][49][50]. The ROC curve is a probability curve constructed from sensibility and 1specificity pairs {(Si,1-Spi)}, obtained using a procedure similar to that of the PB, and the area enclosed by it, the Area Under the ROC Curve (AUC), is interpreted as a measure of the separability between the two classes considered from the selected p parameter. Similarly, a larger area enclosed by the PB implies greater fragmentation of burned areas to that spatial resolution, and therefore a greater contribution to commission and omission errors by the mixed pixels.
The estimate of the area enclosed by the PB was calculated by applying the trapezoidal rule and by considering CE as the independent variable Equation (4) and OE as the dependent variable.
where OE(CE) is the curve representing the PB, defined from the point pairs {(CEi,OEi)}. Finally, to quantify the contribution of burned area fragmentation in the errors of the annual BA maps, linear regression models were constructed that related the annual areas under the PB curve to the total annual errors of each area burned product to its corresponding spatial resolution. Table 3 presents the annual BA percentages detected by the different products based on the AFS reference data. The years in which these products are available within the study period (2000-2017) are indicated. The total values (average percentage) for the entire period under consideration are also included. All the products underestimated the BA total for each of the years analyzed. Fire_CCI 5.1 presented the best overall results (65.89%), almost double that of version 4.1 (34.53%). Similarly, the To assess the impact of the PB on the accuracy of a BA classified map, the area enclosed by it (the Area Under the Pareto Boundary, AUPB) was calculated, by analogy with the area of the ROC (Receiver Operating Characteristic) curve used in other burned area studies using machine learning or data mining [42][43][44][45][46][47][48][49][50]. The ROC curve is a probability curve constructed from sensibility and 1-specificity pairs {(S i ,1-Sp i )}, obtained using a procedure similar to that of the PB, and the area enclosed by it, the Area Under the ROC Curve (AUC), is interpreted as a measure of the separability between the two classes considered from the selected p parameter. Similarly, a larger area enclosed by the PB implies greater fragmentation of burned areas to that spatial resolution, and therefore a greater contribution to commission and omission errors by the mixed pixels.

Temporal Accuracy
The estimate of the area enclosed by the PB was calculated by applying the trapezoidal rule and by considering CE as the independent variable Equation (4) and OE as the dependent variable.
where OE(CE) is the curve representing the PB, defined from the point pairs {(CE i ,OE i )}. Finally, to quantify the contribution of burned area fragmentation in the errors of the annual BA maps, linear regression models were constructed that related the annual areas under the PB curve to the total annual errors of each area burned product to its corresponding spatial resolution. Table 3 presents the annual BA percentages detected by the different products based on the AFS reference data. The years in which these products are available within the study period (2000-2017) are indicated. The total values (average percentage) for the entire period under consideration are also included. All the products underestimated the BA total for each of the years analyzed. Fire_CCI 5.1 presented the best overall results (65.89%), almost double that of version 4.1 (34.53%). Similarly, the MCD64A1 product showed estimates in the same order as Fire_CCI 5.1, although slightly lower (63.22%), an improvement which is more than double that of the previous product, MCD45A1 (28.11%).   Table 3) that there were years with low fire activity even though the estimate percentages in both products were higher. For example, in 2008, when 0.04 million ha were burned, MCD64A1 and Fire_CCI 5.1 gave above-average values for the estimate percentages of 67.49% and 79.43%, respectively. Figure 7 shows the scatterplots for the percentage of annual burned area detected by each of the four products analyzed against the annual burned area recorded by AFS (in millions of hectares). The right panel in Figure 7 represents the same percentages of annual burned area versus the percentage of large fires (above 10,000 ha) in the reference data. For MCD64A1, the highest percentages of burned area (83.36% and 88.14%, respectively) were detected in 2003   To analyze the temporal correlation of the reference data, the temporal distribution of the annual BA for the different products are represented along with the reference data ( Figure 8). All products conformed to the temporal pattern of the AFS reference data, with determination coefficients of 0.991 (Fire_CCI 5.1), 0.987 (MCD64A1), 0.973 (MCD45A1) and 0.817 (Fire_CCI4.1). To analyze the temporal correlation of the reference data, the temporal distribution of the annual BA for the different products are represented along with the reference data (  Table 4 shows the results of the main metrics derived from the confusion matrix (OE, CE and TE) for each year and for each BA product, together with the average for all years. Older versions of Fire_CCI and MCDs had the lowest average CE (5.9% and 6.6%, respectively). However, the lowest average OE (39.0% and 48.0%, respectively) are for the new versions, as are the lower average TE (43.9% and 59.3%, respectively).   Table 4 shows the results of the main metrics derived from the confusion matrix (OE, CE and TE) for each year and for each BA product, together with the average for all years. Older versions of Fire_CCI and MCDs had the lowest average CE (5.9% and 6.6%, respectively). However, the lowest average OE (39.0% and 48.0%, respectively) are for the new versions, as are the lower average TE (43.9% and 59.3%, respectively).   Figure 9 shows the scatterplots of annual omission errors versus annual commission errors for each of the products under consideration. Extreme values stand out, with a TE of around 100%, corresponding to the years 2000 and 2001. A qualitative analysis of these charts shows that, in general, new versions of both products have shifted the point cloud significantly down (less omission error) but have also shifted it slightly to the right (slightly increasing the commission error). The higher point cloud concentration of the Fire_CCI 5.1 product also indicates lower CE and OE variability. Figure 10 shows the annual commission and omission error scatterplots against the annual reference burned area for the latest versions of the MCD64A1 and Fire_CCI products. Both errors fluctuated strongly when the annual reference burned area was below 0.5 million ha, whereas they stabilized around their average for higher values. The fluctuations were higher for the MCD64A1 C6 product than for the Fire_CCI 5.1 product, in line with higher average values. Sensors 2020, 20, x FOR PEER REVIEW 15 of 24 Figure 9. Scatterplots of annual omission errors versus annual commission errors for each of the four burned area products analyzed. Figure 10 shows the annual commission and omission error scatterplots against the annual reference burned area for the latest versions of the MCD64A1 and Fire_CCI products. Both errors fluctuated strongly when the annual reference burned area was below 0.5 million ha, whereas they stabilized around their average for higher values. The fluctuations were higher for the MCD64A1 C6 product than for the Fire_CCI 5.1 product, in line with higher average values.    Figure 10 shows the annual commission and omission error scatterplots against the annual reference burned area for the latest versions of the MCD64A1 and Fire_CCI products. Both errors fluctuated strongly when the annual reference burned area was below 0.5 million ha, whereas they stabilized around their average for higher values. The fluctuations were higher for the MCD64A1 C6 product than for the Fire_CCI 5.1 product, in line with higher average values.

Pareto Boundaries
The Pareto boundaries were constructed at 250/300/500 m for all years in the study period. It can also be seen that the PB are more separated from the Cartesian axes in those years. This indicates a greater fragmentation of the burned areas when the spatial resolution is degraded, which contributes, in part, to the errors of the corresponding annual map.

Pareto Boundaries
The Pareto boundaries were constructed at 250/300/500 m for all years in the study period. Figure  11 shows the Pareto boundaries at 250 and 500 m for 2004, 2006, 2008 and 2015, as well as the annual CE and OE for the MCD64A1 and Fire_CCI 5.1 products. Years 2006 and 2008 were selected because both have the biggest errors in both products. In addition, two years, 2004 and 2015, with good behavior were included, to provide a visual comparison of the PB. One can observe that in 2006 and 2008 the distance from the annual pair of errors (EC, EO) to the PB is much larger than in 2004 and 2015. It can also be seen that the PB are more separated from the Cartesian axes in those years. This indicates a greater fragmentation of the burned areas when the spatial resolution is degraded, which contributes, in part, to the errors of the corresponding annual map. To quantify the impact of the PB on the annual maps of the different BA products, the annual and total AUPB values were calculated at spatial resolutions of 250/300/500 m (Table 5)  To quantify the impact of the PB on the annual maps of the different BA products, the annual and total AUPB values were calculated at spatial resolutions of 250/300/500 m (Table 5) Finally, Figure 12 shows the annual TE scatterplots for each product versus the annual PB area corresponding to their spatial resolution. A linear regression model was constructed for each BA product. The R 2 values are very low for the older versions of the BA products (0.00 and 0.18, respectively), but increase markedly in the new versions. Fire_CCI 5.1 has a slightly higher value than MCD64A1 C6 (0.45 vs. 0.41) and also a linear regression line slope that is almost double (63.1 vs. 32.8).  Finally, Figure 12 shows the annual TE scatterplots for each product versus the annual PB area corresponding to their spatial resolution. A linear regression model was constructed for each BA product. The R 2 values are very low for the older versions of the BA products (0.00 and 0.18, respectively), but increase markedly in the new versions. Fire_CCI 5.1 has a slightly higher value than MCD64A1 C6 (0.45 vs. 0.41) and also a linear regression line slope that is almost double (63.1 vs. 32.8).

Discussion
The evaluation and validation of global BA products derived from satellite imagery require a set of reliable, independent and representative reference fire perimeters that cover as long a period as possible. It is necessary to understand the uncertainty of these products before incorporating them as input data into global carbon, vegetation or climate models, as well as for the management of all the wildfire phases. Numerous prior studies have evaluated the behavior of the MODIS sensor-derived products analyzed in this study, both as part of ESA's Fire Climate Change Initiative Project and in the different versions [13,[51][52][53][54][55][56][57][58] of the MODIS Direct Broadcast Monthly Burned Area Product. However, many of these works construct reference fire perimeters using images from better spatial resolution sensors such as Landsat TM/ETM or Sentinel-2, and which are limited to short time periods, usually one or several years, due to the difficulty in creating larger reference sets. This work presents an independent intercomparison exercise on the accuracy of the two main global BA products derived from the MODIS sensor covering an extensive boreal region over 18 years, containing all the AFS-registered fire perimeters. To the best of our knowledge, there have been no previous intercomparison studies of these four products over such a long period.
In relation to the temporal accuracy (Figure 8), and except for the Fire_CCI 4.1 product, which has the shortest time series, the BA products analyzed conform to the time pattern of the reference data with determination coefficients above 0.97, and with results similar to those found by other authors in different ecosystems. Thus, for example, Turco et al. [56] reported high determination coefficients of 0.96 and 0.97 in the monthly BA estimates for Fire_CCI 5.1 and MCD64A1 C6 when compared to the reference dataset from the European Forest Fire Information System (EFFIS) that includes burned area data for some European countries in the Mediterranean basin (Portugal, Spain, Southern France and Greece). However, all the BA products analyzed underestimated the annual burned area (Table 3) with higher percentages than those obtained by other authors for regions other than Alaska's boreal forest. Turco et al. [56] obtained an underestimation of only 14% for the Fire_CCI 5.1 product with the EFFIS reference set. Campagnolo et al. [57] found a 28% underestimation for the MCD64A1 product using a reference set of more than 100 fire perimeters (1.24 Mkm 2 ) derived from Landsat TM/ETM+ images distributed around the world for 2008 and previously constructed by Padilla et al. [54].
From the results intercomparison determined for the two new versions (Table 3), one can show that MCD64A1 provided better BA estimates than Fire_CCI in some years, despite having a lower spatial resolution (500 vs. 250 m), even in the years with the largest BA, which are the ones that contributed more to the average values (2004, 2005 and 2015); nevertheless, it performed relatively worse in other years (2009, 2010 and 2011). For the MCD64A1 and Fire_CCI 5.1 products, in the years with the highest amount of BA (over 1 million ha), the annual estimate percentages tended to stabilize around the average value with lower variability (Figure 7). One should bear in mind that it is precisely these years that contribute most to the overall average throughout the study period. Conversely, the years with the least amount of BA (below 0.5 million ha) had greater variability; this was less pronounced in the Fire_CCI 5.1 product than in MCD64A1, probably due to its better spatial resolution (250 vs. 500 m).
With regard to spatial accuracy (Table 4), a global analysis showed that the older versions, Fire_CCI 4.1 and MCD45A1 C5.1, had the lowest commission errors (5.9% and 6.6%, respectively) but the omission errors were very high (67.5% and 73.7%, respectively). The new versions, Fire_CCI 5.1 and MCD64A1 C6, had significantly reduced OE (up to 39.0% and 48.0% respectively), even at the cost of slightly worse CE (7.5% and 17.8%, respectively). As the percentage of BA of all products was below 100% (they underestimated the total BA), the weighted sum or TE tended to give a lower weight to the CE than to the OE, favoring the algorithms that reduce the omission error even at the cost of increasing the commission errors. This low imbalance between omission and commission errors (along with the sharp drop in OE) translates into a better BA estimate. These results are comparable to those obtained by other authors using other reference sets [51,58]. In a recent validation of The OE from this latter study, which was specific to the boreal region, differs significantly from that obtained in the present work (48.0%) for the Alaska region; we understand that this may be due to the small reference dataset used by Boschetti et al. [51] for the entire boreal region. In contrast, our study used all the burned area perimeters recorded in the 2000-2017 period.
The PB analysis ( Figure 11) allowed us to partially explain the annual variability in the commission/omission errors of the BA products. The average AUPB values (Table 5) reflect an increase as the spatial resolution worsens, indicating that there is a higher percentage of commission and omission errors attributable to the data's low spatial resolution. Again, there was high variability in the annual values relative to the average values, attributable to the fragmentation of the burned areas. This high accuracy variability attributable to landscape fragmentation was also reported by Rodrigues et al. [59], in the accuracy assessment of MCD64A1 C6 in the Brazilian Cerrado vs. Landsat perimeters for the 2011-2019 period. Rodrigues et al. [59] found that, in the northern Cerrado, which had larger areas affected by fire, the MCD64A1 C6 performance was significantly higher than in the southern Cerrado area, where a more fragmented landscape and smaller patches of fire predominated.
The linear regression analysis between the TE and the AUPB found an upward linear trend for the new versions of the BA products; this partially justifies the high commission and omission errors encountered in some years. However, it is also apparent that in some years (e.g., 2016), the high level of fire fragmentation had little influence on the TE. Conversely, in years with low levels of fragmentation (e.g., 2002 and 2005 for Fire_CCI 5.1 and 2003 and 2014 for MCD64A1 C6), the TEs were significant. For these years, the observations were limited by other factors such as the poor behavior of the detection algorithms, the low severity of burned areas or certain environmental conditions; as indicated by Loboda et al. [60], these factors might explain such errors. The more precise linear fit of the new versions compared to the old ones, with R 2 values close to 0.5 and a higher slope, reflects greater TE sensitivity (CE and OE) to burned area fragmentation.
Of all the years analyzed, and except for the years 2000 and 2001 (years with incomplete data from the MODIS sensor), there are two years in which the latest versions of MCD64A1 and Fire_CCI performed poorly. The year 2006 showed low detection rates (29.6% and 45.83%, respectively), high omission values (83.9% and 60.9%) and high commission values (29.6% and 14.7%). The burned area recorded (0.11 million ha) was below the annual average, as was the same percentage corresponding to large fires (69.82%), which was also below the annual average, while the PB areas at 250 (0.002257) and 500 m (0.006480) were almost double the average values at these resolutions. Similarly, in 2008, the year with the least amount of burned area (0.04 million ha) and the lowest percentage attributable to large fires (36.65%), there were high commission and omission error values (57.3% and 71.2% for MCD64A1 and 29.2% and 43.8% for Fire_CCI) but a large percentage of burned area detected (67.49% and 79.43%, respectively). For that year, the area enclosed by the PB at 250/500 m was the highest, so a significant part of the spatial errors was due to the data's low spatial resolution.

Conclusions
An independent and detailed study was carried out to evaluate the spatiotemporal accuracy of the latest versions (along with the previous versions) of the two main global scale BA products derived from MODIS-sensor satellite images, here restricted to the Alaskan boreal region. As reference data, we used all the polygons of the BA recorded by AFS over the study period. In addition, a detailed 50 m × 50 m pixel analysis of the accuracy of each BA product was performed. Fire_CCI 5.1 and MCD64A1 C6 presented significant improvements over their previous versions in terms of BA estimation. Improvements were achieved by reducing the imbalance between commission and omission errors and, especially, by greatly reducing omission errors even at the expense of worsening commission errors. Both products currently produce similar BA estimation percentages, although the positional accuracy of Fire_CCI is better than that of MCD64A1, which is in line with its higher spatial resolution (250 vs. 500 m). Fire_CCI 5.1 would be the option chosen for users who, through geoinformatics analysis techniques, use BA products for forest fire management. For those users who are studying the increase in greenhouse gas concentrations or the change in the chemical composition of the atmosphere due to fire emissions, any of the latest versions of the products analyzed could be suitable, but always taking into account the errors of omission (over 40%) in both cases.
The high variability in annual results is noteworthy, both in the percentages of BA detection and in the omission and commission errors, which puts into question much accuracy assessment work that uses limited spatiotemporal reference data. It would be recommendable to use all possible reference data when available.
The yearly analysis of burned-area fragmentation across the corresponding Pareto boundaries has allowed us to establish a quantitative measure of the same (AUPB), which relates to the total errors (weighted sum of commission and omission errors) obtained for the latest versions of the Fire_CCI and MCD64A1 products. The results from this work could be extrapolated to other boreal regions that do not have such accurate reference datasets as are available in the northern regions of North America.