Determination of Vegetation Thresholds for Assessing Land Use and Land Use Changes in Cambodia using the Google Earth Engine Cloud-Computing Platform

As more data and technologies become available, it is important that a simple method is developed for the assessment of land use changes because of the global need to understand the potential climate mitigation that could result from a reduction in deforestation and forest degradation in the tropics. Here, we determined the threshold values of vegetation types to classify land use categories in Cambodia through the analysis of phenological behaviors and the development of a robust phenology-based threshold classification (PBTC) method for the mapping and long-term monitoring of land cover changes. We accessed 2199 Landsat collections using Google Earth Engine (GEE) and applied the Enhanced Vegetation Index (EVI) and harmonic regression methods to identify phenological behaviors of land cover categories during the leaf-shedding phenology (LSP) and leaf-flushing phenology (LFS) seasons. We then generated 722 mean phenology EVI profiles for 12 major land cover categories and determined the threshold values for selected land cover categories in the mid-LSP season. The PBTC pixel-based classified map was validated using very high-resolution (VHR) imagery. We obtained a cumulative overall accuracy of more than 88% and a cumulative overall accuracy of the referenced forest cover of almost 85%. These high accuracy values suggest that the very first PBTC map can be useful for estimating the activity data, which are critically needed to assess land use changes and related carbon emissions under the Reducing Emissions from Deforestation and forest Degradation (REDD+) scheme. We found that GEE cloud-computing is an appropriate tool to use to access remote sensing big data at scale and at no cost.


