Monitoring Wildﬁres in the Northeastern Peruvian Amazon Using Landsat-8 and Sentinel-2 Imagery in the GEE Platform

: During the latest decades, the Amazon has experienced a great loss of vegetation cover, in many cases as a direct consequence of wildﬁres, which became a problem at local, national, and global scales, leading to economic, social, and environmental impacts. Hence, this study is committed to developing a routine for monitoring ﬁres in the vegetation cover relying on recent multitemporal data (2017–2019) of Landsat-8 and Sentinel-2 imagery using the cloud-based Google Earth Engine (GEE) platform. In order to assess the burnt areas (BA), spectral indices were employed, such as the Normalized Burn Ratio (NBR), Normalized Burn Ratio 2 (NBR2), and Mid-Infrared Burn Index (MIRBI). All these indices were applied for BA assessment according to appropriate thresholds. Additionally, to reduce confusion between burnt areas and other land cover classes, further indices were used, like those considering the temporal di ﬀ erences between pre and post-ﬁre conditions: di ﬀ erential Mid-Infrared Burn Index (dMIRBI), di ﬀ erential Normalized Burn Ratio (dNBR), di ﬀ erential Normalized Burn Ratio 2 (dNBR2), and di ﬀ erential Near-Infrared (dNIR). The calculated BA by Sentinel-2 was larger during the three-year investigation span (16.55, 78.50, and 67.19 km 2 ) and of greater detail (detected small areas) than the BA extracted by Landsat-8 (16.39, 6.24, and 32.93 km 2 ). The routine for monitoring wildﬁres presented in this work is based on a sequence of decision rules. This enables the detection and monitoring of burnt vegetation cover and has been originally applied to an experiment in the northeastern Peruvian Amazon. The results obtained by the two satellites imagery are compared in terms of accuracy metrics and level of detail (size of BA patches). The accuracy for Landsat-8 and Sentinel-2 in 2017, 2018, and 2019 varied from 82.7–91.4% to 94.5–98.5%, respectively.


Introduction
Wildfires have been a common phenomenon throughout time and are caused by ignition sources, either human or natural, and their interaction with climate factors that foster combustion combustion and propagation. Their extension and frequency are directly related to the temperature, humidity, wind, the type and characteristics of vegetation, relief, among other factors [2]. The locals use the fires to clean up the terrain and prepare the soil for planting [43]. Considering what was previously mentioned, the goal of this research is to develop a decision rules-based routine for monitoring wildfires employing Landsat-8 and Sentinel-2 imagery and using the cloud-based GEE platform. A comparison between the BA products derived from these two satellite imagery collections is also envisaged.

Study Area
The Amazonas department is part of the Peruvian Amazon and extends over a surface of 39,306.47 km 2 (Figure 1, 3 • 0 -7 • 2 S and 77 • 0 -78 • 42 W), with elevations that range from 120 m.a.m.s.l. in the north to 4900 m.a.m.s.l. (meters above mean sea level) in the south. It comprises seven provinces (Bagua, Bongará, Chachapoyas, Condorcanqui, Luya, Rodríguez de Mendoza, and Utcubamba), with a total population of 379,384 inhabitants and a density of 9.6 inh./km 2 [44]. This administrative division is composed of the interandean sector, which occupies 27.0% of its surface, and the Amazon Forest sector, which corresponds to the remaining 73.0% of the area. The Amazonas department contains a great diversity of natural resources, with complex and fragile ecosystems (Andean, lowland forest, highland forest, and tropical dry forest) [45], mostly covered by unexplored tropical forests [46]. The climate changes according to the geographic conditions, although the "warm and wet" climate prevails in Condorcanqui and in part of Rodríguez de Mendoza. There are other dissimilar climates scattered throughout the region, like the "warm and dry" climate in Bagua and Utcubamba provinces, and the "mildly wet and temperately warm" climate in Chachapoyas and Luya. The temperature rises up to 40 • C in the lowland forests (in the north) and the minimum temperature reaches 2 • C in the Andean Mountains (in the south) [45,47]. The land cover and land use in this region are predominantly represented by forests (forest formation and floodplain forest), followed by agriculture, natural formation (floodplain non-forest natural formation and grasslands), and non-vegetated areas [48].
The most outstanding economic activities in this region are agriculture (cropping and cattle raising) and the industrial activity, which respectively account for 51.2% and 14.5% of the regional Gross Domestic Product (GDP) [45]. These activities lead to deforestation, overgrazing, migratory agriculture, soil contamination and erosion, wildfires, and illegal logging and mining. It has been identified that 27.4% (11,533.93 km 2 ) of this department's surface is degraded [41,49]. During the period from 2002 to 2017, there were 130 forest fires [50]. This is a consequence of migratory agriculture and the opening of new accesses and land routes, such as connection axes for departmental and national economic development actions that affect a great variety of flora and fauna species [41,42]. Furthermore, the vegetation cover map of this region reported that only in 2018, 7453 ha of rain forests were lost [51]. Since the Amazonas department contains the greatest share of wildfires events in the Peruvian Amazon (Figure 2), it has been selected as our study area.

