Mapping large area tea plantations using progressive random forest and Google Earth Engine

. A timely and accurate understanding of the spatial distribution of tea plantations is beneficial for agricultural management and regional sustainable development. However, obtaining detailed distribution data on large-area tea plantations remains challenging owing to limitations in computational capabilities, training data, and workflow design. Utilizing the Google Earth Engine, which provides a catalog of multisource data in a cloud-based environ-ment, we developed a methodology to generate a highly accurate tea plantation map, with a 10-m resolution, for Anhui Province, China, by integrating a random forest model with a progressive model. Our major contribution lies in this hybrid approach, which comprises two major com-ponents: (1) an optimal classification band combination derived from Sentinel-2 products and the digital elevation model filtered by the J-M distance model and (2) a progressive random forest method introduced for tea plantation classification. The experimental results show that our proposed workflow achieved an average classification accuracy of 89.27% for the entire Anhui Province. In addition, this approach is semiautomatic and can effectively reduce the labor required during the generation of training data compared with traditional classification approaches. These findings demonstrate the potential of integrating machine learning and progressive models to produce high-precision remote sensing classification maps. © The Authors. Published SPIE Creative Commons License.


Introduction
Tea, which originated in China, has become one of the three most popular beverages worldwide. 1 In China, the global leader in tea production, tea is primarily grown in tropical and subtropical regions, such as Anhui Province, 2 where tea cultivation plays a vital role in the agricultural economy and rural development. 3 According to the Anhui Provincial Bureau of Statistics, the area of tea plantations in Anhui in 2019 was 1870.58 km 2 , an increase of 57.38% from 1990. 4 Although the expansion of tea plantations has promoted local economic development, it has also caused a series of environmental problems, such as decreased soil fertility and soil erosion. 5 Therefore, a timely and accurate understanding of the spatial distribution of tea plantations is conducive to governmental environmental protection and agricultural management. 6 Information on the distribution of tea plantations can generally be obtained through manual reporting and remote sensing image classification. 7 Utilizing remote sensing data to identify tea plantations saves time and labor costs; hence, it is more popular than manual reporting. 8 The biggest difficulty in using remote sensing data is that tea trees are perennial evergreen woody plants, and their spectral characteristics are easily confused with similar woody plants, such as natural forests. 9 Thus, it is difficult to achieve accurate tea plantation classifications using only spectral bands. 1 Hence, previous studies primarily used hyperspectral and high-resolution images, fusions of multiple types of remote sensing data, and classifications based on phenological features to identify tea plantations. 1,2,[9][10][11] The classification accuracies using the former three methods are very high, but often the data are expensive and the study area is small, which makes it difficult to identify tea plantations over large areas. 12 In this situation, phenological feature-based classification using Sentinel-2 data is among best choices for identifying the distribution of tea plantations over large areas. 13 This is because Sentinel-2 data have three distinct advantages over other satellite images: worldwide coverage, high temporal and spatial resolution, and a variety of red-edge bands that are sensitive to vegetation. 14 Nevertheless, there are some challenges in applying this method to identifying tea plantations over a large area. 15 First, determining the optimal combination of bands for classification is difficult. 12 Too few classification bands may lead to low classification accuracy, whereas too many bands may lead to low classification efficiency, over-fitting, and local optima. 16 Second, the classification of tea plantations over a large area requires many accurate training samples. 17 If the number of training samples is small, the classification accuracy is usually low and has a large uncertainty. 18 However, obtaining multiple training samples is costly and time-consuming. 19 Therefore, it is essential to balance the number of training samples with the required labor. Finally, for large-area time-series classifications, a high storage capacity is necessary to store the remote sensing data. 20 Furthermore, strong computing power is required to achieve high classification efficiency. 21 Web-based remote-sensing cloud platforms such as Google Earth Engine (GEE) can greatly reduce computing time and enable the classification of large areas using remote sensing time series data. 22 GEE contains a variety of remote sensing images from different sensors, which reduces personal data storage requirements. 23 In addition, GEE provides high computing power that greatly improves classification efficiency, and it can effectively perform repeated comparison experiments to determine the optimal classification bands. 24 Moreover, GEE provides a visual interactive platform that can dynamically add training samples, providing a balance between training sample acquisition and labor. 25 In this study, we used the Jeffries-Matusita distance (J-M distance) model to select the optimal classification features and utilized the progressive random forest method to classify tea plantations. These processes were conducted on the GEE platform, allowing our workflow to achieve high classification accuracy for tea plantations over a large study area with minimal time and labor costs.
The goals of this study were to classify Sentinel-2 time series data to obtain a high-precision tea plantation distribution map for Anhui Province during 2020 and determine which classification features and methods effectively perform high-precision tea plantation classification over a large area with complex terrain. Specifically, we addressed the following two questions: (1) how can we select the optimal combination of classification features in GEE, so it can classify tea plantations with high classification accuracy and efficiency? (2) How can a suitable classification method be constructed in GEE to ensure that tea plantations can be identified with high accuracy over a large area?
2 Study Area and Datasets

