Mapping Rubber Plantations and Natural Forests in Xishuangbanna (Southwest China) Using Multi-Spectral Phenological Metrics from MODIS Time Series

We developed and evaluated a new approach for mapping rubber plantations and natural forests in one of Southeast Asia’s biodiversity hot spots, Xishuangbanna in China. We used a one-year annual time series of Moderate Resolution Imaging Spectroradiometer (MODIS), Enhanced Vegetation Index (EVI) and short-wave infrared (SWIR) reflectance data to develop phenological metrics. These phenological metrics were used to classify rubber plantations and forests with the Random Forest classification algorithm. We evaluated which key phenological characteristics were important to discriminate rubber plantations and natural forests by estimating the influence of each metric on the classification accuracy. As a benchmark, we compared the best classification with a classification based on the full, fitted time series data. Overall classification accuracies derived from EVI and SWIR time series alone were 64.4% and 67.9%, respectively. Combining the phenological metrics from EVI and SWIR time series improved the accuracy to 73.5%. Using the full, smoothed time series data instead of metrics derived from the time series improved the overall accuracy only slightly (1.3%), indicating that the phenological metrics were sufficient to explain the seasonal changes captured by the MODIS time series. The results demonstrate a promising utility of phenological metrics for mapping and monitoring rubber expansion with MODIS.


Introduction
Mapping land cover is one of the key applications of remote sensing [1].The increased availability of broad scale Earth observation data together with recent developments in multi-temporal analyses techniques have increased the quality of continental to global land cover maps [2], global forest cover maps [3], and global maps of cropland extension [4].The Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), the Medium Resolution Imaging Spectrometer (MERIS), and Satellite Pour l'Observation de la Terre Vegetation (SPOT VEGETATION) capture images of the globe at moderate spatial but very high temporal resolution (2-3 days for MERIS and SPOT VEGETATION sensors; daily for AVHRR and MODIS sensors).This high temporal resolution facilitates monitoring dynamic inter-and intra-annual processes on the Earth surface, which would not be observable using less frequent Earth observation data.This phenological information supports mapping recent land cover and monitoring land cover changes.
The conversion of natural and semi-natural forests to rubber plantations (Hevea Brasiliensis Muell.Arg.) has become a significant land-use change process during the last decade.The conversion occurs throughout the tropical and sub-tropical region, but it is especially prevalent in Southeast Asia [5].Transforming natural forests to rubber plantations has significant ecological impacts on water balance, carbon cycle, and biodiversity [6].In some regions, rubber has replaced traditional subsidence farming [7] and therefore completely changed local economic structure [8].Accurate mapping and monitoring of rubber plantations is thus of great importance to quantify and project the ecological and economic impacts of rubber expansion.
The need to develop methods for improved monitoring of rubber plantations with remote sensing has been recognized already in several studies.Most studies have used high spatial resolution, multi-spectral sensors such as Landsat or hyper-spectral sensors [9][10][11].A high spatial resolution is a clear advantage for capturing the fine spatial detail of many land-use processes.However, extensive cloud cover and limited data availability often diminish the utility of Landsat-like sensors for mapping large tropical areas.In comparison, coarse resolution sensors like MODIS provide data at a higher temporal frequency, i.e., every day, and over larger areas.
MODIS' high temporal resolution not only increases the chance of cloud-free observations but also permits a detailed temporal record (signature) of seasonal vegetation patterns.These temporal signatures can be useful to discriminate vegetation and land cover types that are spectrally distinct only during short periods of time throughout the year.Time series of MODIS Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) have been successfully used to characterize vegetation types in different environmental settings [12][13][14].For example, MODIS Vegetation Index (VI) time series have improved the classification of abandoned farmland, different crop types, and semi-arid vegetation by capturing the specific phenological pattern of each land cover type.Other studies developed phenological metrics from NDVI or EVI time series to describe patterns in vegetation phenology, e.g., length and peak of the vegetation season, and to use this information for mapping different land cover types [15][16][17], including different forest types [18].
To date, only few studies have tried to map rubber plantations with MODIS data [19,20].Li and Fox [19] used MODIS Enhanced Vegetation Index (EVI) time series in combination with sub-national statistical data on rubber growth to map rubber across Southeast Asia.As statistical data they used information on the area of land covered by rubber plantations, on the amount of land in rubber production, and on the total amount of latex production at different sub-national administrative units.Their results proved to identify mature rubber plantations with a producer's/user's accuracy of 67.0%/98.1% and young rubber plantations with a producer's/user's accuracy of 59.4%/97.2%.When using MODIS data alone, the producer's/user's accuracies of mature and young rubber plantations decreased to 60.9%/64.6%and 0%/0%, respectively.The approach from Li and Fox [19] thus heavily depended on statistical data, which are not updated very frequently.Dong et al. [20] mapped three forest types, including rubber plantations, on the Hainan Island using Advanced Land Observing Interferometric Synthetic Aperture Radar (ALOS PALSAR), MODIS NDVI, MODIS EVI, and Land Surface Water Index (LSWI) time series.Their results suggested a good separability of rubber plantations from other forest types using a simple threshold of summer and winter NDVI composites (85% overall accuracy for their binary rubber plantation classification).The method from Dong et al. [20] is an effective application of MODIS time series, but it is likely to work well only in evergreen forests, as they do not shed leaves seasonally in contrast to rubber plantations.To map rubber for the large tropical seasonal forest regions of Southeast Asia based on MODIS time series will require an approach that builds upon an in-depth understanding of the phenological patterns of natural forests and rubber plantations.
Short-wave infrared (SWIR) reflectance and SWIR-based indices have been shown to be important predictors for mapping rice paddies and certain forest types in Asia and Southeast Asia [21][22][23] and some studies have used SWIR to derive phenological metrics [24,25].The importance of SWIR for discriminating forest composition and structure has been known of a while [26].Yet, vegetation indices like the NDVI and EVI are commonly used to characterize vegetation phenology.Recently, Dong et al. [27] used Landsat data acquired during seasonal leave senescence and found that the SWIR-based LSWI in addition to other VIs was important for discriminating rubber plantations and forests.Evaluating the outcome of these studies, it is likely that the mapping of rubber plantations and natural forests can be enhanced by capturing the phenological dynamics of both land cover classes.Furthermore, incorporating the phenological dynamics of the SWIR reflectance might further enhance the mapping of rubber plantations.
In this study, we tested a new approach for mapping rubber plantations and natural forests using phenological metrics derived from MODIS EVI and SWIR reflectance time series.By analysing the importance of each phenological metric on classification accuracy we then explore the seasonal differences of rubber plantations and natural forests based on MODIS time series.Our study is performed in the region of Xishuangbanna, China, where rubber plantations have become one of the major land cover types over the past decades.

