EXPLORING NASA’S HARMONIZED LANDSAT AND SENTINEL-2 (HLS) DATASET TO MONITOR DEFORESTATION IN THE AMAZON RAINFOREST

Deforestation is a threat to biodiversity and the world’s climate. As agriculture and mining areas grow, forest loss becomes unbearable for the environment. Consequently, monitoring deforestation is crucial for decision makers to create polices. The most reliable deforestation data about the Amazon forest is generated by the Brazil’s National Institute for Space Research (INPE) through its PRODES project. This effort is labor and time intensive because it depends on visual interpretation from experts. Additionally, frequent Amazon’s atmospheric phenomena, such as clouds, difficult image analysis which induces alternative approaches such as time series analysis. One way to increase the number of images of an area consists of using images from different satellites. NASA provides the Harmonized Landsat and Sentinel-2 (HLS) dataset solving spectral dissimilarities of satellite sensors. In this paper, the possibilities of HLS for forest monitoring are explored by applying two deforestation detection methods, Break Detection for Additive Season and Trend (BFAST) monitor and Random Forest, over four different vegetation indices, NDVI, EVI, GEMI and SAVI. The SAVI index used as input for BFAST monitor performed the best in this data setup with 95.23% for deforested pixel, 53.69% for non-deforested pixels. Although the HLS data is described as analysis ready, further pre-processing can enhance the outcome of the analysis. Especially, since the cloud and cirrus cover in the Amazon causes gaps in the dataset, a best pixel method is recommended to create patched images and thus a continuous time series as input for any land cover and land use classification.


INTRODUCTION
The Amazon rainforest is the largest rainforest in the world. It is not only one of the largest and more diverse ecosystems but also an important player in the South American hydrological cycle (Lovejoy and Nobre, 2019). The Amazon returns up to 75% of its rainfall to the atmosphere. However, deforestation at large scale can transform it into a tropical Savannah, producing a general loss of moisture and up to 50% of water run-off (Lovejoy and Nobre, 2019). Although, this fact is well-known in the scientific literature, the deforestation rate in 2018-2019 increased by 30% according to data collected by Brazil's National Institute for Space Research (INPE) (INPE, 2019a), and this current trend is to keep increasing during the future (Escobar, 2019).
The same Amazon characteristics of extent, diversity and abundant rain make it a difficult area for monitoring. Cloud coverage is a significant problem when monitoring the forest through satellite images, and cloud cover persist during most of both wet and dry seasons (Collow et al., 2016). Finding clear images for a specific date for the whole Amazon is inconceivable. For that reason, most analysis of Land Use and Land Cover change over the Amazon forest use images of a few dates of the dry season into what is known as multi-temporal analysis. One of the issues with this kind of analysis is the potentially long timespan between the deforestation and its detection, * Corresponding author making it unsuitable for monitoring. The only viable alternative are the time series analysis of images, where the persistent cloud coverage is counteracted by increasing the frequency of the images. However, the revisiting satellite time is measured in weeks. To solve this problem, the National Aeronautics and Space Administration (NASA) has recently developed a workflow that harmonizes Landsat and Sentinel-2 images resulting in a high temporal resolution dataset, called the Harmonized Landsat and Sentinel-2 (HLS) data.
In this paper, we explore the possibilities and restrictions of using HLS to assess deforestation in the Brazilian Amazon. Therefore, we compare four different spectral indices, NDVI, SAVI, GEMI and EVI. After calculating all indices for the given time period, the images were stacked, creating time series as input for the two image classification methods: Break Detection for Additive Season and Trend (BFAST) monitor (Verbesselt et al., 2012) and Random Forest (Breiman, 2001) for Satellite Image Time Series (SITS).
In the following sections, we review scientific literature, then we list the materials and the methods we used to test the HLS dataset. Later, we discuss the outcome of our experiments and finally, we introduce some conclusions and recommendations for future works.

