Mapping the Complex Crop Rotation Systems in Southern China Considering Cropping Intensity, Crop Diversity, and Their Seasonal Dynamics

Crop rotation increases crop yield, improves soil health, and reduces plant disease. Mapping crop rotation is difficult because crop data from a single time point do not sufficiently represent the dynamics of a system. Studies have tried to map crop rotation by sequentially combining crop maps. However, this produced a large number of meaningless crop sequences, hindering the assessment of rotational benefits at regional scales. Here, we propose a crop rotation classification scheme that integrates temporal information into static crop maps and develop an innovative approach to map crop rotation. We chose a typical multiple cropping region in southern China. Given that the landscape is characterized by high crop diversity (e.g., food crops, cash crops, and permanent crops) and variable cropping intensity, our classification scheme first defines three main rotation systems, i.e., paddy, vegetable, and orchard systems, and then further divides the systems into nine subsystems according to their seasonal dynamics. Finally, we apply time series of Sentinel-1 and Sentinel-2 images to identify the systems by a hierarchical rule-based method. The map of crop rotation systems in 2020 had producer, user, and overall accuracies of 81%, 79%, and 81%, respectively. The results indicate that integrating temporal information into the classification scheme is vital to representing complex rotation systems and that remotely sensed temporal dynamics of crops are useful to characterize these systems. It also shows that crop rotation can be mapped directly rather than aggregating multiple crop layers, thus providing a new perspective for mapping and understanding crop rotation systems.


I. INTRODUCTION
C ROP rotation is the practice of growing a series of crop types on the same land across consecutive seasons or years [1]. It has been found to increase crop yield, improve soil health, and reduce plant disease [2], [3], [4], suggesting that crop rotation could be considered an option for ensuring food production without compromising the environment. Crop rotation practices can be found all over the world. Amidst concerns over the impacts of increased chemical inputs on soil and water quality, the United States is encouraging more crop rotation [5]. The Common Agricultural Policy in the European Union also promotes legume-based rotations [6]. In China, crop rotation has traditionally been adopted by smallholder farmers to maintain soil fertility [7].
The adoption of crop rotation provides a variety of benefits, which have been widely reported by field studies conducted in the experimental plots and farms. For example, it was shown that growing velvet bean in rotation with maize was effective in increasing maize yield and improving soil quality [8]. Despite such evidence, crop rotation benefits have not been systematically investigated at regional or larger scales due to the lack of spatial crop rotation data. In developed countries, farm-level cropping decisions might be well documented for subsidy or insurance purposes (e.g., the land parcel identification system [9]), and these data could be used for deriving crop rotation maps [10]. However, farm-level data are scarcely available in developing countries where smallholder farms are predominant. For example, China has a long history of crop rotation practices, and the agricultural landscapes in China are characterized by small fields, high cropping intensity, and diverse crop types [11], [12]. The fragmented agricultural landscapes and the lack of advanced data acquisition channels lead to limited documentation on cropping decisions across Chinese farms.
To derive the geographical distribution of crop rotation patterns, both modeling and remote sensing techniques have been explored [13], [14]. For modeling studies, a two-phase procedure is usually followed. The first step generates possible crop sequences based on crop type information from statistics combined with some auxiliary data (e.g., agronomic criteria, expert knowledge) [15], [16]. And in the second step, the generated crop sequences are spatially allocated according to specific rules to obtain the spatial patterns [17], [18]. However, modeling research articles often lack actual data on crop sequences and, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ thus, need to rely on the analyst defining realistic crop sequences, which may not account for the large variations in farming practices across a region [19].
With the development of remote sensing technology and data availability, remote sensing emerges as a way to map the actually practiced crop rotations [20], [21], [22], [23], [24], [25]. Because crop rotation involves crops grown at different seasons or years, crop rotation mapping is more challenging than crop mapping, and multidate remote sensing images are required. There are usually three steps in crop rotation mapping. First, crops are separately mapped for different seasons or years. Then, crop sequences are generated by sequentially combining multiple crop maps. Finally, representative crop rotations are identified from crop sequences based on the expert knowledge or specific rules. For example, Martinez-Casasnovas and Martin-Montero [21] characterized typical crop rotations bases on the temporal series analysis of crop maps from 1993 to 2000 generated by supervised classifications of Landsat TM images [21]. And Waldhoff et al. [25] obtained diverse crop sequences based on the combined annual crop maps for eight consecutive years and identified dozens of representative crop rotations [25]. Several studies have mapped crop rotation based on multiyear analysis of existing earth observation data, e.g., cropland data layer [23]. Recently, Li et al. [26] reported a phenology-based method to map crop species and rotation types using fused MODIS and Landsat time-series data [26].
However, as most crop maps are not produced for tracking temporal changes, combining these independently produced maps may result in many meaningless sequences due to "salt and pepper" pixels and other mapping errors [24], [27]. These sequences can rarely represent farmers' rotation decisions, limiting their decision-supporting role in agricultural management. In addition, most studies focus on single-cropping regions, while the crop rotation patterns in multicropping areas (e.g., southern China) remain unclear, which hinders the understanding and planning of cropland use. To address the above limitations, we present a case study in Zengcheng district, a typical multicropping region in southern China. And this article would contribute to the pieces of literature in two major respects. 1) We propose a crop rotation classification scheme that integrates temporal information into static crop maps to represent the intra-annual variability of cropland use in Zengcheng. The classification scheme accounts for information on cropland intensity as well as crop diversity, enabling us to understand the dynamics and complexity of local cropping systems. 2) We develop a hierarchical rule-based method, offering an automatic and efficient way to map the crop rotation systems. It is designed to directly map all the systems to reduce error propagation.

