Evaluation of a Bayesian Algorithm to Detect Burned Areas in the Canary Islands ’ Dry Woodlands and Forests Ecoregion Using MODIS Data

Burned Area (BA) is deemed as a primary variable to understand the Earth’s climate system. Satellite remote sensing data have allowed for the development of various burned area detection algorithms that have been globally applied to and assessed in diverse ecosystems, ranging from tropical to boreal. In this paper, we present a Bayesian algorithm (BY-MODIS) that detects burned areas in a time series of Moderate Resolution Imaging Spectroradiometer (MODIS) images from 2002 to 2012 of the Canary Islands’ dry woodlands and forests ecoregion (Spain). Based on daily image products MODIS, MOD09GQ (250 m), and MOD11A1 (1 km), the surface spectral reflectance and the land surface temperature, respectively, 10 day composites were built using the maximum temperature criterion. Variables used in BY-MODIS were the Global Environment Monitoring Index (GEMI) and Burn Boreal Forest Index (BBFI), alongside the NIR spectral band, all of which refer to the previous year and the year the fire took place in. Reference polygons for the 14 fires exceeding 100 hectares and identified within the period under analysis were developed using both post-fire LANDSAT images and official information from the forest fires national database by the Ministry of Agriculture and Fisheries, Food and Environment of Spain (MAPAMA). The results obtained by BY-MODIS can be compared to those by official burned area products, MCD45A1 and MCD64A1. Despite that the best overall results correspond to MCD64A1, BY-MODIS proved to be an alternative for burned area mapping in the Canary Islands, a region with a great topographic complexity and diverse types of ecosystems. The total burned area detected by the BY-MODIS classifier was 64.9% of the MAPAMA reference data, and 78.6% according to data obtained from the LANDSAT images, with the lowest average commission error (11%) out of the three products and a correlation (R2) of 0.82. The Bayesian algorithm—originally developed to detect burned areas in North American boreal forests using AVHRR archival data Long-Term Data Record—can be successfully applied to a lower latitude forest ecosystem totally different from the boreal ecosystem and using daily time series of satellite images from MODIS with a 250 m spatial resolution, as long as a set of training areas adequately characterising the dynamics of the forest canopy affected by the fire is defined.


