Applications of the Google Earth Engine and Phenology-Based Threshold Classiﬁcation Method for Mapping Forest Cover and Carbon Stock Changes in Siem Reap Province, Cambodia

: Digital and scalable technologies are increasingly important for rapid and large-scale assessment and monitoring of land cover change. Until recently, little research has existed on how these technologies can be speciﬁcally applied to the monitoring of Reducing Emissions from Deforestation and Forest Degradation (REDD + ) activities. Using the Google Earth Engine (GEE) cloud computing platform, we applied the recently developed phenology-based threshold classiﬁcation method (PBTC) for detecting and mapping forest cover and carbon stock changes in Siem Reap province, Cambodia, between 1990 and 2018. The obtained PBTC maps were validated using Google Earth high resolution historical imagery and reference land cover maps by creating 3771 systematic 5 × 5 km spatial accuracy points. The overall cumulative accuracy of this study was 92.1% and its cumulative Kappa was 0.9, which are su ﬃ ciently high to apply the PBTC method to detect forest land cover change. Accordingly, we estimated the carbon stock changes over a 28-year period in accordance with the Good Practice Guidelines of the Intergovernmental Panel on Climate Change. We found that 322,694 ha of forest cover was lost in Siem Reap, representing an annual deforestation rate of 1.3% between 1990 and 2018. This loss of forest cover was responsible for carbon emissions of 143,729,440 MgCO 2 over the same period. If REDD + activities are implemented during the implementation period of the Paris Climate Agreement between 2020 and 2030, about 8,256,746 MgCO 2 of carbon emissions could be reduced, equivalent to about USD 6-115 million annually depending on chosen carbon prices. Our case study demonstrates that the GEE and PBTC method can be used to detect and monitor forest cover change and carbon stock changes in the tropics with high accuracy.


