Monitoring Forest Change in the Amazon Using Multi-Temporal Remote Sensing Data and Machine Learning Classification on Google Earth Engine

: Deforestation causes diverse and profound consequences for the environment and species. Direct or indirect effects can be related to climate change, biodiversity loss, soil erosion, floods, landslides, etc. As such a significant process, timely and continuous monitoring of forest dynamics is important, to constantly follow existing policies and develop new mitigation measures. The present work had the aim of mapping and monitoring the forest change from 2000 to 2019 and of simulating the future forest development of a rainforest region located in the Pará state, Brazil. The land cover dynamics were mapped at five-year intervals based on a supervised classification model deployed on the cloud processing platform Google Earth Engine. Besides the benefits of reduced computational time, the service is coupled with a vast data catalogue providing useful access to global products, such as multispectral images of the missions Landsat five, seven, eight and Sentinel-2. The validation procedures were done through photointerpretation of high-resolution panchromatic images obtained from CBERS (China–Brazil Earth Resources Satellite). The more than satisfactory results allowed an estimation of peak deforestation rates for the period 2000– 2006; for the period 2006–2015, a significant decrease and stabilization, followed by a slight increase till 2019. Based on the derived trends a forest dynamics was simulated for the period 2019–2028, estimating a decrease in the deforestation rate. These results demonstrate that such a fusion of satellite observations, machine learning, and cloud processing, benefits the analysis of the forest dynamics and can provide useful information for the development of forest policies.