Data Processing
The adopted methodology is presented in Figure 3, which includes Landsat-8 and Sentinel-2 multispectral data acquisition; masking of clouds, cloud shadow, and water bodies; calculation of spectral indices; classification; spatiotemporal filtering; mapping of BA; and statistical validation.

Data Processing
The adopted methodology is presented in Figure 3, which includes Landsat-8 and Sentinel-2 multispectral data acquisition; masking of clouds, cloud shadow, and water bodies; calculation of spectral indices; classification; spatiotemporal filtering; mapping of BA; and statistical validation. respectively [52]. The surface of burn scars in the Amazonas department totalized 189.87 km 2 [53]. Usually, these fires are caused by human action through agricultural practices, while in some other cases, they are also triggered by ignition sources, both of them with dramatic effects on the Amazon ecosystems [41,42]. In this context, it is extremely important to detect burnt areas so as to assess their spatiotemporal pattern as well as to conceive and implement prevention actions.

Data Processing
The adopted methodology is presented in Figure 3, which includes Landsat-8 and Sentinel-2 multispectral data acquisition; masking of clouds, cloud shadow, and water bodies; calculation of spectral indices; classification; spatiotemporal filtering; mapping of BA; and statistical validation.

Google Earth Engine (GEE)
This platform operates in the logic of "Big Data" and contains not only Landsat-8 and Sentinel-2 data, but further remotely sensed raster data possessing an order of magnitude of petabytes and rendered publicly available by agencies like the National Aeronautics and Space Administration (NASA), the European Space Agency (ESA), and others [14]. GEE also enables access to digital images processing and machine learning algorithms as well as to libraries of Application Programming Interfaces (API) (such as JavaScript, Python, and R) and scripting interfaces, in which users can develop their own codes and process data online. Thus, GEE allows the manipulation, analysis, and visualization of geospatial data without the need to access supercomputers [54].

Satellite Multispectral Data and Fire Products
There are several products that enable to detect BA. In this study, we employed medium spatial resolution satellite data available in the GEE platform, like the collection of Surface Reflectance (SR) multispectral images of Landsat-8 OLI/TIRS Tier 1 (ID: LANDSAT/LC08/C01/T1_SR) and the Top-Of-Atmosphere (TOA) images in the case of Sentinel-2/MIS level-1C (ID: COPERNICUS/S2) [55]. The data acquisition and processing comprised the time span from 2017 to 2019, in view of data availability (Tables 1 and 2). Each imagery collection was individually classified. For Landsat-8, we worked with the original spatial resolution of 30 m, while for Sentinel-2, we adopted a standard spatial resolution of 15 m, obtained by means of resampling using the nearest neighbor interpolation technique.