A. Study Region
The case study region is Zengcheng district, Guangdong Province, southern China (see Fig. 1). The topography is mainly characterized by hills in the north and plains in the south, and  the agricultural landscapes are largely fragmented. The district has a southern subtropical monsoon climate with abundant sunshine and rainfall, which allows crops to be grown all year round, resulting in multiple cropping seasons. Vegetables are widely distributed and intensively managed, and Zengcheng is an important production base in the Vegetable Basket Project in the Guangdong-Hong Kong-Macao Greater Bay Area (GBA). Fruits and rice are the other two major crops. There are also many minor crops, e.g., peanut, corn, and sweet potato. Multiple cropping seasons and diverse crops together lead to highly complex crop rotation systems (see Fig. 2).

B. Remote Sensing Data
Sentinel-1 (S1) images with a 12-day interval from the Level-1 ground range detected product and Sentinel-2 (S2) images with a 5-day interval from the Level-2A surface reflectance product archived in the google earth engine (GEE) platform were used in our study [28]. All the available satellite images acquired during the period from December 1st 2019 to November 30th 2020 were exploited for classification, which completely covers the cropping cycles of typical rotations in Zengcheng.
Specifically, we obtained 289 S2A/B images and 30 S1A images. S1 carries a dual-polarization C-band synthetic aperture radar (SAR) instrument, resulting in vertical transmission/horizontal receiving (VH) and vertical transmission/vertical receiving (VV) bands. S2 carries a multispectral imager sensor that provides high-resolution multispectral images. S1 images have been used to explore the spatiotemporal patterns of paddy rice in cloudy regions due to the ability of the sensor to penetrate clouds [29]. It was previously shown that including S1 data in a study resulted in higher rice mapping accuracy than using data from only a single sensor (i.e., S1 or S2) [30]. Therefore, S1 images were included to improve the accuracy of paddy rice extractions. S2 images have shown promise in crop mapping over fragmented croplands [31], [32], suggesting their potential for crop rotation mapping in southern China. S2 images were gap filled and smoothed to construct the seamless five-day time-series data. More information can be found in Supplementary text Fig. S1 and Table S1.

C. Cropland Mask
The global land-cover product with a fine classification system at 30 m for 2020 (GLC_FCS30-2020) produced by the Chinese Academy of Sciences was applied to extract the cropland extent [33]. It was produced based on GLC_FCS30-2015, Landsat data, S1SAR data, and other auxiliary datasets. The classification system inherited that of the European Space Agency Climate Change Initiative Global Land Cover and cropland comprised rainfed and irritated cropland in which rainfed cropland included herbaceous cover and tree or shrub cover (orchard).
Kang et al. [34] indicated that GLC_FCS-2015 (especially cropland) had higher accuracy than GlobeLand30-2010 and the finer resolution observation and monitoring of global landcover product in 2015 over a complicated tropical region [34]. We further compared the cropland mapping accuracies among GLC_FCS-2020, GlobeLand30-2020, and FROM_GLC-2017 based on 300 cropland samples (randomly selected from ground reference data) and 300 noncropland samples (visually interpreted on google earth). The results revealed that GLC_FCS-2020 and GlobeLand30-2020 had higher overall accuracies. Given that GLC_FCS-2020 had a great advantage in spatial details (see Fig. S2), it was finally used in our study.

