Remote Detection of Large-Area Crop Types: The Role of Plant Phenology and Topography

Sustainable agricultural practices necessitate accurate baseline data of crop types and their detailed spatial distribution. Compared with field surveys, remote sensing has demonstrated superior performance, offering spatially explicit crop distribution in a timely manner. Recent studies have taken advantage of remote sensing time series to capture the variation in plant phenology, inferring major crop types. However, such an approach was rarely used to extract detailed, multiple crop types spanning a large area, and the impact of topography has yet to be well analyzed in mountainous regions. This study aims to answer two questions in crop type extraction: (i) Is it feasible to accurately map multiple crop types over a large mountainous area with phenology-based modeling? (ii) What are the effects of topography in such modeling? To answer the questions, phenological metrics were extracted from MODIS (Moderate Resolution Imaging Spectroradiometer) satellite time series, and the random forests classifier was used to map 12 crop types in South China (236,700 km2), featuring a subtropical monsoon climate and high topographic variation. Our study revealed promising results using MODIS EVI (Enhanced Vegetation Index) and NDVI (Normalized Difference Vegetation Index) time series, although EVI outperformed NDVI (overall accuracy: 85% versus 81%). The spectral and temporal metrics of plant phenology significantly contributed to crop identification, where the spectral information exhibited greater importance. The increase of slope led to a decrease in model accuracy in general. However, uniformly distributed tree plantations (e.g., tea-oil camellia, gum, and tea trees) being cultivated on large slopes (>15 degrees) achieved accuracies greater than 80%.


Introduction
Estimates of crop types and extent can provide essential baseline data for managing agricultural lands, determining food pricing, designing trade policies, and supporting carbon balance research [1][2][3]. Crop mapping is particularly important for the tropical or subtropical mountainous regions. On the one hand, crop types are highly diverse throughout the regions, mainly due to climate and significant agricultural expansion at the expense of intact forests [4,5]. On the other hand, high topographic complexity in mountains can cause significant logistical challenges for land surveying. For example, the Guangxi Zhuang Autonomous Region (hereafter Guangxi) in South China has a subtropical monsoon climate with cool, dry winters and long, hot, wet summers [6]. Enjoying the reputation of "natural greenhouse", Guangxi is the leading fruit cultivation area and the largest winter vegetable base in China [7,8]. However, accurately estimating crop types and extent remains a challenging task, because over half of the region is dominated by karst limestone formations, giving rise to complex topography across much of Guangxi [6].
Crop mapping with single-date remotely sensed imagery has been a traditional practice for over four decades [9][10][11][12][13]. But, it is problematic to use a snapshot of the landscape to distinguish between crop types that have high spectral similarities [14]. Recent studies have taken advantage of satellite image time series and vegetation phenology, i.e., capturing seasonal variations of plant spectral reflectance in its life cycle, to improve the accuracy in crop mapping. The main assumption is that every crop type under investigation has a unique phenology [15]. Applying image time series to characterize phenological phases is similar to change detection, where denser time series data are more likely to capture the detailed crop functioning changes (or trajectories) across seasons. While medium-resolution sensors Landsat and the more recent Sentinel image time series have proven effective in crop mapping [16][17][18][19], MODIS (Moderate Resolution Imaging Spectroradiometer) data with a high-revisiting frequency are more popular in tropical and subtropical studies [20][21][22][23], reducing the challenges of collecting cloud-free satellite images.
Phenology-based crop estimation with dense MODIS time series has successfully extracted single crop types from agriculture-dominated regions, including paddy rice [20,23,24], winter wheat [15,25], maize [26], winter/nonwinter crops [27], and rubber plantations [21]. In addition, previous studies have made extended efforts to concurrently map several crop types in one model or a single framework. For instance, Wardlow and Egbert [28] successfully estimated the extents of alfalfa, summer crops, winter wheat, and fallow in Kansas, one of the topographically flattest states in the United States. Wang et al. [29] identified corn, soybean, and winter wheat in western Missouri, a state adjacent to Kansas. While the results are promising, the performance of phenology-based mapping has yet to be confirmed if multiple crop types (e.g., over 10) are the mapping subjects. On the one hand, the phenology of different crop types may show high similarities. On the other hand, the same type of crop and its phenology may demonstrate high variation caused by the differences in microclimate, e.g., illumination and humidity, [30], which is particularly true in the area featuring complex topography. To the best of our knowledge, few studies have evaluated the effects of topography in phenology-based crop monitoring.
The main goal of this study is to evaluate the feasibility of applying MODIS dense time series and vegetation phenology to extract detailed crop types in Guangxi, a large (236,700 km 2 ) and mountainous province in South China. The task is emphasized on 12 major types of crops cultivated in the region. To achieve the goal, this study asks two specific questions: (1) Is it feasible to accurately map multiple crop types with phenology-based modeling? (2) What are the effects of topography in such modeling?