Masking of Clouds, Cloud Shadow, and Water Bodies
Masking was an important processing step due to the fact that the study area presents clouds during most of the year. This process allowed the masking of clouds, cloud shadow, and water bodies. The algorithm C Function of Mask (CFMASK) was applied to Landsat-8 images, so as to detect clouds and cloud shadow [56]. In the case of the Sentinel-2 images, the algorithm Google CloudScore was applied to mask clouds, and the algorithm Temporal Dark Outlier Mask (TDOM) to mask cloud shadow, with the aid of the Quality Assessment band (QA60) [57]. For the sake of masking pixels presenting noise derived from cloud shadow, a threshold was applied to band B12, where B12 < 0.07 was considered cloud shadow [12]. Considering the often presence of cirrus clouds, a threshold was also applied to band B2 (B2 > 0.15) [58]. Lastly, the Normalized Difference Water Index (NDWI) was applied to detect and mask water bodies [59].

Google Earth Engine (GEE)
This platform operates in the logic of "Big Data" and contains not only Landsat-8 and Sentinel-2 data, but further remotely sensed raster data possessing an order of magnitude of petabytes and rendered publicly available by agencies like the National Aeronautics and Space Administration (NASA), the European Space Agency (ESA), and others [14]. GEE also enables access to digital images processing and machine learning algorithms as well as to libraries of Application Programming Interfaces (API) (such as JavaScript, Python, and R) and scripting interfaces, in which users can develop their own codes and process data online. Thus, GEE allows the manipulation, analysis, and visualization of geospatial data without the need to access supercomputers [54].

Satellite Multispectral Data and Fire Products
There are several products that enable to detect BA. In this study, we employed medium spatial resolution satellite data available in the GEE platform, like the collection of Surface Reflectance (SR) multispectral images of Landsat-8 OLI/TIRS Tier 1 (ID: LANDSAT/LC08/C01/T1_SR) and the Top-Of-Atmosphere (TOA) images in the case of Sentinel-2/MIS level-1C (ID: COPERNICUS/S2) [55]. The data acquisition and processing comprised the time span from 2017 to 2019, in view of data availability (Tables 1 and 2). Each imagery collection was individually classified. For Landsat-8, we worked with the original spatial resolution of 30 m, while for Sentinel-2, we adopted a standard spatial resolution of 15 m, obtained by means of resampling using the nearest neighbor interpolation technique.

Mosaic of Pre and Post-Fire Images
This process initially grouped (mosaiced) the time series data by every two months, which was meant not only to build a mosaic with a minimum amount of clouds, but also to reduce the input data dimensionality. In the sequence, pre-selected filters were applied for excluding pixels with clouds, cloud shadow, and water bodies from the Landsat-8 and Sentinel-2 multitemporal data. For the creation of pre and post-fire mosaics, the NBR and NDVI were applied. The elaboration of the pre-fire mosaic consisted in selecting the pixels with maximum values of the NDVI along the time series, while for the post-fire mosaic, the pixels with maximum values of the NBR were selected.
The NDVI is one of the most known and widely used vegetation indices, which takes into account the amount of absorbed visible red light by the green vegetation. In a general way, the NDVI is sensitive to the content of chlorophyll and other pigments responsible for absorbing solar radiation in the red range of the electromagnetic spectrum, while the Enhanced Vegetation Index (EVI) is more sensitive to the variation in the canopy structure, including the Leaf Area Index (LAI), the plant physiognomy, and the canopy volume [62,63]. The EVI and the Soil Adjusted Vegetation Index (SAVI) make use of atmospheric correction and soil parameters that would lead to a greater uncertainty in the calculation of thresholds to identify burnt areas. The fire-induced damage caused to the vegetation results in a substantial decrease in the NDVI [34,64], and hence, we have selected the NDVI for BA detection.
The time series analyses revealed that the NBR2 has shown to be efficient in discriminating burnt from unburnt areas [12,38]. In this way, the pixels with maximum values of the NDVI and NBR2 were selected for differentiating pre-fire from post-fire areas (BA), respectively.
Based on the generated information, it was then necessary to set appropriate thresholds for the pre and post-fire differential spectral indices, namely dMIRBI [65], dNBR [66], dNBR2, and dNIR [65]. This procedure was executed considering 444 sampled points by means of visual analysis with the aim of defining the BA thresholds [12] (Figure 4). The application of these differential spectral indices to the images led to the classification of "burnt" and "unburnt" areas ( Figure 5) according to the defined thresholds [67], which then framed the decision rules-based routine in the GEE platform [68].

