Identifying Establishment Year and Pre-Conversion Land Cover of Rubber Plantations on Hainan Island, China Using Landsat Data during 1987–2015

Abstract

Knowing the stand age of rubber tree (Hevea brasiliensis) plantations is vitally important for best management practices, estimations of rubber latex yields, and carbon cycle studies (e.g., biomass, carbon pools, and fluxes). However, the stand age (as estimated from the establishment year of rubber plantation) is not available across large regions. In this study, we analyzed Landsat time series images from 1987–2015 and developed algorithms to identify (1) the establishment year of rubber plantations; and (2) the pre-conversion land cover types, such as old rubber plantations, evergreen forests, and cropland. Exposed soil during plantation establishment and linear increases in canopy closure during non-production periods (rubber seedling to mature plantation) were used to identify the establishment year of rubber plantations. Based on the rubber plantation map for 2015 (overall accuracy = 97%), and 1981 Landsat images since 1987, we mapped the establishment year of rubber plantations on Hainan Island (R² = 0.85/0.99, and RMSE = 2.34/0.54 years at pixel/plantation scale). The results show that: (1) significant conversion of croplands and old rubber plantations to new rubber plantations has occurred substantially in the northwest and northern regions of Hainan Island since 2000, while old rubber plantations were mainly distributed in the southeastern inland strip; (2) the pattern of rubber plantation expansion since 1987 consisted of fragmented plantations from smallholders, and there was no tendency to expand towards a higher altitude and steep slope regions; (3) the largest land source for new rubber plantations since 1988 was old rubber plantations (1.26 × 10⁵ ha), followed by cropland (0.95 × 10⁵ ha), and evergreen forests (0.68 × 10⁵ ha). The resultant algorithms and maps of establishment year and pre-conversion land cover types are likely to be useful in plantation management, and ecological assessments of rubber plantation expansion in China.

1. Introduction

The rubber tree (Hevea brasiliensis) is native to the Amazon basin, but historically it has been widely cropped in tropical regions within equatorial zone (between 10°N and 10°S) [1]. To meet the increasing international demand for natural rubber, rubber plantations have expanded rapidly since the 20th century throughout Southeast Asia and China, where more than 80% of the world’s natural rubber is produced [2]. According to the Food and Agriculture Organization (FAO) of the United Nations Global Forest Resources Assessment (FRA) 2010 report, the global area of rubber plantations has steadily increased by 25% during the past two decades [3]. The planting of rubber trees has brought great benefits to the plantation owners and rubber-relevant industries. Knowledge of the spatial extent and stand age (or establishment year) of rubber plantations at regional scales is of great interest to plantation managers who are concern with the yield of natural rubber and rubberwood. However, the rapid expansion of monoculture rubber plantations also created conflicts between forest management and conservation [4]. In China, rubber plantations in Yunnan province have been expanded to marginal growing environments of as far north as 24°N and higher elevations exceeding 1300 m and steeper slopes. Therefore, data on stand age (or establishment year) and pre-conversion land cover of rubber plantations are necessary to correctly assess the impacts of land conversion on biodiversity, carbon sequestration, and water provisioning [5,6,7,8], and are critical to conservation policy development [9].

Stand age or establishment year of rubber plantations is usually obtained through on-site survey (interviewing owners or estimating from the length of tapping line), querying planting record data from farms or statistical yearbooks. Although data collected by these methods has high accuracy, it contains limited spatial distribution information and is difficult to be applied on a large scale. Remote sensing is a viable approach to determine the stand age or establishment year of rubber plantations at large spatial scales. Previous studies have used optical images [10,11,12,13,14,15,16,17,18], synthetic aperture radar (SAR) images [19], or their combination [20] (

Table 1. Selected peer-reviewed references related to remote sensing estimates the stand age of rubber plantations.

Approach	Algorithm	Optical			SAR	Optical+SAR
		VHR/HR	MR	CR
Statistical Models	Regression	Charoenjit and Zuddas (2015)	Chen and Cao (2012); Suratman and Bull (2004); Chen et al. (2012)
	Others *		Razak and Shariff (2017); Chen et al. (2012)
Image Classification	Supervised	Koedsin and Huete (2015)	Li and Fox (2011)	Li and Fox (2012)	Trisasongko (2017)
	Objected	Chen and Yi (2016)	LIU et al. (2012)
	Hybrid		Dibs and Idrees (2017)
LULCC Detection	Decision Tree		Beckschäfer (2017)			Kou and Xiao (2015)

* VHR: very high-resolution imagery; HR: high-resolution imagery; MR: moderate-resolution imagery; CR: coarse-resolution imagery; SAR: synthetic aperture radar; LULCC: Land use and Land cover change; Others: include artificial neural network and leaf area index histograms.

) to map changes in land cover. Historically, moderate-resolution (MR) optical images were the main data sources for estimating the stand age of rubber plantations [10,11,12,13,14,15,16,21], followed by high-resolution or very high-resolution (HR/VHR) optical imagery [4,17,18], and SAR data and optical-SAR integration [20] (Table 1).

The approaches used for estimating stand age of rubber plantations can be grouped into three categories: (1) statistical models; (2) image classification; and (3) land use and land cover change (LULCC) detection (Table 1). Statistical models like multivariate regression [11,15,18,21] and artificial neural networks [21] have been used to establish the relationships between the stand age of rubber plantations and optical-based parameters (spectral reflectance and/or vegetation indices) or SAR-based backscatter coefficients. The image classification approach usually divided the rubber plantations (usually 25–30 years of economic life cycle) into two or three age groups (e.g., ≤6, >6 years), then performed classification according to spectral metrics of the different groups. Different classification methods such as supervised classification [1,17,19], object oriented [16], and hybrid approaches [12] have been employed to identify the stand age of rubber plantations. Both the statistical model and imagery classification approaches assumed a good relationship between stand age and spectral/structural features of the plantation canopy. As tree volume and biomass growth usually have a sigmoid-shaped curve as a function of stand age, large estimation errors (RMSE = 3.7–8.2 years) were found by statistical model approaches [11,15,21]. Image classification approaches were usually unable to classify stands beyond two or three age groups and had wide ranges in accuracy (44–98%, dependent on classification algorithms, input images, and number of age groups) [12,18,19,20].

