Specific Bamboo Forest Extraction and Long-Term Dynamics as Revealed by Landsat Time Series Stacks and Google Earth Engine

Understanding the spatiotemporal dynamics of bamboo forests is of critical importance as it characterizes the interaction between forest and agricultural ecosystems and provides essential information for sustainable ecosystem management and decision-making. Thus far, the specific dynamics of moso bamboo (Phyllostachys edulis) and other bamboo are still unknown. In this study, we used temporal information extracted from Landsat time series stacks with Google Earth Engine (GEE) to characterize the spatiotemporal dynamics of bamboo forests, including moso bamboo and other bamboo, in Lin’an County, China, from 2000 to 2019. The bamboo forests were mapped in four periods: the early 2000s (2000–2004), the late 2000s (2005–2009), the early 2010s (2010–2014), and the late 2010s (2015–2019). The overall accuracy of these maps ranged from 97% to 99%. We then analyzed the spatiotemporal dynamics of the bamboo forests at the county and subdistrict/township scales, and probed the bamboo forest gain and loss with respect to the terrain features. Our findings show that bamboo forests increased by 4% from 2000 to 2014, followed by a sharp decrease of 13% in the late 2010s. The decrease was mainly caused by the loss of other bamboo. Approximately 69% of the bamboo forest gain occurred in non-bamboo forest areas, and the rest occupied non-forest areas. Bamboo forest loss was mainly due to conversion into orchard (59%) and forest plantation (22%). Compared to bamboo forest gain, bamboo forest loss was typically observed in areas with lower elevations and steeper slopes. Our study offers spatially explicit and timely insight into bamboo forest changes at the regional scale. The derived maps can be applied to study the drivers, consequences, and future trends of bamboo forest dynamics, which will contribute to sustainable ecosystem management.


Introduction
Intensive forestry practices in the 21st century have resulted in dramatic changes to the world's forested areas, especially subtropical forests [1]. Bamboo is one of the major types of forest in subtropical and tropical regions. Because of its economic value, bamboo has a long history of being cultivated as a tree crop. According to the National Forest Inventory (NFI) in China, the total bamboo forest area increased by 110.86% (3.04 to 6.41 million ha) from 1973 to 2018 [2]. Recently, a continuous long-term expansion of bamboo forests has generally occurred in mountainous areas during the intense forestry process, influencing forest and agricultural land cover changes as well as causing According to the FRI for Lin'an in 2017, broadleaf, bamboo, conifer, and hickory (Carya cathayensis) plantation areas contributed approximately 31%, 21%, 20% and 16%, respectively, of the entire vegetation coverage. Other contributing vegetation included broadleaf-conifer mixed forest, shrub, and cropland. The spatial distribution of vegetation types is shown in Supplementary Figure S1.