Introduction
As a recurring natural process, fire has been a basic element in preserving the health and biodiversity of numerous terrestrial ecosystems for thousands of years [1][2][3].However, any alteration to this nature-driven dynamic, due either to direct anthropogenic actions or as a result of global warming, can give rise to situations where fire can become a significant disturbance for many ecosystems, causing feedbacks on global warming, regionally and globally, in turn.Prominent amongst said disturbances are the loss of biodiversity due to the decimation of species and biomass [2]; soil deterioration and erosion which, alongside the change in land cover (terrestrial albedo) [4,5], alter water and energy flows; and the emission of large amounts of aerosols and greenhouse gases affecting the carbon cycle and global warming [6][7][8][9].For all the above, several international climate research teams-chiefly members of the Global Climate Observing System (GCOS)-have identified Fire Disturbance as an Essential Climate Variable (ECV) that critically contributes to the characterization of Earth's climate.
The Fire Disturbance ECV consists of three variables: Burned Area, Active Fires and Fire Radiated Power [10].Amongst them, Burned Area is deemed the main variable [11] and, therefore, its tracking becomes imperative, in order to forecast and understand the Earth's climate evolution, as well as to take the necessary steps to mitigate its negative effects.But Burned Area maps are also necessary in other applications of vital interest to our society as an early warning fire alert system; the fire hazard assessment for ecosystems management; the study of changes in atmospheric chemistry; and in biogeochemical modelling and that of the vegetation dynamics at a global scale [12].
It is not easy to gather time-consistent Burned Area data for periods exceeding 30 years for climate-related studies.Only sensors onboard the Earth observation satellites-primarily polar orbit satellites, which boast better spatial resolution (≤1000 m) than the geostationary orbit ones (>1000 m)-could provide images of areas affected by the fire at a local, regional, and global level, allowing for the registry of their distribution and occurrence patterns [13].
Of all remote sensors, the only one that has been taking daily images of our planet from the early 1980s is Advanced Very High Resolution Radiometer (AVHRR) onboard the series of NOAA (National Oceanic and Atmospheric Administration) satellites, whose data have made it possible to develop various algorithms to assess fire effects [14][15][16][17][18][19][20][21][22][23].Other widely used sensors in Burned Area mapping at a regional and global level have been VEGETATION (VGT) onboard SPOT (Systeme Pour I'Observation de la Terre) 4 and 5 satellites [24][25][26][27], and Moderate Resolution Imaging Spectrometer (MODIS) onboard NASA's Terra and Aqua satellites.MODIS, with excellent temporal resolution (12 h) owing to the combination of two satellites, spectral resolution (36 bands) and spatial resolution (250-1000 m), has allowed for the development of various Burned Area detection algorithms that have been globally applied to and assessed in diverse ecosystems, ranging from tropical to boreal [28][29][30][31][32][33][34].
There are several global Burned Area products obtained from the data provided by remote sensors [12,44].These products are free to access, but not all them are operational or cover long temporal series despite the sensors used having had or having a longer lifespan.This is the case of GBS (Global Burned Surfaces) [45], which is limited to the 1982-1999 period.GBS used Pathfinder AVHRR Land (PAL) data at an 8 km spatial resolution.Based on another sensor, SPOT-VGT, we find three Burned Area products: GBA2000 [46], specifically designed for the year 2000; L3JRC [47], available for the April 2000-March 2007 period; and GEOLAND2 [48], based on the previous two and with data available from April 1999 to March 2014.GEOLAND2 has currently been adapted to sensor PROBA-V, with a 300 m spatial resolution.PROBA-V's spectral bands are similar to those of the VGT, thereby preserving an observational consistency with its predecessor.Two other Burned Area products from European satellites were GLOBCARBON, obtained through a combination of the VGT sensor and the ATSR instrument (Along Track Scanning Radiometer) [49] for the April 1998-December 2007 period, and GLOBSCAR, from ATSR-2 images [50].
Two of the most reliable global Burned Area products, MCD45A1 [51] and MCD64A1 [52], have been generated from MODIS Aqua and Terra images, with maps made available from the year 2000 onwards.Both are monthly Level 3 gridded 500 m products containing per-pixel burning and quality information, and tile-level metadata.Another product that has been using MODIS data from the year 2000 is Global Fire Emissions Database (GFED) [31], which estimates Burned Area and monthly emissions of various fire types with a 0.25 • latitude by 0.25 • longitude spatial resolution.GFED have available data from 1997 through to the present, for which reason it is regarded as the longest Burned Area temporal series [31].Active fire data from the Tropical Rainfall Measuring Mission (TRMM) Visible and Infrared Scanner (VIRS), and the Along-Track Scanning Radiometer (ATSR) were used to get GFED where MODIS data was not available yet [53].
Finally, it is worth mentioning the European Project Fire_cci's Burned Area dataset [12].Fire_cci Burned Area has been obtained by combining temporal changes in reflectance from Medium Resolution Imaging Spectrometer (MERIS) onboard the ESA ENVISAT satellite and thermal information on active fires from MODIS, for the years 2005 to 2011, as monthly composites.The main advantage of the Burned Area Fire_cci product could be the improvement in detection of small fires by increasing the spatial resolution to 300 m with respect to MODIS Burned Area products [43].In a second stage, Fire_cci will cover the time series 2000-2017, with a spatial resolution between 250 and 500 m, depending on the sensor being used [54].
Studies comparing the accuracy of global Burned Area products, discussed in the previous paragraphs, prove that there are appreciable differences between them in many regions of our planet [12,23,31,33,43,55,56].This way, using these products in applications at a local or regional level calls for previously validated studies in order to know the exact commission and omission errors and decide whether or not they may be used as input data for climate models and vegetation dynamics in line with the current requirements.It is hard to connect globally developed algorithms for a specific sensor with the various ecosystems affected by fires at different latitudes.
An algorithm for mapping Burned Area in boreal forests was developed as part of a researching project supported in part by the Spanish Ministry of Science and Innovation [23].It used AVHRR archival data Long Term Data Record (LTDR) version 3, with 0.05 • of spatial resolution and was applied to western Canada for the period between 1984 and 1999.Fire event records from Canadian Forest Service National Fire Database (CFSNFD) were used to train the algorithm.The kernel of the algorithm is a Bayesian network classifier, combined with absolute and relative radiometric thresholds, and a post-processing step (neighbourhood analysis) for spatial fire coherence.A chain of software tools was built in order to implement the algorithm and automate the processing from the daily images in their original format to 10 day Burned Area maps.The algorithm was successfully evaluated against reference data and with Burned Area records generated from 1 km AVHRR data by [22].The same Bayesian network algorithm was applied to mapping Burned Area in the North American boreal region for the years 1984-1998 using the AVHRR-LTDR dataset version 3, and two training sites (Western Canada and Alaska).The results showed an overall good agreement compared to reference maps derived from CFSNFD and Alaska Fire Service (AFS), and comparable results for both training sites [57].Later, a Burned Area time series (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) for North America from the MODIS-LTDR dataset version 3 was successfully built as a continuation of the previous time series applying the same methodology, but using the MODIS Brightness Temperature from band 31 T31 (11.03 µm) instead of the AVHRR Top of Atmosphere Brightness Temperature from band 3 (3.7 µm).The results obtained were even better than those derived from other Burned Area products, with higher spatial resolution (MCD45A1 and GEOLAND2) in terms of spatial and temporal accuracy, and comparable with the MCD64A1 Burned Area product [58].
This paper presents the application of the above methodology [23,57,58] to a lower latitude forest ecosystem, totally different from the boreal ecosystem and using daily time series of satellite images from MODIS with a 250 m spatial resolution.The main objectives are: (i) to build a 10 day composite time series of MODIS satellite images using the maximum surface temperature criterion, for the 2002-2012 period in a Mediterranean forest region (the Canary Islands); (ii) to generate Burned Area maps corresponding to the study region using the methodology based on the Bayesian classifier and the same variables, although modifying the values of its parameters by training with reference data extracted from official fire records and post-fire Landsat images; and (iii) to assess the spatial and temporal accuracy of the detected Burned Area with reference data, and to compare the results with those of the official MODIS Burned Area products with analogous spatial resolution (MCD45A1 and MCD64A1).

Study Area
The study area included the five western islands of the Canarian Archipelago (Spain) (29 • N, 18 • W; 27.50 • N, 15 • W): Tenerife, La Palma, La Gomera, El Hierro, and Gran Canaria (Figure 1).These volcanic islands have a particularly rich biodiversity that supports a wide variety of endemic taxa produced by evolutionary processes due to isolation [59].The five islands form a small terrestrial ecoregion [60] known as the "Canary Islands' dry woodlands and forests" [61].Despite its name, in the northern slope of these islands, and favoured by the north-eastern humid wind regime, the trade winds, there is a humid forest known as "Laurisilva" (laurel forest).This forest is comprised of tree species of the laurel family endemic to subtropical climates, which is regarded as one of the most interesting ecosystems in the Canary Islands, and an heirloom of Tertiary forests in the Mediterranean basin.
corresponding to the study region using the methodology based on the Bayesian classifier and the same variables, although modifying the values of its parameters by training with reference data extracted from official fire records and post-fire Landsat images; and (iii) to assess the spatial and temporal accuracy of the detected Burned Area with reference data, and to compare the results with those of the official MODIS Burned Area products with analogous spatial resolution (MCD45A1 and MCD64A1).

Study Area
The study area included the five western islands of the Canarian Archipelago (Spain) (29°N, 18°W; 27.50°N, 15°W): Tenerife, La Palma, La Gomera, El Hierro, and Gran Canaria (Figure 1).These volcanic islands have a particularly rich biodiversity that supports a wide variety of endemic taxa produced by evolutionary processes due to isolation [59].The five islands form a small terrestrial ecoregion [60] known as the "Canary Islands' dry woodlands and forests" [61].Despite its name, in the northern slope of these islands, and favoured by the north-eastern humid wind regime, the trade winds, there is a humid forest known as "Laurisilva" (laurel forest).This forest is comprised of tree species of the laurel family endemic to subtropical climates, which is regarded as one of the most interesting ecosystems in the Canary Islands, and an heirloom of Tertiary forests in the Mediterranean basin.
In the western Canary Islands, as in other zones with Mediterranean climate, forest fires are a constant worry.Summers with very low precipitation and high daytime temperatures reduce humidity and promote the spread of fires resulting from both natural and man-made causes.Canarian endemic pine forests (Pinus canariensis), located in the northern side above 1200 and up to 2100 m of altitude, and between 700 and up to 2300 m in the southern side, are the plant association most affected by those fires.However, fires also affect the endemic Macaronesian heaths (fayal-brezal) and "Laurisilva".Only the pines and heaths have a remarkable ability to regenerate after the fire [62][63][64].As for the rest, being comprised of species not adapted to fire, it could take decades to regain its natural balance.The image is a NDVI composition from the Terra-MODIS sensor at a 500-m spatial resolution on 30 December 2015 retrieved from the online Land, Atmosphere Near Real-time Capability for the EOS (LANCE) system (https://lance-modis.eosdis.nasa.gov/imagery/subsets/?area=af).The red rectangle includes five islands studied.
In the western Canary Islands, as in other zones with Mediterranean climate, forest fires are a constant worry.Summers with very low precipitation and high daytime temperatures reduce humidity and promote the spread of fires resulting from both natural and man-made causes.Canarian endemic pine forests (Pinus canariensis), located in the northern side above 1200 and up to 2100 m of altitude, and between 700 and up to 2300 m in the southern side, are the plant association most affected by those fires.However, fires also affect the endemic Macaronesian heaths (fayal-brezal) and "Laurisilva".Only the pines and heaths have a remarkable ability to regenerate after the fire [62][63][64].As for the rest, being comprised of species not adapted to fire, it could take decades to regain its natural balance.

Burned Area Information from MAPAMA
Given the complex topography of the Canary Islands, with its many ravines and altitudes ranging between sea level up to 3700 m in the island of Tenerife, combined with weather conditions occurring in summer, the fires are very difficult to control, have a great destructive power, and, therefore, have the potential to burn large areas.According to data from the Ministry of Agriculture and Fisheries, Food and Environment of Spain (MAPAMA) [65], large forest fires (LFF) affecting areas over 100 ha have been increasingly frequent in the region over the last two decades.We must emphasize that during the period studied, between 2002 and 2012, there have been only 14 fires of over 100 ha (Table 1), representing less than 1% of the 1319 occurred.Nevertheless, they have burned 97.35% (55,571.20 ha) (Table 2) of the total affected forest area (57,085.04ha) in this period.In addition, there have been several fires in the same year, and at the same time, but on different islands, as seen in 2007 and 2012.The year 2007 holds the annual affected area historical record, with 35,758.62ha in Tenerife and Gran Canaria.On the other hand, in 2012, although the area affected was less than half (12,135.88 ha), there was a notable difference from previous episodes.This year, four LFF burned between 752 and 6512 ha in three of the western islands simultaneously, two in La Palma, one in Tenerife and another one in La Gomera, where the Garajonay National Park-declared a UNESCO World Heritage Site owing to an ecosystem with as high a natural value as Laurisilva-was burning for 86 days (Table 1).

Data and Pre-Processing
Two daily image products, MOD09GQ and MOD11A1, derived from the Terra-MODIS sensor, were downloaded from http://modis.gsfc.nasa.govfor the study area and the 2002-2012 period.The MOD09GQ provides the surface spectral reflectance corrected for atmospheric conditions of band 1/RED (sur_refl_b01_1) and band 2/NIR (sur_refl_b02_1), with a spatial resolution of 250 m, and from MOD11A1 original files land surface temperatures (LST_Day_1km), were extracted.This information was used for the construction of daily files in geoTiff format with a band structure similar to that of the LTDR files [66] (Table 3).The sinusoidal projection of the original data was transformed to a UTM-28N projection with the nearest neighbouring resampling, an output pixel size of 0.0025 • (250 m) and WGS-84 datum.Daily reflectances and temperature bands were combined with their respective quality indicators (QC) to obtain a quality indicator layer according to the LTDR criterion [66].Scaling factors were applied to each of them to obtain the daily files in band sequential format (BSQ) of signed 16 bit integer.Finally, composite images of 10 days were generated with the criterion of maximum surface temperature [67,68] considering the QC for filtering out invalid values in measurements.The Modis Reprojection Tool retrieved from https://lpdaac.usgs.gov/tools/modis_reprojection_tooland IDL ® , a product of Exelis Visual Information Solutions, Inc., (Boulder, Colorado, USA) were used in the preprocessing.The first cloud-free post-fire images collected by the series of Landsat sensors with a spatial resolution of 30 m (Table 4), together with the Burned Area registered by MAPAMA [65], were both used as reference data.MAPAMA Burned Area maps are usually made with GPS by mean of fieldwork and, for the greatest fires, using aerial orthophotos with 50 cm resolution.The polygons of all the fires greater than 100 ha that occurred from 2002 to 2012 were used to produce the annual vector layers of ground truth.Next, these vector layers were reprojected to a UTM-28N projection with a 250 m × 250 m pixel size to generate annual ground truth maps.In order to determine how the pixel is assigned a burned/non-burned value, the maximum area method within the pixel was used [69].In order to compare the results obtained by the methodology proposed in this paper (Section 3.2) with the official Burned Area products derived from the MODIS sensor, the MCD64A1 Version 6 [70] and MCD45A1 Version 5.1 monthly subset were downloaded for the 2002-2012 period.These monthly maps were reprojected to UTM-28N projection with a 500 m × 500 m pixel size and resized to the study area, and finally combined to generate the annual maps of burned/non-burned pixels.MODIS Burned Area Products used are briefly described below.

• MCD64A1 Version 6
This product is generated using an improved version of the MCD64 Burned Area mapping algorithm [52].MCD64A1 is derived from a group of daily active fires with a 1 km spatial resolution and daily land surface reflectances of 1, 5, and 7 spectral bands with 500 m spatial resolution.The algorithm uses a ratio spectral index of channels 5 and 7 [71]) as a primary indicator.This Vegetation Index (VI) detects an immediate decrease in band 5 reflectance after a burn, in comparison with minor changes in band 7, being a good discriminator between burned and non-burned pixels.The synergy algorithm applies dynamic thresholds to composite imagery generated by VI and a measure of temporal texture.Cumulative active fire maps are used to guide the selection of burned and non-burned training samples and to guide the specification of prior probabilities [31].
• MCD45A1 Version 5.1 MCD45A1 has a spatial resolution of 500 m and a time resolution of one month.Unlike MCD64A1, do not use the active fire observations as input data in its algorithm.MCD45A1 employs a change detection algorithm based on a BRDF (bi-directional reflectance distribution function) model [72].The algorithm calculates the difference between the observed BRDF and the predicted BRFD in the near and middle infrared channels at the viewing and illuminating angles of the observation.The decision about a pixel being burned is based on a threshold.This method maps the approximate date of burning by locating pixels showing rapid change in daily surface reflectance time series data over an entire year.

Bayesian Algorithm to Detect Burned Area
The methodology proposed by Moreno-Ruiz et al. [23] based on a Bayesian algorithm, was adapted and applied to the study region to obtain annual maps of burned areas and their temporal distribution.This methodology is formed by a series of stages displayed in Figure 2. In order to carry them out, a software tool allowing for effective implementation and its integration into a process chain has been developed.The stages are as described below: 1.
Pre-processing of MOD09GQ and MOD11A1 time series used for building the daily set of images with the same format as the LTDR dataset (see previous section).

2.
Classification scheme: In order to sense a Burned Area at the pixel level, a dichotomous classification scheme was used, the thematic classes under consideration being Burned and Non-Burned.

3.
Building a set of training cases: The algorithm was trained for the Burned class with 152 supervised burned pixels from a fire of 6512 ha occurred in the south of the island of Tenerife in the year 2012.For the Non-Burned class, we used 450 non-burned pixels from nearby regions.4.
A choice of an optimal set of variables: Variables that best sensed the change in every one of the chosen pixels for the year the fire occurred in, the previous and the following years were determined, with a 2 month period for pre-fire and post-fire each year, based on the hypothetical date of ignition.The start of the fires is determined by the maximum value of the Burn Boreal Forest Index (BBFI) derived from the NIR and LST bands (Table 3), whose expression is [23] The study of fire evolution radiometric behaviour and its connection with the potential indexes eventually led to the choice of 12 variables falling under 3 categories:

•
In order to sense pixels of the Non-Burned class (with vegetation or bare ground) the year before the fire, variables measuring the differences in the GEMI [73], BBFI, and NIR indexes, means between the post-fire period of the year the fire occurred in and the previous year were used, as well as the maximum value of the BBFI index for the year before the fire.GEMI is formulated as GEMI = η(1 − 0.25η) − (RED − 0.125)/(l − RED) where η = [2(NIR 2 − RED 2 ) + 1.5NIR + 0.5RED]/(NIR + RED + 0.5) [73].

•
For the detection of pixels whose Burned class remained in this class in the year after the fire, we used the variables that measure the maximum of BBFI and LST in the year after the fire, together with the difference of the medians of the GEMI index between the year of the fire and the year after.

•
For the detection of pixels whose class was Burned with vegetation in the year of the fire, we used the variables that measured the maxima of BBFI, GEMI, and LST, respectively, and the minimum of NIR in the year of the fire, and the differences of the medians of GEMI before and after the fire in the year of the fire.

5.
Design of the classification algorithm: Based on the selected training cases set, a data file with a format compatible with the WEKA software package from the University of Haikato (New Zealand) [74] was constructed.The value for each variable appears together with the thematic class to which it belongs (Burned or Non-Burned).Resulting from the application of the Bayes-Net classifier in the WEKA set to the training cases set, an algorithm allowing for determination of the likelihood for a pixel to belong to a given class, based on partial probabilities of each of the variables under consideration, is obtained.Based on the minimum commission and omission errors that are obtained from the classification of the training set itself-through a process of refinement of training cases-the classification model is chosen and applied on the dataset in its totality, obtaining a map of probabilities; each variable has an associated probability depending on the class to which it belongs (Burned, Non-Burned).The probability map is obtained by normalizing the joint probability of all variables, taking as positive those belonging to the Burned class and, as negative, those in the Non-Burned class.In order to connect the probability ranges with the thematic classes-Burned and Non-Burned-a discretisation was conducted by adjusting the total number of burned pixels on the grounds of a given probability of the Burned Area estimation obtained in the study area.In this case, we have used the criterion that only pixels with a positive probability are deemed burned.6.
Post-processing: To improve the result, in the post-processing stage we require spatial coherence using a filtering process based on the theory of cellular automata [75].Cellular automata have previously been applied to improve the results of satellite image mapping [76,77].Starting with an initial classification state based on burned probabilities, a series of transition rules are applied repeatedly until a final state is reached.Pixels with probabilities greater than 0.90 were labelled as Burned; pixels with probabilities between 0.00 and 0.90 were labelled as Partially-Burned; otherwise they were labelled as Unburned.The transition rules chosen were: • A pixel whose class is Burned remains in its status if it has at least one neighbour classified as Burned or Partially-Burned.Otherwise, its status becomes Partially-Burned.

•
A pixel whose class was Partially-Burned moves through to the Burned status if it has at least one neighbour classified as Burned, and two or more neighbours classified as Partially-Burned.Otherwise, its status becomes Unburned.

•
A pixel whose class was Unburned cannot change its status.
Finally, all pixels labelled as Burned or Partially-Burned were considered as Burned pixels.
The ten-day Burned Area maps with 250 m spatial resolution in a UTM-28N projection were generated by the Bayesian network algorithm from the daily images.Then, these maps were combined to generate the annual maps of Burned/Non-Burned pixels.Finally, all pixels labelled as Burned or Partially-Burned were considered as Burned pixels.
The ten-day Burned Area maps with 250 m spatial resolution in a UTM-28N projection were generated by the Bayesian network algorithm from the daily images.Then, these maps were combined to generate the annual maps of Burned/Non-Burned pixels.

Accuracy Assessment of the Estimates of Burned Area
All annual maps (BA products and reference) were clipped to the defined study area and masked with coastline boundaries for every island.The maps were resized to a UTM-28N projection with a 250 m × 250 m pixel size 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.To assess the spatial accuracy, the Burned Area proportion of each product map versus the reference map was calculated on a grid of 2.5 km × 2.5 km cells evenly distributed over the entire study area.A simple regression analysis of these proportions was carried out [39].Then, the error matrix and

Accuracy Assessment of the Estimates of Burned Area
All annual maps (BA products and reference) were clipped to the defined study area and masked with coastline boundaries for every island.The maps were resized to a UTM-28N projection with a Remote Sens. 2018, 10, 789 10 of 21 250 m × 250 m pixel size 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.To assess the spatial accuracy, the Burned Area proportion of each product map versus the reference map was calculated on a grid of 2.5 km × 2.5 km cells evenly distributed over the entire study area.A simple regression analysis of these proportions was carried out [39].Then, the error matrix and commission and omission errors [78] were computed considering a 2.5 km pixel size [33] to prevent georeference-derived errors.
Lastly, the efficiency in the detection of fires greater than 2000 ha which occurred during the study period was assessed for each Burned Area product.A set of subscenes measuring 30 km × 30 km that included, at their centre, one of these large burned areas for all of the products and reference data was developed.The total amount of BA was computed for each subscene.The detection of a minimum of 10% of burned pixels inside a subscene, with respect to the reference data, was the criterion considered to identify the presence of a large BA.

Burned Area Distribution in the Period 2002-2012
Table 5 shows the total area sensed by the proposed Bayesian algorithm (hereinafter called BY-MODIS) and official MODIS Burned Area products.For the period under consideration, fourteen registered fires are indicated in the studied area.The total BA estimated over the entire study period for each of these BY-MODIS, MCD45A1, and MCD64A1 products was 36,050, 36,900, 42,850 ha, respectively, and 45,856 ha for LANDSAT.The correlation coefficient between the annual Burned Area estimates for each product in relation to the reference data in the period considered was 0.94, 0.98, 0.98, and 0.99, respectively.The MCD64A1 product sensed the largest Burned Area, despite not being able to sense any of the smaller fires (2003, 2005_1, 2007_2, 2008, and 2012_5).MCD45A1 failed to detect the fire labelled as 2005_1, and barely 25 ha out of a total 306 ha burned in 2012_5.Regarding BY-MODIS, it must be noted it failed to sense the fire occurred in the island of El Hierro in 2006 and, in line with official MODIS products, it had trouble sensing the fire 2005_1, despite estimating a Burned Area close to 20% of the reference total (130 ha).On the other hand, we observed an overestimation in three Burned Area When spatial distribution of the burned areas was analysed by islands (Table 6), the biggest differences in relative percentages found corresponded to the islands of El Hierro and Gran Canaria.In the first case, MCD64A1 estimated close to a 92% of the reference data, whereas other products showed an underestimation around 30%.For Gran Canaria, MCD64A1 recorded a good estimation of the Burned Area (79.6%).Conversely, both BY-MODIS and MCD45A1 only sensed around 40%. Results are comparable to those for the remaining islands.Regarding temporal distribution of fires for the 2002-2012 period, and considering the MAPAMA reference data, a rather uneven distribution was found.The year 2007 exhibited the highest fire activity (35,612.17ha) followed by the year 2012, whereas the years 2002, 2004, 2010, and 2011 did not register any fires.Along with the reference data and BA derived from Landsat imagery, Figure 3 shows the annual distribution of the total BA in the western Canary Islands from 2002 to 2012 for the products under analysis.Regarding temporal distribution of fires for the 2002-2012 period, and considering the MAPAMA reference data, a rather uneven distribution was found.The year 2007 exhibited the highest fire activity (35,612.17ha) followed by the year 2012, whereas the years 2002, 2004, 2010, and 2011 did not register any fires.Along with the reference data and BA derived from Landsat imagery, Figure 3 shows the annual distribution of the total BA in the western Canary Islands from 2002 to 2012 for the products under analysis.

Spatial Accuracy of the Burned Area Products
The results of the accuracy of the BA estimate obtained from the scatter plots of the various products versus the reference map using a 2.5 km × 2.5 km grid are shown in Table 7.The determination coefficient (R 2 ) over the 2002-2012 period for the BY-MODIS, MCD45A1, and MCD64A1 products was 0.82, 0.81, and 0.88, with a slope of 0.76, 0.81, and 0.93, respectively.

Spatial Accuracy of the Burned Area Products
The results of the accuracy of the BA estimate obtained from the scatter plots of the various products versus the reference map using a 2.5 km × 2.5 km grid are shown in Table 7.The determination coefficient (R 2 ) over the 2002-2012 period for the BY-MODIS, MCD45A1, and MCD64A1 products was 0.82, 0.81, and 0.88, with a slope of 0.76, 0.81, and 0.93, respectively.
Regarding the errors derived from the confusion matrix for each year (Table 8), the average commission errors for the BY-MODIS, MCD45A1, and MCD64A1 products were 0.11, 0.13, and 0.12, whereas the omission errors were 0.30, 0.28, and 0.18, respectively.

Detection of Fires Greater than 2000 ha
Table 9 shows the results of the study of the various 30 km × 30 km subscenes containing all fires greater than 2000 ha that were recorded in the 2002-2012 period in the Canary Islands.The total BA encompassed by all these subscenes was 50,175 ha (approximately 90% of the total area burned across the study area in that period).The BA recorded was less than 10% in none of the subscenes of each product, and was greater than 100% in the 2009 subscene of all products.The average relative percentage of Burned Area estimated for Landsat images was 83%.We now proceed to describe the rest of products below.MCD64A1 detected the largest total Burned Area of the six large fires recorded (78.18% over the reference data), followed by MCD45A1 (69.36%), and BY-MODIS (66.73%).MCD64A1 achieved the best results in three of them (2007_1, 2012_1, and 2012_3), BY-MODIS in two (2009 and 2012_4) and MCD45A1 in one (2007_3).On the other hand, BY-MODIS and MCD45A1 got the worst results in the 2007_1 fire, which has the largest amount of Burned Area, detecting around 40% of the Burned Area.For the MCD64A1 product, the result diverging the furthest from the reference value corresponded to fire 2007_3-the second in terms of Burned Area-detecting around 54% of the Burned Area.Also noteworthy is that this product showed the biggest overestimation in the year 2009, about 23% of the Burned Area.
The spatial distribution of the burned pixels for the BY-MODIS, MCD45A1, and MCD64A1 products (in their native resolution) for each of the considered subscenes and the Landsat perimeter are shown in Figure 4.

Discussion
All products, BY-MODIS, MCD45A1, and MCD64A1, including LANDSAT, underestimated the total Burned Area for the set of fourteen fires over 100 ha recorded in the study region for the 2002-

Discussion
All products, BY-MODIS, MCD45A1, and MCD64A1, including LANDSAT, underestimated the total Burned Area for the set of fourteen fires over 100 ha recorded in the study region for the 2002-2012 period.The difference between the total Burned Area according to MAPAMA and the LANDSAT polygons may be due to the existence of islands of unburned vegetation in the official product, but also to the construction of the LANDSAT polygons with images that, in some cases, differ months from the date of the fire.Regardless of the algorithmic strategy used to tell burned and non-burned pixels apart, all methods are normally designed to be as accurate as possible; that it to say, to obtain low commission errors, even if that means omitting a larger number of burned pixels [39,79].In this respect, results of this study are in agreement with most previous published works using MODIS data for various regions and ecosystems [33,34,39,58,[80][81][82][83][84][85][86][87][88].However, a great variability was found in Burned Area percentages detected by the products under consideration for the various fires (Tables 5 and 6).
There are no precedents in the literature about the use of MODIS Burned Area products in the Canary Islands region to analyse a historical series of over ten years, such as the one discussed in this work.Only a few studies have been published on the two large fires that occurred simultaneously in 2007 in Tenerife and Gran Canaria [89,90].For this reason, our results were compared and discussed with works conducted on Mediterranean ecosystems similar to the Canary Islands and those that included fires in regions of topographic complexity.We cannot ignore the fact that the five western islands of the Canarian Archipelago historically affected by forest fires show a highly rugged terrain of a volcanic origin, with numerous narrow ravines scoring them radially which, in many cases, may disguise the correct sensing of burned areas with the spatial resolution of MODIS products [34,91].
Analysis of the temporal accuracy proved that annual distribution of Burned Area for the BY-MODIS, MCD45A1, and MCD64A1 products shows a good correlation with the reference set, resulting in determination coefficients (R 2 ) of 0.82, 0.81, and 0.88, respectively (Table 7).Similarly to results found by Moreno-Ruiz et al. [58] in the North American Boreal Forest, the year 2006 presented the biggest errors, both in terms of accuracy (Table 7) and spatial precision (Table 8), for the BY-MODIS and MCD45A1 products.By contrast, 2006 was found to be one of the years with the best results for the MCD64A1 product.This fact can be observed in the fire that occurred in the island of El Hierro in September 2006 (Tables 5, 7 and 8).In this case, both BY-MODIS and MCD45A1, with omission errors greater than 98%, were hardly able to sense a single burned pixel within the fire polygon, compared to MCD64A1, with only a 2% omission error and a 24% commission error, possibly due to the identification of isolated regions of bare volcanic ground as burned areas.The use of active fires in the MCD64A1 provides very valuable information to correctly estimate the Burned Area in such cases, where the affected vegetation is basically comprised of Canarian pine trees (conifers) with a forest canopy cover ranging from 60% to 90% and under water stress.Three consecutive months without rainfall and with hot temperatures, alongside intense and long-lasting southerly winds (information retrieved from http://www.aemet.es/es/lineas_de_interes/datos_y_estadistica),led to a relative humidity lower than 20% (whereas, under normal conditions, it exceeds 70%), which made it easier for the fire to spread in just six hours.This situation makes it harder for BY-MODIS or MCD45A1 products, which are based on algorithms seeking to sense drastic changes in post-fire reflectances, to be effective for this fire.Tsela et al. [82] obtained comparable results for pine forest plantations in Sabie (South Africa), with omission errors exceeding 80%.
The year 2005 failed to record satisfactory results for either of the products, with omission errors of 39% for MCD64A1 and of 69% for MCD45A1.The fire dominating statistics for this year corresponded to the one in the north of the island of La Palma, with a 1890 ha Burned Area, according to official figures.The fire affected an area of approximately 7 km long, west-to-east, by 3 km long, north-to-south, but with an altitude variation of over 1000 m in the latter-from 1000 to 2000 m above the sea level-and essentially burned the Canarian pine forest in the lower part and brushwood in the higher part.The steep land, with a highly rugged terrain and deep ravines filled with shade areas, made it difficult to map the Burned Area for all products, as other authors previously pointed out in a mountainous area of Northwest Yunnan (China) [34] or for the fynbos, a Mediterranean shrubland of the southwestern area of South Africa [85,92].The same argument would be valid to explain big commission and omission errors (Table 8) for the only fire occurred during the year 2008.The 376 affected hectares were located in the north of the island of La Gomera.The area shows a highly heterogeneous coverage, including numerous farming allotments, brushwood, palm trees and heaths located on a topographically complex terrain, making its modelling in two dimensions-using MODIS sensors' spatial resolution-difficult.
In regard to spatial accuracy of the burned areas estimates, the overall linear regression slope between the results for the BY-MODIS, MCD64A1, and MCD45A1 products, and the reference dataset (Table 7) were similar, with values of 0.76, 0.93, and 0.81, respectively, and are also in line with the results obtained by [79] for the MCD45A1 product.These authors [79] validated MCD45A1 in six different regions with Mediterranean ecosystems in the year 2003.They obtained a commission error of 17.4%, and an omission error of 38.6%, with a 0.79 slope and a 0.84 determination coefficient.However, these results are a long way from the validations conducted by authors [41] for the MCD45 product for Mediterranean Forest biome, with characteristics similar to those of our case study, with commission errors of 69% and omission errors of 95%.
In general, BY-MODIS and MCD45A1 commission and omission errors are similar, but slightly different to those of MCD64A1.The total omission error for the MCD64A1 was around 10% lower than in the other two products.Those results are comparable to those obtained by the authors in the boreal forest region of North America using the same BA algorithm on daily images from the LTDR dataset [33,57].
Another noteworthy difference of the MCD64A1 product with respect to BY-MODIS and MCD45A1 was its difficulty to sense small fires.This product failed to identify a single pixel as burned in the five fires with fewer than 400 ha.Similar findings were reported by Tsela et al. [82], who suggested the combination of the two MODIS BA products (MCD45A1 and MCD64A1) would improve detection of small burned scars in South Africa.Fornacca et al. [34] also report the worse performance of MCD64A1 in the detection of small burned areas in the mountainous region of Yunnan (China) in the years 2006 and 2009, compared to the product that does not use active fires (MCD45A1).
For the largest fires (>2000 ha), although the best total performance corresponded to the MCD64A1 product in both, in the estimated accuracy and spatial precision, we observe a high variability for the Burned Area detected for these fires (Figure 4, Table 9).For the 2007_1 fire-the largest fire, greater than 18,000 ha-MCD64A1 obtained the best Burned Area estimation, twice above the other BA products, which led to an unavoidable bias in favour of this product in global statistics (Table 9).In this case, the Burned Area mainly corresponded to herbaceous and non-wooded vegetation by 63%.The fire took place in the southern slope of the island of Gran Canaria, an area with little forest coverage which is very arid all year round.Once again, the MCD64A1 methodology with active fires allowed for a better total estimation according to the reference set, compared to the other two products.Conversely, for the second largest fire, the 2007_3, which took place mainly in the northern slope of the island of Tenerife, largely affecting a dense and wet Canarian pine tree forest, both BY-MODIS and MCD45A1 improved on MCD64A1 results.The cause could be rooted in the thick smoke and steam produced during biomass combustion, preventing active fires in the central part of the Burned Area from being recorded.There were no remarkable differences between the three products for the remaining big fires.

Conclusions
We have evaluated a Bayesian Algorithm (BY-MODIS) to detect burned areas in the Canary Islands' dry woodlands and forests ecoregion using daily time series of satellite images from MODIS with a 250 m spatial resolution.The results obtained by BY-MODIS were compared with MODIS databased Burned Area products MCD64A1 Version 6 and MCD45A1 Version 5.1 for the 2002-2012 period.This is the first study dedicated to mapping burned areas for such a long period of time in this particular region.Despite these islands' topographic complexity and the diverse type of ecosystems affected by the fire found in such a small and fragmented territory, BY-MODIS and the other two products are capable of offering good temporal accuracy when assessed with reference datasets, resulting in determination coefficients (R 2 ) of 0.82, 0.81, and 0.88, for the BY-MODIS, MCD45A1, and MCD64A1, respectively.In regard to spatial accuracy, the Burned Area was underestimated in all products.However, the results should be separated according to the size of the fires.For fires bigger than 2000 ha, the best total performance corresponded to the MCD64A1 product outperforming the other two products by approximately 10%.The use of active fires seems to provide very valuable information to correctly estimate the Burned Area in such cases.On the other hand, the MCD64A1 could not detect any of the five fires with less than 400 ha, being BY-MODIS, the algorithm with the best estimated accuracy and spatial precision in this case.BY-MODIS and MCD45A1 products, which are based on algorithms seeking to sense drastic changes in post-fire reflectances, are effective for these small fires.Taking into account both spatial accuracy in fires of any size and its applicability to the assorted sizes within the period under consideration, BY-MODIS can be regarded as a reliable alternative to MCD54A1 and MCD64A1 products.
The results of this study prove the capacity to adapt and the versatility of a Bayesian algorithm-previously developed and successfully applied to regions in the boreal forest-in other regions showing different ecosystems and satellite images with diverse spatial resolution, as long as a set of training areas adequately characterising the dynamics of the forest canopy affected by the fire is defined.

Figure 1 .
Figure 1.The Canary Islands Archipelago (Spain) is located off the north-west African coast (Morocco).The image is a NDVI composition from the Terra-MODIS sensor at a 500-m spatial resolution on 30 December 2015 retrieved from the online Land, Atmosphere Near Real-time

Figure 1 .
Figure 1.The Canary Islands Archipelago (Spain) is located off the north-west African coast (Morocco).The image is a NDVI composition from the Terra-MODIS sensor at a 500-m spatial resolution on 30 December 2015 retrieved from the online Land, Atmosphere Near Real-time Capability for the EOS (LANCE) system (https://lance-modis.eosdis.nasa.gov/imagery/subsets/?area=af).The red rectangle includes five islands studied.

22 
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of A pixel whose class is Burned remains in its status if it has at least one neighbour classified as Burned or Partially-Burned.Otherwise, its status becomes Partially-Burned. A pixel whose class was Partially-Burned moves through to the Burned status if it has at least one neighbour classified as Burned, and two or more neighbours classified as Partially-Burned.Otherwise, its status becomes Unburned. A pixel whose class was Unburned cannot change its status.

Figure 2 .
Figure 2. Flowchart of the Bayesian algorithm applied to the study region to obtain annual maps of Burned Area.

Figure 2 .
Figure 2. Flowchart of the Bayesian algorithm applied to the study region to obtain annual maps of Burned Area.

Figure 3 .
Figure 3. Annual distribution of the Burned Area estimate (ha) in the study region for the analysed products and the reference data.

Figure 3 .
Figure 3. Annual distribution of the Burned Area estimate (ha) in the study region for the analysed products and the reference data.

Figure 4 .
Figure 4.The 30 km × 30 km subscenes with fires greater than 2000 ha for the MODIS-based MCD45A1, MCD64A1, and BY-MODIS.The black lines represent the Landsat perimeter of all the fires registered on each subscene.The geographic coordinates indicate the top left and bottom right corners of each of the subscenes.

Figure 4 .
Figure 4.The 30 km × 30 km subscenes with fires greater than 2000 ha for the MODIS-based MCD45A1, MCD64A1, and BY-MODIS.The black lines represent the Landsat perimeter of all the fires registered on each subscene.The geographic coordinates indicate the top left and bottom right corners of each of the subscenes.

Table 1 .
Fires registered for the 2002-2012 period in the Canary Islands: ID (identification code), Region (island name), Detection Date, Days until extinction, geographical coordinates (Latitude and Longitude) and Burned Area (ha) according to Ministry of Agriculture and Fisheries, Food and Environment of Spain (MAPAMA).

Table 2 .
Spatial distribution of the fires registered in the Canary Islands for the 2002-2012 period: Island, Surface (km 2 ), Fires Registered (count and percentage), Burned Area (ha, percentage in relation to the total Burned Area and the Island Surface).

Table 3 .
Bands of the MOD09GQ and MOD11A1 products used for building the daily set of images with the same format than the Long Term Data Record (LTDR) dataset.

Table 5 .
Burned Area distribution (ha) by fire event for the different MODIS-based products in the Canary Islands during the period 2002-2012.
1Reference data are from MAPAMA information.

Table 6 .
Burned Area distribution by Island for the different MODIS-based products in the Canary Islands during the period 2002-2012.

Table 6 .
Burned Area distribution by Island for the different MODIS-based products in the Canary Islands during the period 2002-2012.

Table 7 .
Burned Area estimate accuracy by year (Intercept, Slope and R 2 ) for MCD45A1, MCD64A1, and BY-MODIS.

Table 7 .
Burned Area estimate accuracy by year (Intercept, Slope and R 2 ) for MCD45A1, MCD64A1, and BY-MODIS.

Table 9 .
Detection of fires greater than 2000 ha in the Canary Islands' forestry region for MCD45A1, MCD64A1, and BY-MODIS.