Introduction
Deforestation means a long-term reduction in trees due to natural or anthropogenic activities [1]. It occurs worldwide as a result from complex socio-economic processes including population and inhabitation growth, agricultural expansion and wood extraction in developing countries [2,3]. The economic, political, technological and cultural factors aggravate the occurrences of deforestation, as well [1]. Many severe problems are caused by deforestation, including biodiversity loss, soil erosion, water cycle changes, and potential global effects [4]. Due to deforestation, the land is changing too fast, resulting in the loss of flora and fauna, which impairs ecosystem functioning [5]. One of the typical consequences of soil erosion is that the trees could be washed away through losing the anchoring with the soil [6,7]; landslides can also be triggered by soil erosion; moreover, trees absorb more water through their roots and prevent sediment runoff in cases of a flood. Vegetation plays a crucial role in the hydrological cycle regulation, helping the control of the water level content in the atmosphere. With less forest cover, less water will be returned to the soil, and the inland areas could be more prone to drought [8,9]. When tropical forests are destroyed, an enormous amount of carbon dioxide stored by the vegetation is released to the atmosphere which accelerates global warming [10,11]. After the burning of fossil fuels, deforestation is the second-largest source of carbon emissions [12], making it a top global concern nowadays. Forest change monitoring through remote sensing can be an effective and practical way to reveal the deforestation rate, the correlation between human disturbance and forest clearing, and the effectiveness of regulation-making for forest governance [13]. Moreover, rapid global change urges the understanding of temporal ecosystem dynamics [14], which-with the use of remote sensing techniques-helps to control and prevent further deterioration of forest clearing [15,16].
Most of the tracts in the Brazilian Amazon remained intact until the early 1970s [17]. The deforestation rate started to show an upwards trend since 1991, and a high annual deforestation level in 2004 (27,772 km 2 ) was recorded by Brazil's National Institute for Space Research (INPE) through a survey from 2000 to 2015. The human disturbances have accelerated the deforestation process, including land occupations, cattle ranching, agricultural expansion, fragmentation, etc. [18,19]. Brown et al. [20] conducted a spatial analysis and proved that land occupations in one municipality had a direct influence on deforestation in adjacent areas. Furthermore, cattle ranching and agriculture expansion are also significant drivers for deforestation [21]. In the research done by Aldrich et al. [22], they investigated the role of contentious social processes and the land conflicts in state Pará, and as an outcome developed a regional-scale statistical model to describe the contentious land change. Müller et al. [23] and Cabral et al. [24] concluded that most of the losses were in the continuous tropical forests and along the highway in the Brazilian Amazon. Their results showed that some protected areas located in the Mato Grosso and the Pará states had suffered more significant fragmentation, compared with other regions.
Considerable efforts have been made for estimating and monitoring the changes in Amazon rainforests, for example, the Brazilian government satellite monitoring program PRODES (Program to Calculate Deforestation in the Amazon) [25] and many individual REDD+ (Reducing emissions from deforestation and forest degradation) processes through the United Nations Programme on Reducing Emissions from Deforestation and Forest Degradation [26]. When dealing with time-series monitoring and mapping of deforested areas, Earth Observation programs are commonly used as main data sources [14]. Various methodologies based on remote sensing have been implemented to monitor forest dynamics in the Brazilian Amazon. In the study of Panta et al. [27], an initial unsupervised iso-clustering of Landsat products was applied, to identify spatial clusters, followed by field surveys to generate land use/cover training sets for sequent maximum likelihood classifications, while for the actual classification of the forest canopy, artificial neural networks were used. Shimabukuro et al. [28] mapped deforestation in the state of Mato Grosso, Brazil, using Landsat and RapidEye imagery through linear spectral mixing analysis (LSMA), where it is assumed that each pixel value is a linear combination of the reflectance. A similar study using LSMA [28] also highlighted the importance of pre-processing the data and the quality of the input, while using a knowledge-based decision tree classification method. The applied methods were time-consuming, computationally intensive, and not always appropriate for large areas. The efficiency of the abovementioned approaches for pre-processing, classifying, and analyzing massive datasets will raise problems when the objective is to map a large area through multi-year satellite observations. Three challenges should be considered for the production of multi-temporal land cover maps [29]. Due to varying weather conditions, the first one is to obtain useful images with minimal cloud cover. The second one is to establish a suitable platform for faster image processing. The third one is to satisfy a certain level of accuracy for image classification. Due to the limitation of data storage and computing speed, it is hard for conventional tools to fulfil all requirements. In addition, accessing suitable satellite products for time-series analyses can be a significant issue.
In the presented work Random Forest (RF) was implemented-a state-of-the-art machine learning algorithm-on Google Earth Engine (GEE), a cloud computing service, to evaluate nearly 20 years of forest evolution, through spaceborne multispectral imagery. Random Forest had already proved its application as a reliable and accurate classification technique [30][31][32][33], while GEE has emerged in recent years as a valuable platform where researchers from different domains can deploy their models [34][35][36]. Moreover, GEE gives access to a data catalogue with various global products, as well as to the most recent satellite observations from NASA Landsat and European Space Agency (ESA) Sentinels missions. Except for the almost near real-time data distributions, users can access imagery in a wide time frame allowing them to choose the most suitable ones for temporal analyses. As an outcome of the implementation of the Random Forest algorithm and GEE, those three challenges, previously noted in the literature, have been successfully overcome-GEE with its vast data catalogue and computational power is reducing the time for determining the suitable data and the actual processing time, while Random Forest is ensuring high quality and reliable image classification.
As a case study, an area in the Brazilian Amazon, in the Pará state, was chosen and the forest dynamics were assessed from 2000 to 2019 at five-year steps. With the help of cloud processing, it was possible to easily determine the most suitable cloud-free satellite products and implement the needed land cover classification, which allowed the quantification of the loss and gain processes. Even though internal model validations were performed, the obtained classification outputs were validated using high-resolution satellite images through a surveying open-source software-Collect Earth. In addition to the temporal analyses, a future deforestation simulation based on past trends through the Modules for Land Use Change Evaluation (MOLUSCE) tool, was performed. The findings of this work will overcome the previously mentioned challenges and will introduce a path for local authorities and decision-makers not only to monitor forest dynamics in near real-time, but also to analyze how the trends will develop in the future.
The rest of the paper is organized as follows: in Section 2, the materials and used methodology for generating the land cover maps are presented. Then the results, including accuracy assessments, land cover maps and future deforestation simulations will be displayed in Section 3. After, all the experimental results with a reference to previous research and relevant regulations will be discussed in Section 4. In the end, the framework, the results and the highlights of the study are presented in Section 5.

