Mapping the Natural Distribution of Bamboo and Related Carbon Stocks in the Tropics Using Google Earth Engine, Phenological Behavior, Landsat 8, and Sentinel-2

Although vegetation phenology thresholds have been developed for a wide range of mapping applications, their use for assessing the distribution of natural bamboo and the related carbon stocks is still limited, especially in Southeast Asia. Here, we used Google Earth Engine (GEE) to collect time-series of Landsat 8 Operational Land Imager (OLI) and Sentinel-2 images and employed a phenology-based threshold classification method (PBTC) to map the natural bamboo distribution and estimate carbon stocks in Siem Reap Province, Cambodia. We processed 337 collections of Landsat 8 OLI for phenological assessment and generated 121 phenological profiles of the average vegetation index for three vegetation land cover categories from 2015 to 2018. After determining the minimum and maximum threshold values for bamboo during the leaf-shedding phenology stage, the PBTC method was applied to produce a seasonal composite enhanced vegetation index (EVI) for Landsat collections and assess the bamboo distributions in 2015 and 2018. Bamboo distributions in 2019 were then mapped by applying the EVI phenological threshold values for 10 m resolution Sentinel-2 satellite imagery by accessing 442 tiles. The overall Landsat 8 OLI bamboo maps for 2015 and 2018 had user’s accuracies (UAs) of 86.6% and 87.9% and producer’s accuracies (PAs) of 95.7% and 97.8%, respectively, and a UA of 86.5% and PA of 91.7% were obtained from Sentinel-2 imagery for 2019. Accordingly, carbon stocks of natural bamboo by district in Siem Reap at the province level were estimated. Emission reductions from the protection of natural bamboo can be used to offset 6% of the carbon emissions from tourists who visit this tourism-destination province. It is concluded that a combination of GEE and PBTC and the increasing availability of remote sensing data make it possible to map the natural distribution of bamboo and carbon stocks.


Introduction
Bamboos are evergreen perennial flowering plants that can be found in the tropical and subtropical regions of the world [1]. Bamboos provide important ecosystem goods and services to people: for example, they serve as sources of food and raw materials for construction [2] and play important roles in soil erosion control, water conservation, and land rehabilitation [3]. Bamboo forests also have a potential role in climate change mitigation, acting as sinks or sources of atmospheric carbon [2,4], depending on how they are managed. Despite such importance, studies on the natural distribution of bamboos remain limited, making it difficult to assess their abundance or measure and monitor their carbon stocks. Therefore, it is critically important to map bamboo distributions using fast and timely technologies for transparent, speedy, and monitoring at scale. Understanding the natural distribution of bamboos also allows us to assess the roles of bamboos in climate change mitigation and rural community development.
Bamboos are plants in the subfamily Bambusoideae of Gramineae. More than 1500 species are found in the tropical and subtropical regions of Asia [1]. Bamboos are considered non-timber forest products (NTFPs) with the ability to provide food, raw materials [2], and other ecosystem services to industries and local people [2,5]. Natural bamboo and plantations can be managed to mitigate climate change [3,6]: as fast-growing, high-yield plants they are capable of sequestering between 5 and 12 tCO 2 ha −1 year −1 [3,5,6]. Therefore, the management of bamboo forests can be eligible under the REDD+ scheme (reducing emissions from deforestation and forest degradation, forest conversation, sustainable management of forests, and enhancement of forest carbon stocks) of the United Nations Framework Convention on Climate Change (UNFCCC) [6]. It can also contribute to the achievement of Sustainable Development Goals (SDGs) 7, 13, and 15 [3]. Nevertheless, effective management of bamboo forests requires an understanding of their distribution through accurate mapping [7].
Remote sensing imagery has been used for various phenological applications in land use and land cover classification [8], mapping of rubber plantations [9], cropland [10,11], and bamboo mapping [12]. Researchers have used medium-resolution Landsat satellite data for the assessment of land cover changes over the past 40 years. Such data have become freely available in global archives that offer detail on land use and land cover dynamics [13][14][15]. These Landsat data have 30 m spatial resolution and a 16-day temporal interval. Launched in February 2013, Landsat 8 is the latest data acquired using the Operational Land Imager (OLI) and Thermal Infrared Sensor sensors in addition to the existing Landsat series of NASA. The recent launch of the Sentinel-2 satellite mission of the European Space Agency in 2015 and 2017 provides a new opportunity for land-based mapping and monitoring in the tropics using Sentinel-2 s high-resolution multispectral data and 10-day temporal interval. While Landsat 8 OLI data have been used for phenology-based studies, very few studies have evaluated the potential use of Sentinel-2 imagery data for phenology-based land cover mapping. Previous studies attempted to examine surface phenology for land use and land cover classification [13] and rubber plantation [9,16] and cropland mapping [17] in tropical regions [11]. For example, Schwieder et al. [18] used remote sensing technologies to assess forest carbon stocks [18]. However, challenges remain in determining the phenology-based threshold values for bamboo distribution mapping because of the rapid changes in phenological behavior between seasons [19,20]. Such rapid changes in phenology make it difficult to avoid the misclassification of bamboo. A recent study indicates that the mapping of bamboo distribution needs to employ robust classification methods to avoid misclassification [13]. Freely available moderate-resolution remote sensing data in a cloud-computing platform offers significant potential [12] as a robust methodology, for example in developing countries in Southeast Asia such as Cambodia, where remotely sensed reconstructions of past land cover change have been limited.
Previous studies that employed current mapping techniques for bamboo stands either used very high resolution images without regard to temporal spectral behavior [7,19,21] or time-series Landsat data [12,21] and the enhanced vegetation index (EVI) [22][23][24][25][26][27] with an additional focus on temporal information. The majority of these studies also implemented traditional imagery