Background
The most important monitoring systems for the Brazilian Legal Amazon are the deforestation monitoring program (PRODES) and near real-time deforestation detection system (DETER); both of them are run by INPE. On the one hand, PRODES produces a yearly report of deforestation since 1988 by visually analysing large sets of Landsat, Sentinel, and CBERS images (Souza et al., 2019). This labor-intensive task relies on trained specialists to determine deforested areas based on three criteria: Color, texture and context. On the other hand, the DETER alert system uses images from MODIS and CBERS satellites (Souza et al., 2019). The images provided by the sensors aboard new CBERS-4A satellite will improve deforestation monitoring as they have a finer temporal resolution of two or three days. The DETER system automatically analyses images as soon as they arrive and then proceeds to immediately alert the authorities (the Brazilian Institute of Environment and Renewable Natural Resources) and to release a monthly report of deforestation and forest degradation. Whereas the PRODES project is more accurate, the DETER system is faster.
As the public awareness of deforestation rises, more research is conducted and more projects on deforestation monitoring are initiated (Grinand et al., 2013;Hamunyela et al., 2016;Brown et al., 2016;Cabral et al., 2018;Yang et al., 2019;Schultz et al., 2018). A variety of images from several earth observation projects as well as different combinations of methods and indices were tested in these studies.
For example, Yang et al. (2019) analysed deforestation and its driving factors in Myanmar's rainforest. To determine land use and land cover classifications, the authors analysed Landsat images with 30m spatial resolution from 1988 to 2017. They applied a linear decomposition method on the computed indexes Normalized Difference Soil Index (NDSI), the Normalized Difference Vegetation Index (NDVI), and the Modified Normalized Difference Water Index (MNDWI). Later, they trained a Decision Tree using a sample dataset of land use and land cover collected through visual analysis. The results were used as thresholds for each land use and land cover category. They evaluated the accuracy of their results through a hybrid matrix precision authentication. After successfully extracting the raster values, they analysed the land cover classification and conducted a time change analysis of the forest. In this way, they successfully identified the areas of forest cover reduction and increase (average accuracy of ca. 88%) and the reasons of these change (e.g., cropland expansion).
Another example is Grinand et al. (2013), who identified deforestation in the Madagascar wet forest using Landsat/TM imagery. They calculated the NDVI and the Normalized Infrared Index (NIRI). They did not remove the clouds from the images but rather developed different strategies to handle no-data values. For example, if a pixel was forested at date one and three but cloudy at date two, it was considered as forested throughout the whole time period. Hence, using a cloud-mask was only necessary when a cloud or a shadow was observed at date two, while forest was observed at date one and deforestation at date three. To determine the deforested areas, they applied Random Forest. Their validation revealed that 84.7% was stable land cover 60.7% was land cover change.
Deforestation in the Amazon was researched, among others, by Cabral et al. (2018). They used the Global Forest Change (HD) (Hansen et al., 2013) and PRODES datasets to determine each type of protected area and to sample outside these areas to create forest loss maps from 2002 to 2016. They established a forest fragmentation model using ArcGIS Patch Analyst and Vector-based Landscape Analysis Tools (V-LATE); later, they evaluated the levels of fragmentation using Principal Component Analysis (Cabral et al., 2018). The outcome was an assessment of the deforestation dynamics in the Brazilian Legal Amazon.
Since the HLS dataset is fairly new, only a limited amount of research has been conducted with it (Lei et al., 2018;Pastick et al., 2018;Franch et al., 2019;Zhou et al., 2019;Griffiths et al., 2019;Torbick et al., 2018), especially in the context of Land Use and Land Cover change. Zhou et al. (2019) investigated if the HLS dataset improves the accuracy of the assessment of grassland dynamics compared to the Landsat-8 dataset. Although running into difficulties like varying observation frequency and spatial inconsistencies, the HLS dataset consistently lead to better results than using Landsat-8 data . Others were also able to enhance their frameworks. Pastick et al. (2018) were able to create a continuous time series with HLS data by filling in gaps with additional satellite images. Thus, their modelling framework was better able to monitor land-surface dynamics.

MATERIAL AND METHODS
Two case studies were conducted using the same dataset. The two different approaches and the different use of data in these methods permits making concrete statements about the use of the data and whether or not it is Analysis-Ready Data (ARD) .