Study Area
Anhui Province is in eastern China, between 114°54E and 119°37E and 29°41N and 34°38N. It spans 450 km from east to west and 570 km from north to south, with a total area of 140; 100 km 2 , accounting for 1.45% of the total land area of China. 26,27 Anhui is in the middle to lower reaches of the Yangtze River and the Huaihe River. It borders Jiangsu to the east, Zhejiang and Jiangxi to the south, Henan and Hubei to the west, and Shandong to the north. 28 As of July 2020, there were 16 cities in the administrative regions in Anhui Province: Hefei, Wuhu, Bengbu, Huainan, Ma'anshan, Huaibei, Tongling, Anqing, Huangshan, Fuyang, Suzhou, Chuzhou, Lu'an, Xuancheng, Chizhou, and Bozhou (Fig. 1).
In general, the natural conditions of northern and southern Anhui differ greatly. 29 The terrain is generally high in the west and south and low in the east and north. The north is dominated by plains, whereas the south is dominated by hills and mountains. Anhui is in the mid-latitude zone, which belongs to the transitional region between the warm temperate and subtropical zones. The annual average temperature in Anhui is ∼14°C to 17°C, and the annual average precipitation is 750 to 1750 mm. 26 Precipitation is higher in the south and the mountainous areas and is lesser in the north, the plains, and the hills. Rainfall is abundant in summer, accounting for 40% to 60% of the annual precipitation. 27 The unique climatic and topographical conditions have made Anhui Province one of the major tea-producing provinces in China. 6 In 2019, the tea plantation area reached 187,058 hectares, and the total output of tea was 121,980 tons. 4

Basic Data and Data Preprocessing
The data used in this study include four types of raster data, specifically, Sentinel-2 MSI data, land use and land cover data, topographic data, and Google Earth CNES/Airbus imageries, and two types of vector data, specifically, administrative boundary data and field survey data. All data were open-access datasets, except for the field survey data. See Table 1 for detailed descriptions of the data used in this study.
The Sentinel-2 data were obtained by a satellite cluster composed of two identical satellites, Sentinel-2 A and Sentinel-2 B. 30 The revisit periods of the two satellites were 5 days and were simultaneous. 31 In this study, Sentinel-2 satellite data were selected as the main remote sensing image data for accurately identifying the distribution of tea plantations. Sentinel-2 satellite data products include three levels: level-0, level-1C, and level-2A. Among them, level-2A images are bottom-of-atmosphere reflectance in cartographic geometry after atmospheric correction and, thus, are ready to use. 32 In this study, we used level-2A products; only the identification of clouds and cloud shadows was required prior to use. We used the Fmask 4.0 algorithm by Qiu et al. for cloud detection and achieved good results. 33 Topographic features (elevation, slope, and aspect) are important factors that affect tea tree growth; therefore, they can effectively improve the classification accuracy for tea plantations. 11 To describe the topographic features, we used a digital elevation model (DEM) with a 30-m resolution generated by the Shuttle Radar Topography Mission (SRTM). 34 This DEM is a postprocessed elevation dataset that is widely used because of its high accuracy and extensive coverage. 35 High-precision reference data and an appropriate classification system are prerequisites for classifications. In this study, we combined imagery observations and field investigations to determine the reference points. We collected field survey data acquired with a handheld GPS/GNSS receiver in point-positioning mode in Huangshan and Wuhu City in 2020. The imagery observation data for other cities were incrementally added via Google Earth CNES/Airbus during the progressive random forest classification process. According to Giuseppe Pulighe et al., an error of ∼1 m occurred in visualizing added points. 36 Considering that our classification results have an accuracy of 10 m, these points can be used for classification and accuracy verification.
GlobeLand30 was used for qualitative cross-validation with Sentinel-2 MSI classification results through visual interpretation. Then, according to the distribution of the ground objects in the study area, they were reclassified into six LULC (land use/land cover) types: cropland, forest, grassland, water bodies, built-up land, and unused. 37 After the data were collected, the inputs were resampled to a spatial resolution of 10 m, which corresponds to the fine resolution of the Sentinel-2 MSI. Finally, all data were projected into the WGS84 coordinate system (EPSG: 4326) for tea plantation classification. Figure 2 shows the overall workflow for the classification of tea plantations in the study area in 2020. First, to enhance the separability of the tea plantations from the other LULC types, we used original band features, remote sensing index features, texture features, tasseled cap transformation (TCT) features, and terrain features to generate overall classification features. Then, based on maintaining the separability between the tea plantations and other LULC types, the J-M distance model was applied to obtain the optimal classification features. Finally, we used a progressive random forest classification method for high-precision classification.

