Effect of Textural Features in Remote Sensed Data on Rubber Plantation Extraction at Different Levels of Spatial Resolution

Abstract

The expansion of rubber (Hevea brasiliensis) plantations has been a critical driver for the rapid transformation of tropical forests, especially in Thailand. Rubber plantation mapping provides basic information for surveying resources, updating forest subplot information, logging, and managing the forest. However, due to the diversity of stand structure, complexity of the forest growth environment, and the similarity of spectral characteristics between rubber trees and natural forests, it is difficult to discriminate rubber plantation from natural forest using only spectral information. This study evaluated the validity of textural features for rubber plantation recognition at different spatial resolutions using GaoFen-1 (GF-1), Sentinel-2, and Landsat 8 optical data. C-band Sentinel-1 10 m imagery was first used to map forests (including both rubber plantations and natural forests) and non-forests, then the pixels identified as forests in the Sentinel-1 imagery were compared with GF-1, Sentinel-2, and Landsat 8 images to separate rubber plantations and natural forest using two different approaches: a method based on spectral information characteristics only and a method combining spectral and textural features. In addition, we extracted textural features of different window sizes (3 × 3 to 31 × 31) and analyzed the influence of window size on the separability of rubber plantations and natural forests. Our major findings include: (1) the suitable texture extraction window sizes of GF-1, Sentinel-2, and Landsat 8 are 31 × 31, 11 × 11 to 15 × 15, and 3 × 3 to 7 × 7, respectively; (2) correlation (COR) is a robust textural feature in remote sensing images with different resolutions; and (3) compared with classification by spectral information only, the producer’s accuracy of rubber plantations based on GF-1, Sentinel-2, and Landsat 8 was improved by 8.04%, 9.44%, and 8.74%, respectively, and the user’s accuracy was increased by 4.63%, 4.54%, and 6.75%, respectively, when the textural features were introduced. These results demonstrate that the method combining textural features has great potential in delineating rubber plantations.

1. Introduction

Of the 56.8 million hectares (ha) of plantations in the tropics in 2015, 29.9 million ha are in South and Southeast Asia [1]. These plantations are crucial in global climate change regulation [2] and natural resource protection [1] as they not only increase carbon dioxide fixation and carbon sink with high growth rates, but also provide large quantities of wood and other products and services [3]. Although there are so many benefits, the process of plantation expansion has also caused ecological and environmental problems including soil erosion and natural forest reduction. Natural rubber (latex) is a necessary, strategic material for national development, and the rubber tree is the only source of natural rubber products. Due to its suitable climate and growing conditions, Southeast Asia has the largest proportion of rubber plantations in the world [2]. In recent years, rubber planting areas have dramatically expanded in many Southeast Asian countries due to the increasing global demand for natural rubber products. In order to solve poverty, the Thai government launched a plan to promote rubber plantations in northeastern Thailand from 2003 to 2013. As a result, the planting area of rubber plantations in northeastern Thailand has expanded rapidly in the past decade. While large-scale rubber planting increases the income of local governments and farmers, it also causes many ecological and environmental problems, such as the dramatic reduction of tropical rainforest, soil degradation, and loss of biodiversity in rubber planting areas [4,5]. Timely, accurate maps of rubber plantations are required to manage forests effectively and document the expansion of rubber plantations.

Remote sensing is a powerful and efficient tool to extract rubber plantations from the surrounding landscape. A number of studies have used optical satellite data to detect rubber plantations, mostly relying on spectral signatures like NDVI, EVI, LSWI, and SWIR1 of optical images to delineate rubber plantations [6,7]. However, rubber plantations in Mainland Southeast Asia grow in tropical rainforest areas with complex ecosystems, characterized by high vegetation coverage and less seasonal variation. Rubber plantations have similar spectral characteristics compared to natural forests because the spectral coefficients of rubber plantations tend to be saturated, especially during peak growing seasons or after stands reach a certain age [8,9]. These problems bring considerable uncertainty to rubber plantation extraction and monitoring based on spectral and vegetation indices [9].

Textural features in satellite imagery contain spatial information that is very important for feature extraction because textures do not depend on the tone or brightness of the object surface and can reflect gray information, spatial distribution, and structural information about the image. Textural features are relatively stable compared with spectral characteristics, which are susceptible to the external environment. Some studies have employed textural features for mapping crops, monitoring built-up areas [10], detecting wetlands [11], classifying tree species [12,13], and mapping land cover [14,15]. The effectiveness of textural features on tree species identification has been confirmed by adding them to the classification procedure of high- to relatively low-resolution satellite images. In the study of mapping plantation extent based on Landsat-8 Operational Land Imager (OLI), Synthetic Aperture Radar-2 (PALSAR-2), and Sentinel-1A images, Torbick et al. [9] found textural features could capture the homogeneity of the plantation canopy, spacing, and structure relative to natural forests and were useful for distinguishing oil palm and rubber plantation from other land cover types. Dian et al. [16] confirmed that combining the spectral and spatial information derived from Airborne Hyperspectral Imagery could improve the accuracy of tree species classification. Textures derived from high spatial resolution QuickBird image [17] and the WorldView-3 satellite image [18] were also proved for their potentials in detecting species-specific differences in crown structure, thus for better tree species mapping. Despite many existing efforts, some problems relative to the application of textures to specific tree species like rubber tree mapping require further investigation. Since no single window size could adequately characterize the range of textural conditions in one image [10], what is the effect of window sizes on textures measure and which textural features are more useful and robust for specific tree species discriminating? Furthermore, the textural features strongly depend on spatial resolution of the images used [19], then what is the effect of remote sensing images at different spatial resolutions on specific tree species identification?