Study Area
The Amazon forest is the largest rainforest in the world and spreads across several countries with its biggest area located in north-west Brazil. It is home to over 14,000 species, of which 6,727 are trees (Cardoso et al., 2017). The Mato Grosso State is situated in the Brazilian Legal Amazon region and is the thirdlargest state of Brazil. The study area lies south-west of the The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B3-2020, 2020XXIV ISPRS Congress (2020 capital Cuiabá next to the Xingu Indigenous Park, see Figure 1. From 1988 to 2019, around 33% of the forest has been cut in Mato Grosso State (TerraBrasilis, 2019), recently getting closer to the border of the park.

Orbital data
As deforestation often does not occur on a large scale at one single location, high spatial resolution data is needed to accurately detect it (Kalamandeen et al., 2018). Additionally, high temporal resolution helps when dealing with areas with cloud coverage because it increases the chances of acquiring cloudfree pixels. Thus, in this study, data from the HLS project is used.
The HLS project produces long-term surface reflectance data with high temporal and spatial resolutions so users can explore each pixel through time . Hereto, images from the European Unions Copernicus project from the Sentinel-2 satellite are harmonized with images from the Landsat program and provided through their website. The resulting images have a spatial resolution of 30m, are atmospherically corrected  and have a temporal resolution of two to three days. Additionally, this data provides a quality analysis layer where each pixel is evaluated and characterized by the presence of cloud, cirrus, cloud-shadow, water or with a aerosol quality . For the case study, the tiles 21LYG and 21LYH were used.

Training data
The training data for the Random Forest algorithm belongs to four categories: Agriculture, Pasture, Forest, Deforestation. For the first three categories, samples were collected visually from high resolution satellite images by specialists being careful that the land cover was not changing over the monitoring period. The deforestation samples were collected from the PRODES dataset, where the centroids of the polygons were calculated to avoid border confusion. These points were labeled according to the years that the deforestation happened. For example, in 2018, a total of 266 samples were used for tile 21LYG and 428 for tile 21LYH. With this point layer and the time series for the concerning PRODES year, a matrix with training values was created.

Validation data
The most accurate data concerning deforestation in Brazil can be acquired from the TerraBrasilis platform (INPE, 2019b). Thus, to validate the results of the analysis, the shapefile "Yearly deforestation increments -Shapefile (2008/2018)" was downloaded. It shows the yearly deforestation polygons from 2008 until 2018 with sublayers for each year and was produced by the PRODES project.

Image Pre-processing and calculation of indices
Trying to avoid the problems of different data frequency as described by Zhou et al. (2019), only one image per month is selected for creating a time series. To pick the best suited image for ever month from 2013 to 2018, the images were ranked based on their cloud cover and spatial coverage of the tile. The best ranked image, meaning the least cloud coverage and the highest spatial coverage, was then further processed for calculations. The first processing step was extracting cloud masks from the quality assessment layer to, then, replace the cloud pixels with a no-data value. Secondly, indices for the analysis were calculated with these cloud free images. The indices tested were selected based on the performance evaluation from Schultz et al. (2016a). Accordingly, the normalized difference vegetation index (NDVI) a measurement to evaluate photosynthetic activity Tucker (1979), the enhanced vegetation index (EVI) adjusting the vegetation index in high biomass regions (Huete et al., 1994), the global environment monitoring index (GEMI) enabling large-scale monitoring (Pinty, Verstraete), and the soil-adjusted vegetation index (SAVI) describing dynamic soil-vegetation systems (Huete, 1988). Whilst calculating, all values were compressed to the indices ranges if not a no data value to reduce later confusions when running the algorithms. Lastly, the index images were stacked in temporal order to create a time series as input for the BFAST monitor and the Random Forest algorithms.

Deforestation detection using BFAST monitor
Break Detection for Additive Season and Trend (BFAST) monitor is an algorithm that enables near real-time disturbance monitoring and includes a BFAST seasonal model approach (Verbesselt et al., 2012). It makes use of a history monitoring period to create a stable model, thus having the ability to disregard gaps int the time-series. With this model disturbances in a time-series can be detected during the monitoring period and in newly arrived data. The most important advantages for this study are that it can deal with gaps in the time series and that it takes the full details of the whole time series into account.
BFAST monitor version 0.3 in Python (Gieseke et al., 2019) was used to detect breaks and in the created time series. It is acknowledged that BFAST and BFAST monitor are meant to detect breaks in the established behavior of time series. However, in this case it was used as a satellite image classifier taking into account that PRODES provides masks of primary forest, that is, forest that remains intact since the beginning of the PRODES project. Any change to that forest is considered as deforesta- The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) tion. Even if the forest is removed and it grows again is considered as secondary forest. Under this premise, any change to primary forest is deforestation. Since BFAST is able to detect changes in time series of vegetation indices, it can be assumed that any disturbance to the regular phenological cycle of the PRODES' primary forest is deforestation. Under the aforementioned assumptions, BFAST can be used as a satellite image classifier instead of its traditional use as a time series break detector.
For the initial runs, monthly data from 2016 until 2018 was used. For these runs, the history period was set from January 2016 to December 2017 and the monitoring period, consequently, set from January 2018 until December 2018. The results of these runs were not meaningful since the algorithm detected breaks in nearly all pixels. Therefore, different combinations of parameters have been tried which, neither, led to meaningful results. So, data from 2013 until 2015 was added and the history period was changed to January 2013 until July 2016 and, consequently, the monitoring period from August 2016 to July 2018 to fit the PRODES year (from August to July). After several cycles with adjusted parameters and no significant improvement, the results of the run with the default variables, k = 3, trend = True, hfrac = 0.25 and level = 0.05 are evaluated.