Study Area
Guangxi is situated in South China (104 • 26 -112 • 04 E, 20 • 54 -26 • 24 N). It meets the edge of the Tibetan Plateau in the west, connects to the southwestern mountainous region of Hunan Province in the northeast, and joins the western mountainous region of Guangdong Province in the east. It also shares an approximately 637 km long border with Vietnam in the southwest (Figure 1). Guangxi encompasses a total land area of 236,700 km 2 (2.47% of the total land area of China), making it the ninth largest provincial-level (provinces, autonomous regions, and direct-controlled municipalities) administrative division in China [31].
Topographically, Guangxi is low in the southeast and mostly surrounded by high mountains with low-lying basins, plains, and hills as well as curved mountainous ranges in the center, and thus bears resemblance to a dustpan-shaped basin. As one of the rainy regions in China, most areas of Guangxi have an annual precipitation of 1250-1700 mm, with the majority of precipitation (generally over 75% of the total annual precipitation) occurring in the spring and summer (March through August) and a relatively small amount of precipitation (approximately 25% of the total annual precipitation) occurring in the fall and winter (September through February of the following year) [31]. Major crop Topographically, Guangxi is low in the southeast and mostly surrounded by high mountains with low-lying basins, plains, and hills as well as curved mountainous ranges in the center, and thus bears resemblance to a dustpan-shaped basin. As one of the rainy regions in China, most areas of Guangxi have an annual precipitation of 1250-1700 mm, with the majority of precipitation (generally over 75% of the total annual precipitation) occurring in the spring and summer (March through August) and a relatively small amount of precipitation (approximately 25% of the total annual precipitation) occurring in the fall and winter (September through February of the following year) [31]. Major crop types in Guangxi include rice, tomatoes, sugarcane, mangoes, bananas, grapes, tea trees, tobaccos, mulberries (Morus), tea-oil camellia (Camellia oleifera), and gum trees (Eucalyptus) [31].

Field Data
A modified stratified random sampling strategy was employed to collect a total of 1462 samples circa 2016, covering 12 crop types-gum trees (69 samples), tea trees (82), tomatoes (33), sugarcane (112), mangoes (122), grapes (60), tea-oil camellia (83), mulberries (54), rice (215) ,bananas (116), tobaccos (37), and other vegetables (14)-and four noncrop land-cover types-woodland (159), water (54), impervious surfaces (152), and bare ground (100). Specifically, the study area was divided into strata (i.e., 14 administrative divisions) corresponding to major crop cultivation regions. Our selection of strata was based on the following considerations. First, each of these divisions has relatively homogenous natural environments within the region (e.g., climate, soil type, and topography). Second, local governments often promote certain crop types through agricultural subsidy programs. Within each stratum, field surveys were combined with photo interpretation of Google Earth (Google LLC, Mountain View, CA, USA) high-resolution imagery to extract samples in a random fashion. Here, photo interpretation was used to supplement field surveys. To ensure accurate photo interpretation, field-surveyed polygons were imported to Google Earth, which allowed us to compare crop types and image spectral/spatial patterns. To mitigate the cloud/illumination effects in Google imagery, the images taken from the neighboring years were also compared. Due to logistical