Spectral features of the original bands
The spectral features of the original bands of the remote sensing images record the reflection intensity of different ground objects. 38 By analyzing the different spectral characteristics of these objects, it is possible to distinguish vegetation, water bodies, and other types of ground objects. 39 The 60-m bands in the Sentinel-2 MSI data are mainly used for atmospheric correction; therefore, they were excluded from this study. The other 10 original bands in the Sentinel-2 MSI data were selected for classification, specifically, blue, green, red, red edge1-4, NIR, and SWIR1-2. To analyze the phenological characteristics of the vegetation, we calculated the maximum,

Features of remote sensing indices
Due to the surface environment, the same objects may have different spectral features, and different objects may have the same spectral features, which directly impacts the classification accuracy. 40 To overcome this issue, this study selected five commonly used remote sensing indices to improve classification accuracy. These indices were mainly used to enhance the distinction among other LULC types, rather than tea plantations. Nevertheless, because of the strong spectral similarity between tea trees and forests, eight red-edge vegetation indices were used to further distinguish tea plantations from forests. The five common remote sensing indices include the normalized difference vegetation index (NDVI), soil adjusted vegetation index (SAVI), normalized difference water index (NDWI), modified normalized difference water index (MNDWI), and normalized difference built-up index (NDBI). [41][42][43][44][45] The eight red-edge vegetation indices include the normalized difference vegetation index Red-edge1-3 (NDVI re1 , NDVI re2 , and NDVI re3 ), normalized difference Red-edge1,2 (ND re1 and ND re2 ), inverted red-edge chlorophyll index (IRECI), MERIS terrestrial chlorophyll index (MTCI), and red-edge chlorophyll index (CI re ). 42,46 The equations to calculate these indices are listed in Table 2.

Texture features
Different ground objects have different texture features; therefore, adding texture features helps enable the full capabilities of remote sensing image information. 47,48 To facilitate tea harvests and obtain adequate nutrition (e.g., CO 2 and soil nutrients), tea trees are planted at certain spatial intervals. 1 Due to this unique cultivation method, tea plantations have obvious textural features that differ significantly from natural vegetation such as forests and grass. 9 Currently, the most common method of texture analysis is the gray-level co-occurrence matrix (GLCM). 49 Because the GLCM has demonstrated good adaptability and strong robustness, it was selected to extract texture information. 50 The GLCM method requires an appropriate band when calculating texture features. 49 Because the texture features calculated based on Sentinel-2 MSI images are the most effective for distinguishing tea plantations from forests, the selected band should be highly capable of detecting vegetation changes. Considering that the resolution for our classification results was 10 m, we chose the 10-m resolution NDVI to calculate the GLCM matrix. In addition, to manage the interannual variation in NDVI, the median value of the annual NDVI was utilized in calculating the GLCM matrix. This is calculated in GEE by applying the "glcmTexture()" function; however, three parameters must first be set: size, kernel, and average. 51 Because the individual tea plantations in the study area were generally small and scattered, the size was set to 3, which means that each neighborhood was a 3 × 3 area. The remaining parameters used default values. We used 12 texture features with lower correlations, including angular second moment, contrast, correlation, inverse difference moment, entropy, sum average, sum entropy, difference entropy, information measure of correlation 1, information measure of correlation 2, cluster shade, and cluster prominence.