Study Area
The Brazilian Amazon is one of the most affected areas in the world from deforestation activities. As such, an area of interest (AOI) was defined in the southern Pará state, Brazil. In particular, an area was chosen that is located between Araweté Igarapé Ipixuna and Kayapó Indigenous Territory (see Figure 1), which covers an area of 48,838 km 2 . The final decision to perform the research on this area was supported by the evident variation of the rainforest highlighted in the academic literature [37,38]. The region close to São Félix do Xingu, the main city in the AOI, appears as a deforestation hot spot as it has suffered from pasture expansion, cattle ranching, road construction, land occupations and agrarian development [39][40][41]. According to the data collected by INPE [42], in recent years, the highest level of accumulated deforestation occurred in the Pará state (152,475 km 2 ), Mato Grosso (146,159 km 2 ), Rondônia (61,677 km 2 ) and Amazonas (26,972 km 2 ). Among these regions, because of deforestation and extension, the Pará state appeared to be the most suitable for our research purposes. Therefore, the forest change in the AOI for the period 2000-2019 was modelled, and the possible forest evolution for the years after was simulated.

Data Used
To perform such an extensive multi-temporal analysis, a large amount of spatiotemporal data was needed. For mapping land cover changes multispectral satellite data are extremely valuable especially when it is coupled with medium to high spatial and time resolutions. Another important aspect is that relatively new earth observation missions provide open and free datasets, such as NASA's Landsat and the ESA's Sentinels. Naturally, the products from those missions were first exploited for their availability and suitability for the current task. Due to the predefined timeframe for monitoring forest cover change, Sentinel-2 data (operational from 06/2015) could only be used for the last analyzed period (2019). The Landsat has more past missions that eventually managed to provide suitable datasets (see Table 1). Landsat 5, 7 and 8 have a 16-day cycle at the equator and Sentinel-2 A/B has a 5-day revisit time at the equator. The multispectral data collected from Landsat and Sentinel were used for applying a land cover classification algorithm and for determining the forest change. While Landsat and Sentinel provide sufficient spatial resolution for classification problems, a higher resolution was needed to validate the results. For that purpose, the China-Brazilian Earth Resources Satellite (CBERS) high-resolution panchromatic data for the years 2010, 2015 and 2019 were used, while for the rest of the cases it was not possible to obtain free and open high-resolution images with low cloud cover.

Classification Method
In this study, supervised classification using Random Forest [43] has been applied for binary (forest and non-forest) land cover classification. The method is widely discussed in the academic literature and has been proven as a suitable method for the application as a land cover classification strategy using medium and high resolution satellite data [32,44,45]. The supervised approach for pixel labelling requires the user to select representative training data for each of the predefined number of classes [46]. The training data should be representative and homogeneous for the classes to be modelled [47] since the quality of the output directly depends on the quality of the training samples. Random Forest is constructed by many decision trees systematically operating as an ensemble [48]. Each decision tree is generated by sampling a random vector independently from a training set, then a classifying vote is computed, and finally, the most popular class, among the trees, becomes a classification model [43]. Breiman [43] highlights that the Random Forest algorithm has several advantages compared with other supervised classifiers: the accuracy level of Random Forest is high and can be compared to that obtained with Adaboost, and sometimes it is even better; it is also faster than bagging and boosting; it is relatively powerful in terms of outliers and noise, and it can estimate the intrinsic error and the importance of the variable.
The whole Random Forest classification process was entirely deployed on the GEE-from creating training data, image classification and evaluating model performance. The definition of training data was performed through visual interpretation of the satellite images, where was implemented stratified sampling. The definition of both land cover classes was straightforward since the two classes are easily distinguishable. In a few cases, where some doubt could occur for the actual class the Normalized Difference Vegetation Index was used to support the final decision, or simply the current area was omitted and another better representative was chosen. The training sets were randomly divided into two partitions-one used for the model training (80%) and one used to assess the performance (20%). It should be noted the performance assessment done on GEE was considered just as initial step, where a validation of the classification results was done using high resolution datasets (more in details in Section 2.5).