D. Ground Reference Data
In September 2020, we conducted a field survey in Zengcheng, which involved collecting not only ground samples by recording coordinates and taking pictures but also information from farmers who are the managers of land parcels. In addition to visible land management decisions that were recently made in the field, we asked farmers who were working in the surveyed fields to retrospectively review their land management decisions from the last winter to the current season (December 2019 to September 2020). We also asked them to share their land management plans for the next few months (September 2020 to November 2020). This information about land management decisions and plans was combined to generate crop sequences. Ultimately, we collected data from 738 field points. In the following sections, these samples will be divided into different groups for certain purposes (see Table I). To train the random forest (RF) classifier, 80 samples were used (40 for vegetables and 40 for orchard). We discarded 34 samples due to poor quality, and the remaining 624 samples were exploited for assessing the classification performance of the entire hybrid system (80 + 34 + 624 = 738). The validation sample sizes of rice, vegetables, and fruit trees are 223, 212, and 189, respectively. In addition, we determined the important growing stages of the major crops (see Fig. 3). More information about the field survey can be found in Supplementary Text Figs. S3 and S4.

E. Agricultural Statistics
The Guangzhou Statistics Bureau publishes annual reports on the harvested area of major crops in each district (https: //lwzb.gzstats.gov.cn:20001/datav/admin/home/www_nj/). We acquired the 2021 Guangzhou statistical yearbook that reported the harvested area of crops in 2020.

III. METHODOLOGY
We first provide a crop rotation classification scheme (see Fig. 4) based on the field survey, farmer interviews, expert knowledge, and agricultural statistics to represent the complex crop rotations. The cropping systems in the study region are characterized by both inter-and intra-annual variabilities. Here, we focus on the intra-annual variability rather than the interannual variability because the intra-annual variability is needed to understand the more complex inter-annual variability. Based on the established classification scheme, we further propose a hierarchical rule-based method to generate a 10-m map of crop rotation systems for 2020 (see Fig. 5). Here, the year 2020 refers to the period from December 1st 2019 to November 30th 2020.    Fig. 2. "Rice-vegetables" was the rotation most frequently mentioned by the experts, which shows that the rotations of vegetables with rice are very common (see Table II). In addition, vegetables are intensively managed, with "vegetables-vegetables-vegetables" and "vegetables-vegetables" frequently mentioned.

2) Crop Rotation Classification Scheme:
Based on the field survey, farmer interviews, expert knowledge, and agricultural statistics, we provide a classification scheme that divides crop rotation into three main systems and nine subsystems (see Fig. 4). Rice is the primary grain crop, while vegetables and fruits are the primary cash crops. Three main systems were, thus, identified: paddy, vegetable, and orchard systems. The distinctive management style of fruits is one reason to include orchards in our classification [35]. Other crops, such as maize and peanuts, are not included due to their small planting area (less than 5% total).
We separated subsystems from the three main systems to highlight differences in each group. The paddy system was divided into four subsystems. Single rice and double rice were first defined and were further distinguished by the presence or absence of vegetables, given that it is quite common to rotate vegetables with rice. The vegetable system was classified into three subsystems in which high-intensity and low-intensity vegetables were first separated. Information on cropping intensity can help us better understand cropland use intensity. During the field survey, we found that some farmers tend to plant leafy vegetables year around, while others grow different types of vegetables, indicating that temporal crop diversity varies among vegetable fields. Considering that crop rotational diversity has large impacts on agricultural systems, integrating crop diversity might be useful to understand the characteristics of crop rotation systems. Therefore, high-diversity and low-diversity vegetables were further specified among the high-intensity vegetables to demonstrate differences in temporal crop diversity. For the orchard system, short-term and long-term orchards referring to herbaceous (e.g., banana, papaya) and arbor (e.g., litchi, longan) fruit trees were defined. Short-term orchard usually involves variant management, while long-term orchard involves more stable management. For example, the filed survey indicates that litchi typically extends over several decades with stable management, while banana would be replaced by other crops after growing for three or four years.

B. Mapping Crop Rotation Systems
According to the presented classification scheme, we propose a hierarchical rule-based method (see Fig. 5) to map the crop rotation systems. We initially obtained four key indicators: flooding frequency, cropping intensity, cropping diversity, and coefficient of variation. Then, we mapped the three main crop rotation systems according to the S2 images and the RF algorithm. Finally, the nine crop rotation subsystems were identified by the four indicators. The cropland layer, as described in Section II-B2, was used to mask the nonagricultural classes. All the analyses were performed within the GEE cloud-based public imagery repository and high-performance computing system.