In this study, three satellite data including GaoFen-1 (GF-1), Sentinel-2, and Landsat 8 imagery with different resolutions were used for testing and analyzing the effect of textural features on rubber plantation extraction in tropical Thailand. The specific objectives of this study are to (1) define appropriate window sizes for each individual textural feature of different spatial resolution images; (2) quantify the importance of texture variables and identify the effective textural features; and (3) evaluate the applicability of satellite images at different spatial resolutions on the extraction of rubber plantations.

2. Materials and Methods

2.1. Study Area

Northeastern Thailand is located on the Khorat Plateau with an altitude between 100 m and 200 m. It belongs to a humid subtropical monsoon climate, with obvious dry and wet seasons. Heavy rainfall is concentrated in the rainy season from May to October. The highest temperature, often in April, can reach 40 °C while in December, the night temperature usually falls below 0 °C. Traditionally, northeastern Thailand has been an important cultivation area, with fewer rubber plantations [20]. Stirred by an aggressive government policy for rubber plantations since 2003, dramatic expansions have been seen in the south, the traditional rubber plantation area in Thailand. Then the plantations gradually spread to northeastern Thailand to meet the demand for large areas of land [21]. As a result, many patches of cropland and natural forest have been encroached by rubber plantations. A subset image of GF-1, Sentinel-2, and Landsat 8 image within Northeastern Thailand was derived as the experiment area (Figure 1) in this study, where massive rubber planting activity is prompting a significant encroachment to natural forest and other land cover types.

2.2. Datasets and Pre-Processing

2.2.1. Sentinel-1

A total of 92 scenes of Sentinel-1A and Sentinel-1B Interferometric Wide Swath GRD images of 2018 from the Google Earth Engine (GEE) [22] platform were used for forests mapping at the first step. The Sentinel-1 data was pre-processed with the Sentinel-1 Toolbox [23] using thermal noise removal, radiometric calibration, and terrain correction. A Refined Lee filter was applied to de-speckle each image. The final terrain-corrected values were converted to decibels via log scaling (10 × log10(x)), where x is the digital number value of pixels in vertical transmit and vertical receive (VV) or vertical transmit and horizontal receive (VH). The Ratio and Difference between VV and VH were calculated based on Equations (1) and (2).

$$\mathrm{Ratio}={\mathsf{\gamma }}_{\mathrm{VV}}^0/{\mathsf{\gamma }}_{\mathrm{VH}}^0,$$

(1)

$$\mathrm{Difference}={\mathsf{\gamma }}_{\mathrm{VV}}^0-{\mathsf{\gamma }}_{\mathrm{VH}}^0.$$

(2)

where ${\mathsf{\gamma }}_{\mathrm{VV}}^0$ and ${\mathsf{\gamma }}_{\mathrm{VH}}^0$ are the backscattering coefficients of VV and VH in decibels.

2.2.2. GaoFen-1