Tasseled cap transformation features
TCT features refer to those obtained by TCT, which is a result of special principal component analysis. 52 Unlike general principal component analysis, the conversion coefficient of TCT is a fixed transformation matrix. 53 TCT transforms the original multispectral image into the actual Table 2 Remote sensing indices used in this study and their equations for calculation.
physical meaning of brightness (comprehensive surface reflectance), greenness (integration of the degree of surface vegetation coverage), and wetness (comprehensive surface water conditions) three-dimensional feature space, which can fully reflect surface reflection, vegetation coverage, and moisture information. 32 The TCT process can reduce the feature dimensions and enhance remote sensing image information. The specific transformation equation is as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 6 2 where Y represents the image after TCT, C represents the coefficient matrix corresponding to the TCT, X represents the original multispectral image before transformation, and a is a constant added to avoid negative values. The coefficient matrix C is a fixed transformation matrix that depends on the specific sensor. The TCT matrix coefficients of Sentinel-2 images are shown in Table 3. A total of six bands are required for TCT, including blue, green, red, NIR, SWIR1, and SWIR2.

Topographic features
Topographic features are widely used to distinguish ground objects. 54 Tea trees have extremely strict environmental requirements for growth; thus, the terrain in a tea tree plantation should be neither too flat nor too steep. 9 In flat areas, rainfall will accumulate at the bottom of the tea trees and cause the roots to rot. In steep areas, it is difficult to maintain essential soil nutrients and moisture for tea tree growth, and tea farmers have difficulty picking tea. Therefore, tea trees can only grow in less steep mountainous or hilly areas. Therefore, terrain features can effectively improve the recognition accuracy for tea plantations and other objects, owing to the different terrain requirements for each object. 2 For example, the objects in plains are typically cropland, built-up land, and water bodies, whereas those in mountainous and hilly areas are mainly forests and grasslands with only a small percentage of cropland and built-up land. Three topographic features were considered in this study: elevation, slope, and aspect. Elevation directly affects the temperature in mountainous areas, thereby indirectly affecting the distribution of vegetation. 11 The slope expresses the degree of steepness of a surface, which affects not only the surface energy exchange rate but also human activities and ultimately the spatial distribution of ground objects. 55 The aspect, which is the facing direction of the slope, indirectly affects the growth of vegetation and ultimately affects its spatial distribution. 5 In GEE, these topographic features can be directly calculated through the "ee.Algorithms.Terrain()" function.

Optimization of classification features
Considering a phenological perspective, the maximum, minimum, median, and standard deviation were calculated for the original spectral features, remote sensing indices, and red-edge vegetation indices, resulting in a total of 110 classification features. Notably, strong correlations can exist among these features, resulting in low classification accuracy and efficiency. 56 Therefore, this does not mean that more features are better, and appropriate classification features must be selected.
The J-M distance method is relatively simple and intuitive for directly measuring the distinguishability of different classes; therefore, we utilize it for feature optimization. 57 Essentially, the J-M distance uses the distance between samples of different classes to measure their separability, thereby determining distinguishable features based on the training samples. 58 The larger the J-M distance is, the higher the separability between classes is, and the easier classification it is to them. The calculation method for J-M distance is shown as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 9 9 J-Mðc i ; c j Þ ¼ where x is a time series, c i and c j are the two classes to be compared, and pðxjc i Þ and pðxjc j Þ are the conditional probability density functions for the time series x. When only considering the separability between two classes, Eq. (27) can be simplified to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 6 1 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 5 6 8 where B ij is the Bhattacharyya distance of classes i and j; V i and V j are the sample covariance matrices of i and j, respectively; and M i and M j are the mean vectors of the corresponding samples. The values of the J-M distance range from 0 to ffiffi ffi 2 p . The larger the J-M distance is, the stronger the separability between the two classes is, with ffiffi ffi 2 p corresponding to the largest separation. If JM ij ≥ 1, there is no overlap between classes i and j; if JM ij < 1, overlap exists between classes i and j.
According to Bruzzone et al., 59 when performing multiclass cases, the most important classification features for classification accuracy are determined by the classes with the lowest discrimination. In our classification process, tea plantations and natural forests were frequently misclassified, whereas combinations of other LULC types were easy to classify. The distinguishability between tea plantations and natural forests directly impacted the overall accuracy (OA) of classification results; therefore, the J-M distance between tea plantations and natural forests was calculated to optimize the classification features.