Time Series Landsat Images
In this study, 30-m Landsat surface reflectance (SR) images (path/row: 119/039, 120/039) acquired from January 2000 to December 2019 were used. All available SR images, 1003 in total, derived from the Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) instruments, were selected in GEE. These SR products were developed by the U.S. Geological Survey (USGS, http://earthexplorer.usgs.gov/) and were atmospherically corrected using the Landsat Ecosystem According to the FRI for Lin'an in 2017, broadleaf, bamboo, conifer, and hickory (Carya cathayensis) plantation areas contributed approximately 31%, 21%, 20% and 16%, respectively, of the entire vegetation coverage. Other contributing vegetation included broadleaf-conifer mixed forest, shrub, and cropland. The spatial distribution of vegetation types is shown in Supplementary Figure S1.

Time Series Landsat Images
In this study, 30-m Landsat surface reflectance (SR) images (path/row: 119/039, 120/039) acquired from January 2000 to December 2019 were used. All available SR images, 1003 in total, derived from the Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) instruments, were selected in GEE. These SR products were developed by the U.S. Geological Survey (USGS, http://earthexplorer.usgs.gov/) and were atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Landsat 5 and 7) [22] and the Land Surface Reflectance Code (LaSRC) algorithm (Landsat 8). In addition, these images were accompanied by cloud, shadow and snow masks [23] for individual pixels, leaving the remaining clear observations.
Supplementary Figure S2 summarizes the numbers of Landsat 5/7/8 images used in this study. The distributions of these available scenes by month, by sensor and by path/row were all uneven. The pixels for some months did not have valid observations due to masks, the failure of Landsat 7, and other reasons.

Supplementary Data
The FRI relies mainly on field surveys and is updated at ten-year intervals to provide the dominant vegetation distributions in vector polygon format. The latest FRI map of Lin'an was completed in the summer of 2017. Seventeen major land cover categories are recorded, for instance, arbor, bamboo, shrub, cropland, water body, built-up land, logged land and burned area. This map also provides detailed information on the spatial distributions of various types of bamboo forests, including moso bamboo, Lei bamboo, Gaojie bamboo and so on.
We designed 10 land cover categories, which were extracted from the FRI map and are referred to as the initial classes ( Table 1). The training and validation points of interest (POIs) were obtained from the detailed FRI map for 2017 and Google Earth (GE) images acquired from 2000-2020. For each period, we randomly generated 6500 POIs for different homogenous classes based on a stratified random sampling design [24]. The FRI map served as a critical reference for vegetation types. Each POI in the training and validation datasets was carefully examined by visual interpretation of the continuous high-resolution images from GE. On the whole, all the procedures aimed to ensure that the POIs were unmixed, stable and correctly classified during the four periods. Because relatively few GE images from 2000-2010 were available for our study area, POIs for the 2000s were also collected using Landsat time series images.  [25] in GEE. The terrain features (elevation and slope) were calculated from the DEM. In our study area, the elevation ranges from 5 to 1771 m above sea level (Figure 1), and the slope is below 80 • . Figure 2 presents our framework for analyzing the dynamics of two types of bamboo forests over the past 20 years. First, temporal features were extracted from monthly composites of five spectral indices. Second, the extracted temporal features were processed with a random forest classifier to map moso bamboo and other bamboo for each period. Third, we analyzed and revealed the spatiotemporal dynamic characteristics of the bamboo forests.

Methods
Figure 2 presents our framework for analyzing the dynamics of two types of bamboo forests over the past 20 years. First, temporal features were extracted from monthly composites of five spectral indices. Second, the extracted temporal features were processed with a random forest classifier to map moso bamboo and other bamboo for each period. Third, we analyzed and revealed the spatiotemporal dynamic characteristics of the bamboo forests.

Temporal Feature Extraction
The temporal features of vegetation, regarded as signals of plant growth, usually have a periodic pattern and vary with species. These features can be extracted from time series spectral information.
In this study, a total of five spectral indices were calculated for individual pixels from January 2000 to December 2019 in GEE, namely, the normalized difference vegetation index (NDVI) [26], the normalized difference moisture index (NDMI) [27], the enhanced vegetation index (EVI) [28], the normalized difference water index 1 (NDWI1) [29] and the normalized difference water index 2 (NDWI2) [30]. These indices are generally used for vegetation classification, and some have been shown to be useful for bamboo mapping specifically [3,4,17]. NDVI and EVI are sensitive to vegetation greenness. NDMI has the potential to detect the vegetation canopy water content [27]. NDWI1 and NDWI2 are related to water, bare soil, and vegetation status [29][30][31]. The equations for these indices are as follows:

Temporal Feature Extraction
The temporal features of vegetation, regarded as signals of plant growth, usually have a periodic pattern and vary with species. These features can be extracted from time series spectral information.
In this study, a total of five spectral indices were calculated for individual pixels from January 2000 to December 2019 in GEE, namely, the normalized difference vegetation index (NDVI) [26], the normalized difference moisture index (NDMI) [27], the enhanced vegetation index (EVI) [28], the normalized difference water index 1 (NDWI1) [29] and the normalized difference water index 2 (NDWI2) [30]. These indices are generally used for vegetation classification, and some have been shown to be useful for bamboo mapping specifically [3,4,17]. NDVI and EVI are sensitive to vegetation greenness. NDMI has the potential to detect the vegetation canopy water content [27]. NDWI1 and NDWI2 are related to water, bare soil, and vegetation status [29][30][31]. The equations for these indices are as follows: where B, G, R, NIR, and SWIR1 are the Landsat SR values in the blue, green, red, near infrared, and short-wave infrared 1 bands, respectively. When there were fewer than four valid observations for a given month, a situation referred to as an unclear month, the corresponding index was replaced by the mean value of all observations within the interquartile range for the previous and subsequent months. In some extreme cases where there were 3 consecutive unclear months, we replaced the corresponding index with the mean value of all observations within the interquartile range for the month after next and the month before last. In this way, we were able to ensure that the monthly composite values for each pixel were the best representative values for each month.

Random Forest Classification
Random forest (RF) classifiers [32] were calibrated on the training samples (Supplementary  Table S1) using the extracted temporal features as the input features to create our bamboo forest maps for each of the four periods separately. An RF classifier is a highly efficient classifier for a large amount of data and is robust against imbalanced input features [33]. Conveniently, only two main parameters need to be set for RF classification: the number of classification trees and the number of selected input features to be used to split each node. To achieve a suitable balance between classification accuracy and computation time, we selected an RF model with 200 classification trees and the square root of the total number of input features for splitting after various tests of these parameters for our study area. The initial output of the RF classifier for each period was a 10-class (Table 1) land cover map, including moso bamboo and other bamboo. Then, the initial classes were reclassified into 4 classes, referred to as the reclassified classes (Table 1). We used the maps based on these reclassified classes, i.e., the bamboo forest maps, for further analysis.

Accuracy Assessment and Area Comparison
The accuracy of bamboo forest mapping was then assessed using the confusion matrix based on the validation samples for each period (Supplementary Table S1). The F1-score was calculated for the 10 initial classes as shown in Equation (6). The frequently employed kappa indices have been criticized for being misleading or flawed for map production [34], and therefore, were not used. The Remote Sens. 2020, 12, 3095 7 of 16 overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA), allocation disagreement (AD), and quantity disagreement (QD) [34] were calculated for the 4 reclassified classes.
Furthermore, to evaluate the area accuracy of our classification results, we compared the areas of the two types of bamboo forests derived from our results for the late 2010s against the FRI map from 2017 at the subdistrict/township scales.