The most promising way to identify the establishment year of rubber plantations was to detect LULCC events [10,20]. New plantation establishment involves clearing land and transplanting rubber seedlings. A common practice in land clearing and preparation is to remove all vegetation (e.g., natural forests, old rubber plantations, and crops), and then to construct a gully and dike system in the flat lands (slope < 5°), or contour terraces on hillsides. Then, rubber seedlings are simultaneously transplanted into these patches with a spacing of about 3–4 m apart within rows and 5–8 m between rows with a planting density of ~450 trees/ha [15,18,22]. Rubber plantations have about seven years of non-productivity, of which the first four years or so have an open canopy and can be intercropped with pineapple, banana, sugarcane, and peanuts to increase economic returns [1]. Rubber latex harvesting usually begins at the seventh year when the canopy has totally closed. Latex harvest lasts for about 25 years, then trees are harvested for timber. Rubber trees rapidly defoliate (complete or near-complete defoliation) and foliate (RDF) in areas where have distinct dry seasons [14,23,24]. In the northern hemisphere, rapid defoliation usually occurs in mid-February, and the canopy quickly becomes dense in mid-March. The leaf area index (LAI) during this period (one month or so) can drop from around 3.3 to 0.9 then grow to 3.8 [24]. The rapid canopy change during the RDF period is a key indicator in differentiating rubber plantation from evergreen forest [23], while exposed bare soil during initial plantation establishment stage is valuable in identifying their stand age.

The land surface water index (LSWI) is particularly sensitive to bare soils and less affected by clouds [10], and has been used to estimate the stand age of rubber plantations [10,20]. Kou and Xiao first mapped the stand age of rubber plantations in Xishuangbanna, Yunnan Province, China using the first year when a rubber plantation had minimum LSWI less than zero (LSWI_min < 0) during the defoliation period (from middle January to early February), and then grouped them into ≤5, 5–10, and ≥10 year age groups. Beckschäfer improved Kou and Xiao’s algorithm by detecting the last year when a rubber plantation has LSWI_min < 0 in leaf-on season to estimate the establishment year of rubber plantations annually in Xishuangbanna since 1988. Using LSWI_min during the leaf-on season resolved classification issues associated with defoliation, which degraded LSWI significantly (less than, or close to, zero) [23]. However, both intercropping and weeding would greatly increase the chance of exposing bare soil or producing dry vegetation, which significantly degraded LSWI_min during the leaf-on season to less than or close to zero. Furthermore, severe natural disasters such as typhoons and drought may alter the plantation canopy or lead to leaf shedding [11] and, therefore, lower LSWI significantly. Low LSWI in the leaf-on season leads to an underestimation of plantation establishment years. However, the integration of biophysical metrics like increasing canopy closure in non-production period is possible to significantly reduce these uncertainties.

The objectives of this study were to: (1) develop a new algorithm to map the establishment year of rubber plantation by integrating LULCC events and biophysical metrics with time series Landsat imagery; (2) generate an establishment year map of rubber plantations on Hainan Island, China, and investigate the pattern of rubber plantation expansion since 1987; and (3) explore the sources, time, area, and spatial extent of land conversion during rubber plantation expansion over the past 30 years. The resultant algorithms, maps of plantation establishment year, and pre-conversion land covers are useful for practical land management, ecological assessment, and conservation policy development.

2. Material and Methods

2.1. Study Area

Hainan Island (19°20′N–20°10′N, 108°21′E–111°03′E, Figure 1a) is the second largest island in China and has an area of about 3.4 million ha. The topography of the island is mountainous in the interior with lowlands along the coast. Wuzhi Mountain (1867-m above mean sea level, AMSL) is the highest mountain on the island. The climate varies from tropical to subtropical. The annual mean temperature is approximately 23–25 °C, and monthly temperature varies between ~16 °C in January and ~30 °C from May to July. The average annual precipitation is 1500 to 2000 mm, 80% of which occurs during the rainy season from May to October. However, the spatial distribution of precipitation can be as high as 2400 mm in the central and eastern areas, and as low as 900 mm in the coastal areas of the southwest.

In the 1950s, the island began to plant large areas of rubber trees to meet China’s national defense and economic needs. It is the second largest natural rubber production base in China, with an area of 5.4 × 10⁵ ha in 2015 [25]. According to Hainan Island National Land Cover Dataset in 2010 [26], forest (include rubber plantations), cropland, grassland, water, and built-up account for 63.4%, 25.7%, 3.3%, 3.8%, and 3.6% of the total land area, respectively.

2.2. Landsat Imagery and Processing

Hainan Island is covered by four Worldwide Reference System (WRS-2) path/row (123/46, 123/47, 124/46, and 124/47, respectively). A total of 1981 Collection 1 Tier Top-of-atmosphere reflectance of Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+), and Operational Land Imager (OLI) images from 1987 to 2015 were available in Google Earth Engine (GEE) as image collections. The Landsat Collection 1 includes Level-1 Precision Terrain (L1TP) processed data that have well-characterized radiometry, are inter-calibrated across the different Landsat sensors, and are considered suitable for time-series processing analysis. Poor quality Landsat images, caused by clouds and shadows, were identified using cfmask [27,28] and were flagged in the Quality Assessment (QA) band. We considered observations unaffected by clouds, shadows, and ETM+ scan-line-off strips as good quality observations. The statistical information of Landsat imagery on Hainan Island is presented in Figure 2, and annual cloud-free observations is shown in Figure S1.

For each good quality Landsat TM/ETM+/OLI imagery, we calculated the normalized difference vegetation index (NDVI) [29], enhanced vegetation index (EVI) [30], and LSWI [31,32] using Equations (1)–(3):

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR}+{\rho }_{red}}$$

(1)

$$EVI=2.5\times \frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR}+6\times {\rho }_{red}-7.5\times {\rho }_{blue}+1}$$

(2)