Study Area
The autonomous prefecture Xishuangbanna in the Yunnan Province is the most southern prefecture of China and borders Laos to the east and Myanmar to the west (Figure 1(A)).The prefecture is subdivided into three municipalities: Jinghong, Menghai, and Mengla.The capital of Xishuangbanna, Jinghong, is located in the central low elevation areas (~400 to 1,000 m above sea level) of the Mekong River.The northern parts of the Jinghong municipality are at a higher elevation, reaching 2,000 m above sea level.The western municipality Menghai is characterized by higher elevations reaching up to 2,500 m above sea level, whereas the eastern municipality, Mengla, is dominated by lowlands to the west and highlands to the east.Due to the bordering highlands, Xishuangbanna is the only region in continental China with tropical, monsoon-influenced climate.The wet season starts in April and ends in November; and the precipitation maximum is reached in July/August with an average of 300 mm/month.In the dry season (December to March) heavy fog is occurring frequently.Temperatures do not drop below 15 °C except for the very high elevation areas.Wet-season temperatures are high with peaks around May and August/September.With its junction location between different climate and ecological zones, Xishuangbanna inhabits a diverse flora and fauna, which makes up for 20% of the total species diversity of China [28].However, Xishuangbanna is also home to a massive rubber producing agro-industry [6] established in the late 1990s [29].The rubber trees in Xishuangbanna are deciduous trees that shed their leaves for a relatively short period of two to four weeks during the coldest and driest month (January to March) [30].Throughout the rest of the year, rubber trees stay foliated.Forests in Xishuangbanna can be differentiated into four major forest types: (1) tropical rain forest, (2) tropical seasonal moist forest, (3) tropical montane evergreen broad-leaved forest, and (4) tropical monsoon forest [31].Forest types 1 to 3 are evergreen forests, whereas forest type 4 is a deciduous forest, influenced by annual dryness caused by the monsoon climate.Evergreen forests rely on heavy water deposition from fog in the dry season to overcome water shortages [32].Tropical rain forests and tropical seasonal rain forests present differences in species composition, caused by different soil types.Montane evergreen broad-leaved forests are only found in elevations higher 1,000 m.Monsoon forests mostly occur on the low elevation banks of the Mekong River or in low elevation basins.The Natural Forest Conservation Program (NFCP), which has been established in 1990 [33], strictly protects remaining natural forests in Xishuangbanna.However, before the NFCP was established, forests have largely been converted to shifting cultivation or other cash crops.Many forest areas are therefore segmented by abandoned or fallow lands, which are covered by secondary vegetation such as deciduous monsoon forests, savannah woodlands, bamboo, and grasslands [34].These small patches of deciduous vegetation within the natural evergreen forest can change the phenological response of tropical rain forests to a more seasonal variation [31].Rubber plantations, in turn, can be mixed with other cash crops such as pineapple, which are cultivated year-round.In such cases, rubber plantations might have green vegetation during the dry season.