Processing Platform
GEE was utilized as a cloud processing platform. The platform combines vast amount of satellite datasets and planetary-scale computational capabilities, moreover it is accessible for free (upon registration) for scientists and non-profit users. Users can run geospatial analysis and process satellite imagery or other geospatial data from the database on the cloud [29]. The GEE Engine Code Editor is based on JavaScript application programming interface (API), which can be used for writing and running scripts for complex geospatial analysis [49]. GEE has a large database of geospatial data obtained from numerous institutes and satellites that is available for all the users [50], therefore, there is no need to download the enormous data which is time and storage consuming. Such a platform makes the use of desktop computing resources impractical to exploit and analyze multi-temporal datasets covering large areas [51]. Instead, users can access the public data and deploy their processing model on the cloud, which makes such analyses much more convenient. Nevertheless, GEE gives high flexibility when dealing with issues such as cloud-cover, having a strict time range and other issues that can occur during the selection of images and processing them. As a cloud-based platform for planetary-scale geospatial analyses with massive computational capabilities, GEE can be a valuable tool to analyze a variety of high-impact societal issues such as monitoring forest change, droughts, hydrogeological disasters, water management, climate monitoring and determining related environmental protection measures [52].

Validation
For validating the obtained classification maps, a photointerpretation of high-resolution panchromatic satellite images (HiRes) using the open-source software Collect Earth was performed. Collect Earth is a tool that enables data collection through Google Earth, Bing Maps and GEE. Users can analyze high-resolution satellite imagery on the platform for a variety of purposes including validation of existing maps, quantifying deforestation, monitoring urban and land changes, etc. [53,54]. During the validation phase, the classified land cover maps are verified against the reference high-resolution images to obtain the accuracy level of the resulted maps. An equal set of points was generated for both classes (forest and non-forest) as validation samples. This was followed by crossvalidating each point of the samples with the HiRes reference. A sample can be considered as a subset of area, following a procedure to select its locations, also known as a sample design [55,56]. A required accuracy level and priority design criteria need to be considered for each sample design [57]. Stehman [58] summarized the basic advantages and disadvantages of the basic sampling methods. Stratified land cover random sampling design has advantages when considering the generation of probability samples, the practical level, the precise estimate of class-specific accuracy, the ability to estimate standard errors and the flexibility in changing sample size. In the current study, applied stratified random sampling was used, since it can ensure a precise estimate for each stratum by specifying a sample size [59]. In the current work, an approach following the one by [55] was adopted, where simple random samples of pixels were collected inside each stratum.
The required sample size was estimated according to the formula introduced by Stadelman [60]: where n is the sample size, p is the required accuracy, = 1 − , E means the allowable error which indicates the error range we can make with a sampling strategy during the validation phase, and Z is derived from the standard normal distribution. In the end, 1067 points (534 for the class forest, 533 for the class non-forest) were validated with the survey capabilities of Collect Earth. The distribution of the validation sample for the class non-forest in 2019 is given as an example in Figure 2.