Progressive Random Forest Classification Approach
The progressive random forest classification method used in this study combines prior knowledge and human intervention regarding the use of training samples for machine learning. Progressive random forest classification is similar to classical supervised classification; the goal of both is to obtain the highest possible accuracy through human-computer interaction. The difference is that supervised classification improves classification accuracy by marking a large number of samples before classification, whereas progressive random forest classification uses a strategy based on classification accuracy. After a classification accuracy assessment, which sorts the probability of correct classification for each pixel from low to high, some of the pixels with probability values below a specific value (80% was used in this study) are manually added to the training sample. From the perspective of practical application, progressive random forest classification assumes that the lower accuracy samples are always located in areas that are difficult to discriminate, which can quickly and continuously promote the performance of the classifier. The overall performance of the progressive random forest classification model is controllable and does not fall into a local optimum.
The progressive random forest classification model reflects the human learning process, progressing from easy to difficult. The model is trained using a small number of training samples based on classification accuracy, and complex samples are gradually added to improve the classification accuracy. In our study, a typical tea plantation distribution area was first selected, and a small number of training samples were used for classification. Subsequently, the validation samples were used to calculate the OA and the confidence of each pixel. If the OA was high, the classification area was expanded. If the OA was low, low-confidence pixels were added to the training sample and reclassified. After reclassification, if the OA improved, the newly added training samples were retained; if the OA decreased, the newly added training samples were removed and reselected. The above steps were repeated, expanding the classification area by gradually increasing the number of training samples; thus, the classification accuracy was gradually improved until the entire study area achieved a satisfactory classification accuracy.
As this study aimed to classify tea plantations over a large area, it was preferable to add lowprecision tea plantation sample points before other types when applying the progressive random forest method. This ensured that the classification result for each pixel point within a tea plantation has high confidence. When applying the random forest model in GEE, two parameters must be set: numberOfTrees and minLeafPopulation. After repeated verification, we set numberOfTrees to 100 and minLeafPopulation to 10.

Postprocessing and Accuracy Assessment
There were a small number of blank values in the classification result, which directly impacted the classification accuracy. There are two main reasons that blank values are generated. One reason is that some pixels have classification features outside the normal range, resulting in empty classification results. For these pixels, we rechecked the values for all classification features, eliminated outliers, and then reclassified them. The other reason is that some pixels are recognized as clouds by the cloud recognition method. For these, a multivalue within the 3 × 3 neighborhood of these pixels was assigned.
Tea trees prefer shade; hence, many tea plantations are in the shadows of hills. Thus, some of these plantations were easily misclassified as natural forests; therefore, we also used the 3 × 3 neighboring pixels to eliminate the influence of the hill shadows.
After postprocessing, we used a 10-fold cross-validation to assess the classification accuracy. We calculated the classification error matrix based on the classification results, which were quantified using the OA. The OA is the sum of the diagonal elements divided by the sum of all of the elements of the confusion matrix. Then, we calculated the producer's accuracy (PA), user's accuracy (UA), and F 1 -score. The F 1 -score is particularly useful for class-level accuracy assessment because it assigns equal importance to the PA and UA. Finally, the above process was repeated 10 times, and the average OA, PA, UA, and F 1 -score were calculated. The PA, UA, and F 1 -score are calculated for each class as follows: where r is the number of classes and n ij is the element of the confusion matrix in the i'th column and j'th row, that is, the count of elements of class j classified as class i. PA i , UA i , and ðF 1 Þ i stand for the PA, UA, and F 1 -score for class i, respectively.