Introduction
The development of simple methods for land use and land cover classification is essential to allow researchers and policy-makers to understand land use and land use changes, carbon emissions, and ecosystem dynamics at scale. This understanding is useful for addressing the effects of climate change. World Leaders at the 15th Conference of the Parties (COP15) of the United Nations Framework Convention on Climate Change (UNFCCC) recognized the important roles of Reducing Emissions remote sensing and vegetation phenology approaches can measure heterogeneous seasonality in vegetation, because remote sensing can provide time series data at any scale; therefore, these approaches are essential for understanding the role of vegetation in climate change [17,[37][38][39][40][41][42].
A basic remote sensing phenology approach is the use of leaf-shedding (LSP) and leaf-flushing (LFP) and time series data to determine the start of the season (SOS) and end of the season (EOS) because LSP and LFP can provide accurate vegetation phenology characteristics that are available with single image classification [31,43,44]. Until recently, the derivative-based, threshold-based, and model-fitting phenology approaches were commonly applied for extracting land surface phenology characteristics using spatial and temporal satellite-derived green proxies [45]. Previous studies have demonstrated the effectiveness of using the liner harmonic model to fit the time series remote sensing data [46][47][48][49][50][51][52]. For example, Brookes et al. [47] applied the harmonic regression model to Normalized Difference Vegetation Index (NDVI) time series profiles for Landsat Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) imagery to estimate the biophysical parameters. Vogeler et al. [52] were able to map the Minnesota forest cover harmonizing Landsat TM and ETM+ imagery from 1973-2015 with an overall accuracy of 87%. Simonetti et al. [53] used Landsat 8 operational land imager (OLI) top of atmosphere (TOA) data to map four protected areas and their 20 km buffer zones from different ecoregions in Sub-Saharan Africa using the GEE phenology-based approach. This study had an accuracy level of 90%, indicating that a single image phenology-based approach is the most appropriate method for classifying land cover categories in addition to its ability to avoid the misclassification of land use maps.
Phenology-based threshold methods have increasingly gained public and scientific attention in geospatial science. The Phenology-based classification has shown a promising level of accuracy for mapping vegetation, crops, and land cover using moderate to high spatial resolution satellite data [24,27,[54][55][56] with aid of GEE cloud computing technology [23,53,[57][58][59][60][61]. There has been some research effort to study surface phenology of vegetation, rubber plantation, and cropland mapping in tropical regions [62][63][64]. However, dry tropical regions present challenges in determining the phenology-based threshold values for land use or land cover categories because of the rapid changes in the phenological behavior of the vegetation from one season to another which make it difficult to avoid misclassifying land cover categories. Therefore, this study was designed to determine the vegetation thresholds for individual land cover categories through an analysis of the phenological behaviors of major land cover categories for specific seasons and to develop a robust phenology-based threshold classification (PBTC) approach to be used for a single EVI classification method to map the land cover categories. This could allow continuous monitoring with the aid of cloud-computing GEE classification techniques.
We applied a Landsat EVI time series analysis to identify the phenological behaviors of major land cover categories in different phenological phases in Cambodia. We used Landsat TM (from [2000][2001], Landsat ETM 7 (2002ETM 7 ( -2003, and Landsat 8 OLI imagery (from 2014-2017) to assess the phenological behaviors of the selected land cover categories and applied very first PBTC method to assess the mid-LSP season Landsat composite single EVI classification, which was then used to classify the selected land cover categories in Cambodia using the fast cloud-computing GEE platform.

Study Area
The study area is shown in Figure 1a, and ground reference data are shown in Figure 1b [65]. Cambodia is located between latitudes 10 • and 15 • North and longitudes 102 • and 108 • East.  [66]. Red color indicates high elevation and mountainous regions, while yellow and green colors indicate medium highlands and land surfaces with low elevation, respectively. The Tonle Sap lake and Mekong river are shown as water bodies. (b) Forest cover map in 2002 and geolocations of the field reference data collected in 2018. The star symbols represent the geolocations of land cover categories during field observations in February and July 2018. The dot symbols show the locations of the Permanent Sample Plots (PSPs) in evergreen (EG), semievergreen (SEG), and deciduous forests (DD). The red dots represent semievergreen forest locations where an inventory was conducted (refer to Reference [67] for details.) The photos shown on the maps were taken during the field visits and represent flooded forests (FF), PSPs, mangroves (MG), EG, rice or paddy fields (CR), DD, SEG, and rubber plantations (RB). EG, SEG, DD, and RB aerial photos were taken from the drone Phantom 4 (of DJI company) in July 2018. Reference data for the land cover were obtained from Sasaki et al., 2016 [65].
There are four environmental regions in Cambodia: The plain, Tonle Sap, coastal, and plateau and mountainous areas [68]. Cambodia has a tropical monsoon climate and there are two distinct seasons: The dry season from November-April and the rainy season from May-October [69]. The Northwest monsoons (wet) bring 90% of the rainfall, which generally varies between 1200 and 2000 mm per year across the country. The Northeast monsoons result in dry weather in the period of  [66]. Red color indicates high elevation and mountainous regions, while yellow and green colors indicate medium highlands and land surfaces with low elevation, respectively. The Tonle Sap lake and Mekong river are shown as water bodies. (b) Forest cover map in 2002 and geolocations of the field reference data collected in 2018. The star symbols represent the geolocations of land cover categories during field observations in February and July 2018. The dot symbols show the locations of the Permanent Sample Plots (PSPs) in evergreen (EG), semievergreen (SEG), and deciduous forests (DD). The red dots represent semievergreen forest locations where an inventory was conducted (refer to Reference [67] for details.) The photos shown on the maps were taken during the field visits and represent flooded forests (FF), PSPs, mangroves (MG), EG, rice or paddy fields (CR), DD, SEG, and rubber plantations (RB). EG, SEG, DD, and RB aerial photos were taken from the drone Phantom 4 (of DJI company) in July 2018. Reference data for the land cover were obtained from Sasaki et al., 2016 [65].
There are four environmental regions in Cambodia: The plain, Tonle Sap, coastal, and plateau and mountainous areas [68]. Cambodia has a tropical monsoon climate and there are two distinct seasons: The dry season from November-April and the rainy season from May-October [69]. The Northwest monsoons (wet) bring 90% of the rainfall, which generally varies between 1200 and 2000 mm per year across the country. The Northeast monsoons result in dry weather in the period of November to April,

Collection of Landsat Data and Image Composite
Temporal land cover change monitoring is needed to assess the anthropogenic impacts of climate change on the Earth's ecosystem and to propose appropriate interventions for natural resource management [74,75]. The Landsat satellite acquires over 40 years of data every 16 days. These data have been used for several phenological application purposes in forestry [76,77], cropland mapping [64,[78][79][80], land cover classification [53,54], and carbon accounting [56]. Many land surface phenology studies have found that Landsat images can be used with high accuracy [25,28,37,50,51,56,64,[81][82][83].
Using the freely available Landsat datasets in GEE, we selected existing Landsat  We accessed 2199 Landsat images and used these to extract the phenological behaviors of 12 land cover categories. We selected the Landsat 5 TM TOA collections and OLI TOA collections for single EVI threshold mapping. Accordingly, we obtained 159 collections between December to February. Table A1 (presented in Appendix A) details the open access Landsat datasets in GEE, which we used for assessing phenology and threshold mapping.
Landsat TOA reflectance-based EVI values were generally higher than surface reflectance (SR) based values [78] and therefore TOA data showed greater potential to provide accurate phenology-based classification [53,79,84,85]. Several studies have successfully applied TOA data for mapping the paddy rice [63,78,79,86,87] and mapping the vegetation distribution [17,22,53,84,88] using phenology-based methods. These studies have suggested that TOA data can be used for mapping with higher accuracy. In this study, we used TOA data to estimate phenological behaviors of the respective land cover categories and to determine the threshold values for better classification of land cover categories in the tropical regions.
Google Earth Engine provides a variety of Landsat specific processing approaches. Specifically, methods to compute at-sensor radiance, TOA reflectance, cloud scores, and cloud-free composites are openly accessible from the GEE [20]. The calibration coefficients of the TOA reflectance orthorectified datasets were extracted from the corresponding image metadata [89]. The GEE has already converted digital number (DN) values into TOA data by putting all remote sensing data from different sensors onto a common radiometric scale, thereby minimizing spectral differences caused by acquisition time, Sun elevation, and the Sun-Earth distance. Methods of converting data into TOA reflectance data are described in Reference [89]. In this study, we used Landsat TOA EVI data to estimate the phenology of selected land cover categories and land cover classifications.
We applied GEE algorithms to composite Landsat imagery and reduced the imagery to median pixel values [5]. The image composition defines the target reflectance of the composited TOA image. The median provides the finest balance between oversaturated and lower pixel values. Lower values focus on low vegetation greenness in the dry season composites and avoid oversaturated pixels. Upper values focus on high vegetation greenness and growing season composites; however, these values can include artifacts due to clouds and oversaturated pixels [5,89]. We observed that 60-75% of the study region was covered by cloud from May-November and 30-40% had cloud cover from December-April. The northeastern mountains, northern plains, Mekong alluvial plains, and coastal range have 35-45% cloud over for almost the entire year. To minimalize the cloud problem, we applied the GEE cloud thrash/mask function algorithm to the entire Landsat collection for the purpose of identifying the imagery with less than 40% cloud cover and for avoiding the bad collections.
After that, we applied cloud masking function over the images in order to extract the EVI time series profiles of the selected land use categories between 2000-2003 and 2014-2017. To determine the surface phenology of the land cover categories, we selected Landsat EVI time series data for this study, as previous studies stated that the EVI overcomes the saturation problems of the NDVI, especially in areas with dense vegetation [81]. The formula was used to calculate the EVI for Landsat TOA composite image collections with GEE JavaScript: where the coefficients of the EVI equation are L = 1 (canopy background adjustment factor); C1 = 6 and C2 = 7.5 (aerosol correction factors); and G = 2.5 (gain factor) [23,27,63,90]. NIR represents the near-infrared band (TM band 4 and OLI band 5); RED represents TM band 3 and OLI band 4; BLUE represents TM band 1 and OLI band 2 [85,87].

Reference Data
Ground reference data on land cover categories were randomly collected through field visits in February and July 2018 ( Figure 1b). We geolocated 300 sampling points to ensure the exact locations of each land cover category in different locations in the study region. The locations were as follows: EG (12 locations including two Permanent Sample Plots (PSPs), SEG (11 locations including one PSP), DD (40 locations including 1 PSP), FF (9 locations), paddy fields (120 locations), rubber (90 locations), bamboo (10 locations), and MG (8 locations). A field survey was conducted, which focused on the placement of PSPs and covered the main land cover categories in five provinces, namely Siem Reap, Mondulkiri, Takeo, Koh Kong, and Preah Vihear (Figure 1b). All the field-based geolocations were recorded using WSMaps mobile app (http://swmaps.softwel.com.np/). We flew the DJI Phantom 4 pro drone at five locations: SEG forests in Koh Kong and Mondalkiri provinces, DD forests in Mondulkiri and Takeo provinces, and an RB plantation in Mondulkiri province (Figure 1b). The drone aerial photos were processed using the Pix4D (https://www.pix4d.com/) drone image processing tool for mosaicking the photos and ortho-rectifying the drone images. The orthorectified image geometric corrections were automatically done in Pix4D. The orthorectified images were imported into GEE using the assets manager, and the drone image geolocations were verified by overlaying them on the GEE VHR image before assigning the geometry. After validation of the drone images, we determined the geometry of the forest types.
Prior to utilizing the ground reference data, we improved the reliability of the reference dataset by cross-validating forest cover data from 2002 and 2006 at the commune level ( Figure 1b). A detailed land cover mapping method used in 2002 and 2006 was reported in Reference [65]. We used these data to validate the geometric locations of forest categories in GEE, namely EG, SEG, DD, BB, and Mix WS. These data were also used to assess the accuracy of the resulting threshold spatial maps.
In addition to improving the ground reference data, we collected the locations of PSPs. The forestry administration established the PSP system in 1998 with five or eight clusters with four main plots (each plot 50 × 50 m) and four sub-plots (each 20 × 20 m) for each of the five provinces in Cambodia (32 plots for Kratie and 20 plots each for Karatie, Siem Reap, Kokong, and Ratanakiri) [91]. For this study, we selected 4 PSPs in Siem Reap province followed by EG, SEG, and DD; 3 EG PSPs in Kokong province; and 6 PSPs (5 EG and 1 DD location) in Kampong Thom province. In addition, we geolocated 5 SEG forest inventory locations in Kratie and Ratanakiri provinces by referring to Figure 1b [67]. The forest cover types (2002 and 2006), PSPs, and SEG sampling locations were cross-validated using field reference data and high spatial resolution Google Earth images to produce the reference data for this study.
We imported the reference data (ESRI shapefile) into GEE using the GEE Assets Manager tool and determined the geometry of the reference data using the GEE geometry tool to get the EVI phenology of the individual land cover categories ( Figure 2). Eventually, we validated the EVI profiles and selected Remote Sens. 2019, 11, 1514 8 of 30 10 geometries per land cover category, which resulted in 120 geometries for 12 land cover categories, which were used for the GEE phenological assessment of the selected land cover categories. 8 of 32 Figure 2. shows the 120 geometry locations of 12 land cover categories in GEE which we geolocated based on the reference data and visual interpretation by referring to forest cover data from 2000 and 2006. The temporal EVI chart function presented is in the code editor section. In the main map window, the background image is an OLI TOA composite image (band combination of B6, B5, and B4).

Land Surface Phenology Approach
Specifying the land surface phenology is essential and useful for detailed mapping, and it avoids the frequent misclassification of land cover types. However, having spectral knowledge of the land cover categories is preferable because spectral reflectance values vary over land cover features and time. Therefore, it is advantages to harmonize time-series remote sensing data with the phenological cycle. We used the fast computing abilities of GEE to extract the temporal time series of spectral characteristics of the land cover categories from the GEE harmonic regression model of the TM and OLI collections.
Harmonic regression is a mathematical technique that is used to decompose a complex, static signal into a series of individual sine or cosine waves, each characterized by a specific phase and amplitude [92]. The harmonic fitting method procedure of the GEE was presented by Nick Clinton (GEE relations developer) during the GEE 2016 User's Summit (details on the time series harmonic model algorithms can be found at https://developers.google.com/earth-engine/tutorials#time-seriesanalysis). We adapted Nick Clinton's GEE harmonic regression algorithms for Landsat collections to estimate EVI fitted and original time series values for selected land cover categories. A detailed flowchart of phenology estimations and the threshold mapping approach is shown in Figure 3, and the following steps were carried out to assess the phenological behaviors of major land cover categories in this study: • Landsat TM TOA and Landsat OLI TOA data were collected from the GEE. The underwent preprocessing and the time series of the collections were filtered. The image filter function and the composite, and reducer functions were applied to determine the median values within the study region. The EVI function was used to determine the phases of the composite median imagery.

•
The data were smoothed and the geometry was plotted on reference land cover data to form temporal EVI profiles. The harmonic model function was used to compute the linear regression reducer, and the coefficients were plugged into the linear model to get the phase and amplitude.

•
For phenology estimation, the average and mean EVI profiles of pixels of corresponding land cover types were calculated, and a search process was used to determine the phenology parameters, e.g., high-peak values and low values in time were identified through LSP and LFP as well as SOS and EOS. The threshold values for individual land cover categories were

Land Surface Phenology Approach
Specifying the land surface phenology is essential and useful for detailed mapping, and it avoids the frequent misclassification of land cover types. However, having spectral knowledge of the land cover categories is preferable because spectral reflectance values vary over land cover features and time. Therefore, it is advantages to harmonize time-series remote sensing data with the phenological cycle. We used the fast computing abilities of GEE to extract the temporal time series of spectral characteristics of the land cover categories from the GEE harmonic regression model of the TM and OLI collections.
Harmonic regression is a mathematical technique that is used to decompose a complex, static signal into a series of individual sine or cosine waves, each characterized by a specific phase and amplitude [92]. The harmonic fitting method procedure of the GEE was presented by Nick Clinton (GEE relations developer) during the GEE 2016 User's Summit (details on the time series harmonic model algorithms can be found at https://developers.google.com/earth-engine/tutorials# time-series-analysis). We adapted Nick Clinton's GEE harmonic regression algorithms for Landsat collections to estimate EVI fitted and original time series values for selected land cover categories. A detailed flowchart of phenology estimations and the threshold mapping approach is shown in Figure 3, and the following steps were carried out to assess the phenological behaviors of major land cover categories in this study:

•
Landsat TM TOA and Landsat OLI TOA data were collected from the GEE. The underwent preprocessing and the time series of the collections were filtered. The image filter function and the composite, and reducer functions were applied to determine the median values within the study region. The EVI function was used to determine the phases of the composite median imagery.

•
The data were smoothed and the geometry was plotted on reference land cover data to form temporal EVI profiles. The harmonic model function was used to compute the linear regression reducer, and the coefficients were plugged into the linear model to get the phase and amplitude.
• For phenology estimation, the average and mean EVI profiles of pixels of corresponding land cover types were calculated, and a search process was used to determine the phenology parameters, e.g., high-peak values and low values in time were identified through LSP and LFP as well as SOS and EOS. The threshold values for individual land cover categories were determined at specific phenology phases.

•
To develop the mapping function for single EVI data, the GEE image reducer and median function were used to form composite time series images to give a single median EVI. The PBTC function was applied for median EVI classification. Eventually, the resulting threshold maps were validated using a VHR Geographic Information System (ArcGIS) image and reference forest cover data.

of 32
• To develop the mapping function for single EVI data, the GEE image reducer and median function were used to form composite time series images to give a single median EVI. The PBTC function was applied for median EVI classification. Eventually, the resulting threshold maps were validated using a VHR Geographic Information System (ArcGIS) image and reference forest cover data.

Selection of Land Cover Categories for the Phenology Assessment
To examine the phenological characteristics of the land cover categories, we selected major 12 land cover categories: EG, SEG, DD, FF, MG, MixWS, BB, CR, RB, BL, SN and WA. Due to a lack of reference data and spatial resolution, pine plantations, tree plantations, forest regrowth, and oil palms were merged, where their spectral reflectance was similar to the EG, SEG, DD, and Mix WS categories. The grassland category was merged with the cropland categories as its spectral behaviors are similar to those of croplands during the selected phenological seasons, and the spatial distribution of grass in the study region is minor.
We identified 240 geometries for the 12 land cover categories and 10 geometries for each land cover type ( Figure 4) by referring to the reference data (Figure 1a

Selection of Land Cover Categories for the Phenology Assessment
To examine the phenological characteristics of the land cover categories, we selected major 12 land cover categories: EG, SEG, DD, FF, MG, MixWS, BB, CR, RB, BL, SN and WA. Due to a lack of reference data and spatial resolution, pine plantations, tree plantations, forest regrowth, and oil palms were merged, where their spectral reflectance was similar to the EG, SEG, DD, and Mix WS categories. The grassland category was merged with the cropland categories as its spectral behaviors are similar to those of croplands during the selected phenological seasons, and the spatial distribution of grass in the study region is minor.
We identified 240 geometries for the 12 land cover categories and 10 geometries for each land cover type ( Figure 4) by referring to the reference data (Figure 1a,b) from the TM and OLI collections. BL, SN and WA geometries were marked based on image visual interpretation in GEE.
Google Earth Engine generated fitted and actual EVI profiles were stored in comma-separated value (CSV) files, which we exported outside of GEE and Microsoft Excel was used for further statistical analysis. Based on the spatial resolution and the value of each Landsat pixel (30 m), the harmonic-fitted mean EVI and the original mean EVI values were calculated annually and monthly (Table 2 and Figure 5).
This resulted in the formation of 722 mean EVI profiles for 12 land cover categories, which described the behaviors of the land cover categories vegetation index spectral response at the SOS and EOS as well as during LSP and LFP. palms were merged, where their spectral reflectance was similar to the EG, SEG, DD, and Mix WS categories. The grassland category was merged with the cropland categories as its spectral behaviors are similar to those of croplands during the selected phenological seasons, and the spatial distribution of grass in the study region is minor.
We identified 240 geometries for the 12 land cover categories and 10 geometries for each land cover type ( Figure 4) by referring to the reference data (Figure 1a,b) from the TM and OLI collections. BL, SN and WA geometries were marked based on image visual interpretation in GEE.    Google Earth Engine generated fitted and actual EVI profiles were stored in comma-separated value (CSV) files, which we exported outside of GEE and Microsoft Excel was used for further statistical analysis. Based on the spatial resolution and the value of each Landsat pixel (30 m), the harmonic-fitted mean EVI and the original mean EVI values were calculated annually and monthly (Table 2 and Figure 5). This resulted in the formation of 722 mean EVI profiles for 12 land cover categories, which described the behaviors of the land cover categories vegetation index spectral response at the SOS and EOS as well as during LSP and LFP.  We observed missing original EVI pixel values at some geometry locations of land cover categories on both TM and OLI images between June and September. The missing data could be due to clouds on the imagery. We compared both original and fitted EVI values and found that the fitted mean EVI values had higher accuracy than the original EVI values (Figure 5c,f,i). Therefore, we used the harmonic-fitted EVI values for phenology estimation in this study.

Land Cover Category Phenology and Validation
Previous studies have applied numerous methods to extract land surface phenology metrics from time series of remote sensing images by adopting two definitions for SOS and EOS [45,93]. These are based on the mid-values at the start of the season in the time series and the peak values of the second derivative of the growing season curve [94]. Some studies have suggested that the peakvalues agree more closely with the phenological transitions observed on the ground [93] because of the biophysical meaning of this definition [95], and other studies stated that the mid-values are a more reliable method for determining vegetation phenology because of uncertainty (due to snow cover and cloud contamination) of the remote sensing time series images [96]. (e) OLI monthly average EVI profiles of SEG trees; (f) the difference between annual mean EVI original and fitted EVI profiles of SEG forest produced by OLI and TM; (g) TM deciduous (DD) forest phenology spectral behaviors; (h) OLI monthly EVI phenology behaviors of DD forests; and (i) comparison of the annual mean EVI phenology profiles of DD forest produced by TM and OLI. We generated 722 phenological graphs to validate the phenological behaviors of selected land cover categories at the start-of-season (SOS) and end-of-season (EOS) and during the leaf-shedding (LSP) and leaf-flushing (LFP) periods.
We observed missing original EVI pixel values at some geometry locations of land cover categories on both TM and OLI images between June and September. The missing data could be due to clouds on the imagery. We compared both original and fitted EVI values and found that the fitted mean EVI values had higher accuracy than the original EVI values (Figure 5c,f,i). Therefore, we used the harmonic-fitted EVI values for phenology estimation in this study.

Land Cover Category Phenology and Validation
Previous studies have applied numerous methods to extract land surface phenology metrics from time series of remote sensing images by adopting two definitions for SOS and EOS [45,93]. These are based on the mid-values at the start of the season in the time series and the peak values of the second derivative of the growing season curve [94]. Some studies have suggested that the peak-values agree more closely with the phenological transitions observed on the ground [93] because of the biophysical meaning of this definition [95], and other studies stated that the mid-values are a more reliable method for determining vegetation phenology because of uncertainty (due to snow cover and cloud contamination) of the remote sensing time series images [96].
Here, we observed the seasonal variation in the vegetation indices as a general trend in the EVI time series, where high values were accrued towards the end of the growing season (August-October) and low EVI values were accrued during the dry season (February-April). We identified the seasonal minimum and maximum EVI values using the annual fitted harmonic profiles of each land cover category as indicators of the LSP and following LFP ( Figure 6). We also identified the peak-value phase between SOS and the decaying end of the phenology season (EOS) [27] (Figure 7).

of 32
Here, we observed the seasonal variation in the vegetation indices as a general trend in the EVI time series, where high values were accrued towards the end of the growing season (August-October) and low EVI values were accrued during the dry season (February-April). We identified the seasonal minimum and maximum EVI values using the annual fitted harmonic profiles of each land cover category as indicators of the LSP and following LFP ( Figure 6). We also identified the peak-value phase between SOS and the decaying end of the phenology season (EOS) [27] (Figure 7).  For determination of the EVI phenological mid-values [96], the ratios of the EVI amplitude to the annual mid-value in the early growing season (May-July) and early dry season (December-

of 32
Here, we observed the seasonal variation in the vegetation indices as a general trend in the EVI time series, where high values were accrued towards the end of the growing season (August-October) and low EVI values were accrued during the dry season (February-April). We identified the seasonal minimum and maximum EVI values using the annual fitted harmonic profiles of each land cover category as indicators of the LSP and following LFP ( Figure 6). We also identified the peak-value phase between SOS and the decaying end of the phenology season (EOS) [27] (Figure 7).  For determination of the EVI phenological mid-values [96], the ratios of the EVI amplitude to the annual mid-value in the early growing season (May-July) and early dry season (December- For determination of the EVI phenological mid-values [96], the ratios of the EVI amplitude to the annual mid-value in the early growing season (May-July) and early dry season (December-February) were calculated to determine minimum and peak-values of each land cover category during the mid-phenology vegetation index phase. We observed uncertainties (cloud and shadow) in the Landsat time series images during the growing season (local rainy season) and overlapping profiles of vegetation indices (EG, SEG, DD, BB, MG, and RB). Further, we validated all mean vegetation index values of each category and found a separation gap in the EVI mean profiles during the dry mid-LSP season (December-February).
Subsequently, the mean EVI profiles of the 12 land cover categories during the dry mid-phenology season were identified by separation gaps in the mid-LSP season; however, the spectral characteristics of the FF and MG species were similar to those of other forest categories during the mid-LSP season. To avoid misclassifying FF and MG, we created a boundary for the Tonle Sap flooded forest and the mangrove forest using the GEE geometry tool through visual interpretation, and we separately assigned threshold values to classify them. To determine the threshold values for individual land cover categories, we validated the mid-LSP season EVI mean standard bar values and selected standard bar low values that did not overlap with other categories as minimum threshold values. The mid-LPS season high-peak values were assigned as maximum threshold values.

Phenology-Based Threshold Classification Approach
The GEE cloud computing and image parallel processing capacity make it efficient to run spatial reductions of the image collection. The tile-by-tile processing method of GEE applies a median reducer to each tile of the remote sensing data at a specific season or at a particular time [19].
Here, we applied the GEE image reducer and median function [19,26] to the TM and OLI collections for the dry mid-phenology season (December-February) to get a median EVI for threshold mapping. This involved eight steps. The diagram in Figure 3b illustrates the detailed approach for the single image pixel-based threshold mapping of individual land cover categories: (1) The Landsat collections were accessed using the image collection function in GEE; (2) the reducer function was used to select the seasonal imagery; (3) the tile-by-tile images were converted to the median; (4) the composite function was used to form a median single image for the dry mid-phenology season; (5) the EVI function for the median image was determined; (6) the threshold values for individual land cover categories were assigned; and (7) the PBTC function was used for spatial mapping in GEE.
The GEE basically uses two main programming languages, JavaScript and Python, to compute and analyze the remote sensing imagery. Here, we applied JavaScript to build PBTC algorithms to classify the median EVI data from the single composite Landsat image using predefined knowledge-based fuzzy rules in the GEE. We applied a similar study method developed by Simonetti et al. (2015), this approach which was later implemented in GEE [54].

Accuracy Assessment
The number of reference points collected during the field survey did not cover the entire study region; therefore, the mapping accuracy of the resulting spatial maps was evaluated using VHR imagery and reference forest cover data. Most commonly, stratified random sampling and systematic sampling are implemented to design reference data for validating classification maps. However, stratified random sampling is a practical design that satisfies the basic accuracy assessment objectives of the remote sensing community [97]. Remote sensing time series data with VHR images can be used to determine the accuracy of reference change data, and the process of generating the reference data is more accurate when the established protocol is followed [97].
To assess the accuracy, a total of 1420 spatial points were established for the study region using the stratified random sampling approach [97]. We used the ArcGIS environment to generate random spatial points by creating 355 reference points for each map from 2000 (TM) and 2018 (OLI). Further, to validate our results, we created 355 random points for each forest cover data reference point in 2000 and 2006, and these data were used for a spatial assessment. Here, the sampling size of the reference data (355 sampling points) was used for assessment of the resulted maps. This was to ensure that each land cover category was properly represented by referenced points, thereby reducing bias in the accuracy assessments. We converted the accuracy points to polygons using ArcGIS Future Envelope to Polygon tool to evaluate the classified maps with 30 m pixels. The TM map (2000) reference points and polygons were converted to KML format and exported into Google Earth for validation using Google Earth timescale temporal imagery from 2000. The OLI 2018 map was validated with ArcGIS using VHR imagery.
The reference forest cover maps were classified into 8 categories (Figure 1b): EG, SEG, DD, other forests (OF), wood and shrub evergreen (WE), wood and shrub dry (WD), non-forest (NF), and Bamboo BB. To evaluate the reference forest cover map accuracy, we merged the threshold map categories: CR, BL, SN, and WA were merged into the NF category; FF, MG, and RB were merged into the OF category; and WE and WD were merged into the MixWS categories. The major forest categories, EG, SEG, DD, and BB, were kept the same. The resulting PBTC maps were exported to Google drive using the export function in GEE. We extracted the resulted phenology-based threshold map category labels from the spatial accuracy points using the ArcGIS spatial analysis tool. The reference forest cover maps from 2002 and 2006 and the threshold maps were overlapped in ArcGIS, and they were visually interpreted with the accuracy points.
We analyzed a 30 m Landsat pixel window using accuracy polygon and accuracy sampling points data in ArcGIS and Google Earth, which permitted us to interpret the resulting phenology-based threshold maps to determine the minimum and maximum tree and crown cover in single Landsat pixels (Figure 8).
14 of 32 in the accuracy assessments. We converted the accuracy points to polygons using ArcGIS Future Envelope to Polygon tool to evaluate the classified maps with 30 m pixels. The TM map (2000) reference points and polygons were converted to KML format and exported into Google Earth for validation using Google Earth timescale temporal imagery from 2000. The OLI 2018 map was validated with ArcGIS using VHR imagery.
The reference forest cover maps were classified into 8 categories ( Figure 1b): EG, SEG, DD, other forests (OF), wood and shrub evergreen (WE), wood and shrub dry (WD), non-forest (NF), and Bamboo BB. To evaluate the reference forest cover map accuracy, we merged the threshold map categories: CR, BL, SN, and WA were merged into the NF category; FF, MG, and RB were merged into the OF category; and WE and WD were merged into the MixWS categories. The major forest categories, EG, SEG, DD, and BB, were kept the same. The resulting PBTC maps were exported to Google drive using the export function in GEE. We extracted the resulted phenology-based threshold map category labels from the spatial accuracy points using the ArcGIS spatial analysis tool. The reference forest cover maps from 2002 and 2006 and the threshold maps were overlapped in ArcGIS, and they were visually interpreted with the accuracy points.
We analyzed a 30 m Landsat pixel window using accuracy polygon and accuracy sampling points data in ArcGIS and Google Earth, which permitted us to interpret the resulting phenologybased threshold maps to determine the minimum and maximum tree and crown cover in single Landsat pixels (Figure 8). The pixels with different color codes are the resulting maps of OLI, and the overlaid background image is an ArcGIS VHR image. The labeled names on the map are as follows: EG: Pixels representing evergreen forest with 90%-100% tree cover and closed crown cover; SEG: Semievergreen forest with 50%-80% tree cover and dense crown cover; DD: Deciduous forest with 20%-50% tree cover and sparse crown cover; WS: Mixed wood and shrublands, with minimal percentages of EG and DD tree cover and large areas covered by bush or grass or crops; C: Pixels representing croplands or pixels with no tree cover.
The accuracy points were compared to the corresponding classified pixels in VHR images in ArcGIS and Google Earth and with existing land cover maps. The percentage of agreements was recorded in an accuracy table. Eventually, the assessment datasets were used to assess the classification accuracy using a confusion matrix [97]. We calculated the overall accuracy, the producer's (errors of omission) accuracy, the user's (errors of commission) accuracy, and the kappa coefficients over the confusion matrix for both classified threshold maps. Figure 8. Visualization of the basic concepts of the accuracy assessment applied in this study: The pixels with different color codes are the resulting maps of OLI, and the overlaid background image is an ArcGIS VHR image. The labeled names on the map are as follows: EG: Pixels representing evergreen forest with 90-100% tree cover and closed crown cover; SEG: Semievergreen forest with 50-80% tree cover and dense crown cover; DD: Deciduous forest with 20-50% tree cover and sparse crown cover; WS: Mixed wood and shrublands, with minimal percentages of EG and DD tree cover and large areas covered by bush or grass or crops; C: Pixels representing croplands or pixels with no tree cover.

Phenological Behaviors of Land Cover Categories from Multi-Temporal Landsat Imagery
The accuracy points were compared to the corresponding classified pixels in VHR images in ArcGIS and Google Earth and with existing land cover maps. The percentage of agreements was recorded in an accuracy table. Eventually, the assessment datasets were used to assess the classification accuracy using a confusion matrix [97]. We calculated the overall accuracy, the producer's (errors of omission) accuracy, the user's (errors of commission) accuracy, and the kappa coefficients over the confusion matrix for both classified threshold maps.

Phenological Behaviors of Land Cover Categories from Multi-Temporal Landsat Imagery
Cambodia has a monsoon climate and its spatial phenological behaviors vary greatly depending on the geographical location, land use management, and land cover types. Changes in the phenological behaviors of land cover categories caused by variations in environmental proximate factors may have shaped those behaviors. The results of the EVI profiles of the tropical region land cover categories phenology behaviors are shown in Figure 9a,b.

of 32
Cambodia has a monsoon climate and its spatial phenological behaviors vary greatly depending on the geographical location, land use management, and land cover types. Changes in the phenological behaviors of land cover categories caused by variations in environmental proximate factors may have shaped those behaviors. The results of the EVI profiles of the tropical region land cover categories phenology behaviors are shown in Figure 9a The obtained mean EVI phenological profiles of the forest and cropland categories show the high-peak vegetation index values accrued at the end of the growing season (rainy season) from August-October and the mid EVI values observed during the growing season (May-July) and dry season (November to January) from 2000-2003 (TM) and from 2014-2017 (OLI) (Figure 9a,b). In the dry season, the EVI values of the forest and croplands category rapidly decreased from December-April, and low EVI values were accrued from March-April. The BL and SN profiles did not change much, and their spectral values stayed at a similar amplitude in both seasons. However, the separation of the gap between most of the selected land cover categories in the mid-dry phenology period between December-February is worth noting (in Figure 9a,b), see the three years mean phenology characteristics window; the box indicates the separation of the vegetation index profiles of the major land cover categories.
Comparing the phenology-based EVI profiles of land cover categories (Figure 9a,b), BB and RB showed greater vegetation index differences than EG and SEG over the entire phenological period  (Figure 9a,b), while some deciduous trees may start flowering or leaf-shedding from September; thus, the DD tree profiles dramatically decreased in EOS.
The phenology profiles of the other forest categories, MG and FF, showed a high level of overlap with DD and Mix WS during the growing season and with SE in the dry season (Figure 9a,b). The MG vegetation index high-peak values (TM: 0.381-0.532; OLI: 0.509-0.654) were accrued in the mid-dry phenology period and the FF high-peak values (TM: 0.213-0.518; OLI: 0.471-0.581) were observed at the end of the dry season. The FF profiles dramatically decreased during the months of July to October (rainy season), and they were mixed with croplands from 2000-2003. This is probably due to the increase in the water level in the Tonle Sap lake during the rainy season and the mixing of the reference pixels of FF and WA. In 2014-2017, a separation of the FF profile from that of CR can be seen clearly. This is probably due to a decrease in the WA level in the lake, resulting from the changes in precipitation patterns in the region in the past ten years.
When determining the phenology profile of croplands (paddy fields) per pixel, we clearly observed high EVI values (TM: 0.276-0.291; OLI: 0.417-0.426) during the growing period (October-November) and a decrease in the vegetation index profile was seen right after the harvesting period in dry-mid LSP (December-February). The BL, SN, and WA widely separated from the vegetation index across all seasons studied. It was observed that in the selected land cover categories, the greatest differences in the vegetation index profile among the forest, rubber, croplands, and non-forest categories were found during the mid-dry LSP season (November to February).

Mid-Dry Phenology and Threshold Mapping
As shown in Figure 10, the mid-LSP results indicated TM and OLI threshold values and suitable months for classifying the Landsat EVI data for accurate land use and land use change mapping.
The threshold values were obtained by validating 722 mean EVI profiles during the leaf-flushing and leaf-shedding seasons.  Table A2 in Appendix A for details on the land cover categories mid-dry phenology season threshold values of TM (2000)(2001)(2002)(2003) and OLI (2014-2017) data. Cloud and shadow threshold values were excluded in this study.

Phenology-Based Threshold Map and Accuracy Assessment
The phenology-based threshold maps were shown to be highly accurate according to the confusion matrix estimated from the random reference data from 2000 and 2018, which validated VHR imagery and reference forest cover data. Table 3a,b and Table 4a,b provide an overview of the overall accuracies (OA), Kappa coefficients, producer's accuracies (PA) and user accuracies (UA) of the 2000 and 2018 maps validated by VHR imagery. Table 4a,b provides an overview of the overall accuracies matrix validated by the reference forest cover data.
Validation of the 2018 map with VHR imagery (Table 3a) (Table 3b). The FF, MG, BB, and RB categories were the main sources of confusion in the EG and SEG forests, and MixWS was the source of confusion in the DD forests in both classified years (Table 3a,  We found that the overall accuracy in VHR of the threshold maps varied from 89.58% in 2018 to 87.89% in 2000 with Kappa coefficients from 0.88 and 0.86 (Table 3a,b), and the overall accuracy of the referenced forest cover varied from 87.32-83.38% with Kappa coefficients from 0.85-0.83 over the same period (Table 4a,b). Among the accuracy results, the forest reference land cover data had lower accuracy. The TM and OLI products both had high spatial consistency in the mapping of phenology-based threshold maps in the tropical region. OLI gave more accurate results than that of the TM for both phenology estimation and threshold map accuracy. However, OLI imagery has a higher capability for threshold-based classification using phenology information compared with Landsat TM data.  Figure 11a,b show Cambodia's first mid-dry phenology-based threshold land cover map using the PBTC approach with the aid of the fast cloud-computing GEE tool. The threshold maps were classified into 12 major land cover categories in Cambodia: Evergreen forest, semievergreen forest, deciduous forest, mixed wood and shrub (mixed with evergreen and deciduous trees), bamboo, rubber plantation, flooded forest, mangroves, croplands, built-up land, sand or soil, and water as shown in the spatial map legend. The validation of the results with reference forest cover data found that the EG (PA 84.9-79. We found that the overall accuracy in VHR of the threshold maps varied from 89.58% in 2018 to 87.89% in 2000 with Kappa coefficients from 0.88 and 0.86 (Table 5a,b), and the overall accuracy of the referenced forest cover varied from 87.32%-83.38% with Kappa coefficients from 0.85-0.83 over the same period (Table 6a,b). Among the accuracy results, the forest reference land cover data had lower accuracy. The TM and OLI products both had high spatial consistency in the mapping of phenology-based threshold maps in the tropical region. OLI gave more accurate results than that of the TM for both phenology estimation and threshold map accuracy. However, OLI imagery has a higher capability for threshold-based classification using phenology information compared with Landsat TM data. Figures 11a,b show Cambodia's first mid-dry phenology-based threshold land cover map using the PBTC approach with the aid of the fast cloud-computing GEE tool. The threshold maps were classified into 12 (a) The resulting PBTC maps were identified as being more heterogeneous, distinguished communities of EG, SEG, DD, and Mix WS over the entire study region. The mid-LSP-based threshold maps comparing the TM and OLI images (Figure 11a,b) indicated that the EG forest pixels were well distributed, covering 21.01% of the study region, and this dramatically decreased to 17.56% in 2018. Similarly, SE forest pixels represented 17.40% of the region in 2000 and decreased to 11.67% in 2018. The DD forest pixels were widely distributed (31.08% of the region in 2000) and showed higher accuracy in mapping in the study region, but their distribution declined to 23.45% in 2018. The distribution of the Mix WS pixels decreased from 25.64% to 18.30%. Bamboo was found over the entire region (0.24% in 2000 and 0.17% in 2018), with larger percentages of BB observed in the Northern parts of Ratana Kiri and Stung Treng provinces and in the Western part of the study region, especially in Koh Kong, Pursat, Battambang, Kampot, and Kampong Speu provinces. The distributions of RB (0.14%-0.80%) and cropland (10.03%-19.15%) pixels increased dramatically from 2000-2018, and the distribution of BL increased from 1.18-1.55%. Regarding the other forest categories, FF is largely found around the Tonle Sap lake and MG is found in Koh Kong, Preah Sihanouk and Kampot, and Kep provinces. However, their percentages of pixel cover may vary when the other forest categories, pine trees, pine plantation, tree plantation, and forest regrowth are included in the threshold mapping. The resulting PBTC maps were identified as being more heterogeneous, distinguished communities of EG, SEG, DD, and Mix WS over the entire study region. The mid-LSP-based threshold maps comparing the TM and OLI images (Figure 11a,b) indicated that the EG forest pixels were well distributed, covering 21.01% of the study region, and this dramatically decreased to 17.56% in 2018. Similarly, SE forest pixels represented 17.40% of the region in 2000 and decreased to 11.67% in 2018. The DD forest pixels were widely distributed (31.08% of the region in 2000) and showed higher accuracy in mapping in the study region, but their distribution declined to 23.45% in 2018. The distribution of the Mix WS pixels decreased from 25.64% to 18.30%. Bamboo was found over the entire region (0.24% in 2000 and 0.17% in 2018), with larger percentages of BB observed in the Northern parts of Ratana Kiri and Stung Treng provinces and in the Western part of the study region, especially in Koh Kong, Pursat, Battambang, Kampot, and Kampong Speu provinces. The distributions of RB (0.14-0.80%) and cropland (10.03-19.15%) pixels increased dramatically from 2000-2018, and the distribution of BL increased from 1.18-1.55%. Regarding the other forest categories, FF is largely found around the Tonle Sap lake and MG is found in Koh Kong, Preah Sihanouk and Kampot, and Kep provinces. However, their percentages of pixel cover may vary when the other forest categories, pine trees, pine plantation, tree plantation, and forest regrowth are included in the threshold mapping.

Discussions and Implications
Our study indicates the viability of assessing phenology behaviors and threshold-based mapping of major land cover categories in the tropical region of Cambodia. To our knowledge, this study is the first to assess the phenological behaviors of major land cover categories and to determine threshold values at a 30 m moderate spatial resolution in the tropical region. The GEE freely available fast cloud-computing environment was used to access multitemporal big remote sensing data and computed the massive multitemporal Landsat collections within a limited time period. This study overcame challenges including the effects of cloud cover and missing Landsat observations by determining climatic threshold values and using GEE cloud-computing algorithms.
Firstly, the spatial and temporal resolution of remote sensing data for phenology assessment needs to be selected carefully. In this study, OLI provided high-frequency imagery of 60.35% on average (2014-2017) and TM had average imagery covered of 39.65% (2000)(2001)(2002)(2003) in both the leaf-flushing and leaf-shedding phenology seasons ( Table 2). The harmonic regression fitting process which produced the 16 days composite of the entire Landsat collection was performed to provide considerable smoothing of the EVI time-series profiles for both the TM (Figure 9a) and OLI (Figure 9b) datasets. Accordingly, the high-peak values of the phenology-based vegetation index profiles were deducted in the growing season (rainy season). However, the decline of EG, SEG, and BB in the rainy season from July-August (Figure 9a) in the TM images could be due to the minimal observations and high cloud contamination of TM [62]. A previous study found similar results that vegetation index profiles were drastically reduced to the lowest level during the heavy rainy season [98]. On the other hand, the OLI images had more observations across the entire phenology study and showed a realistic vegetation index frequency in both phenological seasons as compared with the TM images ( Figure 9b).
Secondly, the assessment of the phenology characteristics using TM and OLI showed noticeable separation among the vegetation index profiles of the land cover categories in the mid-LSP season. However, some spectral characteristics of the species in the flooded forest showed similar vegetation index values to SE and EG. Similarly, the MG vegetation index values behaved like those of other natural forest categories (Figure 10a,b). Thus, to avoid misclassification (Tables 3 and 4) of FF and MG with other forest categories, we assigned threshold values separately for FF and MG. However, there were still challenges when it came to determining the exact minimum threshold values for DD and Mix WS and the maximum threshold value of CR in the mid-LSP season ( Figure 10). The Mix WS EVI amplitude was mixed with DD from January to February in TM and with CR in OLI data. Therefore, we took the mean values of two months (December and January) to determine the minimum threshold values for Mix WS in TM and OLI data. Furthermore, the temperature is the major limitation for phenology-based threshold values. It is the climate condition with the biggest effects on the vegetation phenology pattern stages [85,99,100] and can lead to the misclassification of CR [78] and Mix WS in the dry season. The dry season (leaf-shedding period of the SEG, DD, and Mix WS) is when all of the herbaceous vegetation has been dried out and the bare ground soil is exposed. Minor phenological changes in the evergreen trees and shrubs are barely noted in the vegetation index time series data [101]. See the climate-sensitive map and the possibility of misclassification in phenology-based threshold mapping in the study region in Figure 12.
This study suggests that a field-based survey is required in both phenology seasons, and finer resolution spatiotemporal data should be used to find the spectral differences of overlapping categories. Further, we need to be mindful when interpreting the vegetation phenology behaviors while determining the threshold values during the mid-dry phenology season.
Third, GEE provides an up-to-date library of 40 years of spatially referenced remote sensing data and TOA reference data with composite image functionality for the selection of remote sensing data at a particular period and a cloud mask function to exclude the cloud pixels prior to the image threshold-based classification. The GEE online user guide provides helpful examples and its discussion forums are very useful. Users can share their scripts and get support from GEE expertise. The GEE powerful cloud-computing platform allows us to effectively process a large amount of multitemporal Landsat imagery for threshold-based land cover mapping, which speeds up the whole process to less than 2 minutes per map without using any commercial geospatial tools. This study provides a basis for understanding and applying the GEE cloud-computing platform for phenology-based study using remote sensing data, land cover mapping, and monitoring. However, basic programming skills (JavaScript, Python language) and a stable internet connection are required to develop the code in the GEE platform to perform the spatial analysis. This study suggests that a field-based survey is required in both phenology seasons, and finer resolution spatiotemporal data should be used to find the spectral differences of overlapping categories. Further, we need to be mindful when interpreting the vegetation phenology behaviors while determining the threshold values during the mid-dry phenology season.
Third, GEE provides an up-to-date library of 40 years of spatially referenced remote sensing data and TOA reference data with composite image functionality for the selection of remote sensing data at a particular period and a cloud mask function to exclude the cloud pixels prior to the image threshold-based classification. The GEE online user guide provides helpful examples and its discussion forums are very useful. Users can share their scripts and get support from GEE expertise. The GEE powerful cloud-computing platform allows us to effectively process a large amount of multitemporal Landsat imagery for threshold-based land cover mapping, which speeds up the whole process to less than 2 minutes per map without using any commercial geospatial tools. This study provides a basis for understanding and applying the GEE cloud-computing platform for phenologybased study using remote sensing data, land cover mapping, and monitoring. However, basic programming skills (JavaScript, Python language) and a stable internet connection are required to develop the code in the GEE platform to perform the spatial analysis.
Finally, focusing directly on the aspects of land cover classification and forest cover change is seen as an advantage of using the GEE fast cloud-computing platform, and the PBTC approach can be used to report activity data to establish a forest reference emission level (FREL) in the context of the REDD+ scheme under the UNFCCC agreement. On the other hand, access to freely available remote sensing data and the use of the PBTC approach can assist environmental managers to map and monitor natural resources at scale, for example, to monitor forest cover and related carbon stock changes in protected areas and farming areas using a cost-effective and transparent process.  [102] and land surface temperature [103] data were accessed, and the analysis was done in the GEE environment.
Finally, focusing directly on the aspects of land cover classification and forest cover change is seen as an advantage of using the GEE fast cloud-computing platform, and the PBTC approach can be used to report activity data to establish a forest reference emission level (FREL) in the context of the REDD+ scheme under the UNFCCC agreement. On the other hand, access to freely available remote sensing data and the use of the PBTC approach can assist environmental managers to map and monitor natural resources at scale, for example, to monitor forest cover and related carbon stock changes in protected areas and farming areas using a cost-effective and transparent process.
The PBTC method can be adapted to any climatic zones for detailed land use and land cover mapping. The obtained threshold values for land cover categories are potentially useful for monitoring the vegetation changes and land cover mapping in the dry or monsoon tropical regions. However, the tropical rain forests such as that in Indonesia, Malay peninsula may have higher vegetation index values than the monsoon forest in Cambodia because of their favorable climatic and soil conditions in the entire year. Therefore, the threshold values in our study need to be used carefully, or further study is needed to determine the threshold values for tropical rainforests.

Conclusions
Making use of increasingly available data and technologies for assessing the phenology-based land use and land cover change is needed for natural resources management planning and monitoring at scale. We presented a phenology-based threshold classification approach or PBTC to assess the phenological behaviors of 12 land cover categories (evergreen forest, semievergreen forest, deciduous forest, mix wood and shrub, flooded forest, mangrove forest, bamboo, rubber plantation, croplands, built-up area, sand, and water) during the leaf-shedding and leaf-flushing phenology seasons and explored its ability to determine threshold values to accurately map major land cover categories in Cambodia using a set of Landsat median EVI data and Google Earth Engine. The TM and OLI-based harmonic regression mean EVI profiles showed unique phenology behaviors in the major land cover categories in the phenology seasons, and the separation of vegetation index profiles were noticed during the mid-leaf-shedding (mid-dry season) phenology phase. Phenology-based threshold values were then determined for the 12 land cover categories above in the mid-leaf-shedding phase by selecting their high-peak vegetation index values. The application of PBTC and the use of the freely available remote sensing and GEE fast cloud-computing platform makes it possible to map large areas quickly with higher accuracy, and the integration of harmonic regression, EVI profiles can be used to monitor the continuous changes in the vegetation indices. Due to its transparent methodology, the PBTC methods will be of interest for monitoring, reporting, and verifying of the REDD+ activities critically required for performance-based payment scheme under the REDD+ scheme of the United Nations Framework Convention on Climate Change.
A limitation of the obtained phenology-based threshold values is that a finer resolution of spatiotemporal remote sensing data are needed for different climatic regions, where threshold values of land cover categories may be periodically determined or updated in order to improve the accuracy of the land cover classification. Future studies will need to focus on understanding the nexus of climate change and its effects on phenological behaviors of the respective vegetation because vegetation is very sensitive to climate change. Further study on the PBTC methods for applications by climatic zones can further improve the accuracy of phenology-based threshold values, which has important implications for land cover mapping and monitoring at scale.
For further development of the PBTC, the GEE JavaScript will be available at http://bit.ly/2JKsNqL. Acknowledgments: The research benefited from the Asian Institute of Technology, who provided graduate assistance support. We would like to acknowledge the International Tropical Timber Organization (ITTO) fellowship grant for field research and the FRAWAS project team members for assisting during the data collection. In addition, we would like to thank the anonymous reviewers for their constructive comments on earlier versions of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.