1) Acquisition of the Four Indicators:
We obtained the four indicators by integrating the time series of S1 and S2 images. The temporal frequency of S1 images is 12 days. Here, the S2 time series refer to the smoothed enhanced vegetation index (EVI) and land surface water index (LSWI) with a five-day interval, which were produced by gap filling and smoothing the cloud-masked S2 images. The gaps caused by cloud cover were  filled by time-series linear interpolation based on good-quality observations before and after the time step. Then, the Whittaker algorithm was leveraged to smooth the time series [36]. Here, we chose the Whittaker algorithm because it can effectively balance the fidelity and roughness of curves [37] and is computationally efficient and can be conveniently implemented on GEE [38]. Prior to the smoothing, the key parameter λ quantifying the degree of smoothness is determined by the multiple linear regression, as described in [38]. To this end, we randomly selected 50 sample points, and calculated the value of λ for each point based on the cloud-masked EVI and LSWI time series. At last, the mean values were used to determine the final λ for EVI and LSWI. The smoothed time series was shown in Figs. 6 and 7. See Supplementary text for more information.
a) Flooding frequency: This indicator represents the annual frequency of the flooding and transplanting signals that have been widely used to indicate the presence of rice [39]. The flooding and transplanting signals were detected separately for early, middle, and late rice, and then the flooding frequency was computed by overlaying the three layers. We applied data from S1 and S2 images to detect the signals, and the S1-and S2-based outputs were merged separately for early, middle, and late rice. For the S1 images, we used the VH backscatter and method described in [40], which is easy to implement on GEE. For the S2 images, the relationship between the EVI and LSWI, specifically that LSWI values are temporarily greater than EVI values during the flooding and transplanting phases, was used [41]. See Supplementary Text and Fig. S5 for more information.
b) Cropping intensity: This indicator is defined as the number of cropping cycle(s) [42]. We computed cropping intensity from the smoothed EVI time series by counting the number of peaks (see Fig. 6) with a quadratic difference algorithm [43]. To reduce pseudopeaks caused by abnormal EVI values, some constraint conditions were used. Specifically, to filter false cycles, the temporal interval of peaks was required to be greater than 48 days [44]. In addition, the peak values were not allowed to be less than 0.35 [45].
c) Cropping diversity: Here, we define cropping diversity as crop diversity over time and measure it by quantifying differences among crops. Phenological metrics can efficiently capture the phenological characteristics of crops and have been widely applied to distinguish crops [46]. Thus, we applied phenological metrics to estimate cropping diversity. We selected three commonly used phenological metrics: maximum value, growth duration, and seasonal amplitude [47], [48], [49]. We first retrieved the three crop-specific metrics from the smoothed EVI time series (see Fig. 6). Then, we calculated the coefficient of variation for each phenological metric, which was averaged to represent cropping diversity. The coefficient of variation has been proven to be an effective indicator for observing spectral diversity over space to represent plant biodiversity [50], [51].
Specifically, we first obtained the maximum value to represent the peak value of each crop's growth cycle. Then, we estimated the seasonal amplitude by referring to the difference between the maximum value and base value where A i is the seasonal amplitude for crop i in a crop rotation system; EVI max,i represents the maximum value for crop i; and EVI base,i represents the base value for crop i and is extracted as the average of the left and right minimum values for each crop's growth cycle.
In the next step, growth duration was acquired by referring to the interval between the start of seasons (SOS) and the end of seasons (EOS). The SOS and EOS were identified as the day of year when the EVI reached a value using the EVI ratio method where EVI ratio is the EVI ratio; EVI min represents the EVI value of bare soil, which is calculated as the minimum EVI value over one year; and EVI max represents the crop-specific maximum  [52].
Next, we calculated the coefficient of variation for each metric. Finally, cropping diversity was calculated by referring to the average of the coefficient of variation for the maximum value, growth duration, and seasonal amplitude. Larger values reflect higher crop diversity.
where CD is the cropping diversity; CV M is the coefficient of variation for the crop-specific maximum value; CV GD is the coefficient of variation for the crop-specific growth duration; and CV A is the coefficient of variation for the seasonal amplitude. d) Coefficient of variation: This indicator refers to the coefficient of variation in the EVI and is calculated from the smoothed EVI time series to quantify the annual variability in the EVI or the extent of vegetation fluctuation relative to the mean. It was previously applied to separate cropping systems with a lower degree of seasonality (e.g., permanent crops) from other cropping systems [44]. Long-and short-term orchards show different degrees of seasonality [see Fig. 7(d)], and long-term orchards have a lower degree of seasonality than short-term orchards. Therefore, we introduced the coefficient of variation in the EVI to exploit its potential for identifying orchard systems.