Simulating Forest Evolution
After obtaining past and present forest changes for the AOI, the future state of the land cover was analyzed, based on the previously derived maps. For this purpose, a free and open-source plugin for QGIS called Modules for Land Use Change Evaluation (MOLUSCE) was used. It is designed to incorporate already proven models and to analyze and model potential land use and forest transitions [61]. The internal workflow is straightforward passing through the six modules-inputs, evaluating correlation, area changes, transition potential modelling (TPM), cellular automata simulation (CAS), and validation. After defining two land cover maps of the AOI for two different periods, the correlation between the variables was derived, followed by the computation of transition probabilities and the land cover change map. The TPM module based on the probabilities and artificial neural network (ANN) to compute the potential transition maps. In the following phase, cellular automata simulation based on Monte Carlo was used to generate simulated land cover maps. The cellular automata approach first divides the space into a two-dimensional grid of cells, then evolves the initial class distribution into discrete time-steps according to the predefined evolution rules, i.e., if a cell is surrounded by forest, it will be covered by forest in the next step with X% probability [62]. After n iterations, the forest distribution is obtained for the target year. The plugin has an implemented function to validate and compare the simulation result with a reference map. After obtaining past and present forest changes for the AOI, it was decided to analyze the future state of the land cover, based on the previously derived maps. For this purpose, a free and open-source plugin for QGIS called Modules for Land Use Change Evaluation (MOLUSCE) was used. It is designed to incorporate already proven models and to analyze and model potential land use and forest transitions [61]. The internal workflow is straightforward passes through six modules-inputs, evaluating correlation, area changes, transition potential modelling (TPM), cellular automata simulation (CAS), and validation. After defining two land cover maps of the AOI for two different periods, the correlation between the variables was derived, followed by the computation of transition probabilities and the land cover change map. The TPM module based on the probabilities and artificial neural network (ANN) to compute the potential transition maps. In the following phase, cellular automata simulation based on Monte Carlo was used to generate simulated land cover maps. The cellular automata approach first divides the space into a two-dimensional grid of cells, then evolves the initial class distribution into discrete time-steps according to the predefined evolution rules, i.e., if a cell is surrounded by forest, it will be covered by forest in the next step with X% probability [62]. After n iterations, the forest distribution is obtained for the target year. The plugin has an implemented function to validate and compare the simulation result with a reference map. Мore details about the MOLUSCE's inner processing can be found in [61] and the developers' guide [63].
While the final goal of the simulation was to obtain future forest evolution based on the obtained classification maps for the period 2010-2019, a trial was carried for the period 2006-2010 because it could have been verified with the land cover classification of 2015. A few iterations were performed to calibrate and validate the suitability of the algorithm for the current task. As such, the most suitable spatial variables were implemented for land cover prediction from 2019 to 2028, based on the land cover maps for 2010 and 2019. For the final simulation, spatial variables were used for the urban areas, road network and previous forest loss from the period 2006-2010. While the final variables were the same as the trials, the computation of transition probabilities and the ANN trainings were carried out for each iteration. The 2006-2010 loss was derived from the difference between the two classification maps. The road network was obtained from OpenStreetMap [64]. While the urban areas were downloaded from MapBiomas [65,66]. It was considered that possible future deforestation will occur around areas that have been already affected, close to urban areas and close to the road network which gives access to the forest.
It should be noted that due to the different spatial resolution for both observation periods, 2010-30 m/pix and for 2019-10 m/pix, the latter was resampled to match the 2010 dataset using the majority resampling method.

Workflow
The process started with satellite data access and filtering of the images according to their suitability. Two filters were applied: a cloud filter leaving the observations with the least cloud cover and time filter of the observations for a pre-determined time range, from 1 st of May till the 31 st of August for all years in consideration. This time frame was chosen because it coincides with the dry season of the Amazon forest, allowing greater probability to define cloud-free mosaics. For each year training samples were created for the two land cover classes to be modelled, forest and non-forest, followed by the supervised classification using Random Forest. For validating the classified maps for the years 2010, 2015 and 2019, a photointerpretation with a reference to the high-resolution panchromatic images (HiRes) obtained from the CBERS missions was carried out. Once validated, forest change maps were computed for every two periods. Lastly, the potential forest evolution was simulated for the next nine years (until 2028) using MOLUSCE. The overall workflow is displayed in Figure 3.

Classification Results
The land cover classification based on the machine learning algorithm and satellite time-series data, was applied successively for the years 2000, 2006, 2010, 2015 and 2019. The obtained results were analyzed for monitoring deforestation patterns that occurred during the defined time frame. In Figure 4 the resulting land cover maps are presented. From an initial visual interpretation, by comparing the maps in Figure 4a,b, a massive forest change between 2000 and 2005 can be detected. From the next 5-year interval (Figure 4c) it is noted that the forest kept decreasing, but with a lower rate compared to the previous period. The forest change after 2010 is not notable by simple visual comparison, so maps referring to variances in the land cover classes for every five years were computed.