Data and Methods
The objective of our work was to evaluate the utility of phenological metrics derived from MODIS time series to map rubber plantations and natural forests in Xishuangbanna.We therefore derived phenological metrics from MODIS EVI and SWIR time series (see Section 3.3).To evaluate the added value of including SWIR metrics into a classification, we trained and validated three Random Forest classification models (see Section 3.4) based on: (1) phenological metrics from the EVI time series, (2) phenological metrics from the SWIR time series, and (3) phenological metrics from both, EVI and SWIR time series (Figure 2).As a benchmark, we then compared the classification accuracy of the best model based on phenological metrics to the classification accuracy of a model based on the full, fitted time series data (Figure 2).Finally, we calculated the variable importance (see Section 3.4) for the overall best model based on phenological metrics.By doing so, we explore the relative importance of each phenological metric for classifying rubber plantations and forests.

Moderate Resolution Imaging Spectroradiometer (MODIS) Data Preparation
We used the MODIS Terra Vegetation Index product (MOD13Q1, H/V 27/6) from Collection 5 at 250 m spatial resolution.MODIS data were retrieved from the National Aeronautics and Space Agency (NASA) Earth Resource Observation and Science Center via the Global Visualization Viewer (GLOVIS).The time series spans the period from January 2010 to December 2010 with a 16-day interval, resulting in 23 time steps.We extracted the SWIR reflectance and the EVI time series from the MOD13Q1 Vegetation Index product.We used absolute SWIR reflectance instead of a normalized SWIR based index, because we wanted to test direct effects of the SWIR on the mapping accuracy of forests and rubber plantations.The EVI is a vegetation index using the red ( ), NIR ( ) and blue reflectance ( ) (Equation ( 1)): ( We chose the EVI instead of NDVI, because it is less sensitive to atmospheric noise and does not saturate as early as NDVI at high biomass [35].

Reference Data
The collection of reference data in the data scarce and remote areas of Southwest China is a challenging task.In this paper, we used a combination of field data, high resolution Quickbird imagery from Google Earth TM (GE), and RapidEye imagery acquired through the RapidEye Science Archive (RESA) as reference data.With the GE true-colour imagery at very high spatial resolution (2.6 m × 2.6 m) it was possible to visually identify the systematic row structure of rubber plantations.The RapidEye imagery has a lower spatial resolution (interpolated to 5 m × 5 m), but the sensor features a red-edge band and a near-infrared band, which facilitated differentiating between rubber plantations and forests.RapidEye imagery was available for the entire Xishuangbanna prefecture, whereas Quickbird imagery in GE was only partly available (Figure 1(C)).Field mapping was performed during April 2012 in the township of Manlin and South of the city of Jinghong in Xishuangbanna (Figure 1(C)).We used the field data to help with the interpretation of GE and RapidEye imagery.
To construct a representative reference data set, we randomly selected 520 MODIS pixels.The sample size was chosen based on a multinomial distribution as suggested by [36,37], using a confidence level of 95% (precision of 5%) and an estimated probability of 35% of the largest class.This number was selected based on a priori knowledge from a preliminary analysis.As response design, each pixel was labelled according to its primary land cover using the GE Quickbird and RapidEye imagery (Table 1).The primary land cover was defined as the land cover that covered the largest share of the MODIS pixel area.Hence, we included pure and mixed pixels as reference data.

Time Series Analysis
We used the TIMESAT software [38] to extract phenological metrics from the EVI and SWIR time series.First, TIMESAT fits a double logistic function to each pixel's temporal signal using a nonlinear least-square fit [39].The result is a fitted time series for each pixel that is temporally smoothed, i.e., data only shows residual atmospheric noise and is gap-filled (e.g., for clouds).We setup TIMESAT to apply three fitting iterations, an adaption-strength parameter of 2, a seasonality parameter of 1, and a window size of 1 for the median filter.Based on the fitted time series we then extracted nine phenological metrics (exemplified in Figure 3).The metrics can be grouped into temporal metrics that describe the timing of seasonal events (x-direction) and spectral metrics (y-direction) that describe the EVI and SWIR values at seasonal events.The maximum is the maximum value for a given season (C).The base value is the mean of the off-season minima (A).The season amplitude is the difference between the base value and the season maximum.Season start (B) is defined as the day of year (DOY) at which the time series value increased 50% in relation to the season amplitude, and season end is defined as the DOY at which the time series value decreased 50%, respectively (D).The length is the time difference between season start and season end.The left and right derivatives are the first order derivative in point B and C.And finally, the middle is the DOY at which the maximum occurs.For some pixels (n = 938; 0.3%) no season was extractable, because of extreme noise or cloud cover.These pixels were excluded from further processing steps and labelled as "no season found" in the resulting classification maps.

Classification Algorithm, Variable Importance, and Accuracy Assessment
In this study, a Random Forest (RF) classification approach was used [40].RF is ensembles of decision-trees wherein each tree is trained on randomly selected features of a bootstrapped sample of training data.RF are increasingly used in remote sensing and proved to outperform other machine learning algorithms in terms of accuracy and computational resources [41][42][43].Since each tree in the forest is only trained on a bootstrapped sample, the complement of the sample (out-of-bag) is used to estimate the accuracy of each tree, which is then aggregated over all trees.The aggregated out-of-bag accuracy represents an unbiased estimate of map accuracy as long as the reference data were obtained via probability sampling, as given in our case [41,43].We used the out-of-bag predictions to construct a confusion matrix from which we calculated overall, user's and producer's accuracies [44].By using the out-of-bag accuracy in combination with the random feature selection, the importance of each variable for the classification result was estimated.This was done by measuring the mean decrease in accuracy if a feature is left out in building a tree.The RF variable importance measure has proven to be a reliable tool, providing insights into the predictive power of individual variables [45].
The RFs in this study were built using the statistical software R [46] and the randomForest package [47].The number of trees was set to 1000, the sample size of the bootstrapped sample was stratified with 70 pixels per class to avoid effects of imbalanced training data [48], and the number of features used for each tree was set to the square root of the total number of features, being the standard setting in the original RandomForest™ software [40].

Classification Results
The best model resulted in an overall accuracy of 73.5% and included phenological metrics from EVI and SWIR (Table 2).Using SWIR metrics alone achieved an overall accuracy of 67.9%, slightly better than the EVI model with 64.4%.The SWIR model discriminated better between rubber plantations and forests; producer's accuracies were 5.8% and 13.7% higher compared to the EVI model.In comparison, the EVI model better differentiated between non-forest and forests/rubber plantations.Here, user's and producer's accuracies of the non-forest class were 23.9% and 9.1% higher compared to the SWIR model.
The overall best model, trained on metrics from EVI and SWIR time series, resulted in classification accuracies of 63.6% for rubber plantations, 80.2% for forests and 71.6% for non-forest areas.Confusion mostly occurred between rubber plantations and forests, but also between rubber plantations and non-forest areas (Table 3).Forest/non-forest differentiation yielded better results.The classification map resulting from the EVI model (Figure 4(A)) showed many forest pixels interspersed with rubber plantations (Figure 4(A.1)).It reflects the low producer's accuracy of rubber plantations for the EVI model (Table 2).The map resulting from the SWIR model (Figure 4(B)) correctly represented large rubber plantation areas in the central and southern low-lands.However, many mixed pixels at the border of non-forest and forest areas in the northern and north-western parts of Xishuangbanna were falsely classified as rubber plantations (Figure 4(B.2)), which explains the high confusion between non-forest and forest/rubber plantation areas (Table 2).The EVI-SWIR model presented an accurate classification map (Figure 4(C)) that best captured rubber expansion areas known from previous studies on rubber expansion in Xishuangbanna [11] (Figure 4

(C.3)).
The classification accuracy of the SWIR-EVI model based on the full time series data was 74.8%.Thus, the metric-based model was slightly less accurate (1.3% difference).The differences were largest for rubber plantations, which had a producer's accuracy of 63.6% with the metrics-based model and 68.8% with the model based on the full time series (5.2% difference).