2) Identification of Crop Rotation Systems: a) Identification of the main systems:
The three main crop rotation systems were mapped prior to the identification of the subsystems. Compared with vegetables and fruits, rice has more obvious features, and thus, we first distinguished between paddy and nonpaddy systems. Pixels were identified as paddy systems as long as one season of rice was observed, i.e., the flooding frequency was nonzero.
Next, we removed paddy pixels from the cropland extent and adopted RF to differentiate vegetable and orchard systems. RF is more robust than a single decision tree and typically achieves high accuracies in crop mapping [47]. As shown in Table III, we mainly used three types of features for classification: temporal, phenological, and texture metrics. The temporal metrics were composed of seasonal median composites of reflectance bands and spectral indices (see Table S2), which capture seasonal variations in cropland surface spectra.
Pixel-based compositing allows the exploitation of all available images and contributes to overcome the limitations of satellite data quality [53]. Monthly compositing may better capture the cropland dynamics, yet it would make the data volume too large. And we found that monthly compositing might not be appropriate for Zengcheng because incomplete composites would be generated for many months because of severe cloud contamination.
There are 11 spectral indices, including 5 red-edge spectral indices. The red-edge spectral indices have shown great potential to differentiate crops [53], and thus, several commonly used metrics were selected and included in the feature candidates. The phenological metrics were obtained from the smoothed EVI time series to depict spectral means and variances. The texture metrics were calculated from the annual median composite of EVI by a gray-level cooccurrence matrix (GLCM) method in GEE. The implementation of RF was run on GEE. We applied a backward feature elimination method to select the best variable set for RF by progressively eliminating the least important variables and fine-tuned the main hyperparameters (i.e., number of trees, number of variables per split, and minimum size of terminal nodes). More information can be found in Supplementary text and Figs. S6-S10. b) Identification of the subsystems: As shown in Fig. 5, we identified the nine crop rotation subsystems by introducing different indicators. For the paddy subsystems, a two-phase identification procedure was followed (4). The first step distinguished between single and double rice, which differed in rice planting frequency [see Fig. 7(a)]. Pixels of paddy system with a flooding frequency of one were identified as single rice, while pixels with a flooding frequency of two were identified as double rice. In the second step, the presence or absence of vegetables was detected by the relationship between flooding frequency and cropping intensity [see Fig. 7(b)]. Cropping intensity indicates the number of crops harvested in a given period, including rice, while flooding frequency refers to the number of rice harvested. Therefore, when cropping intensity is greater than flooding frequency, it could be inferred that at least there is another harvest of crop besides rice. When the two values are equal, it means that rice is the only crop harvested. Cropping intensity could be underestimated due to cloud cover, so in the case where it is abnormally smaller than flooding frequency, we assume that the two values are equal to reduce the effect of error propagation.
Based on the above logic, we classified pixels as single/double rice rotated with vegetables when cropping intensity was greater than flooding frequency; otherwise, the pixels remained single/double rice. For example, when cropping intensity is 3 and flooding frequency is 2, this pixel would be classified as double rice rotated with vegetables.
where FF represents the flooding frequency, and CI represents the cropping intensity.
For the vegetable subsystems, we first classified high-and low-intensity vegetables and further identified high-and lowdiversity vegetables from high-intensity vegetables. The number of cropping cycles is different between high-and lowintensity vegetables [see Fig. 7(c)], which can be distinguished by cropping intensity. Empirically, we determined pixels of vegetable systems with cropping intensity greater than one as high-intensity vegetables and those with cropping intensity equal to one as low-intensity vegetables. Cropping intensity equal to one indicates that there is only one harvest of vegetables, which aligns well with the concept of "low intensity." Obviously, this threshold could be adjusted flexibly based on regional conditions, and the adjustment would cause changes in the planting area of high-and low-intensity vegetables. Highand low-diversity vegetables exhibited different seasonal dynamics [see Fig. 7(c)], and we leveraged cropping diversity to distinguish them.
Long-and short-term orchards show different degrees of seasonality, allowing the differentiation of the orchard subsystems [see Fig. 7(d)]. Long-term orchards have a lower degree of seasonality than short-term orchards. Therefore, pixels identified as orchard systems with a higher coefficient of variation were identified as short-term orchards; otherwise, they were identified as long-term orchards.
The thresholds of flooding frequency and cropping intensity were set empirically based on the field survey, farmer interviews, and expert knowledge. For cropping diversity and the coefficient of variation, we use Otsu's thresholding algorithm, which was implemented on GEE, to determine the optimal thresholds. This method processes image histograms, segmenting objects by minimizing interclass variance [54]. It has been used in many applications as a satellite image thresholding method [55]. Otsu's thresholding algorithm automatically derives the optimal thresholds and, thus, can improve the universality of the developed algorithm. Ultimately, the thresholds of cropping diversity and coefficient of variation in the EVI were specified as 0.18 and 0.26, respectively.
3) Accuracy Assessment: We collected 738 samples in total from the field survey, 80 of which were used for training the RF classifier (40 for vegetables and 40 for orchard). We discarded 34 samples due to poor quality, and the remaining 624 samples were exploited for assessing the classification performance of the entire hybrid system (80 + 34 + 624 = 738). Specifically, we adopted the ground samples collected in the field survey to assess the accuracies of the four indicators and the crop rotation systems map, excluding the 80 training samples used for mapping vegetable and orchard systems. We first labeled the ground samples with crop rotation information as crop rotation systems. For example, "spinach-romaine lettuce-eggplant" was labeled as a high-diversity vegetable rotation because the pixel was triple cropped, and both leafy and fruiting vegetables were grown. Then, we labeled the ground samples without complete crop rotation information by interpreting the temporal signatures of the original and smoothed EVI and LSWI on GEE. The paddy and orchard systems are easier to recognize, while it is hard for vegetable systems, especially high-and low-diversity vegetables, which may introduce uncertainty. In this case, we have several remote sensing experts (listed as coauthors) who have check all points, and only well-interpreted points with a high level of confidence are kept, leading to 624 sample records. Based on those sample points, we assessed the accuracy of the identification based on four indicators: the overall accuracy (OA), user's accuracy (UA), producer's accuracy (PA), and F1 score (F1).
In addition, we estimated the harvested area of rice, vegetables, and fruits from the map of crop rotation systems and compared them with the statistically determined harvested area. The harvested area refers to the area from which a crop is gathered. If the crop under consideration is harvested more than once during a year, the area was counted as many times as the crop was harvested. It is worth noting that vegetables in paddy systems were included in the harvested area of vegetables.