GaoFen-1 (GF-1) imagery was acquired from the China Center for Resource Satellite Data and Application (http://www.cresda.com). The GF-1 image (path/row = 16/128) was obtained on 5th March 2019. Pre-processing steps for GF-1 imagery using ENVI 5.3 software included radiometric calibration, atmospheric correction, orthorectification, projection definition, and image fusion. The FLAASH model was used for atmospheric correction. The orthorectification was done by Rational Polynomial Coefficient (RPC) files. The spatial resolution of GF-1 varies for the different bands—the panchromatic band has a resolution of 2 m and the multispectral bands are 8 m. We applied the Gram-Schmidt [24] method to fuse panchromatic and multispectral images and obtain multispectral images with a spatial resolution of 2 m.

2.2.3. Sentinel-2

Sentinel-2 data was downloaded from the European Space Agency’s Copernicus Scihub (https://scihub.copernicus.eu). The acquired Sentinel-2 data (T47QQV) was obtained on 16th February 2018. The imagery was pre-processed for atmospheric correction using the Sentinel Application Platform (SNAP) 7.0 open-source software [25]. ESA’s Sen2Cor in SNAP was used to convert the top-of-atmosphere reflectance values (TOA) to corrected surface reflectance values. The spatial resolution of Sentinel-2 data varies from 10 m to 60 m. Bands 1, 9, and 10 were excluded from the dataset due to their sensitivity to aerosol and clouds and their spatial resolution (60 m). Then, the image was resampled at 10 m using a bilinear method.

2.2.4. Landsat 8

The Landsat 8 data was downloaded from USGS (https://earthexplorer.usgs.gov/). The data with orbit number path/row = 129/48 was captured on 13th February 2018. The pre-processing for Landsat 8 using ENVI 5.3 included radiometric calibration, atmospheric correction, and geometric correction.

In order to eliminate the geographic deviation, all acquired data were georeferenced in the WGS84 coordinate system and UTM projection, zone 47. Based on a GF-1 image, relative geometric correction was carried out on Sentinel-2 and Landsat 8 data using a 2nd order polynomial. The characteristics of the datasets used in this study are shown in

Table 1. Summary of the characteristics of the remotely sensed datasets used in the study.

Sensor Name	Sensor Type	Acquisition Date	Band	Resolution (m)
Sentinel-1	C-band Radar	All throughout 2018	VV + VH	10
Sentinel-2	Optical	16th February 2018	490–2190 nm	10–20
GF-1	Optical	5th March 2019	450–890 nm	2–8
Landsat 8	Optical	13th February 2018	450–2290 nm	30

.

2.2.5. Field Survey Data and Sample Datasets

A set of sample points was collected from field survey information and high-resolution satellite imagery from Google Earth (http://earth.google.com/). We conducted the filed surveys in February and August 2017 and collected 284 rubber plantation sample points, 106 natural forest sample points, and 178 other typical land cover types sample points using GARMIN GPSMAP 639sc. Moreover, we obtained 1478 sample points including 679 rubber plantation sample points, 477 natural forest sample points, and 322 other typical land cover types sample points by visual interpretation from Google Earth.

2.3. Texture Processing

The gray-level co-occurrence matrix (GLCM) is currently the most prevalent statistical method to derive textural features [26,27]. We used GLCM to exploit the texture information of GF-1, Sentinel-2, and Landsat 8 images. GLCM is a symmetric matrix with each value representing the probability values of a nearest-neighbor gray tone at a given distance and orientation [26]. This technique first uses a spatial co-occurrence matrix that computes the relationships of pixel values, and then the second-order statistics will be computed by using this matrix. In this study, textural features were retrieved with fifteen window sizes from 3 × 3 to 31 × 31. The following eight textural features including Mean (MEAN), Variance (VAR), Homogeneity (HOM), Contrast (CON), Dissimilarity (DIS), Entropy (ENT), Angular Second Moment (ASM), and Correlation (COR) were computed by the formulas [26] shown in

Table 2. Image texture measure formulas.

Texture Measure	Formula	Description
Mean (MEAN)	$\mathrm{MEAN}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}i{P}_{i,j}$	MEAN represents the average brightness information in the window. It reflects the degree of texture rules. The stronger the rules, the greater the value.
Variance (VAR)	$\mathrm{VAR}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}{\left(P_{i,j}-\mu \right)}^2/\left(N-1\right)$	VAR reflects the contour of each homogeneous region of the image and the change in gray level. When the gray level changes greatly, its value is larger.
Homogeneity (HOM)	$\mathrm{HOM}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}i\frac{P_{i,j}}{1+{\left(i-j\right)}^2}$	HOM is a measure of smoothness of image distribution. The more uniform the image matrix, the larger the value.
Contrast (CON)	$\mathrm{CON}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}i{P}_{i,j}{\left(i-j\right)}^2$	CON reflects texture thickness. The bigger the difference between adjacent pixels and gray value, the bigger the value.
Dissimilarity (DIS)	$\mathrm{DIS}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}i{P}_{i,j}\left|i-j\right|$	DIS is similar to contrast (CON), reflecting the heterogeneity of images.
Entropy (ENT)	$\mathrm{ENT}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}i{P}_{i,j}\left(-\ln {\mathrm{P}}_{i,j}\right)$	ENT is a measure of the amount of information. The more complex the texture in the window, the greater the entropy value.
Angular Second Moment (ASM)	$\mathrm{ASM}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}i{P}_{i,j}{}^2$	ASM is the roughness of image texture. The finer the texture, the smaller the ASM.
Correlation (COR)	$\mathrm{COR}={\displaystyle {\displaystyle \sum }_{i,j=0}^{N-1}}ij{P}_{i,j}-{\mu }_1{\mu }_2$ ${\mu }_1={\displaystyle {\displaystyle \sum }_{i=0}^{N-1}}i{\displaystyle {\displaystyle \sum }_{j=0}^{N-1}}{P}_{i,j}$ ${\mu }_2={\displaystyle {\displaystyle \sum }_{j=0}^{N-1}}j{\displaystyle {\displaystyle \sum }_{i=0}^{N-1}}{P}_{i,j}$	COR reflects the similarity of pixels in row and column directions of the gray level co-occurrence matrix. The higher the correlation, the greater the value.