Deforestation detection using Random Forest for SITS classification
Random Forest is a machine learning algorithm, that randomly selects a set of training data observations and variables and creates decision trees (Breiman, 2001). With aggregation and applying a majority vote rule, the final class for each pixel is established.
For this study, the Random Forest algorithm of the sklearn Python library, version 0.21.3, was used. With the prepared training matrix, described in section 2.3, and an array containing the labels for each pixel a Random Forest classifier was created with 500 trees, default, and using out-of-bag samples to estimate the generalization accuracy. After fitting it to the training data the time series was classified.

RESULTS AND DISCUSSION
The results were evaluated in three steps. Firstly, the deforestation detection layer were visually compared to the PRODES layer. Secondly, the percentage of correctly classified pixels were calculated and lastly the Kappa index was determined.

BFAST monitor
The resulting break detection layer was firstly visually analysed and compared with the results of PRODES from the concerning years in a validation process. As seen in Figure 2 and Figure  3, detected changes intersect with the results from PRODES. Unfortunately, the algorithm not only results in abrupt breaks of forest land cover and land use but in a lot of noise and also false positives.
Additionally, the algorithm performed differently with the different indices. Using the NDVI index as input, barely any breaks were detected, see Figure 2 top left. With the EVI index, BFAST monitor identified breaks nearly everywhere even in the regions where forest covers the land throughout the whole monitoring period. The GEMI and SAVI indices detected, additional to deforestation, changes in most areas that are agricultural used.
Comparing the detected deforestation from PRODES to the detected breaks by calculating the Kappa index and the percentage of correctly classified pixel, one can find additional evidence on how differently the indices performed. Only 0.34% of the deforestation pixels were classified correctly as deforestation when using the NDVI index, as seen in Table 1. In comparison to the NDVI the EVI index performed opposite with a high accuracy in deforested pixels (97.84%) but low accuracy in other pixel (5.07%). Comparing the Kappa index, it shows little difference, with 0.0005 for NDVI and 0.0007 for EVI. Noticeable is, that the GEMI index detects deforestation more accurate than the NDVI and EVI, which contrasts the results of Schultz et al. (2016). The GEMI index showed 80.15% accuracy for deforested pixel, 60.16 % for non-deforested pixel and a Kappa index of 0.0236. But, agreeing with Schultz et al.(2016), the SAVI index resulted in the highest accuracy levels with 95.23% for deforested pixel, 53.69 % for non deforested pixels and a Kappa index of 0.0247. The GEMI index turned out to be a less accurate .
With the SAVI and the GEMI, indices one can find correlations with the temporal extent of their detection. As seen in the detail view of Figure 3, it is also noticeable that the detection time of the BFAST monitor method does not always intersect with the detecting time of PRODES. The difference between the PRODES deforestation layer and BFAST monitor detected breaks is also seen in assessing the correctness of the classified pixel, see Table 1. Thus, detecting deforestation was not successful with these settings and indices.