Images Filtering
A temporal filter based on the NBR was applied for the purpose of eliminating inconsistencies (or changes in BA) and also to correct errors due to cloud plumes, absence of information, or duplicity of fire dates. To this end, the oldest date in the results was considered for classification. Such temporal filter was applied to the time series of each year and consisted of defining rules of permanence for

Images Filtering
A temporal filter based on the NBR was applied for the purpose of eliminating inconsistencies (or changes in BA) and also to correct errors due to cloud plumes, absence of information, or duplicity of fire dates. To this end, the oldest date in the results was considered for classification. Such temporal filter was applied to the time series of each year and consisted of defining rules of permanence for

Images Filtering
A temporal filter based on the NBR was applied for the purpose of eliminating inconsistencies (or changes in BA) and also to correct errors due to cloud plumes, absence of information, or duplicity of fire dates. To this end, the oldest date in the results was considered for classification. Such temporal filter was applied to the time series of each year and consisted of defining rules of permanence for the NBR values. Regarding the Landsat-8 data, an area was assigned as BA if it presented at least one positive value of the dNBR along the whole time series (due to the limited images availability). For each of the Sentinel-2 annual data, at least two positive values of the dNBR in the whole time series were required to ascribe an area as BA (Figure 6). A spatial filter with a minimum mapping unit of 0.5 ha [69] was also applied. Isolated and border pixels were eliminated by means of neighborhood filters, aiming at improving the data spatial consistency. Areas with rice cultivation were also excluded (on sites with slope < 5 • ), particularly in the Utcubamba river valley. the NBR values. Regarding the Landsat-8 data, an area was assigned as BA if it presented at least one positive value of the dNBR along the whole time series (due to the limited images availability). For each of the Sentinel-2 annual data, at least two positive values of the dNBR in the whole time series were required to ascribe an area as BA ( Figure 6). A spatial filter with a minimum mapping unit of 0.5 ha [69] was also applied. Isolated and border pixels were eliminated by means of neighborhood filters, aiming at improving the data spatial consistency. Areas with rice cultivation were also excluded (on sites with slope < 5°), particularly in the Utcubamba river valley.

