Phenology-Based Vegetation Index Differencing for Mapping of Rubber Plantations Using Landsat OLI Data

Accurate and up-to-date mapping and monitoring of rubber plantations is challenging. In this study, we presented a simple method for rapidly and accurately mapping rubber plantations in the Xishuangbanna region of southwest China using phenology-based vegetation index differencing. Temporal profiles of the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Atmospherically Resistant Vegetation Index (ARVI), Normalized Difference Moisture Index (NDMI), and Tasselled Cap Greenness (TCG) for rubber trees, natural forests and croplands were constructed using 11 Landsat 8 OLI images acquired within one year. These vegetation index time series accurately demonstrated the unique seasonal phenological dynamics of rubber trees. Two distinct phenological phases (i.e., defoliation and foliation) of rubber trees were clearly distinguishable from natural forests and croplands. Rubber trees undergo a brief defoliation-foliation process between late December and mid-March. Therefore, vegetation index differencing between the nearly complete defoliation (leaf-off) and full foliation (leaf flushing) phases was used to delineate rubber plantations within fragmented tropical mountainous landscapes. The method presented herein greatly improved rubber plantation classification accuracy. Overall classification accuracies derived from the differences of the five vegetation indices varied from 92% to 96% with corresponding kappa coefficients of 0.84–0.92. These results demonstrate the promising potential of phenology-based vegetation index differencing for mapping and monitoring rubber expansion at the regional scale. OPEN ACCESS Remote Sens. 2015, 7 6042