Remote Sensing Data
Google Earth Engine provides ready-to-use remote sensing data in an up-to-date library with moderate to high-resolution imagery. To characterize vegetation phenology, we used 297 moderateresolution 30 m Landsat 8 OLI top-of-atmosphere (TOA) collections from 2015 to 2018 (Table 1). Forty Landsat 8 OLI collections that cover the study area of Siem Reap Province were sourced between December and February of 2015 and 2018 for bamboo threshold-based classification [13] (Table 2). To further explore vegetation threshold values for Sentinel-2 products, we accessed 422 Sentinel-2 images from January to December 2019 (Table 2).

Remote Sensing Data
Google Earth Engine provides ready-to-use remote sensing data in an up-to-date library with moderate to high-resolution imagery. To characterize vegetation phenology, we used 297 moderate-resolution 30 m Landsat 8 OLI top-of-atmosphere (TOA) collections from 2015 to 2018 (Table 1). Forty Landsat 8 OLI collections that cover the study area of Siem Reap Province were sourced between December and February of 2015 and 2018 for bamboo threshold-based classification [13] (Table 2). To further explore vegetation threshold values for Sentinel-2 products, we accessed 422 Sentinel-2 images from January to December 2019 (Table 2). Table 1. Landsat 8 Operational Land Imager (OLI) collections used for assessing the phenological behavior of selected land cover categories. PATH  ROW  2015  2016  2017  2018   126  50  17  16  16  15  51  21  20  18  16  127  50  20  20  19  18  51  17  22  22  20  Total Collections  75  78  75  69  Table 2. Landsat OLI and Sentinel-2 image collections used in this study for phenology-based threshold classification (PBTC).  Figure 2 illustrates the methodological framework of selected vegetation phenological assessment and phenology-based threshold classification methods for bamboo mapping and bamboo carbon stock assessment.   Figure 2 illustrates the methodological framework of selected vegetation phenological assessment and phenology-based threshold classification methods for bamboo mapping and bamboo carbon stock assessment. Here, we used JavaScript programming in GEE to collect time-series of Landsat 8 OLI image collections from 2015 to 2018 and applied cloud mask functions (less than 45%) to exclude cloud pixels prior to the assessment of land use classification ( Figure 2, Step 1) [10,34,37,45,46]. We then applied an image compositing and reducer function to select the most recent pixel in the collections.