A. Maps of the Four Key Indicators
We mapped the four indicators in Zengcheng for 2020. All four indicators exhibit large spatial variations (see Fig. 8). The pattern of flooding frequency shows that the northern part of the study region has more paddy fields, especially those with two rice harvests [see Fig. 8(a)]. Cropping intensity indicates that most croplands are intensively managed [see Fig. 8(b)]. A total of 40% of the pixels are identified as double cropping, followed by triple cropping (31%) and single cropping (18%). The pixels with a cropping intensity of four occupied nearly 10% of the total pixels. Fig. 8(c) suggests that the croplands with lower cropping diversity are mainly concentrated in the southern part of the study area, especially the southeastern part. For the coefficient of variation, the croplands with higher values are widely distributed in the study region, while the croplands with lower values are mainly located in the southern part of the study area [see Fig. 8(d)].

B. Map of the Crop Rotation Systems
Based on the crop rotation classification system and mapping methods, we mapped the crop rotation systems across Zengcheng in 2020 at a 10 m spatial resolution [see Fig. 9(a)]. Overall, the northern part of the study area is predominantly covered by paddy and orchard systems, with paddy systems spreading in the plains of the valley, while the flat southern part of the study area is primarily covered by orchard and vegetable systems. This result is consistent with the findings of the field survey, which show that the southern part of the study area is closer to markets and is easier to travel to, and hence, farmers there prefer vegetables and fruits with higher profits over rice, despite government rice subsidies. For paddy subsystems, the system consisting of single rice rotated with vegetables is widely distributed throughout the paddy regions. Vegetable subsystems are interspersed throughout the vegetable-growing regions. For the orchard subsystems, short-term orchards are prevalent in the northern part of the study area, where mainly bananas and papayas are grown, while long-term orchards are clustered in the central southern part of the study area, with guava and citrus being the major fruits. The orchard system accounts for approximately 64% of the total cropland area, followed by the vegetable (21%) and paddy (15%) systems. The paddy system is dominated by single rice rotated with vegetables, followed by double rice and double rice rotated with vegetables. Among the vegetable systems, high-intensity vegetable systems are predominant, with low-diversity vegetable systems having slightly larger areas than high-diversity vegetable systems. Within the orchard systems, short-term orchards are predominant over longterm orchards. We further summed the harvested area of the crop rotation systems [see Fig. 9(d) and (e)]. Compared with the proportion of cropland area, the proportion of vegetable systems was nearly double at 43%, confirming the extremely high intensity of local vegetable production.