Optimized Classification Features
In this study, 110 classification features were optimized to ensure classification efficiency and achieve a higher classification accuracy. Both the topographic and TCT features had only three features and weak correlations; hence, they were excluded from the feature optimization process. Therefore, the rest of the classification features, including the original bands, remote sensing indices, and texture features, were optimized using the J-M distance. The feature optimization process was conducted in the following two steps. First, 1000 tea plantation samples and 1000 natural forest samples were randomly selected from the training samples. Then, the J-M distance for each feature was calculated (see Table 8 for details), and the features with J-M distance values lower than a certain standard were eliminated. The J-M distance can determine the discrimination between different classification bands but not the number of optimal features for classification. Here, we selected the optimal feature bands based on Lorenzo Bruzzone et al. 59 As shown in Table 8, the maximum, minimum, and median values of the J-M distance of the original bands were 0.483, 0.085, and 0.229, respectively. Thus, the tea plantations and natural forests were difficult to distinguish. The original bands with J-M distance values greater than 0.2 were selected. Similarly, the characteristics of tea plantations and natural forests were also very similar in the remote sensing indices. Those with J-M distances greater than 0.2 were also selected. However, the distinguishability between tea plantations and natural forests was very strong for the red-edge vegetation indices, which are more sensitive concerning vegetation. Therefore, the red-edge vegetation indices with J-M distances greater than 0.4 were selected. A brief description of the optimized classification features, including the number and J-M distance values, is shown in Table 4 (for details, see Table 9). The total number of classification features was reduced from 116 to 64, a 45% decrease, which dramatically improved the classification efficiency. The overall J-M distance values varied from 1.414 for the original classification features to 1.412 for the optimized classification features, which is a difference of 0.14%; therefore, a comparable result was obtained after the feature optimization process. From the perspective of balancing the classification efficiency and accuracy, the feature optimization process performed excellently.

Classification Results for the Tea Plantations and Other Land
Cover Types Figures 3 and 4 show the results of the tea plantation classification using the progressive random forest classification model from a regional and local perspective, respectively. As shown in Fig. 3, tea plantations were scattered throughout in Anhui Province. Tea plantations were mainly distributed in the southern and western study area, which predominantly have a mountainous terrain. In the middle of the study area, there were a small number of tea plantations in hilly terrain. There were few tea plantations in the northern study area or along the Yangtze River; these places are plains that are not suitable for the growth of tea trees. The city with the most tea plantations was Huangshan City, which was consistent with the statistical data. Figure 4 shows the classification results for three typical regions. The tea plantations identified using the progressive classification model were generally consistent with reality. In hilly areas, tea plantations were well identified, and the cropland and forest classifications were accurate [ Fig. 4(b)]. In the transitional regions between hills and mountains, the overall tea plantation classifications were also relatively accurate; however, a small amount of forest was misclassified as tea plantations [ Fig. 4(d)]. The classification results for forests and water bodies were more accurate; however, some cropland was misclassified as built-up land. In mountainous areas, the classification results for tea plantations and forests were accurate; however, some cropland was again misclassified as built-up land [ Fig. 4(f)]. In general, tea plantations were accurately classified in the main distribution areas.

Classification Accuracy Analysis
To test the classification accuracy for tea plantations, the OA, UA, PA, and F 1 -score were calculated quantitatively using the confusion matrix. In the quantitative analysis, a 10-fold crossvalidation was used to determine the stability of the classification model. The calculation results are presented in Table 5. The average OA reached 89%. The classification accuracy for cropland, forest, and water bodies was greater than 85%, and that for tea plantations reached 82%. With the spatial resolution and the large area of Anhui Province, this result is generally satisfactory. However, the classification accuracies for grassland and unused land were poor, which was mainly caused by their small distribution areas within the study area and the small number of sample points used during classification. Nevertheless, these low accuracies had little influence on the identification of tea plantations because of their small areas and low probability of misclassification as tea plantations.