Random Forest
Running the Random Forest Classification with the same data than BFAST monitor results in less noise and false positives. As seen in Figure 4, the four indices still performed differently. For the analysis the Agriculture, Pasture and Forest classes are summarized as non-deforested. A visual comparison on the whole classified image rules out the GEMI index due to too much noise and foggy outlines. This conjecture gets confirmed when calculating the accuracy. For deforested pixels, only 3.58% were detected. A thorough analysis of the distinct feature showed in the detailed view in Figure 4 leads to the belief that the NDVI index performs the best in combination with the Random Forest algorithm in this area. Assessing the accuracy statistically, it becomes apparent, that the SAVI index created slightly better results than the NDVI index, with a Kappa index of 0.2868 for SAVI and 0.2020 for NDVI, see Table 1 .
Comparing the features detected with Random Forest with the PRODES dataset, the outlines of the deforestation feature are clearer and the whole path of deforestation is detected instead of just pieces of the feature when using NDVI and the SAVI index, see Figures 5 and 6. It is still clear, that depending on the area the indices perform differently. After calculating how many pixel were classified correctly, it shows that the SAVI index accomplished a better outcome, with around 22% for deforested pixels and 98% of non deforested pixels, see Table 1. Although, the SAVI index detects parts of the deforestation correctly, it still fails to detect a great amount of deforested pixels.   The overall outcome of this study contests the result described in other papers, for example by Torbick et al. (2018). In their paper dealing with crop type mapping using HLS and Random Forest they obtained an overall accuracy of 85% in all but one of their study areas with the NDVI and land surface water index (LSWI) indices.

CONCLUSIONS
The HLS project is a fairly new approach to provide a dataset with high temporal and spatial resolution as ARD. The data was tested with the aim to detect deforestation, using two different methods, BFAST monitor and Random Forest algoritm with four vegetation indices NDVI, EVI, GEMI and SAVI.
Masking out the clouds based on the quality assessment layer turned out to be rather inaccurate, which was also described by Pastick et al. (2018). A visual analysis of a random selection of cloud-free images showed some fuzzy edges, where either too much cloud or too little was detected. On other images whole clouds were not detected at all. This inaccuracy might have led to inaccuracies in the later analysis.
In the conducted studies, BFAST monitor and Random Forest were tested under the same conditions. At first glance, the visual analysis showed that the results of the Random Forest algorithm were more precise than the results from BFAST monitor. When calculating the fractions of the correctly calculated deforested pixels, it becomes apparent, that overall the SAVI index lead to the best results in this study area with the HLS dataset with both methods. The NDVI, which resulted in second best results with Random Forest, performed the worst with BFAST monitor. Besides, the GEMI index generated better results using BFAST monitor than NDVI and EVI but produced the most noise using Random Forest. The EVI index resulted in low performance in both algorithms although seeing overall better outcomes with the Random Forest algorithm.
Still, the results of the study using BFAST monitor are overall unsatisfying. This might be due to the holes in the time-series due to high cloud coverage especially during the months of the rain season from December to May, although Verbesselt et al. (2012) describe that data gaps will not effect the outcome. In a newer paper Schultz et al. (2016b) describe the impact of observation frequency when using BFAST monitor. They conclude that data availability is highly correlated to error omission. Additionally, the inaccuracy might also be because of other atmospheric disturbance. The original study of BFAST monitor was conducted in Somalia, an area with semi-arid tropical climate and little seasonal change in temperature (Verbesselt et al., 2012). Thus, it is possible that the constant air humidity in the Amazon can create disturbance that is not accounted for in the BFAST monitor algorithm (Collow et al., 2016). Moreover, the seasonal change in Mato Grosso State is significantly higher than in Somalia, in function of the agricultural dynamics. The agriculture plays not as an important role in Somalia as it does in the Amazon region.

FUTURE WORK
Further investigations of the use of this dataset, especially in the Brazilian Legal Amazon, should take different aspects into account. First of all, a more precise cloud mask needs to be created thus ensuring an accurate time series. An evaluation of different cloud masks was for example done by Beatens et al. (2019). Additionally, to fill the wholes in the time series arising from masking out the clouds, a best pixel method, as for example described by Griffiths et al. (2019), should be established. Evaluating each pixel based on a variety of features, for example the spatial distance to the next cloud or the temporal distance to the target date, facilitates an image aggregation based on the best rated pixel. This results in a high quality patched image, easier processable in time series analysis.
As especially BFAST monitor detected more seasonal change not only in EVI index but also in GEMI and SAVI indices it is further to investigate if removing seasonal variance can enhance the performance of detecting breaks only in the forest region. Hence, methods like Hamunyela et al. (2016) can be applied to this dataset. They developed a method to integrate the spatial context into an early deforestation detection method. The NDVI index derived from Landsat-5/TM and Landsat-7/ETM+ was used in their study case that analyzed dry forest in Bolivia and wet forest in Brazil. They reduced seasonal variance based on their spatial context with a moving window approach. Comparing their method to a season model-based approach, using the spatial context resulted in higher accuracy and less temporal delay.