Introduction
Over the past several decades, massive rubber plantations have become major drivers of deforestation and other types of habitat degradation in Southeast Asia and southern Yunnan of China [1][2][3].The environmental and socioeconomic impacts of rubber expansion have widely been explored at various scales [4][5][6][7][8][9].Understanding and quantifying the large-scale impacts of rubber plantations requires accurate mapping of rubber plantations and timely measurements of the conversion of other land-cover types to rubber plantations [3].
Remote sensing has played an important role in mapping and understanding changes in the areal extent and spatial pattern of rubber plantations.Recently, numerous studies have used satellite imagery for delineating rubber plantations [3,[10][11][12][13][14][15][16].These studies can be categorized into two broad groups: monotemporal image classification [11,[13][14][15][16] and phenology-based classification [3,10,12].Phenology-based classification is effective because it includes the phenological features derived from the time series of satellite images, such as MODIS and China's Feng-Yun-3A [10,12,17].However, the relatively coarse spatial resolutions of these time-series data and frequent cloud cover during the growing seasons of rubber trees make it difficult to delineate and map rubber plantations in fragmented landscapes in tropical regions.To overcome this drawback, Dong et al. [11] used multiyear multitemporal Landsat images for studying the rubber plantation's phenological trajectory, and then combined the PALSAR-based forest mask and single-date Landsat TM image in the foliation stage of rubber trees to map rubber plantations with high classification accuracy.
Multi-seasonal imagery that captures different periods of the growing season is of considerable value for characterizing land cover types [18][19][20][21][22][23].In particular, using paired leaf-on and leaf-off images can result in substantial improvements in the accuracy of classifying forest types [24].Compared with rubber plantations, natural forests have higher vegetation indices during the defoliation stage (dry season) of rubber trees [10][11][12]25].Therefore, it is possible to effectively distinguish rubber plantations from natural forests and other landscapes by using paired defoliation and foliation images.
Vegetation index differencing is often employed to detect changes using remote sensing [26].It is also shown that the inclusion of phenological differences in vegetation indices can greatly improve the accuracy of land cover classifications over a variety of landscapes [22].The objective of this study was to develop and examine a simple and novel approach for delineating and mapping rubber plantations at the regional scale by combining vegetation index differencing and phenological information between the defoliation and foliation stages of rubber trees.This study was conducted in the region of Xishuangbanna in southern Yunnan, China, where rubber trees are presently the most important and prevalent cash crop.The phenology of three major land cover types (i.e., rubber trees, natural forests and croplands) in Xishuangbanna were examined based on the temporal profiles of several vegetation indices derived from Landsat-8 OLI images acquired during the dry season of 2013-2014.Next, the performances of the phenology-based vegetation index differencing method for delineating and mapping rubber plantations were evaluated based on a confusion matrix estimated from random reference data samples.

Study Area
The study area of Xishuangbanna contains the largest rubber plantation areas and accounts for the largest proportion of natural rubber production in China.According to data from the Xishuangbanna Statistical Yearbook in 2013 [27], the total area of rubber trees in the prefecture were 2940.28 km 2 with a total natural rubber production of 31.74 × 10 4 tons, accounting for ~27% and ~37% of the total rubber tree cropping areas and production in China, respectively.Rubber plantations currently cover more than 15% of Xishuangbanna's landscape, and rubber plantations are continuing to expand.Forests and croplands are the other two major land cover types in this area accounting for 75.6% and 5.0% of the total area, respectively.Xishuangbanna is a Dai Nationality Autonomous Prefecture that is located in the most southern region of Yunnan province in Southwest China and borders Laos and Myanmar (Figure 1).It covers an area of 19,125 km 2 and has a mountainous topography.The elevation of this prefecture varies from 2430 m at the top of the mountain in the north to 480 m at the bottom of the lowest valley in the south.Xishuangbanna has a tropical monsoon climate with a seasonal rhythm that is roughly made up of five stages: the fog-cold stage (early December to mid-February), dry-cool stage (late February to late March), dry-hot stage (early April to early May), humid-hot stage (mid-May to mid-October) and humid-cool stage (late October to late November) [28].The annual mean temperature varies from 21.7 °C at lower elevations to 15.1 °C at higher elevations.The annual precipitation ranges from 1100 mm at the river valley to 2500 mm at the top of the mountain, with more than 80% occurring during the humid-hot stage.In the areas of the lower hills and valleys that are dominated by tropical forests, the annual mean temperature is approximately 21 °C and the annual precipitation is greater than 1500 mm [29].Due to its unique location that links different climate and ecological zones, this region harbors diverse flora and fauna that account for ~20% of the total species diversity in China.The primary vegetation of Xishuangbanna can be grouped into four major forest types: tropical rain forest, tropical seasonal moist forest, tropical montane evergreen broad-leaved forest and tropical monsoon forest [29].Natural forests are mainly evergreen, while rubber trees show deciduous characteristics during the fog-cold and dry-cool stages [28].Although some tree species in tropical monsoon forests are also deciduous, they make up a small proportion (<15%).Additionally, tropical monsoon forests are mosaically distributed in seasonal rain forests, and have a greater proportion of understory vegetation.This means they present a phenological rhythm different from rubber trees [28].

Landsat-8 OLI Data and Pre-Processing
Eleven orthorectified and terrain corrected Level 1T Landsat-8 OLI images (path/row 130/45) that were acquired from 5 November 2013 to 16 May 2014 were obtained from the USGS EarthExplorer [30].The specific dates of these Landsat scenes can be found in Figure 2. All images were routinely produced using the Level 1 Product Generation System (LPGS), version 2.2.3.The Landsat LDOPE Toolbelt provided by the MODIS land quality assessment group [31] was employed to detect clouds from the Landsat-8 OLI Quality Assessment (QA) band and to produce a cloud mask for removing cloud-contaminated pixels.Subsequently, Fmask software developed by Zhu and Woodcock [32] was used for automated cloud shadows masking.Cloud and cloud shadow percentages for all the images range from 0.05% to 36.84%.The conversion of scaled Digital Numbers (DN) to the Top Of Atmosphere (TOA) radiance and/or reflectance was subsequently conducted by following the procedure described on the USGS website [33].Atmospheric effects were removed by using the ENVI FLAASH (the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) module based on a MODTRAN4-based atmospheric correction algorithm [34,35].The parameters for removing atmospheric effects with FLAASH were selected following the detailed description in [35].Before differencing the vegetation index, the images were geometrically registered within a mean square error of less than 0.5 pixels.Topographical correction was not considered because no universal topographic correction method applies to all situations [36].Also, the band-ratio vegetation indices could significantly reduce the noise caused by topographical variations [37].Four spectral vegetation indices, the Normalized Difference Vegetation Index (NDVI) [38], Enhanced Vegetation Index (EVI) [39,40], Atmospherically Resistant Vegetation Index (ARVI) [41,42], and Normalized Difference Moisture Index (NDMI) [43,44], were calculated using the surface reflectance.These indices were formulated by using the following equations: where ρ Blue , ρ Red , ρ NIR , and Band 4 (red, 0.64-0.67μm), Band 5 (near-infrared, 0.85-0.88μm) and Band 6 (shortwave-infrared, 1.57-1.65 μm) in the Landsat-8 LOI images, respectively [45].
The Tasseled Cap transformation incorporates six different spectral bands into vegetation indices that can be directly associated with physical characteristics [46].Thus, the Tasseled Cap Greenness (TCG) was also derived from the Landsat-8 OLI at-satellite reflectance following the equation described by Baig et al. [47].

Reference Data
We combined field data, high-resolution imagery from Google Earth™, and a forest inventory vector dataset to produce reference data.Field surveys were conducted between September and November of 2014 using a GPS digital camera with a horizontal position accuracy of less than 10 m as the main tool for capturing information.More than 6000 geo-referenced field photos and records were acquired in plots where the examined land cover areas were larger than 90 m × 90 m.The GPS tracks of the ground survey are shown in Figure 1.We processed all of the geo-referenced photos and records as kml files and geo-linked them with Google Earth™.Subsequently, the polygons of interest (POIs) for the examined land cover types were digitalized by interpreting the high spatial resolution imagery within Google Earth™ and were used as reference data for evaluating the accuracy of mapping the rubber plantations.Numerous studies have suggested that this approach is effective [3,[10][11][12].In this study, furthermore, these reference data were cross-validated with forest inventory vector datasets developed in 2006 by the local government to create a final reference data set (Figure 1), given that forest inventory data recorded stand ages of rubber trees and natural forests.

Image Difference Threshold Selection and Accuracy Assessment
Thresholding is a critical step in detecting change in image differencing [48,49].Different methods were developed for thresholding the difference image for binary change detection.In most cases, a series of discrete threshold values are used to construct 'accuracy assessment curves' or 'calibration curves' that illustrate the variations in accuracy with the threshold value [50][51][52].Then, the optimum threshold value (OTV) that results in the highest accuracy is applied to create a binary change mask.In this study, we used the median and median absolute deviation of the difference images as thresholds for constructing an accuracy assessment curve.A series of discrete threshold values with a narrow interval (e.g., 0.01 median absolute deviations) were tested to generate a relatively continuous curve.The accuracy measures, such as the overall accuracy (OA), producer's accuracy, user's accuracy and kappa coefficient, were calculated for each threshold value by using an error matrix.The OTV with the highest kappa coefficient was applied to build a binary mask for delineating rubber plantations from other land cover types.Subsequently, an Erode filter (3 × 3 kernel sizes) together with a Dilate filter (3 × 3 kernel sizes) was applied to the binary mask for producing the resulting rubber plantation map.For this study, the differences in the vegetation indices were obtained by subtracting the values of the vegetation indices on 9 February 2014 from the values of the vegetation indices on 13 March 2014.Before this, some cloud-contaminated pixels occurred in the north part of the study area on 13 March 2014 and were filled with the cloud-free pixels from 29 March 2014.Accuracy assessment of the rubber plantation maps was performed using the confusion matrix.The random reference data samples consisted of 10,000 pixels was individually extracted from 781 POIs for rubber plantations and 1917 POIs for the non-rubber types of land cover.

Unique Phenological Behaviour of Rubber Trees Observed from Multi-Temporal Landsat 8 OLI Images
As shown in Figure 2, the multi-temporal profiles of five spectral vegetation indices (NDVI, EVI, ARVI, NDMI and TCG) extracted from random reference data samples indicated drastic differences among three examined land cover types.The vegetation indices from natural forests only slightly fluctuated during the examined period, while these indices from croplands consistently decreased from early November 2013 to the middle of April 2014.However, the vegetation indices from rubber plantations substantially decreased during the period from late December 2013 to early February 2014 (defoliation phase) and subsequently experienced a rapid and pronounced rebound until the middle of March 2014 (foliation phase).The lowest vegetation index values of the rubber plantations occurred in early February.The image acquired during the same period indicated that the rubber trees were nearly completely defoliated (Figure 3a).
Compared with the other two examined land cover types, rubber plantations had lower vegetation indices in early February and higher vegetation indices in the middle of March (Figure 2).The largest differences in the vegetation indices between rubber plantations and natural forests were observed in early February when rubber plantations had far lower vegetation indices than natural forests.In mid-March, vegetation indices from rubber plantations rapidly recovered and exhibited higher values than those of natural forests.The images acquired during the same period indicated that rubber plantations translated into light green patches (Figure 3b).With the rebound of the vegetation indices from rubber plantations, the greatest differences in the vegetation indices between rubber plantations and croplands were found in the middle of March.However, the natural forest is always green during the defoliation (a) and foliation (b) stages of rubber trees.Major land cover classes of interest were marked in the zoom-in images, including rubber plantation ("A"), natural forest ("B"), and cropland ("C").

Map of Rubber Plantation Derived from Image Differencing
As shown in Figure 4, the binary classifications of rubber plantations versus non-rubber land cover were obtained by thresholding the difference images of the five spectral vegetation indices between at the nearly complete defoliation and full foliation stages of rubber trees.The resulting rubber plantation maps are highly accurate according to the confusion matrix estimated from random reference data samples (Table 1).The producer's accuracies in mapping rubber plantations and non-rubber land cover were 87%-94% and 96%-99%, respectively, with corresponding user accuracies of 96%-99% and 88%-94%.Consequently, the overall accuracies of the resultant maps varied from 92% to 96% with kappa coefficients from 0.84 to 0.92 (Figure 4).Differencing the ARVI and NDMI images performed better and produced the highest overall accuracy of 96% with a kappa coefficient of 0.92.In comparison, differencing the ARVI image achieved higher user's accuracy of rubber plantations than differencing the NDMI image.Among all the differencing images, the TCG differencing image produced the lowest overall accuracy of 92% with a kappa coefficient of 0.84.The rubber plantation maps represented a highly consistent spatial distribution (Figure 4b).These maps correctly depicted large rubber plantation areas in the central and southern lowlands.However, the rubber plantation maps extracted from the difference images of the NDVI, EVI, NDMI and TCG falsely classified some pixels as non-rubber areas in the western and northwestern regions (i.e., central Menghai County) of the rubber plantations.In contrast, differencing the ARVI images presented fewer misclassified pixels in the western and northwestern portions of Xishuangbanna and achieved the highest user's accuracy, 99%.
Figure 4 showed that the total area of rubber plantations in Xishuangbanna ranged from 3107.31 km 2 to 3841.26 km 2 in 2014.The NDMI differencing image yielded the largest rubber plantation area, while the TCG differencing image yielded the smallest area.The rubber plantation area achieved from the EVI differencing image (3193.00km 2 ) was roughly equivalent to that from the TCG differencing image.The NDVI and ARVI differencing images produced moderate rubber plantation areas of 3514.75 km 2 and 3707.28 km 2 , respectively.

Phenology of Rubber Trees and Its Potential for Mapping Rubber Plantations
The temporal trajectories of the five vegetation indices captured by multitemporal Landsat OLI data in the three major land cover regions showed that the phenology of rubber plantations is distinctly different from that of the other two major land covers.Two specific phenological stages (i.e., defoliation and foliation) can be clearly identified in rubber plantations.The pronounced decrease in the vegetation indices from late December to early February of the next year indicated that substantial defoliation of the rubber trees occurred.The drastic increase in the vegetation indices from early February to mid-March suggested that the rubber trees underwent speedy leaf flushing and canopy recovery.Although the different phenological behaviour of rubber trees from natural forests was also reported in previous studies [3,11], our study further demonstrated that the phenological behaviour of rubber trees was distinguished from that of croplands in the study area.Therefore, it is likely to effectively separate rubber trees from natural forests and croplands based on the differences in phenological information in vegetation indices.This finding is inconsistent with Senf et al. [10], who failed to differentiate effectively between rubber plantations and natural forests based on the MODIS EVI-based phenology information and explained that the phenology of natural forests partly resembled that of rubber trees in Xishuangbanna.However, phenological ground observations showed that the biological rhythms of rubber trees and natural forests were distinctly different in Xishuangbanna [28].This inconsistency may be only explained by the issue of mixed pixels within MODIS data, which confused the MODIS EVI-based phenology information between rubber trees and natural forests [11].
The timing of rubber tree defoliation and foliation varies geographically, and the defoliation-foliation transition process of rubber trees is brief.The findings described here showed that rubber trees defoliated between late December and early February of the next year and subsequently foliated until the middle of March in the Xishuangbanna region.However, later defoliation and foliation time of rubber trees was reported by Dong et al. [11] in the Danzhou Region of Hainan Island, China, where rubber trees defoliate in late February through March and rapidly foliate in late March through April.Li and Fox [3] observed earlier defoliation and later foliation of rubber trees on the mainland of Southeast Asia, where the defoliation of rubber trees occurs between November and early April and the foliation occurs from mid-April to late May.Therefore, temporal inconsistencies in the defoliation and foliation of rubber trees in different regions must be considered for rubber plantation mapping at the trans-regional scale.Even in the same climatic region, climate variability, rubber stand ages and structure, and topography could also cause shifts in phonological phases of rubber trees.Further studies on the effects of these local driving forces on phenological phases of rubber trees should be investigated in the future.
The image differencing of five spectral vegetation indices retrieved in the defoliation and foliation stages of rubber trees more accurately separated rubber plantations from the other land cover types.Overall accuracies of the resultant rubber plantation map vary from 92% to 96% and the kappa coefficients from 0.84 to 0.92 (Figure 4).This result is comparable or even slightly superior to that of Dong et al. [11], who produced an overall accuracy of 92% with a kappa coefficient of 0.88.Likewise, the accuracy achieved in this study is much better than those obtained by Senf et al. [10] and Dong et al. [12].Senf et al. [10] mapped rubber plantations in Xishuangbanna with an overall accuracy of 73.5% using multi-spectral phenological metrics from the MODIS vegetation indices time series.Dong et al. [12] identified rubber plantations on Hainan Island, China, with an overall accuracy of 85% by combining PALSAR and MODIS imagery.Obviously, the phenology-based vegetation index differencing presented in this study can improve rubber plantation classification significantly.Previous studies indicated that combining phenological information with vegetation indices can produce higher classification accuracy [22].The findings described here support the research showing that incorporating phenological information and image differencing can help better differentiate between rubber plantations and non-rubber land cover [3,10,11,22].
Phenology-based vegetation index differencing is highly dependent on the availability of good-quality images in specific phenological phases (i.e., defoliation and foliation).Fortunately, the defoliation-foliation process of rubber trees in Xishuangbanna occurs during the dry season.Thus, the probability of obtaining cloud-free observations from Landsat-like multispectral sensors with high spatial resolution is greater.Additionally, high-frequency temporal satellite imagery (e.g., MODIS and China's Feng-Yun-3A) provide an alternative data source for mapping and monitoring rubber plantations using phenology-based vegetation index differencing.Although ubiquitous mixed pixels within these data with coarse spatial resolution affect the identification of rubber trees in the rapid foliation phase [11], the blending of fine spatial resolution imagery (e.g., Landsat and SPOT) and high-frequency temporal imagery could potentially overcome this obstacle to capture a more specific intra-annual phenology at high resolution in time and space [53,54].

Uncertainties of Rubber Plantation Areas Estimated by Phenology-Based Vegetation Index Differencing
Theoretically, the actual area of rubber plantations in the study area is constant at any given time.However, differencing the different vegetation indices produced different estimated rubber plantation areas in this study.Although high classification accuracy was achieved in this work, overestimates or underestimates of the true rubber plantation areas still existed because of the omission or commission errors.For instance, intercropping young rubber trees with other upland cash crops, such as rice, maize, pineapple, tea or bananas, in tropical regions can bring about the omission of young rubber trees.However, the mosaic distribution of rubber plantations and second-growth forests tend to cause misclassification of second-growth forests [3].There were significant differences among these vegetation indices in the capability to distinguish rubber plantations from other land cover types.Differencing the TCG image produced the largest omission errors of 13%, and differencing the NDVI and TCG images produced the largest omission errors of 4% (Table 1).However, differencing the NDMI image produced the lowest omission errors of 6%, whereas differencing the ARVI image produced the lowest commission errors of 1%.Consequently, the TCG image differencing responded to the lowest classification accuracy of rubber plantations, whereas the ARVI and NDMI image differencing produced the highest one.Under the combined atmospheric and topographic effects, differences in sensitivity to vegetation conditions in the defoliation and foliation stages of rubber trees existed among the examined vegetation indices.The TCG was atmospherically uncorrected and topographically sensitive so that it had less ability to separate rubber plantations from natural forests and croplands (Figure 2).The NDVI and NDMI computed from the normalized combination of spectral bands can reduce atmospheric and topographic interference to some degree [38,43,44].In contrast, the NDMI is more robust to atmospheric effects than the NDVI [43].Therefore, the NDMI presented a larger separability and performed better than the NDVI (Figure 2 and Table 1).The EVI and ARVI were developed to minimize atmospheric aerosol effects [39][40][41][42].However, the EVI is a non-band-ratio based vegetation index, which is relatively sensitive to topographic effects.It exhibited moderate separability and classification accuracy (Figure 2 and Table 1).The ARVI is more robust to topographic and atmospheric aerosol effects than the other vegetation indices, especially in tropical mountainous regions of high atmospheric aerosol content, where are often contaminated by soot from slash-and-burn agriculture [41].The preceding results showed that the ARVI performed the highest separability and classification accuracy in our study area.
Higher overall accuracy and a higher kappa coefficient in the ARVI-based maps indicate that the rubber plantation area estimate of 3707.28 km 2 is more reliable.The estimate is much smaller than the areas reported in several recent works [10,55,56], which estimated that the area of rubber plantations in Xishuangbanna reached 4240 km 2 [55], 5010 km 2 [56] and 5907 ± 788 km 2 [10] in 2010, respectively.However, the reference data obtained from the official statistics of Xishuangbanna showed that the total area of rubber plantations in the prefecture was only 2940.28 km 2 by the end of 2013 [27].The figure reported in the official statistics is far smaller than the rubber plantation areas derived from satellite imagery.In comparison, our results are closer to the official statistics.The difference between our results and the official statistics is likely caused by underreporting of rapidly growing and relatively scattered small rubber plantations in the private sector.It can also be deduced that previous studies overestimated to varying degrees the area and expansion rate of rubber plantations in the study area.Therefore, the underlying environmental and ecological effects of rubber expansion in the study area need urgently to be re-examined in the future through more accurate spatio-temporal trajectories of rubber expansion.The success of the phenology-based vegetation index differencing in delineating rubber plantations showed the use of defoliation-foliation image differencing strengthened the phenological change in the remotely sensed data and greatly improved land cover classification.This method may be transferred to a variety of landscapes where phenological differences always occur.

Conclusions
This study presented a phenology-based vegetation index differencing method for rapidly mapping rubber plantations and explored its ability to delineate rubber plantations in Xishuangbanna in Southwest China, a hotspot of global biodiversity and rubber plantation expansion.The Landsat 8 OLI-based phenological trajectories showed unique phenological characteristics of rubber plantations that were distinguishable from natural forests and croplands.Two specific phenological phases (i.e., defoliation and foliation) of rubber trees were accurately identified from the temporal profiles of several spectral vegetation indices.Substantial defoliation of the rubber trees occurred between late December and early February, and pronounced foliation occurred between early February and mid-March.The bitemporal vegetation index differences that were acquired during the key phenological phases (i.e., almost complete defoliation and full foliation) were used to rapidly and effectively delineate a map of rubber plantations.Our method greatly improved the classification accuracies for rubber plantations and presented more reliable rubber plantation areas.Furthermore, the findings herein indicated that the recent remote sensing-based rubber expansion in Xishuangbanna could have been overestimated in previous studies.Therefore, developing a regular monitoring scheme for evaluating the environmental and ecological impacts of rubber tree expansion in the study area should be urgently considered.

Figure 1 .
Figure 1.Topography of the study area.The black dotted lines show the GPS tracks of the field trips between September and November of 2014.The red polygons denote the polygons of interest (POIs) for the validation of the resulting rubber plantation maps.

Figure 2 .
Figure 2. The temporal profiles of the Landsat 8 OLI NDVI, EVI, ARVI, NDMI, and TCG time series for the rubber plantations, natural forests and croplands.Random samples of 2000 pixels were individually extracted from 134 polygons of interest (POIs) for the rubber plantations, 110 POIs for the natural forests and 136 POIs for the croplands.The boxplot depicts the median (the central horizontal line), the 25th and 75th percentiles (the edges of the box) and the most extreme data points (the whiskers).Rubber plantations present two unique phenology stages: defoliation (the long light yellow column) and foliation (the long light green column).

Figure 3 .
Figure 3.The false colour composition map (R/G/B = Band 7/5/4) of the Landsat 8 OLI images, including (a) the image during the nearly complete defoliation stage of the rubber trees on 9 February 2014 and (b) the image during the full foliation stage of the rubber trees on 13 March 2014.The zoom-in images of a typical area show that the rubber plantations are easily discernible as greyish-purple patches in the defoliation stage in image (a) and as light green patches during the foliation stage in image (b).However, the natural forest is always green during the defoliation (a) and foliation (b) stages of rubber trees.Major land cover classes of interest were marked in the zoom-in images, including rubber plantation ("A"), natural forest ("B"), and cropland ("C").

Figure 4 .
Figure 4. Difference images (a) of the five spectral vegetation indices from the two images acquired on 13 March 2014 (cloud-contaminated pixels were filled with the cloud-free pixels from 29 March 2014) and on 9 February 2014, and the rubber plantation maps (b) obtained by thresholding these difference images.The areas of the rubber plantations (ARP) in these maps were calculated and are shown in (b).

Table 1 .
Error matrices (pixel counts) for five vegetation index differencing methods.Non-rubber includes forest, cropland, built-up land, water body and other land cover types.