Validation Results
Even though Random Forest is considered a highly accurate, with low estimation error, classification algorithm, the output maps had to be validated. As mentioned before, the validation procedure was done using photointerpretation of HiRes CBERS imagery through the Collect Earth software. Due to the reduced availability of the free HiRes data, the validation was done only for the classified maps for the years 2010, 2015 and 2019. Overall accuracy and Cohen's Kappa were implemented as validation metrics, where the related confusion matrices are reported in Table 2. From Table 2 it can be observed that the overall accuracies for those three years are in the range 0.95-0.97. Cohen's Kappa coefficients are 0.91-0.94. Both accuracy and Kappa reached a high level, meaning that the classifications are accurate and reliable.  [67,68] The former represents the percentage of predicted forest classified correctly. While the latter, the percentage of actual forest classified correctly. The actual measure for the quality of the classification is computed by using the area under the curve of the precision-recall curve (AUCPRC), where a perfect classifier is considered as AUCPRC = 1. As it can be noted in Table 2, all three classifications are yielding AUCPRC > 0.95 which can be considered as a more than satisfactory classification.
Due to the lack of suitable reference high-resolution imagery for the years 2000 and 2006, validations were not conducted for them, except model fit estimations during the classification process, which still estimated a high level of precision. However, from the accuracy assessments done for 2010, 2015 and 2019, it can be concluded that the accuracy of the classification process, performed through GEE, can yield more than satisfactory results.

Forest Change
The forest change was computed as the difference between every two time-intervals and each is represented in Figure 5. The areas notated in black represent forest loss, and the white areas-forest gain. The sequence of these four maps revealed the constant evolution of the forested areas, and the estimates of the affected areas can be found in Table 3. The scale of the forest loss inside the AOI can be noted in Figure 5a,d by simple visual inspection. In the period 2000 to 2005, a high rate of forest loss is measured (5081.90 km 2 ) which can be estimated as nearly 10% of the total AOI. During the next five years, the deforestation rate had dropped since the forest loss was estimated as 1942.71 km 2 or 3.93% of the AOI. An even further decrease is observed during the third period, where 1779.41 km 2 of forested areas was cleared. However, from 2015 to 2019, the deforestation rate started to increase again, where an area of 2569.81 km 2 of forest loss was observed.  Table 3 A general trend can be noted from the above-mentioned results: when the forest loss reached its maximum, the forest gain dropped to a minimum level and vice versa. Figures 6 and 7 highlight the cumulative spatial and temporal mapping of forest loss and forest gain modelling from 2000 to 2019. The cumulative forest loss up to 2019 is estimated as 11,373.83 km 2 while the cumulative forest gain in the same period is 4462.79 km 2 . According to the map in Figure 6, the area most affected by deforestation is the area along the urbanization pattern of São Félix do Xingu region. However, the deforestation pattern slowed down along the Iriri river located at the south-east corner of the AOI. From Table 3, it is evident that from 2006 to 2010 the forest loss had dramatically declined by an estimated 61.77%. From 2015 can be seen that the deforestation rate has risen by 44.42%, yet, it stayed at a relatively low level. However, from the spatio-temporal variation trend of forest gain displayed in Figure 7, it should be highlighted that the most recovered areas were located in analogous regions as those that suffered from forest loss. Perhaps this phenomenon could be thanks to the regulations published by the Brazilian government to discourage the illegal cutting of forests.