Statistical Validation of BA
Error matrices were elaborated based on the BA classification and reference data (n sampled observations), which enabled the evaluation of the results' accuracies for each analyzed year [70]. The comparison between the BA product and the reference data was based on the proportion of each pixel classified as burnt or unburnt in the Landsat-8 and Sentinel-2 images. Thus, we calculated the user's accuracy (UA), which corresponds to commission errors (from the user's perspective), and the producer's accuracy (PA), associated with omission errors (from the perspective of the map producer) [71,72]. Furthermore, we calculated the overall accuracy, which indicates the quality of the classification result in terms of quantity; F-Score (F1), which provides relations between data's positive labels and those given by a classifier [73]; the Quantity component (Qt), which reveals the absolute differences in quantity for each category (burnt and unburnt areas) in both observed and classified data; the Exchange component (Et), which assesses the share of changes from burnt to unburnt areas and vice versa with the same proportions in the reference and classified data, without modifying the final quantities [74].

Identification of BA Patterns
With the aim of creating a continuous surface able to represent the wildfires intensity at the Amazonas department level, Kernel Density Estimation (KDE) was used [75]. The KDE [76,77]. Next, the density map was reclassified according to the Jenks natural breaks classification method [78] in five classes related to different levels of BA density: very low, low, moderate, high, and very high [28,41,50]. Finally, the quality of the spatial agreement was evaluated using Yule's Q coefficient [79].

Statistical Validation of BA
Error matrices were elaborated based on the BA classification and reference data (n sampled observations), which enabled the evaluation of the results' accuracies for each analyzed year [70]. The comparison between the BA product and the reference data was based on the proportion of each pixel classified as burnt or unburnt in the Landsat-8 and Sentinel-2 images. Thus, we calculated the user's accuracy (UA), which corresponds to commission errors (from the user's perspective), and the producer's accuracy (PA), associated with omission errors (from the perspective of the map producer) [71,72]. Furthermore, we calculated the overall accuracy, which indicates the quality of the classification result in terms of quantity; F-Score (F1), which provides relations between data's positive labels and those given by a classifier [73]; the Quantity component (Qt), which reveals the absolute differences in quantity for each category (burnt and unburnt areas) in both observed and classified data; the Exchange component (Et), which assesses the share of changes from burnt to unburnt areas and vice versa with the same proportions in the reference and classified data, without modifying the final quantities [74].

Identification of BA Patterns
With the aim of creating a continuous surface able to represent the wildfires intensity at the Amazonas department level, Kernel Density Estimation (KDE) was used [75]. The KDE [76,77]. Next, the density map was reclassified according to the Jenks natural breaks classification method [78] in five classes related to different levels of BA density: very low, low, moderate, high, and very high [28,41,50]. Finally, the quality of the spatial agreement was evaluated using Yule's Q coefficient [79].

Validation of the Achieved Results for Assessing BA
The BA statistical validation was achieved based on the visual interpretation of 456 test sites, so as to allow for a comparative evaluation of the UA, PA, accuracy, F1, Qt, and Et between the Landsat-8 and Sentinel-2 products. The accuracy obtained by Landsat-8 in the years 2017, 2018, and 2019 (

BA in the Amazonas Department
Wildfires in the Amazonas department showed to be of reduced extent in 2017 and 2018 as compared with 2019. An annual comparison of BA per product is presented in Table 5. In 2017, Landsat-8 registered an area of 16.39 km 2 , 0.96% smaller than the Sentinel-2 product, which assessed 16.5 km 2 of surface. In 2018, Landsat-8 recorded an extension of BA of 6.24 km 2 , while the Sentinel-2 product reported an area of 78.50 km 2 (1158.48% greater than Landsat-8). Finally, in 2019 the BA extension was of 32.93 and 67.19 km 2 (the latter one being 104.05% greater), as registered by Landsat-8 and Sentinel-2, respectively.  Figure 7 shows the occurrence of wildfires in the analyzed years. It can be inferred that the wildfires frequency increases between June and December and this behavior tends to repeat throughout the years. ISPRS Int. J. Geo-Inf. 2020, 9,

Spatial Agreement of BA
A spatial comparison of BA patches derived from the Landsat-8 and Sentinel-2 products per analyzed year is shown in Figure 8. The centroids of the BA polygons show a similar spatial distribution pattern in both satellite products, although Sentinel-2 has comparatively produced a greater number and larger patches of wildfires burn scars with respect to those generated by Landsat-8 during the evaluated years, in the face of its higher spatial (15 m) and temporal (5 days) resolutions.

Spatial Agreement of BA
A spatial comparison of BA patches derived from the Landsat-8 and Sentinel-2 products per analyzed year is shown in Figure 8. The centroids of the BA polygons show a similar spatial distribution pattern in both satellite products, although Sentinel-2 has comparatively produced a greater number and larger patches of wildfires burn scars with respect to those generated by Landsat-8 during the evaluated years, in the face of its higher spatial (15 m) and temporal (5 days) resolutions.

Patterns of BA Density
The maps elaborated based on the KDE method considering the three-year time detection of BA using the Landsat-8 and Sentinel-2 products (Figure 9) show "very high", "high", and "moderate"

Patterns of BA Density
The maps elaborated based on the KDE method considering the three-year time detection of BA using the Landsat-8 and Sentinel-2 products (Figure 9) show "very high", "high", and "moderate" wildfires densities in the southwest of the Amazonas department, amid Utcubamba, Luya, Chachapoyas, and Bongará provinces. This indicates that environmental policies and actions should pay greater attention to such vulnerable areas. On the other hand, the provinces Rodríguez de Mendoza, Bagua, and Condorcanqui register "low" and "very low" densities within their territories (Figure 9). wildfires densities in the southwest of the Amazonas department, amid Utcubamba, Luya, Chachapoyas, and Bongará provinces. This indicates that environmental policies and actions should pay greater attention to such vulnerable areas. On the other hand, the provinces Rodríguez de Mendoza, Bagua, and Condorcanqui register "low" and "very low" densities within their territories (Figure 9). It was observed that wildfires take place near human settlements (populated areas) and the road network. Approximately 20% of fires are located within a distance to settlements and roads from 0 to 600 m, 20% between 600 and 1200 m, and the remaining 60% of these events occur beyond 1200 m of distance (Supplementary Materials Figure S1). The wildfires in the study area mostly started on sites with slopes greater than 8° (nearly 90% of all events) and only 10% of them occur on sites with slopes between 0° and 8° ( Figure S2 and Table S8). With respect to the settlements' proximity, the majority of fires (90%) start within distances superior to 1600 m ( Figure S3 and Table S8).
The spatial agreement between BA resulting from Landsat-8 and Sentinel-2 is shown in Table 6. In 2017, the agreeing BA and unburnt areas (UA) between these two products had a surface of 3.51 and 39,277.04 km 2 , respectively, with an omission error of 0.21 and an extension coefficient of 1.00. In 2018, the agreeing BA had a size of 2.35 km 2 and an agreeing UA of 39,224.09 km 2 , presenting an omission error of 0.38 and an extension coefficient of 1.00. In 2019, there was an increase in the BA spatial agreement with regard to previous years (totaling an area of 18.88 km 2 ) and the agreeing UA was limited to 39,225.23 km 2 , with an omission error of 0.57 and an extension coefficient of 1.00. The Q coefficient, in its turn, presented values of 1.00, 0.99, and 1.00 for the years 2017, 2018, and 2019, respectively. It was observed that wildfires take place near human settlements (populated areas) and the road network. Approximately 20% of fires are located within a distance to settlements and roads from 0 to 600 m, 20% between 600 and 1200 m, and the remaining 60% of these events occur beyond 1200 m of distance (Supplementary Materials Figure S1). The wildfires in the study area mostly started on sites with slopes greater than 8 • (nearly 90% of all events) and only 10% of them occur on sites with slopes between 0 • and 8 • ( Figure S2 and Table S8). With respect to the settlements' proximity, the majority of fires (90%) start within distances superior to 1600 m ( Figure S3 and Table S8).
The spatial agreement between BA resulting from Landsat-8 and Sentinel-2 is shown in Table 6. In 2017, the agreeing BA and unburnt areas (UA) between these two products had a surface of 3.51 and 39,277.04 km 2 , respectively, with an omission error of 0.21 and an extension coefficient of 1.00. In 2018, the agreeing BA had a size of 2.35 km 2 and an agreeing UA of 39,224.09 km 2 , presenting an omission error of 0.38 and an extension coefficient of 1.00. In 2019, there was an increase in the BA spatial agreement with regard to previous years (totaling an area of 18.88 km 2 ) and the agreeing UA was limited to 39,225.23 km 2 , with an omission error of 0.57 and an extension coefficient of 1.00. The Q coefficient, in its turn, presented values of 1.00, 0.99, and 1.00 for the years 2017, 2018, and 2019, respectively.

Discussion
In the present study, a comparative performance evaluation of BA products derived from Landsat-8 and Sentinel-2 medium spatial resolution imagery was accomplished. These products were generated for the northeastern Peruvian Amazon by means of an open access decision rules-based routine implemented in the GEE cloud computing platform (Table S1). The accuracy of such products were assessed using error matrices and accuracy indices relying on a common reference dataset, created according to a stratified random sampling design [71,81]. The overall accuracy ranged from 82.7% to 98.5%, which indicates a satisfactory performance of the proposed routine for identifying burnt and unburnt areas. The Sentinel-2 products present greater OA (94.5-98.5%) than Landsat-8 (82.7-91.4%). This can be explained by the fact that the Peruvian Amazon has a wet climate and cloud-free Landsat-8 images are rare (16 days of temporal resolution) [82], while Sentinel-2 presents a better temporal (5 days) and spatial (10-20 m) resolution, with greater availability of images, which tends to improve the detection of BA [67].
The comparison between the Landsat-8 and Sentinel-2 products proved to be of value for a sound detection of BA in the Peruvian Amazon. The literature generally reports that Sentinel-2 offers higher accuracies [12,83], while some authors consider that Landsat-8 is superior to Sentinel-2 when compared with MCD64A1 products [60,84]. According to such authors, Landsat-8 owns better spectral bands than Sentinel-2. Other works integrate both Landsat-8 and Sentinel-2 products to improve BA estimates [85][86][87]. In this study, however, the comparison with further products, such as the series of "MODIS Combined Data" (MCD45, MCD64) and the "Fire Climate Change Initiative" (FIRE CCI), has not been made in the face of their unfeasible application in the Amazon region [4], due to large differences in spatial resolution [60]. Some studies have used the product MCD64A1 (coarse resolution) to detect wildfires at a global scale [88][89][90][91][92][93], but it is unable to detect small and spatially fragmented BA, with which the Landsat-8 and Sentinel-2 products can easily cope [60,94]. On the other hand, several techniques and methodologies have been proposed for mapping BA. One of the most often used procedures regards the bitemporal images that combine algebraic operations [4], such as the calculation of spectral indices (NDVI, NBR, NBR2, and MIRBI), and the pre and post-fire temporal differences between them (dMIRBI, dNBR, dNBR2, and dNIR), which improve the detection of changes in the vegetation cover, and hence, the identification of BA.
The wildfires in the study area generally take place from June to December (Figure 7). Manríquez-Zapata [42] reports that fires in the vegetation cover occur in the third and fourth quarters of the year (July to December). This temporal behavior is associated with the change from the wet to the dry season [95], and with meteorological factors (high temperatures, low rainfall rates, and low relative air humidity) and climate change, which together create favorable conditions for wildfires occurrence and drastically impact the local ecosystem's processes [41,50,96]. All these considerations demonstrate the importance of assessing and taking into account both human and environmental driving factors of wildfires for proper landscape planning and management [76].
The maps of wildfires densities allow the visualization of priority areas to support decision making and resources allocation processes, with the aim to reduce the risk of wildfires and to mitigate their environmental, social, and economic impacts. Moreover, they are useful inputs for landscape planning and management at local, regional, and global scales [76]. Alternately, unraveling the ignition sources of wildfires is also an important task, considering that they are greatly dependent on the relief, natural combustion components (vegetation, soil moisture), and meteorological characteristics [97]. The density of and the distance to roads favor the appearance of susceptible zones to accidental fires, i.e., those unintentionally caused by people, since they enable human activities and displacement [97,98]. The human presence is a major factor that influences the occurrence of wildfires, particularly near populated areas [99]. Unemployment tax and hard living conditions are important socioeconomic variables that lead to the overexploitation of natural resources, which can imminently cause accidental fires [97,98]. Such fires can also be a direct consequence of inappropriate cultural practices of the locals. The terrain slope, in its turn, favors the vertical continuity of vegetation, and hence, the emergence of hillside winds that ease the propagation of wildfires [100,101].
In the Peruvian Amazon, this type of investigation at a national scale is of strategic importance, in face of the fact that the wildfires occur as a direct consequence of agricultural activities (cropping and cattle raising) undertaken in small rural properties (<0.05 km 2 ) [102]. In this sense, this study contributed with a decision rules-based routine for detecting and monitoring vegetation cover areas (including forests, bushlands, grassy shrublands, and pastures) affected by wildfires by means of remote sensing techniques that allowed the quantification of BA in the Amazonas department. The generated information can be used by authorities and decision makers in the design of policies and strategies targeted to prevent, control, and/or mitigate the risk of wildfires. The proposed routine can be further applied to the whole Peruvian Amazon, provided that minor adjustments to the methodology are duly made. Finally, we ought to remind that if one intends to assess the historical series of BA in the last four decades (1984-2020) with free satellite imagery, the data from Landsat-5, 7, and 8 can be used. If one otherwise wishes to identify more recent BA and with a greater level of detail, the data from Sentinel-2, launched in June 2015, are appropriate.
Among the limitations of this study, the first one to be cited is the limited availability of cloud-free Landsat-8 and Sentinel-2 images, due to the presence of clouds throughout the year in the study area. This is the reason why we applied the CFMASK algorithm to the Landsat-8 images and CloudScore and the TDOM to Sentinel-2 images. A possible alternative to overcome the reduced availability of images would be the combination of the two collections. Another solution would be the use of the Sentinel-1 mission radar, which is equipped with a Synthetic Aperture Radar (SAR) that collects data regardless of the weather conditions [103][104][105] and improves the detection of vegetation change in the presence of high cloud cover. Additionally, a further alternative for cloud and cloud shadow detection and masking in Sentinel-2 images would be the use of machine learning [106] and deep learning approaches [107]. The second limitation regards the methodological challenges for discriminating wildfires from fires in rice crops adopted as a rice straw management practice. This differentiation requires great expertise from the interpreter conducting the supervised learning. In the case of collaborative mapping in well-known areas, such discrimination is still possible, but it is seriously limited in the classification of unknown areas.

Conclusions
This study evaluated the potential of a decision rules-based routine relying on Landsat-8 and Sentinel-2 medium spatial resolution imagery to estimate BA in the northeastern Peruvian Amazon.
These estimates were executed through training and supervised classification in the cloud-based GEE platform. The results obtained by the two satellites' imagery were compared in terms of accuracy metrics and level of detail (size of BA patches). This is an unprecedented initiative, given that similar products have been reported for the Peruvian Amazon, especially the Amazonas department, but none of them have used Sentinel-2 data. This study relied on different procedures and products for mapping BA, making use of spectral indices extracted from pre and post-fire images (dMIRBI, dNBR, dNBR2, and dNIR).
The BA resulting from Sentinel-2 were comparatively larger in extension (km 2 ) and presented higher accuracies than the ones derived from Landsat-8. The spectral indices MIRBI, NBR, and NBR2 and the NIR band produced the best results for identifying burnt areas in the Peruvian Amazon. The applied methodology showed promising results for mapping BA, revealing its feasibility to further support the assessment of BA environmental impacts at regional, national, and global levels as well as to subsidize measures for wildfires risk management.
Supplementary Materials: The following are available online at http://www.mdpi.com/2220-9964/9/10/564/s1, Figure S1: Maps of wildfires burn scars proximity to the road network; (a) BA resulting from Landsat-8; (b) BA resulting from Sentinel-2. Figure S2: Maps of wildfires burn scars superimposed on slope ranges; (a) BA resulting from Landsat-8; (b) BA resulting from Sentinel-2. Figure S3. Maps of wildfires burn scars proximity to populated areas; (a) BA resulting from Landsat-8; (b) BA resulting from Sentinel-2. Table S1. Scripts developed for the classification of burnt areas (BA) in the Amazonas department. Table S2. Statistical validation of burnt areas for the Landsat-8 product in 2017. Table S3. Statistical validation of burnt areas for the Landsat-8 product in 2018. Table S4. Statistical validation of burnt areas for the Landsat-8 product in 2019. Table S5. Statistical validation of burnt areas for the Sentinel-2 product in 2017. Table S6. Statistical validation of burnt areas for the Sentinel-2 product in 2018. Table S7. Statistical validation of burnt areas for the Sentinel-2 product in 2019. Table S8. Analyses of the burnt areas' spatial distribution with respect to slope and proximity to roads and populated areas for the Landsat-8 and Sentinel-2 products.