Where i, j are the gray levels of paired pixels; P(i,j) are the probabilities of co-occurrence; and N is the number of distinct gray levels in the quantized image.

.

There is a certain degree of correlation between multispectral bands, resulting in information overlap between different bands [11], so we used the first principal component of the images in the principal component analysis (PCA) [28] to extract the textural features. Bands 1–4 were input into the PCA for GF-1 image, Bands 2–8A, 11 and 12 for Sentinel-2 image, and Bands 2–7 for Landsat 8. The variance contribution of the first principal component of GF-1, Sentinel-2A, and Landsat 8 are 92.41%, 69.53%, and 74.34%, respectively.

2.4. Random Forest Classification

Random Forest (RF), proposed by Leo Breiman [29] in 2001, is an ensemble learning algorithm for supervised classification with a decision tree as the basic classifier. Previous studies have shown RF algorithms can process high-dimensional, massive data in parallel with high accuracy compared with other machine learning algorithms. In this study, Random Forest was employed as a tool to classify and assess the effect of textural features in remote sensed data at different spatial resolutions on rubber plantation extraction (Figure 2).

2.4.1. Forests Mapping

Previous studies have proved that mapping forests using SAR data was feasible and accurate [7,30,31]. In this study, Sentinel-1 SAR data were employed to produce the forest map as the baseline for further assessing the effect of textures on classifying rubber plantations from natural forests. The annual average value images of VV, VH, Ratio, and Difference were input to RF model for classifying forests from other land cover types in reference to [31,32]. In this study, a Python script was used to import the GDAL library to read training sample data and the Random Forest classifier was imported from the Scikit-learn library. Each Random Forest was built using ntree = 800 trees, and the number of variables randomly sampled as candidates at each split (mtry) was set by default (sqrt(p), where p is the number of variables). Stratified random sampling of all the 2046 samples was carried out in ArcGIS 10.5, 70% of which were selected to train RF model and the remaining 30% were used for accuracy assessments of forests mapping using confusion matrix. The 10 m forest/non-forest map was obtained by merging non-forest categories and resampled to 2 m and 30 m to match the GF-1 and Landsat 8 spatial resolutions.

2.4.2. Feature Importance Measure

The RF algorithm is more advantageous for estimating the importance of predictor variables compared with other machine learning algorithms [33,34]. The feature importance is calculated using the Mean Decrease in Gini (MDG), which measures how much a feature reduces the Gini Impurity metric in a class. The Out of Bag (OOB) error is then used for model evaluation to determine the importance of model features. We calculated different window sizes for each texture with 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11, 13 × 13, 15 × 15, 17 × 17, 19 × 19, 21 × 21, 23 × 23, 25 × 25, 27 × 27, 29 × 29, and 31 × 31 to assess the effect of window sizes on textures for rubber plantation extraction. Then the 15 variables were input to RF to determine the optimal window size for each texture using out-of-bag (OOB) estimation. After that, eight textural features which were calculated using individual optimal window size were input to RF to measure their importance for rubber plantation identification.

2.4.3. Rubber Plantation Extraction Using GF-1, Sentinel-2, and Landsat 8 Data

To evaluate the effect of different data sources on rubber plantation extraction, RF models were established using only the original spectral bands and spectral bands combined with textural features selected above to classify the images masked by the forest map. A total of 1528 rubber plantation and natural forest samples were used for training and validation using confusion matrices.

3. Results

3.1. Forest Mapping and Accuracy Assessment

According to the validation results of the confusion matrix (

Table 3. Accuracy validation of forest/non-forest classification result based on 10 m Sentinel-1.

Class	User’s Accuracy	Producer’s Accuracy	Estimation Error 1	Overall Accuracy	Kappa
Forest 2	98.24%	96.34%	3.66 ± 1.71%	0.9593	0.8919
Non-forest 3	89.31%	94.67%	5.33 ± 3.59%

Note: ¹ the error associated to the estimation of a proportion p is z*sqrt(p(1 − p)/n), where n is the number of data points utilized to compute p, and z the tabulated z-score [35]. The number of forest points is 464 and that of non-forest is 150. Z = 1.96 ($\mathsf{\alpha }$ = 0.05) in this study. ² Forest types mainly include natural forest and rubber plantation. ³ Non-forest mainly includes cropland, water, and built-up land.

), the resultant Sentinel-1 forest/non-forest map (Figure 3) has a high accuracy. The overall accuracy is 0.9593, and the Kappa coefficient is 0.8919. The user’s accuracy and producer’s accuracy of the ‘forest’ category are 98.24% and 96.34%, respectively. Therefore, the forest map can be used as a reliable base map for further analysis.

3.2. The Optimal Window Sizes for Textures Calculation