Introduction
The rapid loss of forest cover and acceleration of forest degradation in the tropics have been caused by land clearing, burning, and overexploitation [1,2]. A recent study using remote sensing technology revealed that tropical forests have been degraded drastically since 2000 [3]. Selective logging was a main driver of forest degradation, causing the degradation of approximately 400 million ha of tropical forests [4,5]. In addition to rapid forest degradation, deforestation remains high in the tropics. A recent report by the Food and Agriculture Organization of the United Nations (FAO) stated that global forests were reduced by 178 million ha from 1990 to 2000 [6]. More specifically, annual losses of tropical forests were 7.8 million ha in 1990-2000, 5.2 million ha in 2000-2010, and 4.7 million ha in 2010-2020 [6]. Such losses and forest degradation account for about 20% of total global emissions, and represent the second largest source of global emissions [7].
To reduce these losses, several major international agreements have been reached. Reducing Emissions from Deforestation and Forest Degradation (REDD+) of the United Nations Framework Convention on Climate Change (UNFCCC) is one of these agreements and is a result-based financial scheme for reducing carbon emissions or increasing removals in the tropics in exchange for financial compensation. This compensation can occur only if the REDD+ activities are monitored, measured, and verified, and the actual emissions are below the Forest Reference Emission Level (FREL) [8,9]. Therefore, methods for monitoring, reporting, and verifying (MRV) carbon emissions, reductions, or removals are critically needed so that the results from the implementation of the REDD+ activities can be reported on a regular basis. Forest inventories are the commonly used input for MRV purposes. One of these purposes is the establishment of the FREL because it is a benchmark emission against which emission reductions or removals (results) can be assessed for financial incentives [10,11]. Other field-based measurement approaches use the real weight of all parts of all trees in the targeted sample plots to determine forest biomasses using allometric equations (e.g., diameter at breast height and tree height obtained from field-based forest inventory methods) to convert inventory measurements into biomass estimates [12,13]. Although these methods provide the most accurate assessment of biomasses and carbon stocks, they are time consuming [12,[14][15][16], costly and difficult to scale [17,18].
Satellite remote sensing technologies provide an alternative approach to MRV tasks required for the REDD+ scheme in developing countries in the tropics, where the availability of field measurements is sparse [19]. In recent years, phenology-based classification methods using spectral remote sensing data and the enhanced vegetation index (EVI) have shown satisfactory accuracy [20][21][22][23][24]. Previous research [24,25] has shown promising results with high reliability and accuracy for vegetation [26], cropland [27] and land cover mapping [24,28]. The most recent studies classified land cover categories and changes using the moderate to high spatial resolution of the Landsat remote sensing data [24,29,30]. Based on moderate-solution satellite remote sensing data, other studies used basic image analysis techniques such as the maximum likelihood, artificial neural network, decision tree, support vector machine and random forest techniques to classify land cover categories [31,32] using various commercial image processing tools [33,34]. However, these methods are not without challenges such as the costs of acquiring remote sensing time series data, limitation of the spatial extent, variation of spectral properties, acquisition conditions, atmospheric perturbations, data storage, and the limit of image processing speed [35,36]. Many of these challenges can be overcome using features inherent available in the Google Earth Engine (GEE) and its phenology-based threshold classification (PBTC) methods [24]. GEE's cloud computing capability has the potential to quantify vegetation and land cover change, and assess ecosystem dynamics, in a short time at no cost [37]. Furthermore, a larger quantity of training data may not be required for image classification when the GEE PBTC method is used.
The use of PBTC methods in the remote sensing field (including GEE) has raised public and scientific attention because these methods have shown promising results with high accuracy for land cover mapping using moderate to high spatial resolution satellite data [26,38]. Previous research has attempted to study the surface vegetation phenology such as the phenology of natural vegetation, rubber plantations, and rice for land cover mapping in tropical regions [21,23,39,40]. For example, Venkatappa et al. [24] determined threshold values for 12 land cover categories by assessing the phenological characteristics of each land cover category during the mid-dry phenological season. They used GEE's harmonic model function for Landsat 5 Thematic Mapper (TM) and Landsat Operational Land Imager (OLI) EVI time series data to develop a PBTC method by applying the appropriate algorithms for classifying Landsat imagery [24]. Although a combination of GEE and phenology-based approaches can be used to estimate forest carbon stocks, emissions and removals from deforestation and forest degradation activities in the tropics, only a handful of studies exist on the applications of these approaches despite their importance for performing the MRV tasks of the REDD+ scheme [41].
This study aims to apply the PBTC method for detecting and mapping the forest cover and carbon stock changes by different land cover categories using the available Landsat 5 and Landsat 8 satellite remote sensing data in GEE over the last 28 years from 1990 to 2018. We chose Siem Reap province, Cambodia, as our case study location because of its strategic location, in which all major forest types in Cambodia are found. This study was designed to provide a review of the PBTC method and the required data and description of various steps and methods for establishing the subnational FREL and project emission level. Subsequently, carbon emission reductions and potential carbon revenues were estimated. These are discussed and policy recommendations are presented to achieve the result-based payments from the REDD+ scheme.

Study Area
Siem Reap province has 12 districts as shown in Figure 1. It ranks as one of the 10 largest provinces in Cambodia. Siem Reap is located in a tropical area from 14.8 • N to 9.9 • N, and from 102.2 • E to 107.9 • E. The elevation ranges from 6 m above sea level at Tonle Sap lake (the largest freshwater lake in Southeast Asia) to 469 m on Phnom Kulen mountain National Park. The landscape is a mosaic of dryland and edaphic forests, rice fields, shifting cultivation, and urban areas ( Figure 1). Landscape patterns and particularly vegetation phenology are influenced by inter and intra annual precipitation patterns [42] as seasonal monsoons bring wet, moisture-rich air from the southwest from May to November, whereas December to April is characterized by drier, cooler air that flows from the northeast. Most of the rainfall occurs during the wet season with an annual precipitation range from 1050 to 1800 mm [43]. The total population in the Siem Reap province increased from 896,443 in 2008 to 1,006,512 in 2019, with an annual growth rate of 1.1% [44]. Siem Reap receives a large number of tourists each year because it hosts the World Heritage Site, Angkor Archeologic Park. It hosted 2.2 million tourists in 2018, increasing from 2.1 million in 2015 [45].  [46]. The land cover map of 2018 was produced from this study.
Siem Reap was chosen for this study because it has all land cover categories except mangrove forest. The major forest types are Evergreen, Semi-Evergreen, Deciduous, Mixed Woods and Shrubs, Bamboo, and Flooded forest [24]. The provincial economy is dominated by agriculture and tourism [47]. Forest resources are important to the local poor because of their high dependency on forest ecosystem services for their livelihoods [47]. The majority of the rural residents depend on the forest products for a living [48].

Google Earth Engine Remote Sensing Data
The Google Earth Engine offers georeferenced and atmospherically corrected real-time remote sensing data with calibrated top-of-atmosphere reference (TOA) imagery of the entire Landsat satellite image collections, comprising more than 40 years of data. These collections are free of charge for education and natural and environmental research applications [37,49]. In this study, we collected Landsat 5 TOA data from 1990 to 2010 and Landsat 8 OLI TOA data from 2015 to 2018 covering the entire Siem Reap province. Table 1 illustrates the open-access Landsat datasets in GEE from 1990 to 2018 at five-year intervals.  [46]. The land cover map of 2018 was produced from this study.
Siem Reap was chosen for this study because it has all land cover categories except mangrove forest. The major forest types are Evergreen, Semi-Evergreen, Deciduous, Mixed Woods and Shrubs, Bamboo, and Flooded forest [24]. The provincial economy is dominated by agriculture and tourism [47]. Forest resources are important to the local poor because of their high dependency on forest ecosystem services for their livelihoods [47]. The majority of the rural residents depend on the forest products for a living [48].

Google Earth Engine Remote Sensing Data
The Google Earth Engine offers georeferenced and atmospherically corrected real-time remote sensing data with calibrated top-of-atmosphere reference (TOA) imagery of the entire Landsat satellite image collections, comprising more than 40 years of data. These collections are free of charge for education and natural and environmental research applications [37,49]. In this study, we collected Landsat 5 TOA data from 1990 to 2010 and Landsat 8 OLI TOA data from 2015 to 2018 covering the entire Siem Reap province. Table 1 illustrates the open-access Landsat datasets in GEE from 1990 to 2018 at five-year intervals.

Forest Land Cover Category Threshold Values
In our previous study [24], we were able to determine the major land cover category threshold values, which were derived by assessing phenological characteristics of land cover categories during the mid-dry phenological season (December to March) using the GEE harmonic model function for Landsat TM and Landsat OLI EVI time series data as described in our previous paper [24]. The phenology-based threshold values for individual land cover categories were set for the Landsat 5 (1990Landsat 5 ( -2010 and Landsat 8 (2015Landsat 8 ( -2018 median EVI data. These land cover category threshold values can be applied to the Landsat 5 and Landsat OLI EVI data to classify the forest land cover categories in tropical regions, such as Cambodia. A detailed description of EVI land cover category phenology assessment can be found in our previous paper [24]. Phenology-based threshold values were applied to seven forest categories and three non-forest categories (Croplands, Water and other land cover) in the mid-leaf-shedding phase by selecting their high-peak vegetation index values in this study [24]. The selected forest land cover category threshold values are presented in Table 2 and the JavaScript used for PBTC is available at [24].
In our previous study [24], the EVI vegetation index values of the Flooded forest were difficult to determine because they were similar to the those of the other forests, specifically Evergreen and Semi-Evergreen forest ( Table 2) during the mid-dry phenological season. To avoid this confusion, we created a boundary around the Tonle Sap lake for the Flooded forest area in GEE using the geometry tool, and then separately assigned the threshold values (TM: 0.381-0.519 and OLI: 0.382-0.581) in both Landsat TM and Landsat OLI collections to avoid misclassification and improve the map accuracy.

Phenology-based Threshold Classification
The Google Earth Engine offers a range of image processing approaches such as the compute images at-sensor radiance, TOA reflectance, cloud score and cloud-free composites [37]. We applied JavaScript algorithms to the entire Landsat collections, converted these to median values, and then applied a cloud thrash/mask function to obtain a cloud-free image. Finally, we applied a filter collection function to limit the image to those pixels within the Siem Reap province between December and March in each year. We applied the GEE PBTC algorithms for a composite EVI for forest land cover and carbon stock changes in the study region, which involved five steps ( Figure 2):

1.
Landsat TM and Landsat OLI TOA collections were accessed using an image collection function in GEE and then a filter function was applied to obtain the collections for a specific season. We used the Siem Reap province area boundary to filter the collections within the study region.

2.
We used a cloud mask function to minimize the cloud cover on the image to less than 60%, and applied reducer functions to reduce the median values per pixel [50]. The reducer functions decrease the dimensionality of image collections by calculating simple statistics, such as the median value for each pixel. The output reducer median image object (single raster layer) characterizes the quality of the complete image collection.

3.
The EVI function for the median Landsat collections were applied to entire image collections from 1990 to 2018 [24,51,52].

4.
We assigned the phenology-based threshold values for individual land cover categories by referring to our previous study [24] and then applied the PBTC function in GEE for the forest land cover classification. The resulting maps were validated using very high-resolution images (VHR) in Google Earth-Pro time-lapse.

5.
We assessed the forest cover changes from the PBTC maps and calculated the carbon stock changes and emissions by applying equations adapted from Good Practice Guidelines of the period. Finally, we established the subnational FREL (using the retrospective approach) up to 2030, an upcoming milestone under the Paris Agreement to which Cambodia is a signatory country.

Accuracy Assessment and Data for Validation
Google Earth is highly detailed digital representation of the globe and has recently been recognized for its potential to significantly improve the visualization of land use land cover change science. Higher resolution imagery from Google Earth is free of charge and can be used directly in land use land cover map validation at any scale [54]. The images are constantly updated when new data become available. Depending on the sensor, imagery resolution ranges from 30 m to 15 cm. Utilizing the time-lapse feature in Google Earth provides access to zoomable images as far back as 30 years, which are ideal for monitoring land use land cover changes [55]. Furthermore, it is powerful source of location intelligence data that can be used for investigation and preliminary studies with suitable accuracy [56,57].
Google Earth VHR time series imagery can be used to assess the accuracy of reference data. This accuracy can be enhanced by applying recognized protocols as recommended by Olofsson et al. [58]. In this study, systematic centric sampling techniques were employed to obtain the reference data for validating the PBTC maps (

Accuracy Assessment and Data for Validation
Google Earth is highly detailed digital representation of the globe and has recently been recognized for its potential to significantly improve the visualization of land use land cover change science. Higher resolution imagery from Google Earth is free of charge and can be used directly in land use land cover map validation at any scale [54]. The images are constantly updated when new data become available. Depending on the sensor, imagery resolution ranges from 30 m to 15 cm. Utilizing the time-lapse feature in Google Earth provides access to zoomable images as far back as 30 years, which are ideal for monitoring land use land cover changes [55]. Furthermore, it is powerful source of location intelligence data that can be used for investigation and preliminary studies with suitable accuracy [56,57].
Google Earth VHR time series imagery can be used to assess the accuracy of reference data. This accuracy can be enhanced by applying recognized protocols as recommended by Olofsson et al. [58]. In this study, systematic centric sampling techniques were employed to obtain the reference data for validating the PBTC maps (  Using the ArcMapFishnet tool, the systematic centric points were allocated to the respective points for accuracy assessment on the forest land cover maps within the study region. For any point in the sample, we manually assessed whether the point fell on a certain land cover category however, consideration of the size of the forest land cover category was not required. A total of 3771 reference points were spatially established, of which 2933 were for seven PBTC maps (1990-2018) and 838 were for the two reference land use land cover maps (2000 and 2005).
Subsequently, the PBTC maps were exported to Google Drive using the export function in GEE for accuracy validation. We used these maps in ArcGIS to transform land use labels into 3771 points for validation using ArcGIS spatial-join and related tools. The generated 5 × 5 km systematic centric sampling points were then converted into a "kml" file and imported into Google Earth for validating the PBTC maps.
All of the 3771 sampling points were visually analyzed using the Google Earth "time slider" tool to compare and validate the PBTC land cover classes with historic images from 1990 to 2018. This procedure enabled the verification and validation of the proximate PBTC land cover class with great accuracy and reliability because computation could be executed on a point-by-point basis ( Figure 4). Using the ArcMap Fishnet tool, the systematic centric points were allocated to the respective points for accuracy assessment on the forest land cover maps within the study region. For any point in the sample, we manually assessed whether the point fell on a certain land cover category however, consideration of the size of the forest land cover category was not required. A total of 3771 reference points were spatially established, of which 2933 were for seven PBTC maps (1990-2018) and 838 were for the two reference land use land cover maps (2000 and 2005).
Subsequently, the PBTC maps were exported to Google Drive using the export function in GEE for accuracy validation. We used these maps in ArcMap to transform land use labels into 3771 points for validation using ArcMap spatial-join and related tools. The generated 5 × 5 km systematic centric sampling points were then converted into a "kml" file and imported into Google Earth for validating the PBTC maps.
All of the 3771 sampling points were visually analyzed using the Google Earth "time slider" tool to compare and validate the PBTC land cover categories with historic images from 1990 to 2018. This procedure enabled the verification and validation of the proximate PBTC land cover category with great accuracy and reliability because computation could be executed on a point-by-point basis (  We then calculated the producer's accuracy (PA) by dividing the number of 5 × 5 km accuracy points in an individual land cover class identified accurately by the respective reference points total. The user's accuracy (UA) was computed by dividing the number of reference points in an individual land cover class identified accurately by the classified total. Finally, the overall classification accuracy (OA) was derived by dividing the total number of accurately classified land cover class by the total number of reference points, and the Kappa coefficient (K) of agreements was calculated [23,24,58,60] for the PBTC maps and reference land cover maps using a confusion matrix. We then calculated the producer's accuracy (PA) by dividing the number of 5 × 5 km accuracy points in an individual land cover category identified accurately by the respective reference points total. The user's accuracy (UA) was computed by dividing the number of reference points in an individual land cover category identified accurately by the classified total. Finally, the overall classification accuracy (OA) was derived by dividing the total number of accurately classified land cover categories by the total number of reference points, and the Kappa coefficient (K) of agreements was calculated [23,24,58,60] for the PBTC maps and reference land cover maps using a confusion matrix.

Carbon Stocks and Emission Reductions
To implement the REDD+ scheme, Cambodia employs a national definition of forest consistent with the Global Forest Resources Assessment [59], which is part of the Intergovernmental Panel on Climate Change (IPCC) "FOREST" land use category (Table 3).
Sasaki et al. [10], estimated above-ground biomass, below-ground biomass, deadwood, litter, and total carbon stocks in seven forest categories (Evergreen Forest, Semi-Evergreen, Deciduous Forest, Bamboo, Other forest, Wood Shrubland Dry and Wood Shrubland Evergreen) and estimated FREL at a provincial level using a retrospective approach. Here, we focused only on major forest categories and rubber plantations. For Mixed woods and shrubs carbon density, we used the average total carbon stocks (39.45 MgC/ha) of Wood Shrubland Dry and Wood Shrubland Evergreen. The weighted average (MgC/ha) was calculated only for Siem Reap Province using 1990 forest cover carbon stocks as a baseline (Table 3).
We used the formulas below for calculating the carbon stocks and emission reduction in seven forest categories, namely: Evergreen, Semi-evergreen, Deciduous, Mixed Woods and Shrubs, Bamboo, Flooded Forest, and Rubber Plantation. We estimated total forest area carbon stocks (CS) and summed carbon stocks from all four pools using forest categories as reported in [10]. The retrospective approach was used to develop forest cover changes based on past deforestation trends for the province by applying Equation (1) [10].
where FA i (t) is the total forest cover (ha) of each forest land category (i) at time t (year), FA i (0) is the forest cover (ha) of forest land cover category (i) at the start of the model (0) (i.e., in 1990), a i is the rate of change of each forest land cover category i between 1990 and 2018. Using data in Table 4 where TCS(t) is total above ground carbon stocks in Siem Reap province (MgC); FA i is the area of the 7 forest categories i in the province of 12 districts; CS i is the carbon stocks (MgC/ha/year) in 7 forest categories i .

Annual Carbon Emissions Due to Forest Cover Change
The Gain-Loss method was applied to obtain carbon emissions by forest category using Equation (3): where, CE (t) = annual carbon emissions (MgCO 2 /year (t)) due to deforestation in Siem Reap province. CS i is the average of carbon stocks for forest category i (MgC/ha) in Table 4. We excluded the Bamboo and Rubber plantation categories in the projected carbon emissions because their combined total forest land cover area and carbon stocks are less than 10% of total carbon stocks. Therefore, these two forest categories can be ignored as per IPCC guidelines [53]. The ratio 44/12 is the molecular weight ratio of carbon dioxide to carbon [10,53].

Forest Reference Emission Level between 2020 and 2030
The FREL is an essential benchmark for carbon emissions. The more countries agree to reduce emissions, the more financial support they receive according to REDD+ schemes for performance-based payments [61]. However, it is important to have a good understanding of the historical trend of land use and land cover change and forest deforestation at a moderate to high resolution scale because this helps to establish FREL for measuring the performance when the REDD+ activities are implemented. FREL can be established by: where FREL(t) = Forest reference emission level in Siem Reap at time t (MgCO 2 /year). Because we used the retrospective approach, FREL for Siem Reap province during the Paris Agreement (2020-2030) is the same as the projected CE(t) during the same period between 2020 and 2030.

Project Emissions
If REDD+ project is implemented to reduce the deforestation and forest degradation in the Siem Reap province, the emissions from such project can be estimated by: where PE(t) = Project emission (MgCO 2 /year); RPI(t) is relative project impact taken from Ty et al. [62].

Emission Reductions
Accordingly, carbon emission reductions can be obtained by: where ER(t) is annual emission reduction (MgCO 2 /year) at time t between 2020 and 2030 2.6.7. Carbon Revenues where: CR is carbon revenues from emission reduction (MgCO 2 /year) in Siem Reap province between 2020 and 2030; CP is the carbon price of USD 7 per MgCO 2 reported in World Bank and Ecofys [63].

PBTC Land Cover Map Accuracy
The GEE PBTC maps were shown to be more accurate according to the confusion matrix estimated in this study. The details of the forest land cover category's overall accuracies (OA), Kappa coefficients (K), producer's accuracies (PA) and user accuracies (UA) of the maps are presented in Tables 5 and 6. Appendix A (Tables A1-A7) shows the PBTC map accuracy matrix and reference map accuracy matrix.   (Table 5). The overall cumulative accuracy of this study was 92.1% and its K statistic was 0.9. Generally, the results of PBTC map accuracy are acceptable as suggested by the United States geological survey satellite imagery classification schema [64]. K statistics results greater than 0.85 represent stronger agreements between the classification made and the Google Earth imagery validation information [65]. By comparison, the overall accuracy of the reference maps was 81.38% in 2000 and 83.05% in 2005, and the K statistic for 2000 was 0.77, and that for 2005 was 0.78 (Table 6), which is acceptable.

Siem Reap Forest Cover Change
The PBTC maps between 1990, 1995, 2000, 2005 and 2010 ( Figure 5A-E) were derived from Landsat TM and Landsat OLI from 2015 and 2018 ( Figure 5F,G), and show the geographic distribution of forest category land cover in Siem Reap Province. The major types of forest are found in Siem Reap province: Evergreen forest, Semi-Evergreen forest, Deciduous forest, Mixed woods and shrubs, Flooded forest, and Bamboo. Our findings show that the Deciduous forest covers a larger area than other forest, followed by Semi-Evergreen and then Evergreen forest. The Evergreen and Semi-Evergreen are forests largely found in north and northeastern districts, protected areas and the Angkor wat temple surroundings. The Flooded forest was found adjacent to the Tonle Sap lake and the southern part of the region. Bamboo distribution was found in elevated regions in the most northern part and was mixed with Evergreen and Semi-Evergreen forests in the Banteay Srei, Varin, Svay Leu, Angkor Chum, and Chi Kraeng districts. The expansion of Rubber plantations can be seen in the Chi Kraeng and Varin districts ( Figure 5G). The mixed woods and shrubs category was distributed with Deciduous forest and close to Croplands ( Figure 5).  The Siem Reap province forest area gradually decreased from 1990 to 2000. A significant reduction can be seen from 2010 to 2015 ( Figure 6). See Appendix A, The Siem Reap province forest area gradually decreased from 1990 to 2000. A significant reduction can be seen from 2010 to 2015 ( Figure 6). See Appendix A, Table A10 for forest category land cover change in Siem Reap province from 1990 to 2018. Among the seven main forest categories, the annual change in Evergreen forest was 3467 ha. The rate of decrease was 1.8% from 1990 to 2018; Semi-Evergreen forest decreased at the rate of 1.3% and Deciduous forest decreased at 1.0% over the 28-year period. The Mixed woods and shrubs forest area showed the second highest rate of decrease compared to other forest categories (1.7%). Flooded forest areas experienced substantial losses of 61,420 ha annually at a rate of 1.6%. The province's forest area decreased by 322,694 ha in all forest types with an average annual rate decline of 1.3% over a 28-year period (1990-2018). A rapid decline in forest cover was observed between 2010 and 2018 and a slight increase in the Rubber plantation and Bamboo forest cover can be noted during the same period ( Figure 6 and Appendix A, Table A10).

Changes in Carbon Stocks by Forest Category
Using the GEE PBTC maps, we were able to calculate the forest carbon stocks (MgC) from 1990 and 2018 (Figure 7). Among the seven main forest categories, the annual change in Evergreen forest was 3467 ha. The rate of decrease was 1.8% from 1990 to 2018; Semi-Evergreen forest decreased at the rate of 1.3% and Deciduous forest decreased at 1.0% over the 28-year period. The Mixed woods and shrubs forest area showed the second highest rate of decrease compared to other forest categories (1.7%). Flooded forest areas experienced substantial losses of 61,420 ha annually at a rate of 1.6%. The province's forest area decreased by 322,694 ha in all forest types with an average annual rate decline of 1.3% over a 28-year period (1990-2018). A rapid decline in forest cover was observed between 2010 and 2018 and a slight increase in the Rubber plantation and Bamboo forest cover can be noted during the same period ( Figure 6 and Appendix A, Table A10).

Changes in Carbon Stocks by Forest Category
Using the GEE PBTC maps, we were able to calculate the forest carbon stocks (MgC) from 1990 and 2018 (Figure 7).

Forest Cover and Carbon Stock Changes and Carbon Loss
Chi Kraeng District had the most forest area in the province; however, 35% was deforested between 1990 and 2018 (81,548 ha). Svay Leu District was the second most forested district (189,826 ha in 1990) and it also lost 35% by 2018. Over the same period, Krong Siem Reap lost 38%, Angkor Chum 37%, Srei Snam 35%, Soutr Nikom 35% and Angkor Thum 30%. These 12 districts accounted for 31% of the total deforestation in Siem Reap province between 1990 and 2018. Fifty-two percent of the province remains forested, mostly in Chi Kraeng, Svay Leu, Puok and Varin Districts. Table 8 presents the district level forest cover changes during the past 28 years.  (Table 9). In all districts from 1990 to 2018, the total carbon stocks decreased by 39,198,938 MgC and were responsible for annual carbon emissions of 143,729,440 MgCO 2 over the 28-year period. Over the same period, Svay Leu District emitted 39,719,781 MgCO 2 , followed by Chi Kraeng District with 33,792,193 MgCO 2 . These figures are for four IPCC total carbon pools, including above ground, below ground, litter, and deadwood. However, estimated overall emissions vary depending on the methods used to calculate them (Table 9).

Carbon Emission, FREL and Emission Reductions
The total FREL for Siem Reap province was estimated to be 35,303,937 MgCO 2 and the project emission reduction to be 8,256,746 MgCO 2 for the 10-year project period from 2020 to 2030. If REDD+ activities are initiated in Siem Reap province to maintain the existing forest, the province has the potential to receive carbon revenues of about USD 57.8 million during the project period from 2020 to 2030 (Appendix A, Table A11). However, the carbon price (USD7/MgCO 2 ) can vary depending on international negotiations with emission trading and carbon offset schemes [63].

Discussions and Implications
Our study found that the PBTC mapping accuracy (Tables A1-A7) is higher than that of previous studies in the region (FREL, 2016) and our results show a significant decrease in forested land in Siem Reap province from 1990 to 2018 (Table A10). Although nearly half of the province is still covered in forest, the quality and value of forest ecosystems has generally declined by about 31% or 322,694 ha in 28 years (Table 8). This decline can be attributed to forest clearing to accommodate an increase in population and new households [66], expansion of Croplands [67], and the government's Economic Land Concession program to promote plantations [68]. Other major drivers of deforestation and degradation were the expansion of agriculture to grow cassava, orchard plantations, and paddy fields [16,47]. Evergreen, semi-evergreen, and deciduous forests have been subjected to logging [59]. The flooded forests and floodplains that surround the Tonle Sap fresh water lake provide shelter for fish to breed in addition to important feeding areas. The flooded forest area declined from 141,083 ha to 76,662 ha from 1990 to 2018 (Table A10). The main reason could be the use of smoke to harvest honey, burning firewood, and leaving cooking fires unattended. In addition, farmers may have burned flooded forests for conversion to rice growing, hunting animals, or setting long fishing nets across river channels [69]. Between 2010 and 2018, most districts lost forest cover as a result of rapid population growth, tourism, and the need for land for agricultural production [66]. The main factors affecting carbon stock changes were deforestation, expansion of croplands, and economic land concessions in the region [66]. Chi Kraeng and Svay Leu district forest lands have been converted to plantations, croplands, and economic land concessions, which may be the cause of high concentrations of emissions (MgCO2) in both districts. Due to forest deforestation, the province lost 39,198,938 (MgC) carbon stocks over all forest types between 1990 and 2018 or about 143,729,440 MgCO 2 of total carbon emissions over the same period (Table 9).
Previous studies showed the rate of deforestation in the province to be 2.6% between 2002 and 2006 [10], which is higher than our result. This could be due to a lack of long-term temporal data, as the previous study used activity data from only two time periods (2002 and 2006), during which forest cover sharply declined during the same period [66]. It is generally the case that the more time series activity data available, the better the resulting relationship for estimating FREL [61]. Nevertheless, both of the predicted linear trends show similar fitted values (R 2 = 0.97) of projected FREL from 2020 to 2030. As seen in Figure 8, total FRELs were estimated to be 20,431,376 MgCO 2 for 10 years from 2020 to 2030, which is higher than that projected by Sasaki et al. [10] (Figure 9). Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 34 and new households [66], expansion of Croplands [67], and the government's Economic Land Concession program to promote plantations [68]. Other major drivers of deforestation and degradation were the expansion of agriculture to grow cassava, orchard plantations, and paddy fields [16,47]. Evergreen, semi-evergreen, and deciduous forests have been subjected to logging [59]. The flooded forests and floodplains that surround the Tonle Sap fresh water lake provide shelter for fish to breed in addition to important feeding areas. The flooded forest area declined from 7,957,072 ha to 4,492,959 ha from 1990 to 2018. The main reason could be the use of smoke to harvest honey, burning firewood, and leaving cooking fires unattended. In addition, farmers may have burned flooded forests for conversion to rice growing, hunting animals, or setting long fishing nets across river channels [69]. Between 2010 and 2018, most districts lost forest cover as a result of rapid population growth, tourism, and the need for land for agricultural production [66]. The main factors affecting carbon stock changes were deforestation, expansion of croplands, and economic land concessions in the region (Table 7) [66]. Chi Kraeng and Svay Leu district forest lands have been converted to plantations, croplands, and economic land concessions, which may be the cause of high concentrations of emissions (MgCO2) in both districts. Due to forest deforestation, the province lost 39,198,938 (MgC) carbon stocks over all forest types between 1990 and 2018 or about 143,729,440 MgCO2 of total carbon emissions over the same period (Table 9). Previous studies showed the rate of deforestation in the province to be 2.6% between 2002 and 2006 [10], which is higher than our result. This could be due to a lack of long-term temporal data, as the previous study used activity data from only two time periods (2002 and 2006), during which forest cover sharply declined during the same period [66]. It is generally the case that the more time series activity data available, the better the resulting relationship for estimating FREL [61]. Nevertheless, both of the predicted linear trends show similar fitted values (R² = 0.99) of projected FREL from 2020 to 2030. As seen in Figure 8, total FRELs were estimated to be 20,431,376 MgCO2 for 10 years from 2020 to 2030, which is higher than that projected by Sasaki et al. [10] (Figure 9).  Table A12 for total Forest reference emission levels (FREL) for Siem Reap province from 1990 to 2030.  REDD+ monetary incentive payments result from the amount of accountable carbon credits and the price paid per ton of carbon. Carbon emission reduction prices rose between 2017 and 2018; for example, the European Union Allowance fee rose from USD7/MgCO2 to USD16/ MgCO2 between 2017 and 2018. However, carbon emission reduction prices range widely between initiatives, from less than USD1/MgCO2 (Mexico, Poland and Ukraine carbon tax) to a maximum USD 139/MgCO2 (Sweden carbon tax) [63]. The REDD+ scheme offers financial incentives for reducing emissions or increasing removals or carbon sequestration in forest and agricultural lands [70][71][72][73]. Carbon revenue is dependent on carbon prices. If the carbon price of USD 16 is used, the total carbon revenue is USD132.11 million. On the other hand, if the price used is the same as that of Sweden, then the total revenue would be USD1147.69 million at the time of the Paris Climate Agreement between 2020 and 2030. Because the number of sustainability travelers has increased by about 48% [74], the REDD+ project site could become a tourist destination, thereby generating more revenue for local people.
Our study map accuracy shows that the Mixed woods and shrubs category is the main source of confusion in the Deciduous forest and Semi-Evergreen forest. The Semi-Evergreen forest is the main source of confusion in the Evergreen forest category in all the classified years (1990 to 2018). Flooded forest improved the accuracy of PA and UA because the Flooded forest was mapped separately by creating the Flooded forest boundary. In contrast, the Water, Bamboo, Rubber and other classes were found to have a higher accuracy of PA and UA compared to other land cover classes, as the accuracy of samples were minimized for these classes (Appendix A, Tables A1-A7).
The validation of the reference forest cover maps found that the Evergreen and Semi-Evergreen PA and US land cover classes had higher accuracy in both 2000 and 2005 (Table 5). In contrast, the Mixed woods and shrubs, Other forest and Non-forest classes yielded low UA and PA in both assessment years. Most of the confusion regarding Non-forest, Other forest, and Mixed woods and shrubs may have resulted from the lower accuracy or pixel difference of the refence maps. The distribution of Mixed woods and shrubs in most of the study area was considered as part of the Nonforest class, whereas Bamboo, Rubber plantation, and Flooded forest were considered as part of the other forest class. This could be the reason for the confusion of the Other forest categories with the Evergreen Semi-Evergreen and Deciduous classes. However, the overall accuracy of Evergreen and Semi-Evergreen improved in both validation maps (Appendix A, Tables A8 and A9).
There are possible errors in our study due to the resolution of remote sensing data, limited IPCC land cover category threshold values and the distribution of forest and non-forest land cover REDD+ monetary incentive payments result from the amount of accountable carbon credits and the price paid per ton of carbon. Carbon emission reduction prices rose between 2017 and 2018; for example, the European Union Allowance fee rose from USD7/MgCO 2 to USD16/ MgCO 2 between 2017 and 2018. However, carbon emission reduction prices range widely between initiatives, from less than USD1/MgCO 2 (Mexico, Poland and Ukraine carbon tax) to a maximum USD 139/MgCO 2 (Sweden carbon tax) [63]. The REDD+ scheme offers financial incentives for reducing emissions or increasing removals or carbon sequestration in forest and agricultural lands [70][71][72][73]. Carbon revenue is dependent on carbon prices. If the carbon price of USD16 is used, the total carbon revenue is USD132.11 million. On the other hand, if the price used is the same as that of Sweden, then the total revenue would be USD1147.69 million at the time of the Paris Climate Agreement between 2020 and 2030. Because the number of sustainability travelers has increased by about 48% [74], the REDD+ project site could become a tourist destination, thereby generating more revenue for local people.
Our study map accuracy shows that the Mixed woods and shrubs category is the main source of confusion in the Deciduous forest and Semi-Evergreen forest. The Semi-Evergreen forest is the main source of confusion in the Evergreen forest category in all the classified years (1990 to 2018). Flooded forest improved the accuracy of PA and UA because the Flooded forest was mapped separately by creating the Flooded forest boundary. In contrast, the Water, Bamboo, Rubber and other categories were found to have a higher accuracy of PA and UA compared to other land cover categories, as the accuracy of samples were minimized for these categories (Appendix A, Tables A1-A7).
The validation of the reference forest cover maps found that the Evergreen and Semi-Evergreen PA and US land cover categories had higher accuracy in both 2000 and 2005 (Table 5). In contrast, the Mixed woods and shrubs, Other forest and Non-forest categories yielded low UA and PA in both assessment years. Most of the confusion regarding Non-forest, Other forest, and Mixed woods and shrubs may have resulted from the lower accuracy or pixel difference of the refence maps. The distribution of Mixed woods and shrubs in most of the study area was considered as part of the Non-forest category, whereas Bamboo, Rubber plantation, and Flooded forest were considered as part of the other forest category. This could be the reason for the confusion of the Other forest categories with the Evergreen Semi-Evergreen and Deciduous categories. However, the overall accuracy of Evergreen and Semi-Evergreen improved in both validation maps (Appendix A, Tables A8 and A9).
There are possible errors in our study due to the resolution of remote sensing data, limited IPCC land cover category threshold values and the distribution of forest and non-forest land cover categories. One possible source of error is limited ground reference data. Different vegetation types at various growing stages show different phenological behavior [75] and their spectral values vary during phenological stages [75,76]. More field-based studies are needed to improve the accuracy of the PBTC method for the assessment of land use, land-use change, and forestry canopy disturbance, and carbon stock changes. If carbon emissions from all land use sectors are included in the assessment, the change detection approach would be appropriate [26,41].
Since we focused mainly on IPCC forest land cover category (i.e., the seven forest categories as shown in Table 7), more categories of land cover should be added to make our methods and results consistent with the land cover categories currently used by the Cambodian government for its forest resource assessment at the national level [59]. If such categories are added, the level of map accuracy might change [59], and so the forest carbon stock would also change. Therefore, future studies will need to focus on determining the threshold values for any new cover categories.

Conclusions
Using a combination of PBTC, GEE, and Landsat 5 and Landsat 8 OLI, we were able to classify seven forest cover categories in Siem Reap province, Cambodia over a 28-year period from 1990 to 2018. These categories are bamboo, plantation (Rubber), evergreen, semi-evergreen, deciduous, mixed woods and shrubs, and flooded forests. Evergreen forest suffered rapid loss (−1.8%) while rubber plantation increased by 109% annually. On average, 11,525 ha (1.3%) of forest cover were lost annually. This loss was responsible for 14,372,944 MgCO 2 /year of the carbon emissions, of which 28% occurred in Svay Leu district alone. Using the retrospective approach, FREL for Siem Reap was estimated at 35,303,937 MgCO 2 MgCO 2 /year during the Paris Agreement between 2020 and 2030. If REDD+ activities are implemented, 825,674.6 MgCO 2 of carbon emissions could be reduced annually over the same period. Depending on chosen carbon prices, result-based carbon revenues could be USD6 to 115 million annually.
The overall cumulative accuracy of this study was 92.1% and its cumulative Kappa was 0.9, which are sufficiently high to apply the PBTC method to detect forest land cover change. With these levels of overall accuracy, it is possible to conclude that PBTC and GEE fast cloud computing can be used to assess forest cover and carbon stock changes by category at any scale, quickly and without any costs. Such technologies become increasingly important for monitoring, measuring, reporting, and verifying the performance of the REDD+ activities in developing countries as required under the UNFCCC REDD+ scheme before they can claim for the result-based payments for their efforts in emission reductions or removals.
Our study demonstrated the capability of using the GEE cloud computing platform in combination with the PBTC method, and with the use of Landsat moderate resolution (30 m) remote sensing data, for estimating forest cover and carbon stock changes in a tropical developing country with limited retrospective data or resources to acquire direct or repeated measurements of forest carbon stocks in the field. Our methods, along with GEE and increasingly available satellite imagery, make it possible to perform the MRV tasks required under the REDD+ scheme of the UNFCCC and/or to detect land cover changes anywhere, at speed and scale.

Acknowledgments:
We would like to acknowledge the Second Century Fund (C2F) of Chulalongkorn University, Thailand for their post-doctoral fellowship. We would also like to thank the LEET intelligence Co., Ltd. for research support and the FRAWASA project team for their invaluable help during field data collection.

Conflicts of Interest:
The authors declare no conflict of interest.    Note: Forest reference emission levels (FREL) are equivalent to carbon emissions starting from 2020 through 2030. Emission reductions and carbon revenues were estimated only for the 2020-2030 period.