Methodology
We applied a composite median reducer function to calculate the median value of each image collection. The median reducer function removes clouds in the image, which have high values, and shadows, which have low values [35,47]. The output composite value is the median in each band over time [34,46]. We applied a clip function to subset the image collections within the study region and then calculated the EVI using Equation (1) (Figure 2, Step 1) (Appendix A): Here, we used JavaScript programming in GEE to collect time-series of Landsat 8 OLI image collections from 2015 to 2018 and applied cloud mask functions (less than 45%) to exclude cloud pixels prior to the assessment of land use classification ( Figure 2, Step 1) [10,34,37,45,46]. We then applied an image compositing and reducer function to select the most recent pixel in the collections.
We applied a composite median reducer function to calculate the median value of each image collection. The median reducer function removes clouds in the image, which have high values, and shadows, which have low values [35,47]. The output composite value is the median in each band over time [34,46]. We applied a clip function to subset the image collections within the study region and then calculated the EVI using Equation (1) (Figure 2, Step 1) (Appendix A): where G = 2.5 (gain factor) [48][49][50]; C1 = 6 and C2 = 7.5 (aerosol correction factors); L = 1 (canopy background adjustment factor) [48][49][50]; and NIR (near infrared) is Landsat band 5 and Sentinel-2 band 8, RED is band 4 in both Landsat and Sentinel-2, and BLUE is band 2 in both Landsat and Sentinel-2 [49].
In assessing the vegetation phenological profiles of selected land cover classes, we applied a time-series analysis for EVI by adapting harmonic and regression model JavaScript algorithms from our recent study [13] (Figure 2, Step 2). We selected 65 geometry locations (Figures 1 and 3D), including 30 for bamboo forest, 30 for evergreen forest, and 5 for rubber plantations, where no large land cover changes had occurred during the phenological analysis period (2015 to 2018). As we observed during the field survey in 2019, the larger portion of bamboo distribution was found in and around the Kulen Mountain area, including Varin, Banteay Seri, and Svay Leu Districts, and in the northern part of Chi Kreng District ( Figure 1). Thus, we selected 10 sites with a large homogeneous extent of bamboo around the Kulen Mountain for creating the geometry for analyzing the phenological behavior of bamboo.
In assessing the vegetation phenological profiles of selected land cover classes, we applied a time-series analysis for EVI by adapting harmonic and regression model JavaScript algorithms from our recent study [13] (Figure 2, Step 2). We selected 65 geometry locations (Figures 1 and 3D), including 30 for bamboo forest, 30 for evergreen forest, and 5 for rubber plantations, where no large land cover changes had occurred during the phenological analysis period (2015 to 2018). As we observed during the field survey in 2019, the larger portion of bamboo distribution was found in and around the Kulen Mountain area, including Varin, Banteay Seri, and Svay Leu Districts, and in the northern part of Chi Kreng District ( Figure 1). Thus, we selected 10 sites with a large homogeneous extent of bamboo around the Kulen Mountain for creating the geometry for analyzing the phenological behavior of bamboo.
For the rubber plantation samples, we created four geometry locations in Chi Kraeng District, two in the northern part of Varin, two in Angkor Thom, and two locations in Banteay Srei Districts. For evergreen forest geometries, we created three in the Angkor Wat region, where there are permanent sampling plots for evergreen forest [51], and seven geometries were created within the study area ( Figure 3). Eventually, we generated the time-series and fitted and original EVI profiles for each selected vegetation category to assess phenological behaviors [52] and determine the threshold values [13] from Landsat 8 OLI collections during the leaf-shedding phenology stage (December to March) from 2015 to 2018 ( Figure 2, Step 2). In Figure 3A, is the code window for the GEE time-series indices and harmonic regression JavaScript code developed in GEE to assess the phenological behaviors of selected vegetation categories; Figure 3B is the graphical window showing the harmonic and time-series phenological behaviors of evergreen, bamboo, and rubber vegetation's fitted (blue line) and original (red dots) vegetation index values from 2015 to 2018; Figure 3C shows geometry inputs of selected vegetation categories; Figure 3D is the map window showing the geometry of evergreen (green), bamboo (yellow), and rubber plantation (purple) in Siem Reap Province. The background image is Landsat 8 OLI composite EVI data in GEE. For the rubber plantation samples, we created four geometry locations in Chi Kraeng District, two in the northern part of Varin, two in Angkor Thom, and two locations in Banteay Srei Districts. For evergreen forest geometries, we created three in the Angkor Wat region, where there are permanent sampling plots for evergreen forest [51], and seven geometries were created within the study area ( Figure 3). Eventually, we generated the time-series and fitted and original EVI profiles for each selected vegetation category to assess phenological behaviors [52] and determine the threshold values [13] from Landsat 8 OLI collections during the leaf-shedding phenology stage (December to March) from 2015 to 2018 ( Figure 2, Step 2).
In Figure 3A, is the code window for the GEE time-series indices and harmonic regression JavaScript code developed in GEE to assess the phenological behaviors of selected vegetation categories; Figure 3B is the graphical window showing the harmonic and time-series phenological behaviors of evergreen, bamboo, and rubber vegetation's fitted (blue line) and original (red dots) vegetation index values from 2015 to 2018; Figure 3C shows geometry inputs of selected vegetation categories; Figure 3D is the map window showing the geometry of evergreen (green), bamboo (yellow), and rubber plantation (purple) in Siem Reap Province. The background image is Landsat 8 OLI composite EVI data in GEE.
Google Earth Engine generates EVI time-series profile datasets in a comma-separated values file format (CSV). Using the GEE export option, we exported 25 (10 evergreen, 10 bamboo, and five rubber The time-series analysis of fitted vegetation indices allowed us to quantify the phenological behaviors of bamboo, rubber, and evergreen forest categories. To do this, we first assessed the fitted vegetation index profiles derived from the time-series harmonic regression model using compressed (15- Previous phenological studies showed that fitted index values have higher accuracy in determining the vegetation threshold values for mapping during the LSP stage, locally called EOS [13].
We selected the fitted EVI profiles and assessed the phenology-based vegetation index values from November to April (EOS) or during the LSP season to determine the separation of bamboo index profiles from evergreen forests and rubber plantations from 2015 to 2018. We then created 12 average-fitted EVI phenological profiles for the three land use categories by year to estimate the average EVI profiles and determine the threshold values for bamboo forest, evergreen forest, and rubber plantation and to classify bamboo distribution in the region during the mid-dry season ( Figure 2, Step 3). Because the purpose of this study was to determine the natural distribution of bamboo forest, we did not consider the vegetation index values for other land cover categories in this study, which fell below the evergreen forest threshold values.
Since the resulting Landsat 8 OLI phenology threshold values did not match the Sentinel-2 composite EVI for the selected land use categories, to classify the Sentinel-2 composite EVI data by using the PBTC method, we assigned Min = 0.6712 and Max = 0.7760 for bamboo, 0.6695-0.6623 for rubber, and 0.555-0.660 for the evergreen forest category by modifying Landsat Thematic Mapper (TM) threshold values from our previous phenological study [13]. We used JavaScript to build the PBTC algorithm [13] for classifying the median EVI data from composite Landsat OLI and Sentinel-2 images. As depicted in Figure 2 (Step 4), the steps used for the composite image and pixel-based threshold mapping of individual land cover categories were as follows: • Imagery preprocessing and time-series collections were filtered for the dry mid-phenology season (from December to February). • Image filter functions were applied to select image collections in the particular phenological period within the study region [13].

•
To develop the mapping function for EVI data, the GEE image reducer and median function was used to form composite time-series images to obtain a median EVI [35,53].

•
The PBTC function was applied to the bamboo forest, evergreen forest, and rubber plantation categories.
The final PBTC maps were used to map the bamboo distribution and calculate the bamboo carbon stocks in Siem Reap Province ( Figure 4).

Accuracy Assessment
Google Earth higher-resolution imagery (VHR) is free of cost and can be used directly as a base source for validating land use, land cover maps [54] because it provides the most accurate data for selecting the most appropriate reference points that are very important for validation. The images are updated whenever new images become available. Depending on the sensor, the image resolution ranges from 30 m to 15 cm. Utilizing the time-lapse feature in Google Earth Pro allows the user to view zoomable images as far back as 30 years, which is ideal for validating land use, land cover maps [55] and performing investigation and preliminary studies with suitable accuracy [56,57] by applying recognized protocols, as recommended by [58]. Here, a stratified random sampling technique was employed for the validation of phenology-based threshold classified maps. From the number of random points generated in ArcGIS, a total of 1500 accuracy points were selected within the evergreen forest, bamboo, and rubber plantation land cover categories. We excluded the random points that fell within the water and other land cover categories. Because the spatial distribution of bamboo and rubber plantation areas is smaller than that of the evergreen forest and other land cover classes, we selected random reference points within the three forest categories to acquire a considerable number of reference points for the bamboo and rubber plantation category areas to better assess mapping accuracy [59].
PBTC maps were exported to Google Drive using the export function in GEE, and we used these maps in ArcGIS to transform land cover labels into 1500 accuracy points by applying ArcGIS spatialjoin and related tools. We then used these accuracy points in Google Earth Pro to verify and validate the 1500 points for Landsat OLI (2015 and 2018) and Sentinel-2 (2019) maps (

Accuracy Assessment
Google Earth higher-resolution imagery (VHR) is free of cost and can be used directly as a base source for validating land use, land cover maps [54] because it provides the most accurate data for selecting the most appropriate reference points that are very important for validation. The images are updated whenever new images become available. Depending on the sensor, the image resolution ranges from 30 m to 15 cm. Utilizing the time-lapse feature in Google Earth Pro allows the user to view zoomable images as far back as 30 years, which is ideal for validating land use, land cover maps [55] and performing investigation and preliminary studies with suitable accuracy [56,57] by applying recognized protocols, as recommended by [58]. Here, a stratified random sampling technique was employed for the validation of phenology-based threshold classified maps. From the number of random points generated in ArcGIS, a total of 1500 accuracy points were selected within the evergreen forest, bamboo, and rubber plantation land cover categories. We excluded the random points that fell within the water and other land cover categories. Because the spatial distribution of bamboo and rubber plantation areas is smaller than that of the evergreen forest and other land cover classes, we selected random reference points within the three forest categories to acquire a considerable number of reference points for the bamboo and rubber plantation category areas to better assess mapping accuracy [59].
PBTC maps were exported to Google Drive using the export function in GEE, and we used these maps in ArcGIS to transform land cover labels into 1500 accuracy points by applying ArcGIS spatial-join and related tools. We then used these accuracy points in Google Earth Pro to verify and validate the 1500 points for Landsat OLI (2015 and 2018) and Sentinel-2 (2019) maps ( Figure 5).
These accuracy points were compared to the matching land cover labels using the VHR imagery in Google Earth Pro. The percentage of agreement was recorded in an accuracy table, and classification accuracy was assessed using a confusion matrix in Microsoft Excel. We then calculated overall accuracy, producer's accuracy, user's accuracy, and kappa coefficients from the confusion matrix [13,58] for PBTC maps. These accuracy points were compared to the matching land cover labels using the VHR imagery in Google Earth Pro. The percentage of agreement was recorded in an accuracy table, and classification accuracy was assessed using a confusion matrix in Microsoft Excel. We then calculated overall accuracy, producer's accuracy, user's accuracy, and kappa coefficients from the confusion matrix [13,58] for PBTC maps.

Bamboo Carbon Stock Assessment
To implement the REDD+ scheme, Cambodia employs a national definition of forest that is consistent with the Global Forest Resources Assessment [60] and considers bamboo to be one of the main forest categories for net carbon sequestration from land use, land use change, and forestry [60]. Conventionally, forest inventories are used for measuring, reporting, and verification (MRV) of deforestation and forest degradation in the tropics to provide the needed information for establishing baseline emissions against which mitigation measures and performance can be assessed [61][62][63].
Equation (2) was applied for calculating carbon stocks in the existing bamboo areas and is presented below [30,60] where TCS is the total carbon stocks of bamboo forest in 2015, 2018, and 2019; BFi, is the area of bamboo forest in district i in Siem Reap Province at any given time; and CS is the carbon stocks in bamboo forest in the respective years (i.e., 2015, 2018, or 2019). Sasaki et al. [30] estimated the bamboo carbon stocks in Cambodia to be 57.4 Mg C ha −1 [64] for all pools except organic soil carbon. Here, we used the same value to calculate the bamboo carbon stocks in 12 districts in Siem Reap Province in 2015, 2018, and 2019 ( Figure 2, Step 5) [30,60].
If natural bamboo forests are completely cleared, total emissions from such clearing were estimated using where CE is carbon emission (Mg CO2), TCS is the total bamboo carbon stock at any given time of clearing (in 2015, 2018, and 2019) in Equation (2), and 44/12 is the molecular weight ratio of carbon dioxide to carbon [30,65].

Bamboo Carbon Stock Assessment
To implement the REDD+ scheme, Cambodia employs a national definition of forest that is consistent with the Global Forest Resources Assessment [60] and considers bamboo to be one of the main forest categories for net carbon sequestration from land use, land use change, and forestry [60]. Conventionally, forest inventories are used for measuring, reporting, and verification (MRV) of deforestation and forest degradation in the tropics to provide the needed information for establishing baseline emissions against which mitigation measures and performance can be assessed [61][62][63].
Equation (2) was applied for calculating carbon stocks in the existing bamboo areas and is presented below [30,60]: where TCS is the total carbon stocks of If natural bamboo forests are completely cleared, total emissions from such clearing were estimated using where CE is carbon emission (Mg CO 2 ), TCS is the total bamboo carbon stock at any given time of clearing (in 2015, 2018, and 2019) in Equation (2), and 44/12 is the molecular weight ratio of carbon dioxide to carbon [30,65].
On the other hand, if all these natural bamboo forests are fully protected and the avoided carbon emissions (ACE) are used to offset the carbon emissions from tourists visiting the province, the contribution of bamboo forest to climate change mitigation was estimated using where COR is carbon offsetting rate (%), T is the total number of tourists visiting Siem Reap in 2015 (2,100,000) and 2018 (2,200,000) [66]. Due to the lack of data in 2019, tourist-based emissions for 2019 were not included in our carbon offsetting calculation, EP is the tourist's per capita emissions for traveling to and staying in Siem Reap (4.14 Mg CO 2 /person) [67]; ACE is the avoided carbon emissions. If 100% is avoided, ACE = CE.

Phenological Profiles and Threshold Values
The time-series analysis of fitted vegetation indices allowed us to quantify the differences in phenological behavior between bamboo, rubber, and evergreen forest land cover categories, as shown in Figures 6-8.
On the other hand, if all these natural bamboo forests are fully protected and the avoided carbon emissions (ACE) are used to offset the carbon emissions from tourists visiting the province, the contribution of bamboo forest to climate change mitigation was estimated using = × 100 (4) where COR is carbon offsetting rate (%), T is the total number of tourists visiting Siem Reap in 2015 (2,100,000) and 2018 (2,200,000) [66]. Due to the lack of data in 2019, tourist-based emissions for 2019 were not included in our carbon offsetting calculation, EP is the tourist's per capita emissions for traveling to and staying in Siem Reap (4.14 Mg CO2/person) [67]; ACE is the avoided carbon emissions. If 100% is avoided, ACE = CE.

Phenological Profiles and Threshold Values
The time-series analysis of fitted vegetation indices allowed us to quantify the differences in phenological behavior between bamboo, rubber, and evergreen forest land cover categories, as shown in Figures 6-8.       Figure 9.
The cumulative average phenology profile of land use categories (Figure 9), namely, bamboo, rubber, and evergreen, were assessed by validating 120 index profiles during the LSP and LFP stages. The rubber plantation vegetation index values (0.784-0.798) were greater than the bamboo values (0.707-0.716) in the SOS (monsoon/rainy season) from May to September, and they were slightly lower (0.741-0.748) right after the start of the mid-LSP stage (January-April). However, the bamboo and rubber vegetation index profiles overlapped in the EOS season between September and October, as shown in Figure 9. The evergreen phenological values and the phase of amplitudes are smaller than those for bamboo and rubber in all phenological stages. Therefore, classification confusion of evergreen forest and that of bamboo or rubber is not likely to occur.
Bamboo was found to have the highest vegetation index in the mid-dry phenological phase. Amplitudes of the minimum and maximum threshold values were determined (0.854-0.882) and are shown in Figure 10. Similarly, rubber plantations showed the highest vegetation index, but the range of the threshold amplitude (0.815-0.841) between bamboo and evergreen forest was smaller (0.581-0.648) ( Figure 10). The cumulative average phenology profile of land use categories (Figure 9), namely, bamboo, rubber, and evergreen, were assessed by validating 120 index profiles during the LSP and LFP stages. Comparing the phenology-based index profiles revealed that bamboo had greater index values over the entire LSP period from 2015 to 2018. The amplitude of bamboo minimum and maximum peak values (0.860-0.863) occurred at the end of the LSP stage (December-February), and the lowest index values (0.707-0.716) were found at the start of the rainy season (SOS) in May-July.
The rubber plantation vegetation index values (0.784-0.798) were greater than the bamboo values (0.707-0.716) in the SOS (monsoon/rainy season) from May to September, and they were slightly lower (0.741-0.748) right after the start of the mid-LSP stage (January-April). However, the bamboo and rubber vegetation index profiles overlapped in the EOS season between September and October, as shown in Figure 9. The evergreen phenological values and the phase of amplitudes are smaller than those for bamboo and rubber in all phenological stages. Therefore, classification confusion of evergreen forest and that of bamboo or rubber is not likely to occur.
Bamboo was found to have the highest vegetation index in the mid-dry phenological phase. Amplitudes of the minimum and maximum threshold values were determined (0.854-0.882) and are shown in Figure 10. Similarly, rubber plantations showed the highest vegetation index, but the range of the threshold amplitude (0.815-0.841) between bamboo and evergreen forest was smaller (0.581-0.648) (Figure 10).

Natural Distribution of Bamboo Forest
Using the obtained threshold values of selected land cover categories, we were able to derive the natural bamboo distribution for the entire Siam Reap Province in 2015 (Figure 11), 2018 ( Figure  12), and 2019 ( Figure 13).
The PBTC map of land use in Siem Reap in 2015 and 2018 (Figures 11 and 12) represents the bamboo, rubber, evergreen, and other land cover distributions. Evergreen forest is widely distributed

Natural Distribution of Bamboo Forest
Using the obtained threshold values of selected land cover categories, we were able to derive the natural bamboo distribution for the entire Siam Reap Province in 2015 (Figure 11), 2018 (Figure 12), and 2019 ( Figure 13).
The PBTC map of land use in Siem Reap in 2015 and 2018 (Figures 11 and 12) represents the bamboo, rubber, evergreen, and other land cover distributions. Evergreen forest is widely distributed around Kulen Mountain and the Angkor Wat temple. Bamboo distribution was found largely in the mountain region, and some patches can be observed around the Angkor Wat temple. Rubber distribution can be seen in only a small area ( Figure 11C) in 2015, whereas the increase in rubber is clearly seen in 2018 ( Figure 12C).
In Figure 11A, the map represents the distribution of bamboo around Kulen Mountain and Kulen National Park. Figure 11B (Boeng Peae Wildlife Sanctuary) shows the distribution of bamboo in natural evergreen forests in the northwestern part of the study region, including Chi Kraeng District. Figure 11C shows a mix of vegetation of evergreen forest (green) and rubber (purple) in Chi Kraeng District. Figure 11D shows the rich evergreen forest cover and small bamboo distribution around the Angkor Wat temple. Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 24 In Figure 11A, the map represents the distribution of bamboo around Kulen Mountain and Kulen National Park. Figure 11B (Boeng Peae Wildlife Sanctuary) shows the distribution of bamboo in natural evergreen forests in the northwestern part of the study region, including Chi Kraeng District. Figure 11C shows a mix of vegetation of evergreen forest (green) and rubber (purple) in Chi Kraeng District. Figure 11D shows the rich evergreen forest cover and small bamboo distribution around the Angkor Wat temple. As shown in Figure 13, bamboo distribution was higher in 2019 compared with that in 2015, and the distribution of rubber plantations can be seen clearly in Varin and Chi Kraeng Districts (Figures 12C and 13C). Figure 13 shows a higher bamboo distribution area compared with that in the Landsat OLI maps (2015 and 2018) in most of the study region. On the other hand, rubber plantation distribution can be seen clearly in Varin and Chi Kraeng Districts ( Figure 13C). Our analysis found that the bamboo cover was 2674 ha in 2015, 2711 ha in 2018, and 2672 ha in 2019. A larger area of bamboo was in upland areas and along riverbanks [42], especially in Phnom Kulen National Park and the northern part of Chi Kraeng District.

Area of Bamboo Forest and Its Carbon Stocks by Districts
We

Area of Bamboo Forest and Its Carbon Stocks by Districts
We estimated a total bamboo area of 2674 ha in 2015, 2711 ha in 2018, and 2672 ha in 2019. Carbon stocks in bamboo forest were estimated at 153,511; 155,583; and 153,367 Mg C in 2015, 2018, and 2019, respectively. Specifically, among the 12 districts, the Banteay Srei District had a higher amount of bamboo carbon stocks (47,500 Mg C) between 2015 and 2019, and Svay Leu and Varin Districts had more than 41,600 Mg C in the same period. Other districts covered smaller portions of bamboo, and their carbon stocks were less than 3000 Mg C (Table 4). Siem Reap is a popular tourist destination in Cambodia because of its rich historical temple attraction, Angkor Wat heritage, and natural attractions. Siem Reap welcomed 2.

Discussion
Our method has limitations in identifying bamboo in dense and tall vegetation areas such as evergreen forest, as well as in smaller bamboo patches in other forests. There are possible errors in our study due to data limitation and the occurrence of natural stands of bamboo. Limited groundbased data also constrains potential accuracy. Different stages of bamboo growth in the study areas were obtained through remotely sensed temporal data. Different species at various growing stages

Discussion
Our method has limitations in identifying bamboo in dense and tall vegetation areas such as evergreen forest, as well as in smaller bamboo patches in other forests. There are possible errors in our study due to data limitation and the occurrence of natural stands of bamboo. Limited ground-based data also constrains potential accuracy. Different stages of bamboo growth in the study areas were obtained through remotely sensed temporal data. Different species at various growing stages show different phenological behaviors [20], and their spectral values vary during phenological stages [19,20]. Nath et al. [20] found that ground-based data are needed to reduce the possible errors while increasing the map accuracy. For example, in Sichuan Province, China, where small patch of bamboo could be identified through the use of ground-based data and k-nearest neighbor classifier for the Worldview-2 imagery [21]. Even in such a difficult situation, they achieved a bamboo map producer's accuracy of 82.65% and a user's accuracy of 93.10%.
The second source of error could be the distribution of bamboo mixed with other vegetation types. During the field survey, we found that other lowland areas in Siem Reap have small bamboo patches, usually along riverbanks and intermingled with tall, dense areas of evergreen and deciduous forests. This makes it difficult to identify bamboo from moderate-resolution satellite imagery, which might account for the omission errors in Table 3 and Figure 16. When observing Figures 9 and 10, the rubber category is the main source of confusion for the bamboo category because the ranges of the minimum and maximum threshold values are close to each other. To avoid this problem, big data [35] and deep learning-based technology [68,69] may be employed. In addition, because bamboo is usually mixed with evergreen forest in most of the study region, selection of ground-reference points for accuracy assessment without ground data can create confusion between bamboo and evergreen forests. If such incidence occurs, evergreen or bamboo can be mistakenly classified into just one category. One possible solution would be to use very high resolution imagery such as Planet's daily imagery for bamboo phenological assessment and IKONOS, QuickBird, WorldView [21], GeoEye, and Pleiades, with 0.5 meter or finer spatial resolutions, or to acquire ground-based inventory data [30] or imagery from advanced unmanned aerial vehicles (UAV) [70,71]. Apart from the problems stated earlier, other possible errors when classifying the bamboo forests could occur if we want to classify the bamboo distribution by ages. Chen et al. [28] and Langner et al. [72] noted the difficulty in obtaining an integrated phenological cycle for young and mature bamboo using moderate-resolution remote sensing data.
Errors in estimating carbon stocks in bamboo forest could arise from the use of generalized data of carbon stocks per unit area from Sasaki et al. [30] which was based on average data for the whole When observing Figures 9 and 10, the rubber category is the main source of confusion for the bamboo category because the ranges of the minimum and maximum threshold values are close to each other. To avoid this problem, big data [35] and deep learning-based technology [68,69] may be employed. In addition, because bamboo is usually mixed with evergreen forest in most of the study region, selection of ground-reference points for accuracy assessment without ground data can create confusion between bamboo and evergreen forests. If such incidence occurs, evergreen or bamboo can be mistakenly classified into just one category. One possible solution would be to use very high resolution imagery such as Planet's daily imagery for bamboo phenological assessment and IKONOS, QuickBird, WorldView [21], GeoEye, and Pleiades, with 0.5 meter or finer spatial resolutions, or to acquire ground-based inventory data [30] or imagery from advanced unmanned aerial vehicles (UAV) [70,71]. Apart from the problems stated earlier, other possible errors when classifying the bamboo forests could occur if we want to classify the bamboo distribution by ages. Chen et al. [28] and Langner et al. [72] noted the difficulty in obtaining an integrated phenological cycle for young and mature bamboo using moderate-resolution remote sensing data.
Errors in estimating carbon stocks in bamboo forest could arise from the use of generalized data of carbon stocks per unit area from Sasaki et al. [30] which was based on average data for the whole country. To improve the accuracy of carbon stocks assessment, field measurement of the targeted bamboo forest in the study region could be a solution, but this is costly and time consuming.

Conclusions
Using PBTC method and GEE with remote sensing imagery from Landsat 8 OLI in 2015, 2016, 2017, and 2018, we were able to determine the minimum and maximum vegetation thresholds for bamboo forest, evergreen forest, and rubber plantation in Siem Reap province, Cambodia. With these thresholds, we could classify and map the distribution of natural bamboo forest in 2015 and 2018. We achieved the overall Landsat 8 OLI bamboo maps for 2015 and 2018 UAs of 86.6% and 87.9% and PAs of 95.7% and 97.8%, respectively. Similarly, we achieved UAs of 86.5% and PAs of 91.7% for Sentinel-2 imagery for 2019. Accordingly, spatial distribution of carbon stocks by districts in the Siem Reap was estimated and mapped across the province. We found that area of bamboo was 2674 ha in 2015, 2711 ha in 2018, and 2672 ha in 2019. Total carbon stocks in the respective years were 153,511, 155,583, and 153,367 Mg C.
It can be concluded that GEE can be used to classify and map the natural distribution of the natural bamboo forest and related carbon stocks with the aid of PBTC method, multi-spectral Sentinal-2 and Landsat 8 OLI imagery. However, using very high-resolution imagery, UAV, big data technology, and deep learning could improve the misclassification during the closer ranger of threshold values of various land cover categories. To avoid confusion between natural bamboo and evergreen forests, it is important to select the ground reference points through field visit, forest measurement, and/or using UAV.
If bamboo forest is fully protected, the avoided emissions could offset about 6% of carbon emissions by tourists visiting the province, indicating that bamboo can also play an important role in climate change mitigation in addition to providing non-timber forest products to the local people. As bamboo forest has multiple roles in local livelihood improvement as well as climate change mitigation, management of this forest through proper mapping and zoning could increase awareness of the benefits of the bamboo forests, which could help promote effective management for long-term sustainable use.