Figure 4 shows how the importance of textural features, derived from GF-1, Sentinel-2, and Landsat 8, varies with the window size in distinguishing rubber plantations from natural forest. With the increase in the window size, the importance of the textural features exhibited different variation modes for the three data sources.

The importance of GF-1 textural features increases with small fluctuations as the window size enlarges, and the increasing trends of the MEAN and CON are especially obvious. The importance of the eight textural features of Sentinel-2 follows up-down-up modes with the increase in window size and reaches the maximum mainly at 11 × 11, 13 × 13, and 15 × 15. Among the eight textural features, the increasing and decreasing trend of the importance of VAR and COR is most obvious before and after the peak. In the case of Landsat 8, high importance can be observed in smaller window sizes (3 × 3 to 7 × 7). Except for COR, the rest of the textural features are of relatively high importance at window sizes ≥25 × 25. Considering the 30 m resolution of Landsat 8, too large of a window size will over-smooth the image, introducing erroneous spatial information, and thus is not suitable for calculating Landsat 8 textural features. To define the appropriate window size is a very important step when textural features are considered for classification because the window sizes selected have had a great effect on feature contribution.

The optimal window size was determined where the peak importance was located for rubber extraction. The resultant window size for extracting eight textural features from GF-1 is 31 × 31. For Sentinel-2, the window sizes are 15 × 15 for MEAN and VAR, 13 × 13 for HOM, CON, ENT, ASM, and COR, and 11 × 11 for DIS. For Landsat 8, the textural features calculation window sizes of MEAN, ENT, and ASM are 3 × 3; CON and DIS are 5 × 5; and HOM, VAR and COR are 7 × 7.

3.3. The Importance Sorting of Textural Features for GF-1, Sentinel-2, and Landsat 8

We calculated the importance of eight textural features derived from the optimal window size defined in Section 3.2. The results are shown in Figure 5. Among the textural features of GF-1, ENT is the most important, accounting for 32.02%, followed by MEAN, COR, VAR, accounting for 25.48%, 19.30%, and 10.96%, respectively. The textural features of Sentinel-2 with relatively high importance are HOM, ENT, COR, and ASM, with values of 20.22%, 19.52%, 17.84%, and 13.78%, respectively. For Landsat 8, COR is the most important textural feature, accounting for 17.87%, followed by MEAN, CON, and HOM, accounting for 17.49%, 15.82%, and 13.13%, respectively. Among the top 3 important textural features, COR is the only robust feature in rubber plantations delineation concurrently for all three data sources.

The cumulative contribution rate of the top 3 textural features is 76.80% for GF-1, 57.58% for Sentinel-2, but only 51.18% for Landsat 8. Similarly, the cumulative contribution rate of textural features ranked top 4 is 87.76% for GF-1, 71.36% for Sentinel-2, but only 64.31% for Landsat 8. So, we can find the contribution is more concentrated on a small number of contextual features for GF-1, and less concentrated for Sentinel-2. For Landsat 8, the contribution disperses among all textural features, with no one more than 18%.

According to the ranking results, MEAN, VAR, ENT, and COR were selected as important textural features to combine with spectral information for Random Forest classification of GF-1. Similarly, HOM, ENT, ASM, and COR of Sentinel-2 imagery and MEAN, HOM, CON, and COR of Landsat 8 were chosen to build the Random Forest model.

3.4. Classification Results Using GF-1, Sentinel-2, and Landsat 8 Data

The resultant rubber plantation and natural forest maps produced from GF-1, Sentinel-2, and Landsat 8 imagery are shown in Figure 6. The classification accuracy was validated using a confusion matrix for each data source (

Table 4. Accuracy validation of RF classification using only spectral information.

	Class	User’s Accuracy	Producer’s Accuracy	Estimation Error	Overall Accuracy	Kappa
GF-1	Rubber plantation	92.65%	79.37%	20.63 ± 6.27%	0.8271	0.6448
	Natural forest	70.65%	88.75%	11.25 ± 3.66%
Sentinel-2	Rubber plantation	92.13%	81.82%	18.18 ± 5.98%	0.8386	0.6639
	Natural forest	72.92%	87.50%	12.50 ± 3.83%
Landsat 8	Rubber plantation	89.84%	80.42%	19.58 ± 6.15%	0.8161	0.6162
	Natural forest	70.53%	83.75%	16.25 ± 4.28%

Note: The number of rubber plantation points is 286 and that of natural forest is 160. z = 1.96 (α = 0.05) in this study.

and

Table 5. Accuracy validation of RF classification using spectral and texture information.

	Class	User’s Accuracy	Producer’s Accuracy	Estimation Error	Overall Accuracy	Kappa
GF-1	Rubber plantation	97.28%	87.41%	12.59 ± 5.14%	0.9036	0.7985
	Natural forest	80.95%	95.63%	4.37 ± 2.37%
Sentinel-2	Rubber plantation	96.67%	91.26%	8.74 ± 4.38%	0.9238	0.8379
	Natural forest	85.80%	94.38%	5.62 ± 2.67%
Landsat 8	Rubber plantation	96.59%	89.16%	10.84 ± 4.82%	0.9103	0.8108
	Natural forest	82.97%	94.38%	5.62 ± 2.67%

). As expected, rubber plantations were confused with natural forest in the classification that only used spectral bands. The producer’s accuracy of rubber delineation results for Sentinel-2 and Landsat 8 data was less than 82%, and for GF-1 the rubber plantations producer’s accuracy was the lowest with an accuracy of only 79.37%, which implied that the resultant rubber plantation map produced by GF-1 had the largest number of rubber plantations that was misclassified as natural forests. Regarding the user’s accuracy, GF-1 had the highest accuracy (92.65%), followed by Sentinel-2 (92.13%) and Landsat 8 (89.84%), indicating that it had the largest incorrect classification of the rubber plantation objects in Figure 6e. The overall accuracy of rubber plantations and natural forests classifications using GF-1, Sentinel-2, and Landsat 8 were 0.8271, 0.8386, and 0.8161, respectively, and the Kappa coefficients were 0.6448, 0.6639, and 0.6162, respectively.