Impacts of Different Feature Selection Schemes
To compare the impacts of different classification features on the classification accuracy for tea plantations, this study designed 10 different classification feature combination schemes. Because the red-edge bands had the largest J-M distance, comparisons with and without these bands were carried out to analyze whether they have prominent effects when designing the feature combination method. Due to the large study area and the number of comparative experiments, only Huangshan City, which has the most tea plantations, was selected for comparative experiments. The experimental results are presented in Table 6. Table 6 shows that the combination of the original bands and the remote sensing indices obtained a high OA value of 83.99% because the constructed classification feature space could extract the phenological features of each class. The red-edge indices had the most significant impacts on improving the classification accuracies for tea plantations due to their high sensibility for vegetation. The topographical features had the second largest positive effect on classification accuracy. Tea plantations have strict environmental requirements; therefore, elevation, slope, and aspect all have direct effects on the growth of tea trees. The TCT features also showed positive effects because they comprehensively reflect the brightness, greenness, and wetness of ground objects. For example, the water in the roots of tea plantations and natural forests differs, resulting in a difference in TCT features. Tea trees are planted at certain spatial intervals to facilitate both growth and harvesting; thus, their greenness is generally lower than that of natural forests. The  smallest improvement in the accuracy was from the texture features. This was mainly because the spatial resolution of the Sentinel-2 data was 10 m, which does not fully reflect the textural features of the tea plantations, resulting in a minor improvement. When using all features, the UA for tea plantations was 89.12%, the PA was 95.25%, and the F 1 -score was 0.9208. In addition, these accuracy indices obtained using the optimized features were 88.67%, 95.33%, and 0.9188, respectively. Although the UA is slightly decreased when using the optimized features, the PA improved, and the F 1 -score, which measures the comprehensive classification accuracy, dropped by 0.002. The OA using all features was 86.50%, whereas that using the optimized features was 86.41%, which is a decrease of 0.09%. Although the use of optimized features resulted in a slight loss of accuracy, the number of classified features was dramatically reduced from 116 to 64. A reduction in the number of features can significantly improve the classification efficiency. When using the progressive random forest model, it is necessary to continuously add training samples and then reclassify them. Therefore, compared with the slight loss of accuracy, the improvement in the classification efficiency is more important.

Improvement in Classification Accuracy by Progressive Random Forest
In this study, we designed a progressive random forest model that combines random forest and progressive models. The progressive random forest model improves classification accuracy by gradually increasing the number of training samples on GEE. These newly added training samples were based on the probable classification result for each pixel. A total of 10 iterative classifications were performed, and the results are shown in Fig. 5. Figure 5 shows that, in the process of iterative classification, the OA, UA, and PA of the tea plantations constantly improved. This shows that the proposed progressive random forest model was effective. For the first iteration, the samples were derived from the results of our field surveys in Huangshan and Wuhu. Huangshan and Wuhu are typical mountainous and hilly tea-producing cities in the study area, respectively. For the first iteration, the samples were concentrated in mountainous and hilly areas, resulting in low OA, UA, and PA for the tea plantations. Then, the visually selected training samples for other cities were gradually added according to the accuracies of the classification results. With additional training samples, the OA, UA, and PA greatly improved. On the 10th iteration, the OA and PA exceeded 85%, which met our objectives; hence, we stopped the classification process.
The first visual samples were from the areas surrounding Huangshan and Wuhu; subsequently, the area was gradually expanded. The results of the visually selected additions and land survey samples are shown in Fig. 6. The figure shows that there were more samples in the southern study area but few samples in the northern region. This is because tea plantations and natural forests were initially misclassified; thus, we added samples of tea plantations and natural forests. With the increase in natural forest samples, tea plantations were mainly misclassified as croplands. Therefore, we added samples of cropland in the middle of the study area. When the UA and PA of the tea plantations exceeded 80%, we added samples from other classes to improve the OA.