Analyzing the Spatiotemporal Dynamics of Bamboo Forests
To further analyze the dynamics of bamboo forests, we reclassified all land cover categories into non-bamboo forest, bamboo forest and non-forest. Based on the observed trends, land cover changes were divided into five types: (1) stable: the area of non-bamboo forest/bamboo forest/non-forest was unchanged, or only intra-class changes (e.g., from cropland to orchard) occurred; (2) bamboo forest gain: the area of non-bamboo forest or non-forest was changed to bamboo forest, or the area of bamboo forest was changed to non-forest and then back to bamboo forest (e.g., from other bamboo to cropland to other bamboo); (3) non-forest gain: the area of non-bamboo forest or bamboo forest was changed to non-forest, or the area of non-bamboo forest was changed to bamboo forest and then back to non-forest; (4) forest plantation: the area of bamboo forest or non-forest was changed to non-bamboo forest, or the area of non-bamboo forest was changed to bamboo forest/non-forest and then back to non-bamboo forest; and (5) unknown: areas were changed in an uncertain manner as a result of classification errors.

Accuracy Assessment and Area Comparison
The assessment results revealed that the accuracy and robustness were reasonably satisfactory for most initial classes (Supplementary Table S2). The F1-scores for moso bamboo and other bamboo were all greater than 90% in the four periods ( Table 2). The major classification errors were commission errors for other bamboo and mixed forest and for moso bamboo, mixed forest, and shrub (Supplementary  Table S2). Furthermore, Table 2 also shows excellent classification performance in terms of the UA and PA for the reclassified classes in the different periods. High OAs of 97%, 99%, 98%, and 97%, low ADs of 0.02, 0.00, 0.00, and 0.02, and QDs of 0.01, 0.01, 0.02, and 0.01were achieved for the reclassified classes in the early 2000s, late 2000s, early 2010s and late 2010s, respectively. In most periods, the UAs for moso bamboo were slightly higher than the PAs, suggesting larger errors of omission than commission. By contrast, the UAs for other bamboo were lower than the PAs. These results imply that these bamboo forest maps were comparable and could satisfactorily reflect the dynamics of the bamboo forests in the study area over the past 20 years.
We compared our late 2010s bamboo forest map with the bamboo forest map obtained from the FRI for 2017 (Figure 3). The spatial distributions of the four reclassified classes in these two products were in visually satisfactory agreement (Figure 3a,b). At the subdistrict/township level, our classification results showed good agreement with the FRI, with acceptable R 2 values for both moso bamboo (0.72) and other bamboo (0.68) (Figure 3d,e). Nevertheless, the slopes of 0.91 and 1.36 revealed a slight difference between these two products. Therefore, two zoom-in regions (black boxes in Figure 3a,b) were selected for analysis. The comparison showed some discrepancies between the two types of bamboo forests in the case regions. Therefore, we selected another two zoom-in views (white boxes in Figure 3g,i,k,m) for further analysis based on GE images. Figure 3h,l,j,n show the distributions of moso bamboo and other bamboo in the zoom-in views extracted from our map and the FRI map, respectively. For the background of the figure panels, we used the same images from GE, dated 3 November 2017 (Figure 3h,j) and 9 November 2017 (Figure 3l,n), respectively. This comparison implies that the FRI map is missing/containing some other bamboo areas in cropland areas, whereas our map contains some classification errors for moso bamboo caused by mixed pixels. These may result in some differences in area estimates between these two maps.  Figure 4 presents the spatial distributions of the two types of bamboo forests (moso bamboo in red and other bamboo in yellow), non-bamboo forest, and non-forest in the late 2010s. Recently, bamboo forests were widely distributed in eastern Lin'an, such as in the BQ, THY and TMS townships. Moso bamboo was present mainly in four clusters, in northeastern, southeastern, northwestern and southwestern Lin'an. The northeastern cluster mainly covered GH and QSH. The southeastern cluster included LL, JN and BQ. The southwestern cluster mainly covered TK and QC, while the northwestern cluster mainly covered LG and CH. By contrast, other bamboo was located in almost all subdistricts/townships. Critical percentages of moso bamboo (52%) and other bamboo (59%) were concentrated at medium elevations (Figure 4c). Compared with other bamboo clusters,  Figure 4 presents the spatial distributions of the two types of bamboo forests (moso bamboo in red and other bamboo in yellow), non-bamboo forest, and non-forest in the late 2010s. Recently, bamboo forests were widely distributed in eastern Lin'an, such as in the BQ, THY and TMS townships. Moso bamboo was present mainly in four clusters, in northeastern, southeastern, northwestern and southwestern Lin'an. The northeastern cluster mainly covered GH and QSH. The southeastern cluster included LL, JN and BQ. The southwestern cluster mainly covered TK and QC, while the northwestern cluster mainly covered LG and CH. By contrast, other bamboo was located in almost all subdistricts/townships. Critical percentages of moso bamboo (52%) and other bamboo (59%) were concentrated at medium elevations (Figure 4c). Compared with other bamboo clusters, moso bamboo clusters were observed on relatively steep slopes (Figure 4d). Most of the bamboo forests were adjacent to rural settlements, with other bamboo being generally planted on cropland.