Field Data
A modified stratified random sampling strategy was employed to collect a total of 1462 samples circa 2016, covering 12 crop types-gum trees (69 samples), tea trees (82), tomatoes (33), sugarcane (112), mangoes (122), grapes (60), tea-oil camellia (83), mulberries (54), rice (215), bananas (116), tobaccos (37), and other vegetables (14)-and four noncrop land-cover types-woodland (159), water (54), impervious surfaces (152), and bare ground (100). Specifically, the study area was divided into strata (i.e., 14 administrative divisions) corresponding to major crop cultivation regions. Our selection of strata was based on the following considerations. First, each of these divisions has relatively homogenous natural environments within the region (e.g., climate, soil type, and topography). Second, local governments often promote certain crop types through agricultural subsidy programs. Within each stratum, field surveys were combined with photo interpretation of Google Earth (Google LLC, Mountain View, CA, USA) high-resolution imagery to extract samples in a random fashion. Here, photo interpretation was used to supplement field surveys. To ensure accurate photo interpretation, field-surveyed polygons were imported to Google Earth, which allowed us to compare crop types and image spectral/spatial patterns. To mitigate the cloud/illumination effects in Google imagery, the images taken from the neighboring years were also compared. Due to logistical challenges, some sample locations were adjusted, making them close to roads for easy access. Each sample was chosen from a relatively homogeneous site (within a 300 m cell), reducing the spectral mixing effects in MODIS pixels. Since identifying crop types in certain seasons may face challenges (e.g., due to land clearing for new seeding), the high accuracy of cropland recognition was ensured on the basis of our extensive experience in conducting land management projects in the region and communicating with local farmers to understand planting schedules and history.

Satellite MODIS Time Series
The MODIS MOD13Q1 level-3 product, consisting of 250 m resolution, 16 day composites was used as remote sensing input data. Both EVI (enhanced vegetation index) and NDVI (normalized difference vegetation index) have been used in crop mapping [21,28]. In this study, both indices were used, with their performance compared in the estimation of multiple crop types. The data spanned three years from 2015 to 2017. Although our study was focused on 2016, the three years of data facilitated an accurate derivation of crop phenology, as compared with the use of only one year of data ( [32]; see detail below).

Digital Elevation Model (DEM)
A Digital Elevation Model (DEM) at the 30 m resolution was acquired from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM V2) [33]. This DEM product encompasses 99 percent of Earth's landmass. It has proven effective in our previous remote sensing vegetation studies (e.g., [34]).

Derivation of Phenological Metrics
Crop phenology was extracted using the TIMESAT program package [32]. In the current version of TIMESAT, multiple key phenological metrics can be computed for each full season of a time series. Here, 11 metrics that have been considered as crucial seasonality parameters [32] were calculated, including time for the start of the season (SOS), time for the end of the season (EOS), length of the season (LOS), base level (BL), time for the mid of the season (MOS), largest data value for the fitted function during the season (MaxV), seasonal amplitude (Amp), rate of increase at the beginning of the season (RIBOS), rate of decrease at the end of the season (RDEOS), large seasonal integral (LSI), and small seasonal integral (SSI). See Table 1 for explanations of the parameters. Those parameters can capture various aspects of the seasonal variation in the life cycle of a crop [35].
Prior to the calculation of phenological metrics, fitting functions were applied to reduce the effects of atmospheric noises in EVI (Enhanced Vegetation Index) and NDVI (Normalized Difference Vegetation Index) time series. Three popular fitting algorithms provided by TIMESAT were compared: Savitzky-Golay filter (SG), asymmetric Gaussian (AG), and double logistic function (DL). These three algorithms have proven to be highly competitive for filtering time series data [36]. SG applies a moving window to preserve the area and mean position of a seasonal peak [37], with window width defined by the user. After a trial-and-error test, the window width was set to 4.0. AG and DL fit local functions to data in intervals around maxima and minima in the time series, with the difference being the function type, i.e., Gaussian function for AG and double logistic function for DL [35]. For all the three fitting functions, the season start and season stop value was consistently defined as 0.1 and 0.9, respectively. The values were determined to maximize the value ranges of phenological metrics, facilitating the succeeding classification using the metrics. The default values were used for all the other input parameters. Figure 2 shows the original EVI observations for sampled pixel areas and the corresponding fitted curves using SG, AG, and DL, while Figure 3 shows the same type of curves using NDVI time series as input. Here, three years (2015-2017) of data were used to simulate crop cycles. This was mainly due to that data from surrounding years can reduce the risk for erroneous determinations in noisy time series data (Jönsson and Eklundh, 2004). Such strategy has proven effective in previous studies, e.g., Hird et al. (2009), Senf et al. (2013). The performance of the fitting algorithms was evaluated with two criteria. First, the root of the mean squared difference (RMSE) between the original EVI (or NDVI) observations and the modeled values is minimized. To do so, 70% of the 1462 sampled crop pixels taken separately from each group were used to calculate RMSEs for all the three algorithms. Second, the fitted result of the algorithm has minimized null-value areas. Based on the two criteria, SG was found to be the optimal, balancing the needs for low RMSE and small null-value areas. It was used for fitting both the NDVI and the EVI time series in our study. Once the trajectories had been simulated, the 11 metrics were calculated using the fitted curves in 2016, leading to crop classification in the same year. Table 1. Explanations of the extracted phenological metrics [32]. EVI: enhanced vegetation index, NDVI: normalized difference vegetation index.