Advantages and Disadvantages
Because the spectral bands of tea plantations and natural forests are highly similar, obtaining high-precision classification of tea plantations over a large area remains challenging. With readymade remote-sensing cloud platforms (e.g., GEE) and the use of nonparametric classifiers, the main challenge becomes selecting suitable classification bands and training samples.
The high classification accuracy and F 1 -score for tea plantations and other main land use classes demonstrate the practicality of the progressive random forest method. The effectiveness of this method depends on accessing sufficient computing power. In this study, this was not a problem because the GEE provided sufficient computing power.
In this study, a new method for the classification of sentinel-2 MSI data, which was combined with a random forest and progressive model running on GEE, was proposed. This is advantageous for the following reasons: (1) it retains the advantages of the random forest and the progressive model, which can effectively and robustly run high-dimensional data for a large Fig. 6 The results of the field survey and the addition of visually selected training samples.
area. (2) It can classify sparse objects (such as the tea plantations in this study) with a high classification accuracy. (3) The entire classification process is carried out on GEE to ensure highly efficient overall workflow.
We used classification and regression trees (CART), support vector machine (SVM), and deep neural networks (DNN) to verify the generality and reliability of our conclusions. In the comparative experiments, the samples used for the CART, SVM, and DNN models were the same as for the RF models. Default parameters were set for CART and SVM in GEE. The Keras library was utilized for the DNN hyperparameter optimization in Google Colab, and the key parameters of DNN were set as follows: 64-128-256 for the hidden layers, 10 for the epochs, ReLU (the rectified linear units) for the activation function, Adam for the optimizer, and 20% for the dropout ratio. The OA values of the comparative experiment classifications are shown in Table 7. In general, when a small number of samples was used, a low OA value was obtained. In addition, the OA values of the compared classification models increased with the increase in the number of samples, indicating that increasing the number of samples could significantly improve classification accuracy, which was similar to the RF model. However, the contributions of increasing samples to the OA were found to be different among the compared classification models. The OA values of the CART model were moderately improved with the increase of samples, but the highest OA of the CART model (79.52%) was significantly lower than that of the RF model (89.27%). The SVM model only needs a small number of samples to obtain high classification accuracy, but it easily overfits and its highest OA (83.89%) was also lower than that of the RF model. Increasing the number of samples can dramatically improve the classification accuracy of the DNN model, which needs a large number of samples to obtain high classification accuracy. Although the highest OA of DNN (75.96%) was significantly lower than that of the RF model, it is worth noting that the classification accuracy value of the DNN model continuously improved as the samples increased, which means that it may be close to or exceed the highest classification accuracy of the RF model if the number of samples is more than 6193. For balancing the classification accuracy and the limit number of samples, the combination of the progressive model and the RF model is the optimal choice.
Although our progressive random forest model has many advantages, it also has some disadvantages. First, because the GEE visualization interactive platform only provides current remote sensing images, the training data added by manual visual inspection only utilizes the current data, and the generation of sample points using historical data requires other platforms. Second, progressive random forest is semiautomated. Manual assignment is required each time that new samples are added to the training sample dataset. In future studies, we will consider introducing models (such as LandTrender) to generate long-term sample data. We are also considering adding other existing land use/land cover products to reduce manual assignment and improve the automation of the classification workflow.

Conclusions
A methodology was proposed for the classification of tea plantations in Anhui Province in 2020.
The J-M model distance was applied to the original bands, remote sensing indices, texture features, TCT features, and terrain features to select the optimal combination of classification features. A progressive random forest model was used to classify 3554 scenes from Sentinel-2 MSI images. Using these methods, we obtained a map of tea plantations in Anhui Province in 2020 and analyzed its classification accuracy. The progressive random forest classification method used in this study combines the advantages of the random forest and progressive models. First, remote sensing classification for a large area was carried out quickly by utilizing random forest. Moreover, owing to the advantages of the progressive model, the OA of classification remained high even with a small number of training samples. Furthermore, the progressive random forest classification model accurately identified the spatial distribution of tea plantations.
By comparing the J-M distance, the number of classification features was reduced dramatically, which improves the classification efficiency. Moreover, high classification accuracy was obtained using the reduced classification features.
In addition, when using only the original bands and remote sensing indices, good classification results were still obtained. This is because the calculations of statistical values, such as the maximum, minimum, and standard deviation, reflect the phenological characteristics of the vegetation. Furthermore, the red-edge indices are the most effective for the identification of tea plantations, followed by topographical features. The TCT and texture features have minimal positive effects on improving the classification accuracy for tea plantations.
In future studies, we will explore the effects of applying other machine learning models to the classification process along with the progressive model to obtain a more precise map of tea plantations. Meanwhile, we will consider extending our research workflow using Landsat time-series data to obtain long-term spatial distributions and changes in tea plantations.  Table 8 shows the details of the calculated J-M distance for each feature. The low value of the J-M distance represents low separability and the features with J-M distance values lower than a certain standard (i.e., 0.2) were further eliminated. Table 9 shows the details of the optimized classification features, including the number and the J-M distance values. The total number of classification features was reduced from 116 to 64, a 45% decrease, which dramatically improved the classification efficiency.