Bamboo Forest Map for the Late 2010s
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 16 moso bamboo clusters were observed on relatively steep slopes (Figure 4d). Most of the bamboo forests were adjacent to rural settlements, with other bamboo being generally planted on cropland.      Figure 6 provides an insight into the characteristics of the dynamics of bamboo forests in terms of the types of land cover changes occurring over the past two decades. Most areas remained unchanged, or only intra-class changes occurred (stable, 70.2%), while the least amount of area was changed to non-bamboo forest (forest plantation, 2.9%). In addition, bamboo forest gain, non-forest gain, and unknown accounted for 7.3%, 16.0%, and 3.5% of the total area, respectively. For bamboo forest gain, changes from non-bamboo forest (69%) were over twice as prevalent as changes from non-forest (31%). Note that approximately 91% of these non-forest areas were changed into areas of other bamboo. Changes from bamboo forest to non-forest and from non-bamboo forest/non-forest to bamboo forest and then back to non-forest contributed approximately 63% of the non-forest gain. Moreover, changes from bamboo forest to non-bamboo forest and from non-bamboo forest to bamboo forest and then back to non-bamboo forest accounted for approximately 79% of forest  Figure 6 provides an insight into the characteristics of the dynamics of bamboo forests in terms of the types of land cover changes occurring over the past two decades. Most areas remained unchanged, or only intra-class changes occurred (stable, 70.2%), while the least amount of area was changed to non-bamboo forest (forest plantation, 2.9%). In addition, bamboo forest gain, non-forest gain, and unknown accounted for 7.3%, 16.0%, and 3.5% of the total area, respectively. For bamboo forest gain, changes from non-bamboo forest (69%) were over twice as prevalent as changes from non-forest (31%). Note that approximately 91% of these non-forest areas were changed into areas of other bamboo. Changes from bamboo forest to non-forest and from non-bamboo forest/non-forest to bamboo forest and then back to non-forest contributed approximately 63% of the non-forest gain. Moreover, changes from bamboo forest to non-bamboo forest and from non-bamboo forest to bamboo forest and then back to non-bamboo forest accounted for approximately 79% of forest plantation. Figure 6c shows the areas of the different types of changes (bamboo forest gain, non-forest gain and forest plantation) in the different periods. There was a declining trend in the bamboo forest gain throughout the whole period.