Simulation Results
As mentioned before, a few iterations were performed to validate the reliability and usability of the algorithm. As a trial, a forest cover simulation was obtained, based on the trend between 2006 and 2010. The accuracies of the trial models were estimated as the percentage of correctness in the range 90.09-92.55, and Cohen's Kappa in the range of 0.65-0.72. The validation metrics were not the only criteria that had to be fulfilled. The spatial distribution of the simulated deforested areas was also examined; maps containing mainly outlying clusters were discarded since an expansion of already present hotspots was expected. Considering the period of simulation, the obtained results should output a trend for 2014 which can be compared to the actual values from 2015. The forest loss computed between the simulated result and 2010 was estimated as an area of 1564.01 km 2 which is following the trend of the actual loss between 2010 and 2015 (1779.41 km 2 ), considering that there was one-year difference between the actual and simulated result (2014). The simulation result is considered as satisfactory and it is displayed in Figure 8. However, it should be noted that this verification has some limits because there was still a lack of suitable high-resolution data to carry out validation for the current period. The final simulation was based on the trend derived from the classification maps from 2010 and 2019. The simulated forest evolution is shown by the map in Figure 9 and the forest loss computed between it and the land cover map from 2019 is 5056.13 km 2 . According to this result, the deforestation rate has increased compared to the forest loss of 3057.07 km 2 in the previous ten years. This is also coincident with the increasing trend of the deforestation rate, estimated in the period 2015-2019. From visual interpretation, the simulated deforestation pattern for 2019-2028 was not randomly distributed. The predicted deforestation areas were mostly located around existing deforested areas in south-east parts of the AOI, focused around the existing urban areas and the road network. The tendency of the future deforestation in 2019-2028 is following the trend observed in 2000-2019 along the expansion of the region of São Félix do Xingu. This region is one of the deforestation hot spots as it is a target for land speculation, cattle expansion and road construction  [39]. The simulation done with MOLUSCE for 2019-2028 can be considered as valid and reliable from both the quantitative and spatial aspects mentioned above. However, it is based only on the historical trend, depicting a possible scenario only in the case that it will not change. Some possible future interferences in the observed trend were not considered, such as possible restrictions enforced by the local government or, oppositely, effects of forest fires of different magnitude. More accurate simulations which consider external factors, as regulations or hazards, or which take into account shorter time intervals, deserve more attention and analysis however this is beyond the aim of this paper.