The accuracy of combining spectral and textural features is higher than that of only using spectral features (Table 5). The improvement in accuracy varies under different resolutions. The overall accuracy of Sentinel-2 is the highest at 0.9238, followed by Landsat 8 at 0.9103 and GF-1 at 0.9036. The accuracy changes before and after adding textural features are shown in Figure 7. The producer’s accuracy of rubber plantations based on GF-1, Sentinel-2, and Landsat 8 was improved by 8.04%, 9.44%, and 8.74%, respectively, and the estimation error was reduced correspondingly. The user’s accuracy was increased by 4.63%, 4.54%, and 6.75%, respectively. The overall accuracy and Kappa coefficient have been increased to varying degrees. In addition, the “salt and pepper” phenomenon shows an obvious improvement, and the uniformity effect of rubber forests is better after adding textural features. The GF-1, Sentinel-2, and Landsat 8 classifications after adding textural features could accurately extract the rubber plantation areas, with an area of 179.44 km², 204.83 km², and 221.01 km², respectively.

4. Discussion

4.1. The Optimal Window Sizes for Textural Features Calculation

The importance of textural features for land cover classification has been studied by several researchers [27], and window size plays a critical role in texture calculation [14]. The standard method of texture analysis is to process spectral information in a single band with a fixed window size. However, this method has limitations because a single window size is not able to adequately describe multi-scale characteristics of different objects [10], nor can it fully characterize all the textures for the image. In previous studies [36,37], all textural features were defined by only one fixed window size; we further defined the optimal window size for each individual textural feature.

Window size selection depends on the spatial resolution of the image and the relationship between features (such as canopy size) [16] and texture indices to be calculated [38]. In fact, the spatial characteristics of specific land cover types cannot be exploited sufficiently if the window is too small. Too large of a window may cause mixed pixels at the boundary of two categories to be misclassified, and even pure pixels may be misclassified into another category due to the edge effect [39,40]. Therefore, selecting the texture window is a very important step in calculating textural features. In this study, we used data in three different spatial resolutions to study the multi-scale effect. With the decrease in spatial resolution, the edge effect became more and more obvious for the “rubber plantation” category. The edge effect was apparent for rubber plantations in Sentinel-2 and Landsat 8 images, especially in small rubber plantation patches, since there were more non-rubber plantation pixels involved in the texture measuring area. The lesser edge effect in GF-1 may be due to the higher spatial resolution and higher spectral variation of the pixels at the edges, which leads to insufficient spatial characteristics of rubber plantations, even though the window size had already been set to 31 × 31.

As window size for texture analysis is related to image resolution, it is meaningful to choose the texture calculation window size for images with different resolutions. At finer spatial resolutions, the suitable window size used for extracting textural features should be larger when extracting rubber plantations. In this study, the appropriate texture window sizes for GF-1 images is 31 × 31, then 11 × 11 to 15 × 15 for Sentinel-2, and 3 × 3 to 7 × 7 for Landsat 8 as the textural features of rubber plantations and natural forests could achieve the maximum differentiation among such windows. The results for window size are partially consistent with previous, similar studies. Aguera et al. [17] used QuickBird imagery (2.5 × 2.5 m) to detect plastic-covered greenhouses. Their results revealed that the optimal window size was around 15 × 15. Zhou et al. [38] used 37 window sizes from 3 × 3 to 75 × 75 to test the optimal window size of Sentinel-1 imagery for urban land cover classification. They revealed that for each individual textural feature (MEAN, VAR, ENT, DIS, HOM, COR, CON, and ASM), the optimal window sizes for the best classification results were: 5 × 5, 23 × 23, 25 × 25, 13 × 13, 19 × 19, 49 × 49, 51 × 51, and 9 × 9, respectively. Gao et al. [41] suggested that a smaller window size (5 × 5) may be more representative to obtain multi-seasonal textures of Landsat 8 to map the spatial distribution of the six forest types. Compared to rubber plantations, plastic greenhouses have a smoother texture and larger contrast with surrounding objects, so a small window size can adequately describe the difference. When selecting the window size of a texture, the special conditions of the study area, the bands used to calculate textures, the contents within the images, and the main objective of the classification and other factors should be taken into consideration [16,27].