Phenological Metric Explanation
Time for the start of the season (SOS) Time for which the value (EVI/NDVI) has increased to a user-defined level (i.e., 10% of the distance between the left minimum level and the maximum).
Time for the end of the season (EOS) Time for which the right edge has decreased to a user-defined level measured from the right minimum level (i.e., the 10% of the distance between the right minimum level and the maximum). Computed as the mean value of the times for which, respectively, the left edge has increased to the 80% level and the right edge has decreased to the 80% level.
Largest data value for the fitted function during the season (MaxV) It may occur at a different time compared with MOS.
Seasonal amplitude (Amp) Difference between the maximum value and the base level.
Rate of increase at the beginning of the season (RIBOS) Calculated as the ratio of the difference between the left 20% and 80% levels and the corresponding time difference.
Rate of decrease at the end of the season (RDEOS) Calculated as the absolute value of the ratio of the difference between the right 20% and 80% levels and the corresponding time difference. The rate of decrease is thus given as a positive quantity.
Large seasonal integral (LSI) Integral of the function describing the season from the season start to the season end.

Small seasonal integral (SSI)
Integral of the difference between the function describing the season and the base level from season start to season end.

Random Forests Classification and Accuracy Assessment
Random forests (RF) are ensembles of decision trees wherein each tree is independently determined using a bootstrap sample of the data set, and a simple majority vote is taken for final prediction [38,39]. In this study, RF was applied as the classifier to estimate crop types and map their extents. The 11 tested phenological metrics were used as input. The metrics were calculated for individual crop types using their crop cycles (see Section 4.1). By taking advantage of the RF's unique characteristic of ranking the relative importance for the predictor variables, the research team aimed to understand the importance of phenological metrics in estimating different types of crops. The calculation of importance was based on the optimal classification performance splitting the training data, which simultaneously returned an array of each metric's importance. According to our previous experience in vegetation mapping [40] and the exploratory trials with the current data, two major parameters were specified in RF: 100 for the number of trees and 3 for the number of predictor variables for each tree. 70% of the field samples were randomly selected for training RF. The remainder was used for validating the model's performance with overall, and user's and producer's accuracies reported. The classification was completed in the EnMAP-Box environment, which is a free and open source software [41].
The topographic effects on crop mapping accuracy were evaluated using slope and aspect. Slopes were classified into five groups: ≤2 • , 2-6 • , 6-15 • , 15-25 • , and ≥25 • . The grouping follows the suggestion by the "Technical Specification for the Survey of Land Use Status" issued by the China's Agricultural Regional Planning Commission. Then, classification accuracies for each of the slope groups were summarized. Similarly, aspects were categorized into north, northeast, east, southeast, south, west, southwest, and northwest, with classification accuracies summarized for each aspect group.