$$LSWI=\frac{{\rho }_{NIR}-{\rho }_{SWIR}}{{\rho }_{NIR}+{\rho }_{SWIR}}$$

(3)

where ${\rho }_{blue}$, ${\rho }_{green}$, ${\rho }_{red}$, ${\rho }_{NIR}$, and ${\rho }_{SWIR}$ are blue (450–520 nm), green (520–600 nm), red (630–690 nm), near-infrared (NIR, 760–900 nm), and shortwave infrared (SWIR, 1550–1750 nm) bands of Landsat TM/ETM+/OLI imagery, respectively.

Annual minimum and standard deviation of LSWI in the RDF period (January–March, hereafter referred to as LSWI_{RDF_min} and LSWI_{RDF_Std.Dev}, respectively) [23], annual minimum and average LSWI during the leaf-on season (April–December, hereafter referred as LSWI_{leafon_min} and LSWI_{leafon_avg}, respectively), annual maximum NDVI (NDVI_max) [23], and annual greenness frequency [33] (F_Greenness, calculated using both EVI and LSWI, see Section 2.4.3), were generated for 1987–2015 to map rubber plantations and identify the establishment year of rubber plantations and pre-conversion land covers (see Section 2.4). Annual LSWI_{RDF_min} and LSWI_{RDF_Std.Dev} were composited using two consecutive years (e.g., 2015 composite image was generated using images acquired in RDF period of 2014 and 2015), while annual LSWI_{leafon_min}, LSWI_{leafon_avg}, NDVI_max, and F_Greenness were generated using images acquired in the leaf-on period of the corresponding year.

2.3. Ground Reference Data

2.3.1. Data for Rubber Plantations Map in 2015

The region of interest (ROI) samples for mapping rubber plantations in 2015 initially came from our previous study, which mapped forest and rubber plantations in 2010 on Hainan Island [23]. These ROIs data were updated with geotagged landscape photos taken in February 2015 and October 2016, and Google Earth VHR imagery acquired circa 2015, by adding, removing, or adjusting polygons to fit the rubber plantation boundaries. The non-rubber ROI data consisted of forest, water, cropland, and built-up areas, and they were come from our recent work which studies the spatial-temporal changes of forest on Hainan Island between 2007 and 2015 [34]. In total, we had 804 ROIs (27,539 pixels at 30-m spatial resolution) for rubber plantations and 469 ROIs (81,858 pixels at 30-m spatial resolution) for non-rubber. These ROI data (Figure 1a) were randomly divided into training (30%) and validation (70%) datasets using NOAA/NOS/NCCOS/CCMA Biogeography Branch’s Design Tool for ArcGIS.

2.3.2. Data for Establishment Year Map of Rubber Plantation

ROI samples for identifying establishment year of rubber plantation were obtained in three ways. First, we obtained six plantation vector maps of Hainan State Farm, which were produced using HR imagery (SPOT-5 acquired around 2005), field interviews, and querying historical land-use records from Hainan State Farm. These maps included detailed information for each plantation including establishment year, clones, farmer identification, and tapping scheme. Second, we obtained soil sampling data from the Fertilizer Application department of the Rubber Research Institute, Chinese Agriculture of Tropical Academy Science (RRI, CATAS). Establishment year and clone information of these rubber plantations were collected by reviewing the production departments of the Hainan State Farm or interviewing the private owners. Third, we visually interpreted the establishment year of rubber plantation with Google Earth VHR imagery and time series Landsat imagery.

ROI samples from the first and second ways (vector maps and soil sampling) were used to find rubber plantations established before 2007, while the third way (visual interpretation) was used to find rubber plantations established during 2007 and 2011. Rubber plantations established during 2012 and 2015 were seedlings, which cannot be mapped accurately by either optical and SAR imagery due to complex understory conditions and, therefore, were excluded. The vector maps from the first source was updated using the 2015 rubber plantation map, but only 10% of them were randomly selected as they were spatially clustered. Finally, we obtained 522 ROI samples with establishment year information. The spatial distribution of these ROI samples was presented in Figure 1b and plantation statistics by year was presented in Figure S2. The clustered ROIs in Figure 1b were selected from the vector maps. The spatial distribution of field data was relatively consistent with the spatial distribution of rubber plantations on Hainan Island [23]. We randomly selected five rubber plantations that were established in 1990, 1995, 2000, 2005, and 2010 from the 522 ROI samples for algorithm training, and the rest were used for validation.

2.3.3. Data for Pre-Conversion Land Cover Map