C. Accuracy Assessment
As the four indicators are essential for the research, we first evaluated the accuracies of the four indicators. The validation results showed that flooding frequency and coefficient of variation had higher overall accuracies than cropping intensity and cropping diversity (see Table S7). Then, we evaluated the accuracy of the map of crop rotation systems (see Table IV). The producer, user, and overall accuracies are 81%, 79%, and 81%, respectively. The vegetable system has a lower accuracy (75%) than the orchard (87%) and paddy (80%) systems. The paddy system has more errors of commission than errors of omission (PA > UA). The higher commission errors are partly due to the misclassification of vegetables and fruits as paddy systems. Moreover, there are misclassifications among the paddy subsystems. The misclassification of single rice as double rice subsystems is comparatively severe. Among the vegetable systems, high-and low-diversity vegetables are easily confused, which is the main source of errors resulting in lower F1 scores.
The harvested areas of rice, vegetables, and fruits estimated from the map were compared with the statistically determined harvested areas in 2020 (see Fig. 10). The estimated harvested area of rice is very close to the statistically determined area, with a slight difference. However, the harvested area of fruits is overestimated by 5%, while that of vegetables is underestimated by 5%.

A. Uncertainties
Multisourced data were used in our study and may introduce uncertainties. For example, the errors in the cropland layer are likely to be propagated to the final map [44]. The area overestimation of orchard systems (see Fig. 10) could be due to the inclusion of forests as cropland in GLC_FCS30-2020. The developed mapping method relies much on the time-series satellite data, which can be a source of uncertainty [53]. It was found that the accuracies of different rice types were related to the satellite data availability, and higher accuracies came from more valid observations [30]. The four key indicators estimated with the time-series images are more likely to be affected, which could cause direct impacts on the final mapping accuracy, e.g., an incomplete counting of cropping cycles could underestimate the harvested crop area [42].
Specifically, the mapping accuracies of paddy subsystems are jointly controlled by flooding frequency and cropping intensity. Flooding frequency first determines whether paddy cropland is precisely mapped and then determines whether single and double rice are correctly distinguished from paddy cropland. Afterward, flooding frequency and cropping intensity together determine whether vegetables exist in single and double rice. The mapping accuracies of vegetable subsystems are jointly controlled by flooding frequency, cropping intensity, and   cropping diversity. Flooding frequency first determines whether nonpaddy cropland is precisely mapped. Then, cropping intensity determines whether high-and low-intensity vegetables are correctly distinguished. Afterward, cropping diversity determines whether high-and low-diversity vegetables are correctly distinguished. The mapping accuracies of orchard subsystems are jointly controlled by flooding frequency and coefficient of variation. Flooding frequency first determines whether nonpaddy cropland is precisely mapped, and then the coefficient of variation determines whether long-and short-term orchards are correctly distinguished.
Therefore, we measured the coverage of good-quality observations for each acquisition day (see Fig. 11) and counted the numbers of good-quality observations for individual pixels (see Fig. 12) to explore this potential impact of including S2 data on the results. Temporally, the coverage of good-quality observations was lower during the period from May to September. This would create additional uncertainties in the four indicators. The use of sparser good-quality S2 images from February to April and from July to September covering the transplanting phase will involve risks of losing the flooding signals in paddy fields. Spatially, the quality of the S2 data is better in the northern part of the study area than in the southern part, and the average good-quality observation frequency is 26.8 for the entire year, which is sufficient to identify crop characteristics.
We also found that the higher commission errors for the paddy system could arise from the S1-based ARM-SARFS algorithm despite the reduction in omission errors. This could induce false flooding and transplanting signals due to the presence of diverse vegetables and fruits, suggesting that the application of ARM-SARFS in complex agricultural landscapes is slightly In addition, the training sample size might influence the binary RF classifier's performance [56], thereby affecting the final results. Therefore, here we carried an experiment to examine how training sample size (20-160 per class) affects classification performance of the binary RF classifier. To enable comparable results, a fixed set of stratified random sampling points was selected to validate the results. This validation set is about 30% of the total sample size of the two categories, comprising 144 samples (76 for vegetables and 68 for orchard). We found that the classification results were greatly affected by training sample size, showing a trend of higher accuracies with more samples (see Table V). But notably, when training sample size per class exceeded 40, the increase rate of accuracy became smaller. It indicated that satisfactory results could be obtained with an appropriate sample size, and too many samples may lead to overfitting of the results.