Classification Accuracy Assessment
Crop classification using EVI and NDVI dense time series both achieved an overall accuracy greater than 80%. Using EVI data as input (overall accuracy: 85.27%) revealed generally better performance than using the NDVI data (overall accuracy: 81.01%). Figure 4 further shows the producer's and user's accuracies for individual crop types. For most of the crops, the classification results were comparable using EVI versus NDVI, although EVI led to better performance for estimating tomatoes, tea-oil camellia, and gum trees. While both NDVI and EVI have been used to characterize vegetation on croplands [15,20,21,27], our classification and comparison using 12 crop types suggest a preference for using EVI over NDVI in crop mapping. This finding can be explained by EVI's better responses to vegetation biophysical parameters through reducing background and atmospheric noises [42], which are common in subtropical low-growing and sparsely distributed crops. It is also consistent with the conclusions made by several other researchers [21][22][23]. Hence, the following discussion is emphasized on the results using the EVI time series.
The comparison among crop types revealed the highest accuracies in mapping gum trees, tea trees, tomatoes, tea-oil camellia, and mulberries, followed by mangoes and grapes, with both the user's and producer's accuracies being greater than 90% (Figure 4). The mapping accuracies for sugarcane, rice, and bananas were greater than 70%, with the lowest performance found in mixed vegetables (Figure 4). The results are comparable with those from earlier studies using phenology-based approaches and MODIS time series, although crop types varied across regions and studies. In general, mapping accuracies are higher than 70% for all major crops. For example, a recent study by Boschetti et al. (2017) [22] achieved an 80% accuracy in rice mapping, while winter crops (e.g., wheat) showed superior performance between 85% and 95% [15,27]. The cash crop rubber trees demonstrated a slightly lower accuracy 74% [21]. In our study, the majority of the crop types had reached accuracies higher than 70% (Figure 4). Rice estimation was relatively weak, which was lower than the 80% accuracy as obtained by Boschetti et al. [22]. The main reason is possibly due to not distinguishing between once-a-year and twice-a-year harvested rice, which have different EVI trajectories. In addition, rice was cultivated Agriculture 2019, 9, 150 8 of 14 extensively in highly fragmented areas, causing spectral mixing in MODIS pixels. Apart from the similar spectral mixing issue, the areas with mixed vegetables often have shifted plant types across a year. The variation in crop-shifting practices may have also contributed to the lowest accuracy in mapping mixed vegetables. Because the collected field samples were from homogenous areas, they may have boosted the classification accuracies. In areas of mixed land cover and small farms, the 250 m resolution imagery was unable to capture fine-scale distributions of multiple crop types. The comparison among crop types revealed the highest accuracies in mapping gum trees, tea trees, tomatoes, tea-oil camellia, and mulberries, followed by mangoes and grapes, with both the user's and producer's accuracies being greater than 90% (Figure 4). The mapping accuracies for sugarcane, rice, and bananas were greater than 70%, with the lowest performance found in mixed vegetables (Figure 4). The results are comparable with those from earlier studies using phenologybased approaches and MODIS time series, although crop types varied across regions and studies. In general, mapping accuracies are higher than 70% for all major crops. For example, a recent study by Boschetti et al. (2017) [22] achieved an 80% accuracy in rice mapping, while winter crops (e.g., wheat) As demonstrated in the EVI-based crop classification map (Figure 5), the top five cultivated crop types in Guangxi were rice, gum trees, sugarcane, mangoes, and bananas (see their areal sizes in Table 2). Geographically, rice was cultivated mostly in the plains in central Guangxi, in northern Guangxi, and sporadically in other areas. The total cultivation area for rice reached 21,560.2 km 2 . Gum trees were cultivated extensively as well (20,200 km 2 ) and were distributed in the central, eastern, and southeastern areas of the region. They are followed by sugarcane, which was mainly cultivated in the city of Chongzuo in southern Guangxi, the cities of Nanning and Laibin in central Guangxi, the city of Liuzhou in northern Guangxi, and the coastal areas with an areal size of 19,400 km 2 . Mangoes were found in the Youjiang River Valley in the city of Baise as well as in the hilly and mountainous areas on both sides of the Youjiang River Valley. In addition, the city of Guilin in northern Guangxi and some small regions in southeastern Guangxi also had mango farms. Bananas are spread mostly in the cities of Nanning, Baise, and Chongzuo and in the southern area of the city of Hechi. The spatial distribution of the main crops was consistent with our local experience and observations via field surveys. Each of the other seven crop types covers less than 10,000 km 2 ( Table 2). According to recent statistical data (from field surveys) published by Guangxi government in the book Guangxi Almanac 2017, the cultivation area for rice in 2016 was reported at 19,598 km 2 , which is 9% less than our estimate 21,560.2 km 2 . While errors exist in both estimates, our work is able to demonstrate the detailed spatial distribution of rice.
samples were from homogenous areas, they may have boosted the classification accuracies. In areas of mixed land cover and small farms, the 250 m resolution imagery was unable to capture fine-scale distributions of multiple crop types.
As demonstrated in the EVI-based crop classification map (Figure 5), the top five cultivated crop types in Guangxi were rice, gum trees, sugarcane, mangoes, and bananas (see their areal sizes in Table 2). Geographically, rice was cultivated mostly in the plains in central Guangxi, in northern Guangxi, and sporadically in other areas. The total cultivation area for rice reached 21,560.2 km 2 . Gum trees were cultivated extensively as well (20,200 km 2 ) and were distributed in the central, eastern, and southeastern areas of the region. They are followed by sugarcane, which was mainly cultivated in the city of Chongzuo in southern Guangxi, the cities of Nanning and Laibin in central Guangxi, the city of Liuzhou in northern Guangxi, and the coastal areas with an areal size of 19,400 km 2 . Mangoes were found in the Youjiang River Valley in the city of Baise as well as in the hilly and mountainous areas on both sides of the Youjiang River Valley. In addition, the city of Guilin in northern Guangxi and some small regions in southeastern Guangxi also had mango farms. Bananas are spread mostly in the cities of Nanning, Baise, and Chongzuo and in the southern area of the city of Hechi. The spatial distribution of the main crops was consistent with our local experience and observations via field surveys. Each of the other seven crop types covers less than 10,000 km 2 ( Table  2). According to recent statistical data (from field surveys) published by Guangxi government in the book Guangxi Almanac 2017, the cultivation area for rice in 2016 was reported at 19,598 km 2 , which is 9% less than our estimate 21,560.2 km 2 . While errors exist in both estimates, our work is able to demonstrate the detailed spatial distribution of rice.