A stratified sampling was conducted to generate random ROI samples based on mask layers generated by the maps of pre-conversion land cover and establishment year for accuracy assessment. First, old rubber plantations that were established before 1990 and small rubber plantations with connected pixel fewer than 10 were masked out in advance. Pre-conversion land cover for old rubber plantations established before 1990 was difficult to identify accurately with limited amount Landsat TM images which starts in 1987 here, and small rubber plantations usually contain mixed pixels. Then, 50 random square ROI samples with areas of 0.25 ha for evergreen forests, old rubber plantations, and cropland were respectively created using GEE’s random function and the mask layers. We first screened and excluded mixed ROI, which have multiple land cover, multiple conversion sources, and different establishment years with Google Earth VHR images. We determined pre-conversion land covers for these ROIs samples by respectively observing annual minimum/average NDVI/EVI/LSWI composited from Landsat images in GEE platform, and double checked in TimeSync, a Landsat series visualization tool and web application (Version 3) (http://timesync.forestry.oregonstate.edu/). After these processes, 28, 30, and 37 relatively pure ROIs were identified that had been converted from evergreen forests, old rubber plantations, and cropland, respectively, to be used for map accuracy assessment (Figure 1c).

2.4. Workflow and Algorithms

2.4.1. Mapping Workflow

The workflow is presented in Figure 3. We generated rubber plantation map from a forest base map of 2015 (derived from Phased Array type L-band Synthetic Aperture Radar-2 (PALSAR-2) and Landsat images) and time series Landsat data, then identified the establishment year of rubber plantations using a LULCC- and biophysical-based algorithm. Finally, we tracked the pre-conversion land cover types, which were divided into old rubber plantations, evergreen forests, cropland, or unknown.

2.4.2. Algorithm for Mapping Rubber Plantations in 2015

To generate the rubber plantation map of 2015, we first identified forest, then distinguished deciduous rubber plantations from the forest base map using phenological metrics [23,35,36]. The 2015 forest base map came from our previous study on the annual spatial-temporal of forest on Hainan Island during 2007–2015 using PALSAR/PALSAR-2 and time series Landsat TM/ETM+/OLI images [34,37]. This forest base map had an overall accuracy (OA) of ~96% when validated with ground reference data [34] and, therefore, can serve as a reliable base map for mapping rubber plantation. Then, deciduous rubber plantations were identified from this forest base map using the criteria of dense canopy in the growing season (NDVI_max > 0.85) and the criteria of rapid defoliation and foliation (LSWI_{RDF_min} < 0.1, and LSWI_{RDF_Std.Dev} > 0.08, See Figure S3 for justification). The LSWI_{RDF_min} and LSWI_{RDF_Std.Dev} were composited with Landsat ETM+/OLI images acquired in the RDF period of 2014 and 2015. A 3 × 3 median filter was performed on the resultant map to reduce the salt-and-pepper noise.

2.4.3. Algorithm for Identifying Establishment Year of Rubber Plantations

We used two features to identify the establishment year of rubber plantations: (1) exposed soils during initial plantation establishment, which have low vegetation indices (NDVI, EVI, and LSWI); and (2) increasing canopy closure in non-production period of rubber plantations (from rubber seedling to maturity plantation), which has a linear increase in vegetation indices. Figure 4 shows the temporal profile of Landsat-based annual mean and minimum NDVI, EVI, and LSWI during leaf-on season (April to December) in a rubber plantation, where old rubber trees were removed, and rubber seedlings were transplanted in 2001. All vegetation indices suddenly dropped in 2001, but gradually recovered within about seven years (non-production period). LSWI was found to be more sensitive to bare soil [10,20,38], and was thus selected to identify establishment year of rubber plantations (Figure 4). Specifically, we selected LSWI_{Leafon_min} to track bare soil and the slope of LSWI_{Leafon_avg} against year during the non-production period (hereinafter referred as S_LSWI) to track the change of canopy closure. We marked the last year when a rubber plantation has LSWI_{Leafon_min} < T₀ as Y_i, and then search ahead for N_s years of seedling period to find a proper year as the establishment year of rubber plantation. Within the period of (Y_i-N_s, Y_i), we set the year that held the lowest LSWI_{Leafon_min} and met LSWI_{Leafon_min} < T₀ and S_LSWI > T_Slope as the establishment year of rubber plantations. Here, T₀ and T_Slope are the thresholds of LSWI_{Leafon_min} and S_LSWI, respectively.

The thresholds of T₀ and T_Slope were determined by the statistical value of the training ROI samples. The LSWI_{Leafon_min} was always less than zero in the year that rubber plantation was established, and remained relatively small in the first few years (Figure 5). The S_LSWI ranged from 0.03 to 0.05, with a mean value (µ) of about 0.04 (Table S1). We set T₀ as 0 and the rounded value of mean µ-σ from the five rubber plantations as the T_Slope (0.03 here), where σ is the standard deviation. However, the rapid recovery of understory in young rubber plantations and frequent presence of clouds in the tropics made it difficult to acquire ideal observations (un-vegetated during plantation establishment and linear increase of LSWI_{Leafon_avg}). Thus, we employed the algorithm again with slightly higher thresholds of 0.1 for T₀ and 0.02 for T_Slope to address the less optimal conditions (Table S1). After executing the algorithm twice, we still had many rubber plantations with an unidentified establishment year because they were established earlier than 1987 when Landsat TM imagery was lacking. We set the establishment year of these rubber plantations as 1986 if they had high greenness frequency in leaf on season for all years between 1987 and 2015. Otherwise, the establishment year was marked as unknown. We first calculated greenness in a year using Equation (4) [33].

$$Greenness=\{\begin{array}{l}1\hspace{1em}LSWI>0andEVI>0.2\\ 0\hspace{1em}\hspace{1em}\hspace{1em}Othervalues\end{array}$$

(4)

Then we calculated greenness frequency using Equation (5):

$$F_{Greenness}=\frac{{\displaystyle \sum N_{Greenness}}}{{\displaystyle \sum N_{Total}-N_{Bad}}}\times 100$$

(5)

where F_Greenness is the greenness frequency scaled between 0 and 100, N_Greenness is the number of observations with LSWI > 0 and EVI > 0.2, N_Total is the total number of observations, and N_Bad is the number of bad observations (e.g., clouds, shadows, and ETM+ scan-line-off strips) [39]. We used F_Greenness > 90% as the criteria (see Figure S4 for justification) because rubber plantations on Hainan Island frequently suffered canopy damage from severe typhoons, which may cause rubber plantations to not meet a higher greenness criterion.

2.4.4. Algorithm for Tracking Pre-Conversion Land Cover of Rubber Plantations

After we identified the establishment year of rubber plantations, we tracked pre-conversion land cover by analyzing the time-series vegetation indices. The conversion of grassland to rubber plantations was limited because grassland accounts for only 3.3% of the total land area [26], and the conversion of water and built-up areas to rubber plantations were highly unlikely, so our study focused on the conversion of evergreen forests, old rubber plantations, and cropland to rubber plantations. First, we identified the pre-conversion land cover as cropland or forest by observing LSWI_{leafon_min} of N_x years before plantation was established. A cropland usually has some observations with LSWI_{leafon_min} < 0 in a year due to exposed soils after harvest and during the fallow period. Here we set N_x as 5 to avoid classifying eucalyptus industrial plantation (usually having ca. six-year rotation) as cropland. Second, we divided the pre-conversion forest cover into old rubber plantations and evergreen forests. We used the annual NDVI_max, LSWI_{RDF_min}, and LSWI_{RDF_std.dev} of the five years prior to plantation establishment and classified the pre-conversion forest cover as old rubber plantation if it met NDVI_max > 0.85, LSWI_{RDF_min} < 0.1, and LSWI_{RDF_std.dev} > 0.08 (Figure S3). If it did not meet that condition, then we classified the pre-conversion forest cover as evergreen forests. Landsat TM imagery for Hainan Island is not available before 1987, so pre-conversion land cover types were labeled as unknown for rubber plantations established before 1988.

2.5. Accuracy Assessment

The error-adjusted matrix proposed by Olofsson et al. in terms of proportion of area and estimates of producer’s accuracy (PA), user’s accuracy (UA), and overall accuracy (OA) with 95% confidence intervals was used to evaluate the accuracy of our rubber plantation map in 2015 [23,34]. The error-adjusted matrix was a statistically robust and transparent approach for assessing accuracy, a method, which has recently become common practice [40,41,42].

The accuracy of the establishment year map of rubber plantations was evaluated using ROIs at both the pixel and plantation scale. Each ROI (polygon boundary) was assigned an independent establishment year as rubber trees were planted simultaneously. However, each ROI contained many pixels (30 m spatial resolution), and those pixels may have the same, or very closely estimated establishment year when the algorithm performs well. We compared estimated and observed establishment year for all the ROI pixels (pixel scale, ignoring plantation boundary). We also calculated mean estimated establishment year for each ROI (plantation scale, processing unit was plantation), and then calculated the mean and standard deviation of the mean estimated establishment year group by year. The group mean and standard deviation of observed and estimated establishment year of rubber plantations were plotted.

We assessed the accuracy of the pre-conversion land cover map using all random ROI data, which were visually interpreted with Google Earth VHR imagery and time series Landsat imagery. Olofsson’s error-adjusted matrix was built and PA, UA, and OA for the pre-conversion land cover map was reported.

2.6. Spatial and Areal Analysis of Plantation Establishment Year and Pre-Conversion Land Covers

We generated annual maps of the establishment year of rubber plantations and plotted a histogram of the observed and estimated rubber plantation area since 1987. We analyzed the age structure of existing rubber plantations in 2015 by dividing the stand age into groups of <5, 5–10, 10–15, 15–20, 20–25, and >25 years. The area and proportion of rubber plantations in each age group were explored. We also divided the island elevation into groups of 0–50, 50–100, 100–200, 200–300, 300–400, 400–600, >600 m, and plotted a histogram to explore the elevation distribution of plantation establishment at five-year intervals (<1990, 1991–1995, 1996–2000, 2000–2005, and >2005). Furthermore, the elevation and slope variation of plantation establishment since 1986 were explored annually. We calculated and plotted the annual and total area of rubber plantations converted from evergreen forests, old rubber plantations, and cropland since 1988.

3. Results

3.1. Rubber Plantation Map and Accuracy Assessment

The spatial distribution of rubber plantations in 2015 on Hainan Island is presented in Figure 6a. Rubber plantations are mainly located in the western, northern, and east-central regions. The densest rubber plantations are in Danzhou City and its surrounding regions, which comprise the largest natural rubber production base on Hainan Island.

The PA, UA, and OA for the resultant map of rubber plantations were 95.18% (±0.16%), 88.20% (±0.29%), and 96.99% (±0.07%) when validated with ground reference data (Table S2). The total mapped area of rubber plantations in 2015 was 5.83 × 10⁵ ha, and the adjusted area was (6.30 ± 0.04) × 10⁵ ha.

3.2. Establishment Year Map of Rubber Plantations and Accuracy Assessment

The establishment year of rubber plantations on Hainan Island is shown in Figure 6b. Many rubber plantations have been established across the island since 2000, especially after 2003. New rubber plantations were mainly concentrated in the northwest and northern regions such as Baisha County, Danzhou City, Lingao County, and Chengmai County. Old rubber plantations established before 1995 were mainly distributed in the central western regions, like Qiongzhong and Tunchang County, and along the coastal strip from south to northeast.

Estimated and observed establishment year of rubber plantations with all ROI pixels is presented in Figure 7a. Pixel density (plotting based on the number of pixels which have same observed and estimated establishment year) is shown in different colors, and contours are pixel density boundary (jointed points with equal pixel density). Pixels clustered along the 1:1 line have hotter colors, and pixel density contours are mostly symmetrical to the 1:1 line. The R² and RMSE of the observed against estimated establishment year is 0.85 and 2.34 years, respectively. However, a few pixels are scattered far from the 1:1 line. At the plantation scale (Figure 7b), the mean estimated establishment year of rubber plantations against the observed year are closely distributed along the 1:1 line. The R² of the fitted line is 0.99, and the RMSE is 0.54 years. Slight underestimation was observed for the old rubber plantations that was established prior to 2000.

At the annual scale, the estimated area of newly established rubber plantations by our map is close to the estimates provided by agricultural census data (Figure 8a). However, there were a few years (1987–1990, 1992, 1997, and 2009) in which there were large differences between our area estimate and the agricultural census.

3.3. Spatial-Temporal Distribution of Rubber Plantation Expansion on Hainan Island

At the end of 2015, rubber plantations on Hainan Island with stand ages of 4–5, 6–10, 11–15, 16–20, 21–25, and >25 years accounted for 9.83% (55,867 ha), 17.29% (98,301 ha), 14.82% (84,261 ha), 5.61% (31,909 ha), 4.99% (28,389 ha), and 47.36% (269,174 ha) of the total area, respectively (Figure 8b). About 0.09% (497 ha) of rubber plantations failed to estimate their establishment years.

On Hainan Island, rubber plantations were mainly distributed at elevations less than 400 m AMSL and most were established before 1990 or after 2000 (Figure 9a). The average elevation of rubber plantations established after 1986 was about 150 m AMSL, and the standard deviation of elevation was about 100 m (Figure 9b). The average slope of rubber plantations was about 7°, and the standard deviation of slope was about 5°. In the past 30 years, no significant variation in elevation and slope were found during rubber plantation expansion (Figure 9b).

3.4. Pre-Conversion Land Covers of Rubber Plantations on Hainan Island

Before the 1990s, evergreen forests and old rubber plantations were the main land sources for new rubber plantations, with annual conversion area of >2000 and ~1000 ha, respectively (Figure 9c). The area of rubber plantations converted from old rubber plantations increased during the 1990s (>2000 ha), followed by evergreen forests and cropland. Since 2001, the area of newly-established rubber plantations increased substantially. Old rubber plantations continued to be the predominant pre-conversion land cover type for 2003–2007, and in 2010 were about 10,000 ha. During this period, cropland was the second largest pre-conversion land cover type, with average annual conversion area of ~8000 ha. There was a boom in the amount of evergreen forests converted to rubber plantations between 2002 and 2007, with an annual conversion area of ~6000 ha. However, the conversion area from evergreen forests decreased rapidly after 2007 (Figure 9c). The total mapped area of rubber plantations converted from old rubber plantations, cropland, and evergreen forests since 1988 were 1.26 × 10⁵, 0.95 × 10⁵, 0.68 × 10⁵ ha, and adjusted area were (1.32 ± 0.1) × 10⁵, (0.96 ± 0.1) × 10⁵, (0.61 ± 0.1) × 10⁵ ha, respectively (Figure 9d).

The PA, UA, and OA of the pre-conversion land cover map were >86.60% (±0.05%), >84.21% (±0.07%), and >88.60% (±0.03%), respectively (Table S3). Rubber plantations converted from cropland were densely distributed in the northwest and northern regions, specifically in Danzhou City and Changjiang, Baisha, Lingao, and Chengmai Counties (Figure 10). In addition, we found that many scattered rubber plantations in the western region (Dongfang City and Ledong County) were previously cropland. Rubber plantations converted from old rubber plantations were mainly located in the northwest (e.g., Baisha County and Danzhou City) and northern-central regions (e.g., Qiongzhong, Tunchang, and Chengmai Counties) (Figure 10). In Baoting County, Qionghai City, and Wanning City we also found that many rubber plantations were established on old rubber plantations (Figure 10). Rubber plantations converted from evergreen forests were mainly distributed in Chengmai County, near the boundary between Dingan and Tunchang Counties, and in the interior mountainous area. However, these evergreen forests were relatively sparsely distributed. We did not observe any large hotspots where evergreen forests were converted to rubber plantations.

4. Discussion

4.1. Data and Algorithm for Identifying the Establishment Year of Rubber Plantations

Establishment year (stand age) of rubber plantations can be estimated by field interview, measuring the height of the tapping line, or reviewing the records/statistical data of local and national forestry departments. Compared with mapping using remote sensing, these approaches are difficult to frequently execute at large scales [25]. Remote sensing is a suitable technique to estimate the establishment year of rubber plantations at large scale, fine spatial resolution, and with an annual frequency. When mapped through remote sensing, the amount of satellite imagery was a very important factor. Most previous studies attempted to estimate the stand age (establishment year) of rubber plantations with a limited number of satellite images, such as single/mosaic optical images [11,15,19] or multi-temporal images [1] (

Table 2. Comparison of peer-reviewed studies for estimating establishment year (stand age) of rubber plantations.

Sensor	Method	Nimages/NPathRow	Age Resolution	Accuracy	Location	References
PALSAR-2	Image classification	1/1	3–8, >8	OA = 44–67%	West Java, Indonesia	Trisasongko (2017)
SPOT-5	Image classification	1/1	≤7, 7–25, >25	OA = 81–97%	Selangor, Malaysia	Dibs and Idrees (2017)
RapidEye	Image classificaiton	48/48	≤6, >6	-	Xishuangbanna, China	Chen and Yi (2016)
TM	Statistical model	1/1	-	RMSE = 3.7–5.3 years	Hainan, China	Chen et al. (2012)
TM	Statistical model	1/1	-	RMSE = 5.96 years	Hainan, China	Chen and Cao (2012)
TM	Statistical model	1/1	-	RMSE = 6.4–8.2 years	Selangor, Malaysia	Suratman and Bull (2004)
Pléiades	Image classification	1/1	≤7, 7–15, >15	OA = 94–97%	Phuket, Thailand	Koedsin and Huete (2015)
Thaichote	Statistical model	1/1	≤12, 12–18, >18	-	Rayong, Thailand	Charoenjit and Zuddas (2015)
TM	Image classification	11/11	≤2, 2–4, >4	-	Phuket, Thailand	Li and Fox (2011)
TM/ETM+, MODIS	Image classification	6/3(46) #	≤10, >10	OA = 85.2%	Xishuangbanna, China	LIU et al. (2012)
TM/ETM+	Image classification	29/1	≤6, >6	OA = 79%	Selangor, Malaysia	Razak and Shariff (2017)
MODIS	Image classification	29/4	≤4, >4	OA = 97–98%	Southeast Asia	Li and Fox (2012)
TM/ETM+, PALSAR	LULCC detection	226/1(1) #	≤5, 5–10, >10	OA = 85%	Xishuangbanna, China	Kou and Xiao (2015)
TM/ETM+	LULCC detection	276/1	Annual	RMSE = 2.5 years	Xishuangbanna, China	Beckschäfer (2017)
TM/ETM+/OLI	LULCC-and biophysical-based	1981/4	Annual	R2 = 0.85/0.99, RMSE = 2.34/0.54 years	Hainan Island, China	This study

#. Number in bracket is the number of images for MODIS and PALSAR, respectively.

). The accuracy of estimated stand age was susceptible to the image quality, acquisition time, and the information contained in the images. Here we used all Landsat TM/ETM+/OLI images that were available (1981 scenes; ~500 scenes per path/row), substantially more than the number of Landsat images used in the studies by Kou and Xiao (226 images per path/row) and Beckschäfer (276 images per path/row), (Table 2). Our establishment year map of rubber plantation in 2015 achieved high accuracy (R² = 0.85/0.99, RMSE = 2.34/0.54 years at the pixel/plantation scale) by using time series Landsat images, which were also illustrated by several previous studies [10,20,43,44,45,46].

Unlike statistical models and imagery classification, which assumed a good relationship between spectral/structural features of the plantation canopy and stand age, our algorithm was built upon the understanding of physical land cover change with dense time series Landsat images, including exposed bare soil and the linear increase of canopy closure from rubber seedlings to maturity. We used LSWI < 0/0.1 to identify bare soil, which has been proposed to identify the establishment year of rubber plantations [10,20]. Kou and Xiao used the first year when LSWI_{Leafoff_min} < 0 as the establishment year, but this condition may occur during litter fall (defoliation), a period which usually has negative LSWI [23] and may, thus, introduce uncertainty for rubber plantations converted from cropland (Figure 5d,e). Beckschäfer improved Kou and Xiao’s algorithm by detecting the last year in which LSWI_{Leafon_min} < 0 as the plantation establishment year, and reached a RMSE of 2.5 years in Xishuangbanna, Yunnan Province, China. Beckschäfer’s algorithm reduced the uncertainties raised from defoliation around February and for plantations which were converted from cropland.

However, we found that solely using the last year in which LSWI_{Leafon_min} < 0 may underestimate the stand age of rubber plantations by 1–4 years (Figure 11), mainly due to intercropping (e.g., banana and peanut) and weed-eradication management practices (clearing or using herbicide) (Figure 12a,b). The exposed bare soil or dry vegetation by intensive intercropping or weeding would affect LSWI significantly. Beckschäfer did not observe this uncertainty because validation plantations came from visual interpretation via Landsat or Google Earth VHR imagery, while most of our validation data came from historical records from Hainan State Farm or field interviews. In addition, establishment years obtained from the last year in which LSWI_{Leafon_min} < 0 may be affected by severe disasters, such as typhoon and drought, which are prevalent in non-traditional rubber cultivation regions, such as China (Figure 12c,d) [11]. We additionally introduced a biophysical parameter, specifically, the slope of LSWI_{leafon_avg} against the year during the non-production period, to capture the canopy change from seedling to mature plantation and, therefore, reduced the uncertainties from severe natural disasters, intercropping, and weeding. The mean estimated establishment year of rubber plantations was closely distributed along the 1:1 line (Figure 7b) and mean RMSE was reduced to 2.34/0.54 year at the pixel/plantation scale. Compared with Beckschäfer, we improved the algorithm in three ways: (1) reduced underestimation caused by intercropping and weeding; (2) used slope of LSWI_{leafon_avg} against the year as a canopy change factor to reduce the commission error, which had been increased by severe disasters, such as typhoon, drought, and cold injury; and (3) excluded newly-established rubber plantations (here we excluded plantations <4 years). Some studies try to map rubber plantations at very early stages (e.g., ≤2, 2–4; 1–3 year age groups) [1]. However, these rubber plantations are difficult to identify accurately, even with Google Earth VHR imagery. Both spectral and structural features of these early rubber plantations are determined by weeds, shrubs, intercropping crops, or bare soil rather than rubber seedlings themselves (Figure 12a).

4.2. Uncertainties for Mapping Plantation Establishment Year and Pre-Conversion Land Covers

Fragmented landscapes: The fragmented landscape and complex topography in the tropics makes classification challenging [47]. About 25% of Hainan Island is mountainous with elevation greater than 500 m, and the remainder includes lowlands, hills, terraces, and plains which are favorable for agricultural development. The landscape in the low elevation regions is often fragmented into small patches of different land use and land cover types (Figure 12e). Severely fragmented landscapes are not conducive to the precise identification of rubber plantation and affect our ability to determine the establishment year of rubber plantations. We failed to estimate the establishment year of a few rubber plantations (about 1%), mainly due to uncertainties in classification.

Limitation in cloud-free satellite imagery and level of data product: Although we used the entire Landsat archive in our study, cloud-free imagery was severely limited in some years, especially before 2000 when only Landsat TM was available. Cloud-free satellite images during the RDF period (January to March) were important for generating the rubber plantation base map [23] (later used for identifying the establishment year) and tracking pre-conversion land covers. The short RDF period and frequent presence of cirrus cloud cover in tropical regions [27,48] made it difficult to generate an accurate rubber base map and to identify the pre-conversion land cover. In addition, we used Collection 1 TOA reflectance data instead of Collection 1 SR data, as the latter was not yet ready in the GEE platform at the beginning of this study. The TOA reflectance data is affected by atmospheric scattering and absorption [49], and the usage of higher level product of Collection 1 SR data may improve mapping accuracy.

Weeding management and intercropping: Bare soil was exposed during plantation establishment, but vegetation, like grasses and shrubs, recover quickly because of the favorable tropical climate (Figure 12a). Therefore, weeding (manual or herbicidal removal) is a mandatory management practice during rubber seedling establishment. In addition, intercropping (e.g., banana, pineapple, sugarcane, and peanuts) in the open canopy of early rubber plantations (usually less than four years) on relatively flat terrain is an important source of income for farmers during the non-production period (Figure 12b). The LSWI_{Leafon_min}, S_LSWI, and the ability to estimate establishment year of rubber plantations were influenced by management intensity, crops, and the availability of cloud-free Landsat data around the date of human-related disturbances.

Severe natural disasters: Rubber trees always experience several natural disasters during their economic life cycle (~30 years) in China. Typhoon, cold injury, and droughts often hit Hainan Island and Yunnan Province [11]. For example, Hainan Island experienced one of its most severe droughts between September 2004 and July 2005, which was followed by the strongest super typhoon since 1974 on 25 September 2005 (Darry) [50] and a severe cold spell in 2008 [51]. Severe typhoons will harm the canopy (Figure 12c) while drought and disease can delay phenology or cause leaves to shed (Figure 12d). Such disasters will affect LSWI significantly (Figure S5). Therefore, it is better to avoid using satellite images acquired during these periods. Figure S6 shows a serious overestimation of new rubber plantations established in 2004 and 2005 when images acquired during severe drought and post-typhoon were used. However, overestimation was relieved when excluding these images (by only using images on DOY of 90–304 for 2004, and DOY of 120–268 for 2005) (Figure 6d).

4.3. The Spatial-Temporal Dynamic of Rubber Plantations on Hainan Island

The second rapid expansion of rubber plantations on Hainan Island began in the early 1980s, when the old rubber plantations established in the 1950s need to be updated and China started its Reform and Opening policy [52]. That was the main reason for about half of the old rubber plantations on Hainan Island by the end of 2015 (Figure 6b). We found these old rubber plantations were mainly distributed in the central western regions and along the coastal strip from south to northeast (Figure 6b). Unfortunately, we could not accurately obtain their stand age because the earliest Landsat TM imagery on Hainan Island available starts in 1987. During the 1990s, rubber plantations on Hainan Island grew at a very low and stable speed [53] (Figure 8a). During the early 2000s, rubber plantations expanded sharply because of the rapid rise of the natural rubber price, where the average annual price increased from ~$1200/ton in 2002 to the highest of ~$5500/ton in 2011 [54]. Since the southeastern coast is often hit by typhoons [11,55], newly-established rubber plantations were no longer expanding in these areas. Instead, they were concentrated in the northwestern region with good weather conditions and light typhoon disasters (Figure 6b).

In the past 30 years, the pre-conversion land cover for new rubber plantations changed greatly. Before the early 1990s, evergreen forests were the predominant land cover for new rubber plantations. The forest inventory in 1987 showed that the tropical natural forest fell from 863,000 ha in 1954 to 385,000 ha. Vast demand for timber, rubber, and tropical crops were the main factors leading to forest loss [56]. Therefore, the unknown pre-conversion land cover shown in Figure 9d was most likely evergreen forest. With the foundation of Hainan Province in 1988 (previously an administrative district of Guangdong Province) and strengthened governmental protection of tropical forest [56], the conversion of evergreen forests to rubber plantations fell significantly in the 1990s (Figure 9c). The area of rubber plantations converted from evergreen forests increased sharply from 2001 to 2007, which can be explained by the rapid increase in the price of natural rubber (Figure 9c) and development of privately-owned rubber plantations [52]. Stimulated by the natural rubber price, smallholders are more likely to convert evergreen natural forests and secondary forests into rubber plantations. However, forest conversion decreased significantly after 2007 because of conservation policies, such as the ‘forest eco-compensation mechanism’ in 2005 [57].

Old rubber plantations have always been an important land source for new rubber plantations on Hainan Island, even before 2000 (Figure 9c). Large-scale commercial rubber tree plantations began in the 1950s, so some of them were replanted with rubber trees after a 25–30 years life cycle. The sharp increase in the area of rubber plantations converted from old rubber plantations in the 2000s (Figure 9c, increased about 2 × 10⁵ ha, one third of the total area) can be explained by the rapid rise of the natural rubber price and the replacement of old, less-productive rubber plantations which were established after China’s Reform and Opening up in 1978.

The large area of rubber plantations converted from cropland since 2000 can be explained by the rapid development of privately-owned rubber plantations. Statistical data indicated that the area of privately-owned rubber plantations on Hainan Island exceeded the amount of state-owned rubber plantations for the first time in 2009 [54]. Smallholders were more likely to convert their cropland to rubber plantations because of strengthening ecological protection policies and forested lands were unsuitable for conversion to rubber plantations. A direct consequence for this conversion was the increasing complexity and fragmentation of the rubber plantation landscape (Figure 12e).

There was no significant variation in altitude and slope during rubber plantation expansion on Hainan Island in the past 30 years (Figure 9b). The average altitude of rubber plantations was about 150-m AMSL and the mean slope was about 7°. This phenomenon can be explained by the fact that most mountainous regions where rubber plantations are suitable (recommended less than 350-m AMSL) have already been extensively converted to rubber plantations prior to 1990 [56] (Figure 9c). The large area of rubber plantations converted from cropland may offset the trend of having more rubber plantations at higher elevation and steeper slope regions (Figure 9c), since cropland usually have relative flat terrain and are distributed in the low elevation regions. In addition, new rubber plantations converted from old rubber plantations are unlikely to change the pattern of elevation and the slope of rubber plantations.

The expansion pattern of rubber plantations on Hainan Island was quite different from that of Xishuangbanna in Yunnan Province, another important natural rubber production base in China. The area of rubber plantations in Xishuangbanna doubled from 2002 to 2010 [4]. The common practice for rubber plantation expansion in Xishuangbanna was to clear-cut evergreen forests on entire mountainsides and then contour ledges were built for rubber trees (Figure S7). Previous studies indicated that the mean elevation of rubber plantations established in Xishuangbanna before ~2005 was around 800 m, but after ~2005 mean elevation increased to a peak of 980 m in 2013 [10] and the plantations were being established at high elevations of >1200 m AMSL and on steeper terrain [4]. Unlike the complex, fragmented landscape on Hainan Island, it is easier to detect and estimate the establishment year of rubber plantations when mountainsides are clear-cut.

5. Conclusions

Establishment year (or stand age) of rubber plantations is an important parameter for practical management and ecological assessment. In this study, we developed new algorithms to estimate the establishment year of rubber plantations and identify the pre-conversion land cover types with high frequency time-series Landsat data. By applying these algorithms on Hainan Island, we produced a rubber plantation map in 2015 (OA > 97%), and identified the plantation establishment year (R² = 0.85/0.99, and RMSE = 2.34/0.54 years at pixel/plantation scale) and pre-conversion land cover types (evergreen forests, old rubber plantations, and cropland) since 1987 at annual temporal resolution (OA > 88%). The resultant map indicates that rubber plantations on Hainan Island have increased significantly since 2000. The increased plantations were mainly concentrated in the northwest and northern regions and were predominately converted from old rubber plantations and cropland. From a longer time scale, the expansion of rubber plantations on Hainan Island since 1987 predominantly occurred on the fragmented lowlands owned by smallholders, and not towards regions of high elevation and steep slopes. The resultant maps of rubber plantations, establishment year, and pre-conversion land cover are useful for sustainable plantation management and ecological assessment, and the methodology has the potential to be applied to other regions using the GEE platform.