Discussion
As a relatively new tool, the efficiency and accuracy of the GEE platform when processing multitemporal satellite imagery on a large scale was proven. Moreover, the simultaneous access to remote sensing datasets from different satellite missions and the smooth execution of complex algorithms on the cloud makes the platform very convenient, especially in the domain of land cover monitoring [69,70].
In this work, the deforestation dynamics were mapped for an area in the Pará state in the Brazilian Amazon. Using the machine learning algorithm Random Forest, multi-temporal satellite datasets from the Landsat and Sentinel missions were classified into two land cover classes. The variation in forest cover was derived as the difference between two time-intervals. The resulting land cover maps for the years 2010, 2015 and 2019, were manually validated with the help of highresolution panchromatic images (CBERS missions) through the Collect Earth platform. The validation metrics depicted more than satisfactory results.
The tendency of forest change identified from our results ( Figure 10) coincides with the report of INPE for Brazil [42], for the same study period (2000 to 2019). The descending trend of the deforested area revealed in Figure 11 confirmed the conclusion that the forests in this region experienced the greatest loss at the beginning of the 21st century, and it lasted for about eight years until 2008. Moreover, Hansen et al. [   As for the reasons for the high deforestation rate before 2004 which is displayed in our results, Fearnside [73] concluded that 80% of the deforested areas are under cattle pasture or secondary forests abandoned or degraded as pastures. In addition, in eastern and southern Pará, multiple large cattle-ranching enterprises have been established along the Belém-Brasilia Highway and its branch roads [74]. The regional development policies published by the Brazilian government in the early stages, pushed the deforestation rate growth and the frontier expansion [75]. A more specific study was carried out by Mertens et al. [39] for São Félix do Xingu in south Pará. They have concluded that land speculation, credits and fiscal policies for livestock and crops, road construction and the investments in electrical energy are the main driving factors for deforestation in this region. An explanation for the forest gain in the results could be due to the recently applied regulations for modulating the deforestation and the forest recovery of the legal deforestation of the Amazon in Brazil. For example, the Law of Public Forest Management published in 2006 [76], the Terra Legal Program established in 2009 [77], and several individual REDD+ projects such as the Juma project in Amazonas and the Surui project in Acre [78,79]. These regulations partially influence the deforestation rates. Assunção et al. [80] outlined that the Brazilian government adopted effective policy after the deforestation rates in Amazon rainforests had peaked between 2003 and 2004. In 2004, the Brazilian Federal Government started the Action Plan for Prevention and Control of Legal Amazon Deforestation (PPCDAm). It was a governmental effort which significantly contributed to the decrease in deforestation rates; a reduction of 59% has been achieved over the 2005-2007 period. Several activities were implemented: increase in the number and coverage of protected areas, enhancement of the monitoring of the environment, establishment of more environmental enforcements, etc. [81]. In order to promote continuous and consistent deforestation reduction results, the PPCDAm has been renewed for the periods 2012-2015 and 2016-2020 to keep it up to date. Moreover, there have been other efforts that the Brazilian government has made to regulate illegal deforestation in the Amazon forest, such as, the frameworks for the management of Public Forests from 2006 [82]. These measures taken by the Brazilian government promoted the design of public policies to keep the deforestation rate at bay through actions like monitoring, control, and inspection [83]. Nevertheless, after the deforestation rate stayed at a low level for six years from 2009 to 2014, it started to rise again gradually. Part of the reason for this ascending trend was due to failures in deforestation control and the fact that stakeholders like soybean planters, cattle ranchers and timber merchants found ways to avoid agreements and legislation [84].
Based on the obtained classification results and the computed time-series forest dynamic trends, the land cover change prediction done with MOLUSCE proved as reliable from two aspects. Firstly, from the quantitative analysis: the test model which produced simulations till 2014, based on the period 2006-2010, yielded a forest loss that follows the trend margins depicted in the supervised classification for the period 2010-2015. The second aspect is that instead of being randomly distributed in the AOI, the simulated deforestation followed patterns from 2010-2019 especially due to the urban expansion around the region of São Félix do Xingu. Moreover, in the simulation, the marginalized forest cover stayed mostly intact. The estimated forest loss is highlighting a continuous increase in the deforestation rate, which is alarming for two reasons. Firstly, the simulated trend follows the deforestation dynamics that have re-initiated after 2013 (Figure 11), and secondly new regulations and tools for following them are needed.
GEE has proved as a valuable computational tool with a great data catalogue. Its use is not limited just to land cover monitoring, nor to the Pará state, rather it can be used for problems with various area sizes. Whether the implementation would be for the whole Amazon forest or even for a whole continent, the most significant limitations could be defined as data availability and programming skills.
Such a data availability issue appeared during the current work-historical high resolution free datasets were scarce and for that reason the initial five-year time-step had to be adjusted for some cases. Another improvement would have been to use higher resolution images for the classifications, since lower resolution could overestimate the forest loss/gain. This can be obtained using Sentinel-2 instead of Landsat. However, there were some limitations: Sentinel-2 is a relatively new mission and could not cover the whole period of interest; secondly, there was no free alternative to Landsat that provides such a revisit time.
Shorter time-steps (a year or two, even on the basis of a few months) could depict more detailed forest dynamics. Such a short time-step, for time-series classification or scenario prediction, would be of greater benefit to forest preservation policies which are relying on short-term results. On the other hand, it is possible that regulations that are difficult and slow to be implemented could be favored from the longer time interval. Yet, the focus of the current work was to highlight the free and open tools that allow for the performance of precise time-series analyses and to simulate future development scenarios, through the Google Earth Engine computing platform, datasets from Landsat and Sentinel, and software packages, such as MOLUSCE and Collect Earth. Such outputs can be a reference material when considering the implementation of sustainable policies and can benefit policymakers who can follow and adjust the efficiency of their regulations.

Conclusions
Time-series land cover analyses were presented to map the deforestation dynamics in the Brazilian Amazon for the period 2000-2019. A machine learning classification of multispectral satellite imagery on cloud computing service was implemented. The output classification maps were validated through photointerpretation of high-resolution images, the variation in forest loss/gain trends was accurately computed and compared to the results obtained from similar studies. Moreover, a model for simulating future dynamics of the land cover entirely based on the previously derived historical trends, was implemented. The forest dynamics agree with the previous reports from the INPE and other studies.
The experimental results show that with the implemented approach, high classification accuracies can be obtained. Moreover, cloud processing solutions such as GEE can greatly help in the acquisition of cloud-free satellite images and improving the processing time of the deployed models. Pairing a fast time-series analysis with simulating future forest change can be of great benefit to decision makers in designing and implementing proper regulating actions to protect forests.