Importance of Phenological Metrics
Normalized variable importance for all the tested phenological metrics were derived from the RF classification ( Figure 6). Using EVI data as input, we found that the importance of SSI (1.20) was the highest in the model, which was comparable with that of BL (1.19), and was followed by SOS (1.09). The other metrics had relatively low contributions, with the weakest metric being RIBOS ( Figure 6). According to the metrics' definitions, SSI represents the integral of the difference between the function describing the season and the base level from season start to season end, and BL is the average of the left and right minimum values [32]. Both SSI and BL are mainly determined by the spectral values from EVI time series, while SOS (start of the season) and part of SSI are a representation of the crop's temporal characteristic. This suggests a critical role of both spectral and temporal information in crop-type classification, although the spectral information demonstrated slightly higher importance in the study. Our result is in an agreement with some of the earlier findings. For instance, in a project extracting rubber plantations, Senf et al. [21] discovered higher importance of the EVI's spectral metrics (e.g., EL) than their timing (e.g., LOS).

Importance of Phenological Metrics
Normalized variable importance for all the tested phenological metrics were derived from the RF classification ( Figure 6). Using EVI data as input, we found that the importance of SSI (1.20) was the highest in the model, which was comparable with that of BL (1.19), and was followed by SOS (1.09). The other metrics had relatively low contributions, with the weakest metric being RIBOS ( Figure 6). According to the metrics' definitions, SSI represents the integral of the difference between the function describing the season and the base level from season start to season end, and BL is the average of the left and right minimum values [32]. Both SSI and BL are mainly determined by the spectral values from EVI time series, while SOS (start of the season) and part of SSI are a representation of the crop's temporal characteristic. This suggests a critical role of both spectral and temporal information in crop-type classification, although the spectral information demonstrated slightly higher importance in the study. Our result is in an agreement with some of the earlier findings. For instance, in a project extracting rubber plantations, Senf et al. [21] discovered higher importance of the EVI's spectral metrics (e.g., EL) than their timing (e.g., LOS). Figure 6. Importance of phenological metrics extracted from the EVI time series. Figure 6. Importance of phenological metrics extracted from the EVI time series.