Variable Importance
We estimated the relative importance for each phenological metric based on the combined EVI-SWIR model.For better visibility we analysed importance of EVI and SWIR metrics in the combined EVI-SWIR model separately (Figure 5).The by far most important phenological metric of the EVI was the base value followed by the length, left derivative, and amplitude.Within the SWIR, the base value was also the most important metric followed by the start, the middle, and the maximum of the season.Hence, the spectral values of the EVI metrics were more important than their timing, whereas with the SWIR, temporal (start and middle) and spectral (base value and maximum) metrics were important predictor variables in the combined EVI-SWIR model.

Discussion
Using phenological metrics for mapping natural forests and rubber plantations resulted in an overall accuracy of 73.5%.The finding that phenological metrics are reliable predictors of land cover is in agreement with other studies that used phenological metrics to map other land cover types [16].We also identified that overall accuracies only slightly increased when using the full time series instead of the phenological metrics (1.3% increase).However, differences were largest for rubber plantations (5.2% increase in producer's accuracy), which was our main class of interest.This suggests that the phenological metrics were overall sufficient to capture the spectral-temporal differences between the different land cover classes, but using the full time series generally achieved higher accuracies.Nonetheless, using phenological metrics has several advantages: The feature space is substantially reduced (i.e., 46 vs. 18 features for the EVI-SWIR model), which reduces redundancies and multi-collinearity in the explanatory variable set [49].Second, phenological differences in land cover can be broken down into single key metrics describing a specific land cover.These more simple metrics facilitate the interpretation of phenological differences and their importance for land cover classifications.These more generalized metrics might also increase the transferability of the results to other study areas [49].Third, metrics that have been identified as relevant for differentiating between natural forests and rubber plantations can be monitored over time more easily than tracking changes in multi-dimensional annual time series.These three arguments support the use of phenological metrics as predictor for land cover.
The study also demonstrated the importance of SWIR for mapping rubber plantations.While phenological metrics based on EVI differentiated better between non-forest and forests/rubber plantations, phenological metrics based on SWIR were important to differentiate between rubber plantations and forests.Based on the EVI, the temporal-spectral profile of rubber plantations and forests were very similar (Figure 6(A,B)).In comparison, the SWIR profile showed distinct differences between rubber plantations and forests (Figure 6(D,E)).In the SWIR, the peak of the season occurred on average much earlier over rubber plantations (day of year, DOY: 128 ± 80) than over forests (DOY: 176 ± 48).Similarly, the season start was much earlier in rubber plantations (DOY: 48 ± 80), compared to forests (DOY: 96 ± 48).The differences in the timing of the season start and peak in the SWIR were also reflected by the high importance of these metrics as predictor variables in the RF classification (Section 4.2).
The temporal differences in the SWIR time series between rubber plantations and forests may be explained by several reasons.In the dry season, green vegetation cover of rubber plantations is low and the spectral signal is dominated by open soils, which leads to an increase in SWIR reflectance.The decrease in soil water content in the dry season may further increase SWIR reflectance [50].The SWIR maximum of rubber plantations therefore represents the coldest and driest time of the year.Though, the actual timing may also vary with site conditions, understory vegetation, and rubber tree density, visible in the high standard deviation of the SWIR maximum in rubber plantations (Figure 6(D)).In comparison, the SWIR signal of forests in the dry season is less influenced by soil reflectance, probably because seasonal forests include mixtures of evergreen and deciduous trees, and a greater proportion of understory vegetation.These winter/dry season differences between seasonal forests and rubber plantations are also supported by the significantly higher EVI base value of forests compared to rubber plantations (t = −2.56,df = 346, p-value < 0.05).This finding is in agreement with Dong et al. (2012) [20] who used the difference in NDVI during the winter months to differentiate between evergreen forests and rubber plantations, and with Dong et al. (2013) [27] who highlighted the importance of NDVI, EVI, and LSWI reflectance in the defoliation period for differentiating between rubber plantations and forests.However, in our study, the differences in EVI base value alone were not sufficient to accurately separate forests and rubber plantations (lower accuracies with EVI model).This may be explained by the varying phenology of forests in Xishuangbanna, which partly resembles the phenology of rubber trees (Section 2).
Interestingly, the SWIR peaked similarly to the EVI during the summer months over forests but not over rubber plantations.This summer increase in SWIR for forests may be explained by structural changes in the forest canopy, i.e., less shadows and higher reflectance caused by higher leave area.In rubber plantations we would expect the same effect, because rubber plantations also increase their leave area in summer (higher EVI).However, the winter SWIR reflectance peak in rubber plantations caused by the soil signal (discussed above) may simply superimpose the summer peak, because of the higher SWIR reflectance of soils.This explanation is supported by the fact that SWIR reflectance of rubber plantations and forests in the summer months are equally high (around 10%, Figure 6(D,E)).Therefore, EVI and SWIR time series of forests have coincident seasonal peaks (Figure 6(B,E)), whereas rubber plantations have their SWIR peak in the winter/dry season (Figure 6(D)) and the EVI peak in the summer (Figure 6(A)).
The importance of the EVI base value can also be attributed to the separability of forest areas (including rubber plantations) and non-forest areas (i.e., cropping, urban, and water areas).All non-forest pixels presented a significant lower base value in the EVI time series (mean of 0.27 ± 0.08) compared to rubber plantations (mean of 0.39 ± 0.07; t = 11.90, df = 201, p-value < 0.01) and forests (mean of 0.41 ± 0.07; t = 14.75, df = 168, p-value < 0.01).Furthermore, the importance of the SWIR base value can be attributed to differences between forest and non-forest areas.The SWIR base value of non-forest areas was significantly higher (mean of 0.09 ± 0.02) compared to rubber plantations (mean of 0.07 ± 0.02; t = −4.69,df = 183, p-value < 0.01) and forests (mean of 0.05 ± 0.01; t = −13.84,df = 136, p-value < 0.01).This may be explained by the majority of crop pixels within the non-forest class, where the soil in the dry season causes a high SWIR reflectance (Figure 6(F)).In comparison to the study from Li et al. [19], we achieved similar producer's accuracy for rubber plantations if compared to their classification based on MODIS and statistical data (63.6% compared to 67.8%/59.4%),but lower user's accuracies (64.9% compared to 98.1%/97.2%).However, we achieved higher accuracies if compared to their product solely relying on MODIS data (producer's accuracy: 63.6% compared to 60.9%/0%; user's accuracy 64.9% compared to 64.6%/0%).Since their study mainly relies on the temporal profile of the EVI we can highlight the importance of the SWIR temporal profile for mapping rubber plantations.Nevertheless, their study covered a larger region and accuracy measures that were derived from different validation protocols should always be compared with caution.In comparison to Dong et al. [20], we achieved a lower overall accuracy compared to their binary rubber plantation map (73.5% compared to 85%).However, they only mapped rubber plantations/non rubber plantation areas and did use a different validation protocol (i.e., no class accuracies were reported).It is therefore difficult to compare their classification result with our results.Nevertheless, their assumption of using winter vegetation index differences for differentiation between (evergreen) forests and rubber plantations is supported by our study.One side aspect, which must be mentioned at this point, is that we relied on MODIS data that has been recorded at 500 m spatial resolution (the blue band included in the EVI and the SWIR band) and subsequently resampled to 250 m resolution.Even rubber plantations in Xishuangbanna are very homogeneous; we might have an error introduced by mixed pixels.Since we selected mixed pixels in the reference data, our reported accuracies account for this error.However, if one wants to transfer the method to a more heterogeneous area, this might be an issue to be considered.
Considering that we solely relied on MODIS data and that the forest phenology in Xishuangbanna closely resembles the seasonal patterns of rubber plantations, we could present a novel method for mapping rubber plantation that potentially can be extended to a monitoring framework for rubber expansion in Southeast Asia.Especially the temporally monitoring of expansion "hot spots" may be a major task.These "hot spots" can then be analysed in detail using finer resolution, but temporally infrequent data sources such as Landsat or RapidEye.Hence, a combined monitoring scheme could help understanding recent trends in rubber expansion and the associated ecological and economic consequences.
In 2011, rubber plantations made up approximately 30% ± 4% of the total land cover in Xishuangbanna, whereas natural forests covered 49% ± 3% of Xishuangbanna.Our results present the newest estimate of rubber expansion in Xishuangbanna.Previously, Li et al. [11] estimated that rubber plantations made up 11% and natural forests 50% of the land cover in 2003.Comparing our results and those of Li et al. [11] (which are based on Landsat data) suggests nearly a tripling of the rubber plantations since 2003.At the same time, total forest area was stable during that period, most likely because of the NFCP (see Section 2).Our results also suggest that rubber plantations were established mainly on former fallows and arable areas.

Conclusion
We presented a new approach for mapping rubber plantations and forests in Xishuangbanna that combined phenological metrics derived from Moderate Resolution Imaging Spectroradiometer (MODIS) Enhanced Vegetation Index (EVI) and short wave infrared (SWIR) reflectance time series.Based on our results we can draw three general conclusions: (1) Phenological metrics derived from MODIS time series were useful predictors for mapping rubber plantations and natural forests (overall accuracy of 73.5%), but did not achieve results as good as using the full time series data (1.3% difference in overall accuracy).(2) Phenological metrics derived from SWIR substantially improved the classification accuracy of rubber plantations compared to classifications based on EVI time series alone (9% improvement in overall accuracy).(3) The key phenological metrics discriminating rubber plantations and natural forests were: timing of seasonal events in the SWIR time series and EVI/SWIR reflectance during the dry season.Our study achieved classification accuracies comparable to other studies that developed methods for mapping rubber plantations, but we did not rely on data from other sensors or statistical data on rubber growth and production.Hence, our approach could potentially be extended to map rubber expansion annually over tropical regions.Further research on the influence of site conditions and climate variability on rubber and forest phenology could help improve the accuracy and generalizability of the phenology-based classifications.Finally, results highlighted the recent expansion of rubber plantations in Xishuangbanna, transforming areas formerly used for subsidence or cash cropping into rubber plantations.These rapid changes underpin the need for developing a regular monitoring scheme for rubber expansion in Southeast Asia.

Figure 1 .
Figure 1.(A) Topography of the study area.(B) Location of the study area.(C) Available Quickbird imagery.

Figure 2 .
Figure 2. Processing scheme showing the comparison of the models based on phenological metrics to each other and to the full time series model.RF: random forest classification; Acc.: accuracy assessment; TS: time series; Var.Imp.: variable importance.

Figure 3 .
Figure 3. Raw and fitted Moderate Resolution Imaging Spectroradiometer (MODIS) Enhanced Vegetation Index (EVI) (A) and short-wave infrared (SWIR) (B) time series of a natural forest pixel with phenological metrics.A: Base values of season; B: Start of season; C: Maximum of season; D: End of season.

Figure 5 .
Figure 5. Variable importance measure for the overall best model based on Enhanced Vegetation Index (EVI) and short-wave infrared (SWIR) metrics, estimated by the mean decrease in accuracy if a feature is left out in building a decision tree (Section 3.4).

Figure 6 .
Figure 6.Smoothed Enhanced Vegetation Index (EVI) (A-C) and short-wave infrared (SWIR) (D-F) time series extracted and averaged for all reference pixels of the rubber plantation, forest, and non-forest classes.Dotted lines indicate one standard deviation.

Table 1 .
Reference classes and number of reference pixels.

Table 2 .
Overall producer's and user's accuracies for the three models based on phenological metrics.

Table 3 .
Confusion matrix for the overall best model based on Enhanced Vegetation Index (EVI) and short-wave infrared (SWIR) metrics.