4.2. Contribution of Textural Features on Rubber Plantation Delineation

This study showed that texture information is a valuable feature for rubber plantation delineation as the addition of textural features makes it easier to accurately separate rubber plantations from natural forests. After adding textural features, overall accuracies of the resultant rubber plantation map were increased by 0.0765 to 0.0942 and the Kappa coefficients were increased by 0.1537 to 0.1946. Several studies have found that individual textural feature contributed differently to specific land cover classification [9]. Pei et al. [42] chose CON as the best textural feature to participate in classification. Berberoğlu et al. [27] also found that CON was a useful discriminator for any individual class. Mongus et al. [43] found that texture measures including MEAN, CON, HOM, and DIS extracted from Sentinel-2 data had increased the classification method’s total F1 score by over 7%. Moreover, CON introduced the largest improvements (approximately 3%). Torbick et al. [9] found that the Landsat textural features based on NIR had substantial separation between plantations and forests as forests have higher values of VAR, CON, DIS, while plantations have higher ENT values. The four highest information entropy values of the textural features, including MEAN, CON, HOM, and VAR, were selected to classify the image using the information entropy theory in the study by Zhang et al. [44]. After adding textural features, the overall accuracy was improved by 10.99–24.10% compared with spectral information alone. Similar to these studies, we evaluated the importance of eight textural features when distinguishing between rubber plantations and natural forests. Our results also showed that MEAN, HOM, CON, VAR, and ENT have important effects on improving the separation of rubber plantations and natural forests.

Sorting the results of textural feature importance (Figure 5) showed that COR was the most robust feature when distinguishing rubber plantations from natural forests as its importance has always been in the top three in images of all three spatial resolutions, although not always of the highest importance. COR may provide additional insight into the characteristics of rubber plantations.

4.3. Applicability of Remote Sensing Data at Different Resolutions to Rubber Plantation Extraction

Many studies have shown that the specific land cover type area is closely related to the sensor data used. We also found that there were significant differences in the rubber plantation areas produced by GF-1, Sentinel-2, and Landsat 8. With the decrease in spatial resolution, the rubber plantation area that was extracted increased gradually. The outcome of the study by Paul et al. [45] also showed that glacier extents derived from the 30 m Landsat 8 imagery are 4–5% larger than those from 10 m Sentinel-2. The results are related to the generous interpretation of mixed pixels, especially along the target category perimeter. In this study, the rubber plantation area estimated with GF-1 was less than the actual area as the ability of GF-1 to recognize edges for rubber plantations is still insufficient (Figure 8), although GF-1 rubber plantation extraction results were rich in detail.

According to the classification accuracies shown in Table 5, after adding textural features, the moderate spatial resolution Sentinel-2 and Landsat 8 data outperformed GF-1 with high spatial resolution in extracting rubber plantations. This may be because of the absence of a SWIR band, which was also demonstrated by Shapiro et al. [46] and Wang et al. [47] when mapping mangrove extent. Wang et al. [47] found that the moderate Sentinel-2 and Landsat 8 imagery could outperform high-resolution Pléiades-1 imagery with the help of more refined spectral divisions when using different sensors to map larger mangrove zonation and extent. In this study, the ability of the GF-1 sensor to discriminate rubber plantations with low canopy density was lower than that of Sentinel-2 and Landsat 8 (Figure 8a). This may be because GF-1 imagery has only four multispectral bands, lacking a SWIR band compared with Landsat 8 and Sentinel-2. In previous studies [9,48], LSWI and NBR indices (the function of SWIR) played important roles in distinguishing rubber plantations from natural forests. Consequently, the absence of a SWIR band may limit the accuracy of GF-1 imagery in identifying rubber plantations. In addition, the rubber plantations with low canopy density were more affected by the vegetation and soil under the plantations, which will enlarge the variability of its interior and increase the difficulty of identifying other rubber plantations.