Topographic Effects on Map Accuracy
The correlation between the overall map accuracy and the change in slope or aspect was moderate, in general. However, the classification accuracies showed a decreasing trend as the slope gradient increased (Figure 7). Accuracies were the highest for crops on flat lands or slopes with gradients less than 2 • . The results also indicated a slightly higher accuracy for crops on north-facing than on south-facing slopes (Figure 7), although the difference was not significant.
The cultivation areas with relatively flat terrain showed high accuracies in crop monitoring. This was mainly due to the low level of spatial fragmentation and small variation in elevation change, thus reducing spectral variation and the probability of spectral mixing. The relationship between classification accuracy and slope for individual crop types is shown in Figure 8. The accuracies were calculated following the accuracy assessment steps described in Section 4.2, and the slopes were averaged for the areas covered by individual crop types. Plant structure, cultivation patterns, and topography were found to exert joint effects on mapping. For example, the low-growing crops (e.g., tomatoes, grapes, and mulberries) are mainly cultivated in river valleys and basins. While background soil adds noises to EVI data, the flat agricultural fields may have contributed to the high classification accuracies. For some tree plantations cultivated on slopes with gradients higher than 15 degrees (e.g., tea-oil camellia, gum, and tea trees), the modeling performance remained relatively good, with accuracies greater than 80% (Figure 8). The topographic effects may have been offset by the plantation uniform cultivation patterns.

Topographic Effects on Map Accuracy
The correlation between the overall map accuracy and the change in slope or aspect was moderate, in general. However, the classification accuracies showed a decreasing trend as the slope gradient increased (Figure 7). Accuracies were the highest for crops on flat lands or slopes with gradients less than 2°. The results also indicated a slightly higher accuracy for crops on north-facing than on south-facing slopes (Figure 7), although the difference was not significant. The cultivation areas with relatively flat terrain showed high accuracies in crop monitoring. This was mainly due to the low level of spatial fragmentation and small variation in elevation change, thus reducing spectral variation and the probability of spectral mixing. The relationship between classification accuracy and slope for individual crop types is shown in Figure 8. The accuracies were calculated following the accuracy assessment steps described in Section 4.2, and the slopes were averaged for the areas covered by individual crop types. Plant structure, cultivation patterns, and topography were found to exert joint effects on mapping. For example, the low-growing crops (e.g., tomatoes, grapes, and mulberries) are mainly cultivated in river valleys and basins. While background soil adds noises to EVI data, the flat agricultural fields may have contributed to the high classification accuracies. For some tree plantations cultivated on slopes with gradients higher than 15 degrees (e.g., tea-oil camellia, gum, and tea trees), the modeling performance remained relatively good, with accuracies greater than 80% (Figure 8). The topographic effects may have been offset by the plantation uniform cultivation patterns.   The cultivation areas with relatively flat terrain showed high accuracies in crop monitoring. This was mainly due to the low level of spatial fragmentation and small variation in elevation change, thus reducing spectral variation and the probability of spectral mixing. The relationship between classification accuracy and slope for individual crop types is shown in Figure 8. The accuracies were calculated following the accuracy assessment steps described in Section 4.2, and the slopes were averaged for the areas covered by individual crop types. Plant structure, cultivation patterns, and topography were found to exert joint effects on mapping. For example, the low-growing crops (e.g., tomatoes, grapes, and mulberries) are mainly cultivated in river valleys and basins. While background soil adds noises to EVI data, the flat agricultural fields may have contributed to the high classification accuracies. For some tree plantations cultivated on slopes with gradients higher than 15 degrees (e.g., tea-oil camellia, gum, and tea trees), the modeling performance remained relatively good, with accuracies greater than 80% (Figure 8). The topographic effects may have been offset by the plantation uniform cultivation patterns.

Conclusions
This study investigates the feasibility of using MODIS satellite time series and phenology-based modeling to extract multiple crop types in a subtropical mountainous region in South China. Our results show that, (i) crop classification using EVI time series outperformed NDVI time series for the majority of the 12 evaluated crop types. (ii) The application of RF to classify EVI-derived phenological metrics achieved an overall accuracy at 85.27%. (iii) The crops' spectral metrics showed slightly higher importance than those of their temporal counterpart, although both types of phenological metrics made significant contributions to crop classification. (iv) Slope revealed a generally negative impact on model accuracy. However, uniformly distributed tree plantations (e.g., tea-oil camellia, gum, and tea trees) that were cultivated on large slopes (>15 degrees) were found to achieve accuracies greater than 80%.