Spatiotemporal Dynamics of Bamboo Forests
gain and forest plantation) in the different periods. There was a declining trend in the bamboo forest gain throughout the whole period.
The terrain features of the bamboo forest dynamics from 2000 to 2019 were analyzed along both elevation and slope gradients (Figure 6d,e). The results reveal that bamboo forest gain was concentrated in areas with medium elevations or moderate slopes. Specifically, the newly planted bamboo forests in non-bamboo forest areas tended to be at higher elevations and on steeper slopes than those in non-forest areas (Supplementary Figure S3). By contrast, an obvious bamboo forest loss was observed in non-bamboo forest and non-forest areas at low and medium elevations (Supplementary Figure S3a). Non-bamboo forest areas with moderate slopes and non-forest areas with gentle slopes widely experienced bamboo forest loss during 2000-2019 (Supplementary Figure  S3b).  Figure 6a (the same as in Figure 5a). (c) shows the areas of different types of changes (bamboo forest gain, non-forest gain and forest plantation) in the different periods. (d) and (e) show the area dynamics along elevation and slope gradients throughout the whole period.  Figure 6a (the same as in Figure 5a). (c) shows the areas of different types of changes (bamboo forest gain, non-forest gain and forest plantation) in the different periods. (d,e) show the area dynamics along elevation and slope gradients throughout the whole period.
The terrain features of the bamboo forest dynamics from 2000 to 2019 were analyzed along both elevation and slope gradients (Figure 6d,e). The results reveal that bamboo forest gain was concentrated in areas with medium elevations or moderate slopes. Specifically, the newly planted bamboo forests in non-bamboo forest areas tended to be at higher elevations and on steeper slopes than those in non-forest areas (Supplementary Figure S3). By contrast, an obvious bamboo forest loss was observed in non-bamboo forest and non-forest areas at low and medium elevations (Supplementary Figure S3a). Non-bamboo forest areas with moderate slopes and non-forest areas with gentle slopes widely experienced bamboo forest loss during 2000-2019 (Supplementary Figure S3b).