B. Improvements
The uncertainties identified in this study point to future research directions for further improvements. To reduce the uncertainties caused by earth observation data, images with high temporal and spatial resolutions should be considered. For example, PlanetScope will benefit crop rotation monitoring over rapidly changing cropland surfaces [57]. To improve the accuracy of flooding frequency, supervised classification methods that do not rely on dense images could be considered. A comparison between phenology-based approaches and supervised classification methods could be conducted to determine which are more reliable for mapping rice over complex agricultural landscapes.
In addition, better frameworks to retrieve cropping intensity need to be adopted [44], [52]. The suitability and efficiency of multiple phenological metrics for identifying temporal crop diversity could be evaluated to improve the reliability of the proposed cropping diversity estimation algorithm.

C. Implications
Earth observation provides an effective way to make crop maps [14], [47]. And a few prior studies have managed to map crop rotation by combining sequential crop maps produced for specific seasons or years [20], [21], [22], [23], [24], [25]. This approach could be easily upscaled and extended to other regions where continuous crop monitoring is possible. However, it may not be suitable in regions with intensified agriculture and diversified rotation because identifying crops for each season is challenging. And combining multiple crop maps could produce many sequential crop combinations, which are hard to interpret [24].
Our study provides a timely solution to the issue of crop rotation mapping over complex regions. Especially, the proposed hierarchical rule-based pipeline is vital to the good classification performance, which has reference value for mapping cropping systems in other regions. It breaks the classification task into a two-step process, simplifying the complex issues. In this case, the RF classifier only needs to distinguish two categories, and thus, its efficiency and accuracy can be better guaranteed. To justify this perspective, we carried an experiment in which RF was applied to classify all the nine subsystems directly, including the four key indicators as features. To enable comparable results from the two approaches, we generated a new validation dataset containing 211 samples (about 30% of the total sample size) by stratified random sampling. The old classification result was validated again based on this new dataset. All the remaining samples were used for training the pure classifier-based approach to ensure good performance. The results showed that the overall accuracy reduced by 0.17 compared with the hierarchical rule-based method. The two double rice systems and orchard systems had acceptable results. However, the pure classifier-based approach appeared hard to distinguish between the three vegetable systems and the two single rice systems. This might be because the crop rotation systems were defined manually to represent farmers' diverse decisions, and machine learning methods, such as RF, have limitations when it comes to the highly complex man-made systems. Based on our model, the binary RF classifier (vegetables or orchard) in combination with the manually assigned subcategorization rules could well identify the complex crop rotation systems.
The proposed crop rotation classification system and mapping framework could provide implications for mapping crop rotation in other regions. First, crop types have been the basic units for crop rotation mapping [20], [21], [22], [23], [24], [25], while we update the classification scheme from crop types to crop rotation systems. Thus, the crop rotations could be mapped directly without identifying crops for each season, combining crop maps, and interpreting the confusing sequential crop combinations [24]. And because the systems are predefined based on the actual farming situations, the resultant map could have more important practical benefits. Second, previous studies only focused on the crop sequences themselves [58], while this study suggests that temporal information (e.g., temporal crop diversity and cropping intensity) could be integrated into crop rotation maps, which can produce more informative crop rotation maps with implications for multiple objectives, e.g., crop diversity analysis [59], [60] and soil organic carbon mapping [61], [62]. This may be more consistent with the concept of land system science [63]. In addition, several novel ideas proposed in our study might have some reference values. For example, the algorithm for measuring temporal crop diversity by phenological metrics may be leveraged to estimate spatial crop diversity efficiently without detailed information on crop types and their distribution [64], [65].

VI. CONCLUSION
We developed an innovative approach for mapping highly diversified crop rotation systems based on earth observation data and applied it to a case study in southern China. The resulting map suggests that the local vegetable production is characterized by high intensity and low diversity. The results of this study indicate that an appropriate classification scheme that integrates seasonal information is vital to understanding complex crop rotation systems. They also show that crop rotation can be mapped directly rather than aggregating multiple seasonal crop layers, providing a new perspective for mapping and understanding crop rotation.