When using different sensors to map small-sized features and the details of rubber plantations, spatial resolution might be more important than the spectral bands. All three rubber plantation maps from GF-1, Sentinel-2, and Landsat 8 show generally consistent spatial distribution but differ in geometric details such as boundaries and gaps between two rubber plantation patches. The higher the spatial resolution, the richer the detailed information, especially in scattered, small plantation plots (Figure 8b). The Landsat 8 sensor resulted in under classification of line shape (Figure 8b). This may be because the spatial extent of rubber plantation is less than the pixel size, which leads to the rubber plantations being dissolved by surrounding objects. Plantations along the gaps and boundaries were underestimated, even though the map produced by the high spatial resolution GF-1 data could better distinguished the gaps between two rubber plantation patches and boundaries between rubber plantations and surrounding objects better than Sentinel-2 and Landsat 8. This may be due to the high spatial resolution of GF-1, which resulted in large spectral variation at the edges of the rubber plantation, and the rubber tree shadow at the edges (Figure 8c) and boundaries could be misclassified as natural forest (Figure 8d).

When we decide which data source to use, a variety of factors needs to be taken into account, such as spatial resolution, spectral resolution [47], and the patch size of rubber plantation. The 2 m GF-1 imagery presents more abundant details of the rubber plantations with its smaller pixel size, which can provide excellent imagery for visual interpretation. However, fine spatial resolution will bring more variability than moderate resolutions. In addition, texture calculation is based on the spectral relationship among pixels within a certain window range; the low spectral resolution of GF-1 may limit the distinction between rubber plantations and natural forests in texture characteristics. Owing to those reasons, GF-1 has limited advantages in automatic extraction of rubber plantations. Regarding the 10 m Sentinel-2 imagery with the most refined spectral division, the highest rubber plantation accuracy and accurate boundary of rubber plantation were achieved with a 10 m Sentinel-2 image. Furthermore, as frequent cloud cover is a major obstacle for creating precise images in tropical regions, the shorter repetition interval of Sentinel-2 might allow an increased availability of appropriate images. Compared with GF-1, Landsat 8 has more spectral bands. However, a 30 × 30 m Landsat 8 pixel tends to contain several types of land cover with a fragmented context, which limits its suitability for mapping smaller rubber plantations, as Landsat 8 over- or under-classified rubber plantations to some extent. On the whole, for continuous or large areas of rubber plantations, Sentinel-2 and Landsat 8 are good choices; for scattered rubber plantations and small plot sizes, GF-1 may be better.

5. Conclusions

Land use in northeastern Thailand has undergone tremendous changes in the past two decades because the Thai government launched a plan to promote rubber plantations from 2003 to 2013. Due to the limitation of similar spectral characteristics between rubber plantations and natural forests, textural features were introduced in this study to improve the accuracy of rubber plantations recognition. We explored the performance of textural features for three different image resolutions images including 2 m GF-1, 10 m Sentinel-2, and 30 m Landsat 8. Several conclusions can be drawn from this study:

(1)

The higher the resolution, the larger the window size needed. The resultant window size for extracting eight textural features from GF-1 images is 31 × 31. For Sentinel-2, the window sizes are 15 × 15 for MEAN and VAR, 13 × 13 for HOM, CON, ENT, ASM, and COR, and 11 × 11 for DIS. For Landsat 8, the textural feature calculation window sizes of MEAN, ENT, and ASM are 3 × 3, 5 × 5 for CON and DIS, and 7 × 7 for HOM, VAR and COR.

(2)

The Random Forest importance measures show that MEAN, VAR, ENT, and COR were important for GF-1 classification, while HOM, ENT, ASM, and COR were chosen for Sentinel-2 classification. MEAN, HOM, CON, and COR had more influence in the separation of rubber plantations and natural forests for Landsat 8.

(3)

The importance of COR was always in the top three among the eight textural features in GF-1, Sentinel-2, and Landsat 8, so COR is a robust textural feature in three different resolutions when distinguishing rubber plantations from natural forests.

(4)

Adding textural features as additional inputs improved the overall accuracy and the producer’s and user’s accuracies compared to using spectral features only for rubber plantation classification. With the help of textural features, the 10 m Sentinel-2 imagery could accurately extract rubber plantations with the producer’s accuracy reaching 91.26% and the user’s accuracy reaching 96.67%. Meanwhile, the 2 m GF-1 imagery underestimated rubber plantations with a producer’s accuracy of 87.41% and a user’s accuracy of 97.28%, and the rubber plantation area was overestimated by the 30 m Landsat 8 imagery with a user’s accuracy of 96.59% and a producer’s accuracy of 89.16%. Compared with GF-1, the higher accuracies of the Sentinel-2 and Landsat 8 sensors may be attributed to the SWIR band.

Textural features can capture the homogeneity of plantation canopy, spacing, and structure. Textural features are more stable than spectral characteristics as textural features reflect the varying frequency of tone in a certain window size, rather than depending on the tone and brightness of the object surface. Due to such advantages, the textural-based approach can be duplicated or transferred to other geographical regions and other plantation species such as oil palm (Elaeis guineensis) and teak (Tectona grandis). We finally contend that these achievements are of great significance to guide the sustainable management of forest resources, soil erosion control, and sustainable development of society in Southeast Asia and other regions.