Discussion
Bamboo forest mapping is vital to forestry and agricultural resource management and biodiversity. The rapid forestry process of bamboo forests in China's subtropical regions has led to frequent conversion on agriculture-forest frontiers. This study has attempted to reveal the spatiotemporal dynamic characteristics of two typical types of bamboo forests at a regional scale over the past two decades, based on a combination of Landsat time series stacks, the GEE platform and a machine learning method. The adopted methodology exhibits the potential to be repeatedly used for bamboo forest time series analyses and to be extended to the extraction of other plant species by adjusting the features input into the classifier. The generated maps clearly show that bamboo forest gain and loss, indicating intensive forestry practices, were observed in all subdistricts/townships in the study area during 2000-2019.
The dynamics of bamboo forests represent a good example of the changes in cash crops occurring in subtropical areas. In our study area, one of the core tasks of local forestry production in the early 2000s was to promote the development of both economic forests and ecological forests [35]. By 2019, the total value of the agricultural output of Lin'an had increased from 1.86 billion yuan in 2003 to 6.69 billion yuan, of which the bamboo shoot, hickory and fruit industries accounted for 33.3%, 18.6%, and 8.5%, respectively [35,36]. These observations suggest that local policy is a fundamental determinant of cash crop growth and distribution. As our results show, approximately 159 km 2 of non-bamboo forest area and 50 km 2 of non-forest area were changed into bamboo forest area over the past two decades. These phenomena are consistent with previous studies, which have reported a clear impact of cash crop expansion on both forestry and agricultural areas in subtropical China [3,37,38]. Bamboo forest gain and distribution mainly occurred in relatively low-elevation regions with gentle and moderate slopes (Figures 4 and 6). Approximately 11% of the cropland in Lin'an experienced bamboo forest gain from 2000 to 2019. This finding indicates that some natural habitats and croplands with high soil quality were probably occupied by newly planted bamboo forests, challenging the sustainability of biodiversity conservation and grain production [39]. Moreover, a considerable proportion of bamboo forests were converted into orchard (constituting approximately 59% of bamboo forest loss) since orchard had a higher economic value than bamboo forests. Economic gain plays a key role in land conversion [38].
Remarkably, approximately 30% of bamboo forest gain was observed in non-bamboo forest areas with steep slopes. According to our results, in the late 2010s, a great number of bamboo forests (27%) and orchard (41%) were planted on steep slopes. Because few understory plants grow in these regions, such expansion of these cash crops will lead to soil erosion and landscape fragmentation [40,41]. This implies that the terrain features of cash crop areas should be considered in forest management. Although an increasing attention has been paid to ecosystem health, the increased extent of forest plantation is still not comparable to the expansion of cash crops. Our study reveals a rapid expansion of these cash crops in Lin'an during 2000-2019. Our results therefore call for a better balance between promoting the forestry production of cash crops and protecting natural resources.
In this study, only spectral indices were used in the RF classifier. Ideally, the spectral information used for classification should have spatial and temporal consistency. However, the uneven quantity and quality of Landsat images pose great challenges for land monitoring applications [42]. In this study, approximately 50% of the pixels were affected by at least one unclear month over the four periods (Supplementary Figure S4). Although SR products have shown small differences between the comparable sensor bands and indices [43], no transformation was applied to integrate the data in our study. Instead, to alleviate the uncertainties, we calculated the mean value of all observations within the interquartile range for each month. The errors in our maps can also be attributed to the limited availability of reference data. Due to the lack of available high-resolution images acquired during 2000-2010, some POIs for this time period were collected using historical Landsat images. It was difficult to select unmixed pixels for reference across the complex landscape at a 30 m spatial resolution. Nevertheless, even with these multiple sources of uncertainties, high classification accuracies were achieved in the bamboo forest maps. The results also illustrate the robustness of the classifier against some uncertainties in the composite images. At present, limited at both the spatial and temporal scales, this study has produced bamboo forest maps of Lin'an for 2000 to 2019. The Landsat SR data can serve as an improved data source for long-term monitoring, and the GEE platform offers the flexibility to collect and process related data from the 1980s to the present. Thus, maps of the long-term bamboo forest dynamics since the 1980s could be produced at large scales in the future to conduct relevant studies.

Conclusions
Accurate spatiotemporal distribution data for bamboo forests can be pivotal for supporting sustainable ecosystem management and biodiversity conservation. To date, long-term dynamic maps of moso and other bamboo forests at the regional scale have been unavailable, limiting our ability to understand the spatiotemporal dynamic characteristics of bamboo forests and to predict their dynamic trends. In our study, we have addressed this problem by using a pixel-based machine learning method to map these two types of bamboo forests in a subtropical area from 2000 to 2019 on the basis of 1003 Landsat 5/7/8 images and the GEE platform. The bamboo forest maps were generated for five-year intervals over the past two decades. We then analyzed the spatiotemporal dynamics of the bamboo forests at the county and subdistrict/township levels and probed the bamboo forest gain and loss along both elevation and slope gradients. In future studies, it will be necessary to trace the bamboo forest dynamics in subtropical areas back to the 1980s and to investigate the causes and consequences of the bamboo forest dynamics inferred from the derived maps.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/18/3095/s1, Figure S1: Spatial distribution of vegetation types in Lin'an, which was developed from the Forest Resources Inventory (FRI) product in 2017, Figure S2: The annual numbers of Landsat surface reflectance images used in this study by (a) month, (b) sensor, and (c) path/row, Figure S3: The area dynamics of the different types of changes along elevation (a) and slope (b) gradients throughout the whole period, Figure S4: Percentages of valid observations in the four periods, Table S1: The numbers of POIs in the training and validation datasets for different periods, Table S2: Error matrices of F1-score values (%) for the initial classes in the four periods.

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