Patterns and causes of winter wheat and summer maize rotation area change over the North China Plain

Cropland area and cropping frequency play very crucial roles in determining regional food production. However, rapid urbanization accompanied by declining surplus-agricultural labor force has greatly altered patterns of agriculture land use and cropping frequency. Due to lack of continuous cropland and cropland-use intensity maps, our knowledge is still limited to understand whether the urbanization process must have a negative effect for changes in cropland-use intensity. Herein, we took the North China Plain (NCP), both the largest winter wheat and summer maize rotation area and rapidly urbanized area in China, as the study area, and used 250 m moderate resolution imaging spectroradiometer NDVI anomaly data, the correlation of NDVI time series in two neighboring years and machine learning algorithms to investigate spatiotemporal patterns and trends of cropland area and cropping frequency change over the NCP from 2000 to 2019. Results showed a significantly decreased cropland area observed since 2004 (slope = 783.8 km2 a−1, p < 0.01), while area of double-season cropping presented a relatively steady trend (slope = 446.9 km2 a−1, p = 0.335). As expected, decreased croplands were mainly occupied by urban and built-up land expansion, however, existing cropland-use intensity was yet improved. Patterns and trends of double-season cropping types were varied spatially. Particularly, the area of winter wheat and summer maize rotation presented a significantly increasing trend (slope = 3423.3 km2 a−1, p < 0.01). Furthermore, the respective area of winter wheat and summer maize both displayed significantly increasing trends with slope of 2953.8 and 2874.9 km2 a −1(p < 0.01) in entire period. Land-use and grain subsidy policies are considered as largely responsible for this phenomenon. These satellite-observed findings highlight that positive land-use policies and managements will be helpful for profitably keeping/improving the harvest area.


Introduction
Population growth and rapid socio-economic development pose a tremendous challenge for global agricultural systems and food supply (Godfray et al 2010, Bryan et al 2018, Liu 2018, Abu Hatab et al 2019. Improving cropland-use intensity and expanding cropland area are regarded as the traditional and widely-used means to maximize crop yields in most regions of the globe (Lobell et al 2008, Ray and Foley 2013, Spera et al 2014, Cohn et al 2016. Double cropping plays a crucial role in meeting food demands of the temperate zones and it is also an important evaluation indicator of land-use intensification. However, rapid urbanization associated with declining surplusagricultural labor force has greatly altered patterns of agriculture land use and cropping frequency. This process was usually regarded as having a negative effect for changes of cropland area and/or croplanduse intensity (Cui et al 2019). For example, some rural household-based surveys found that some farmers passively or actively give up farming due to rapid urbanization and sagging grain prices (Spera et al 2014, Liu et al 2014b, Liu 2018. Croplands in some regions were readily swallowed by urban and built-up lands (Liu et al 2014a, Ning et al 2018 or appeared abnormally with low-intensity use, such as singleseason or double-season agricultural land abandonment instead of intensive cultivation (Prishchepov et al 2012, Spera et al 2014. Yet, it is worth noting that lack of available large-area maps could largely limit our systematic understanding for changes of cropland area and cropping frequency in response to rapid urbanization. Therefore, lots of communities, including remote sensing, agriculture, land system science and climate change, have paid more attention to studies of cropland area and cropping frequency changes in recent years (Ray and Foley 2013, Spera et al 2014, Cohn et al 2016, Wu et al 2018. Satellite-based remote sensing data have provided great opportunities to monitor large-area changes of land-use types, given their features of open access, continuous time series and wide-coverage observations (Mueller-Warrant et al 2016, Leroux et al 2017, Waldhoff et al 2017, Lasko et al 2018, Weiss et al 2020. Among the satellite data, 30 m Landsat, 10 m Sentinel and 250 m moderate resolution imaging spectroradiometer (MODIS) are the most widely applied for the related studies of crop mapping and monitoring (Wardlow et al 2007, Spera et al 2014, Dong et al 2015, Veloso et al 2017, Hao et al 2018. In comparison with the data of medium-high resolution Landsat and Sentinel satellites, MODIS data have a slightly coarse spatial resolution but have much better-quality time series information owing to relatively weaker impacts by the noise, e.g. cloud contamination. By contrast, MODIS data thus have greater potentials to map croplands and cropping frequency in view of the importance of time series information for agriculture (Wardlow et al 2007, Cohn et al 2016. Some studies employed 250 m MODIS time series data to investigate the changes of cropping frequency in Mato Grosso state of Brazil and they also provided very valuable information for us to map cropland area and crop types by using time series data (Spera et al 2014, Cohn et al 2016. The automatic and continuous mapping, by directly using time series vegetation index, could be effective in relatively shortterm mapping. Moreover, CO 2 fertilization and land management practices reportedly increase the vegetation biomass (Liu et al 2020). Our preliminary tests also found that the above-mentioned effects cannot be ignored for the long-term continuous mapping of crop types. Fortunately, annual vegetation index anomaly could minimize these impacts of CO 2 fertilization and land management as they are usually calculated by using time series values subtracting averaged value each year. To do so, we can better reserve the temporal behavior of each crop, which is important for improving specific-crop classification. On the other hand, proper approaches are also very crucial for cropland-use and cropping frequency mapping. Previously, machine learning approaches and Twi-difference approach are usually used to map cropping frequency (Panigrahy and Sharma 1997, Wardlow and Egbert 2008, Yang et al 2013, Spera et al 2014. Compared to Twi-difference approach, machine learning approaches are more promising. When we map land-use types in a specific region, a core limitation is that there are rarely fine-resolution maps of entire study area available to assess the robustness of inferred land-use maps (Lunetta et al 2010). One of the methods being practiced by researchers to tackle this challenge is to use results from multiple machine learning algorithms (or termed as ensemble learning). Using the base map derived from multiple machine learning algorithms as well as the image change detection analysis method, we can produce the continuous maps.
In this study, our aims are: (a) to produce annual map of cropland and double cropping area across the North China Plain (NCP) by using machine learning algorithms and corresponding time series of MODIS-NDVI anomaly data, (b) to investigate the spatiotemporal dynamics of cropland and double cropping in the context of rapid urbanization, and (c) to clarify the reasons behind changes in croplands and double cropping practices in 2000-2019.

Study area
The NCP (32 • 00 ′ N-40 • 24 ′ N, and 112 • 48 ′ E-122 • 45 ′ E) is located in the eastern part of China (figure 1), covering southern Beijing municipality, Tianjin municipality, southeastern Hebei province, eastern Henan province, Shandong province, northern Anhui province and northern Jiangsu province. The overall area of this study is approximately 40 × 10 4 km 2 with the average elevation of ∼50 m above the sea level. It has a typical temperate monsoonal climate with annual precipitation of 500-1000 mm and annual averaged temperature of 8 • C-15 • C. The prevailing double cropping system is winter wheat and summer maize rotation. It is one of the most important granaries of China .

Remote sensing data and preprocessing
The 250 m MODIS vegetation index collection 6 (C6) data with 16 d interval (MOD13Q1), covering the period of 2000-2019, were collected from the Land Processes Distributed Active Archive Center. An earlier study has evidently reported that multitemporal vegetation index data track similar seasonal responses for the crop growth, and are highly correlated across the growing season (Wardlow et al 2007). Therefore, this study used NDVI data in view of its higher values and more widely-used characteristics (Wardlow and Egbert 2008, Estel et al 2015. The noises induced by poor atmospheric conditions and cloud contamination are inevitable in satellite-based data. To minimize the impacts of noises in time series satellite-based data, previous studies suggested the advantages of using modified Savitzky-Golay (mSG) filter for reconstructing a high-quality 16 d NDVI data (Chen et al 2004. Thus, this study first employed mSG filter with smoothing window size of 2m + 1 (m = 4) and quartic polynomial to smooth all NDVI data.
Besides, some auxiliary data, e.g. existed 1 km land-use data and 300 m land cover data, were used for the cross-validation of typical-year cropland distributions. The 1 km spatial resolution land-use data in 2000, 2010 and 2018, derived from Landsat TM/ETM+ images with the help of human-computer interaction mapping technique, were collected from Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (CAS), thus were also known as CAS land use data (Liu et al 2014a, Ning et al 2018, Cui et al 2019. The first-class landuse types include forests, grasslands, croplands, urban and built-up land, water body and unused land. The dataset has widely been used in previous studies of China's land-use management and change monitoring (Liu et al 2014a, Ning et al 2018, Cui et al 2019. The 300 m climate change initiative (CCI) land cover products were collected from the European Space Agency. Croplands (label codes of 10, 11, 12, 20 and 30) and urban areas (label codes of 190) of CCI data were used to help investigate the impacts of urbanization on cropland area.

Croplands and specific crops mapping 2.3.1. Cropland mapping
The primary land-use types in this study area included forests, grasslands, croplands, built-up lands, barren lands and water bodies (Liu et al 2014a, Ning et al 2018. In this study, forests, grasslands, built-up land, barren land and water body were also regarded as non-croplands. A random sample selection was conducted to collect training pixels for each land-use type. The total number of 2993 samples were manually labeled with the help of 2013-2019 fine-resolution images of Google Earth by considering the dominated land-use types of each 250 m pixel (i.e. over 90% homogeneity), including 1234 cropland samples and 1759 non-croplands samples. More than 45% of samples were from the year of 2017. This study used the latitude, longitude and time information of samples to extract the corresponding time series of NDVI anomaly data (1) where ANDVI i and NDVI i are the ith NDVI anomaly and NDVI values, and NDVI is the averaged NDVI during Julian days 65-289 of each year. A hierarchical mapping strategy was used in this study, i.e. mapping annual cropland mask and then mapping annual crop types with the help of annual cropland mask (figure 2). To classify croplands and non-croplands, three simple but efficient machine learning algorithms, including a feedforward neural network (NN) algorithm, a decision tree (DT) algorithm and a support vector machine (SVM) algorithm, were used in this study. Specifically, the NN algorithm included 15 input nodes corresponding to 15 time series MODIS-NDVI anomaly images from Julian days 65-289 (i.e. the growing season from March to October). The remaining data were discarded in view of low information contents (Lunetta et al 2010). Reducing redundancy input nodes is also helpful for improving the operating efficiency of NN algorithm. In the hidden layer, this study set seven nodes. The output layer included two nodes, namely croplands and noncroplands. The NN classifier was trained by using 2/3 data of all samples collected from Google Earth images. The remaining 1/3 data were used for validating the performance of the NN classifier. Meanwhile, the DT and SVM algorithms were also employed to extract the cropland mask. Similarly, 2/3 data of all samples were trained to obtain the threshold values, and the remaining 1/3 data were used for validating the performance of the classifier. These trained NNbased, DT-based and SVM-based classifiers were subsequently employed to classify croplands and noncroplands of the entire study area in 2017 (considering the most training data from the year of 2017), respectively. A decision rule of cropland mask is defined as two croplands out of three predictions at least. A final cropland mask map in 2017 was obtained according to this decision rule of the ensemble learning.
Furthermore, based on the derived cropland map in 2017, we updated the cropland mask map in 2016 and 2018, respectively. Taking the cropland mapping in 2016 as a case, the specific updated method was detailed explained. This study firstly calculated the correlation and significance of NDVI of Julian days 65-289 and the bias values of maximum NDVI (NDVImax) between 2016 and 2017. For a specific pixel, if time-series NDVI data between 2016 and 2017 had a significant correlation (p < 0.05) and the bias absolute value of 2 year NDVImax was less than 0.1 (∼1 standard deviation of cropland NDVImax, figure S1 (available online at stacks.iop.org/ERL/17/044056/mmedia), land-use type of the specific pixel was not updated. Otherwise, the pixel would be reclassified by using the best performance algorithm out of the three algorithms (e.g. NN algorithm, see the following section 3.1) as well as corresponding NDVI anomaly data in 2016. Using the similar method, we updated the cropland mask map in 2015 based on the derived cropland map in 2016. Similarly, we finally produced all cropland maps in 2000-2019.

Peak information extraction
A 6 • polynomial fitting algorithm was used to fit the curve of vegetation growth for Julian days 65-289 and reconstruct daily NDVI time series year by year  y (DOY) = p 0 + p 1 · X + p 2 · X 2 + · · · + +p n · X n .
( 2) where DOY is an abbreviation of day of year; y (DOY) is the daily NDVI value at time DOY; X is the normalized value of corresponding Julian (e.g. when DOY is 1, X is 1/365), p 0 , p 1 ,…, p n are the fitting parameters; and n is 6 in this study. The corrected Akaike Information Criterion (AICc) (Akaike 1974, Hurvich andTsai 1989), root mean square error (RMSE) and correlation coefficient (R) were used to test potential overfitting risk of the six-degree polynomial fitting algorithm. To obtain the DOY and corresponding NDVI values of peaks and valleys in the curve of vegetation growth, the fist derivative of function y was computed by using the following function where, y ′ (DOY) represents the first derivative of function y. When y ′ (DOY) is zero, the corresponding y (DOY) values are the peak or valley values. When y ′ (DOY − 1) > 0 and y ′ (DOY) ⩽ 0, we can obtain the y (DOY) values and corresponding DOY in peaks. Note that this study also introduces a threshold to remove fake peaks with relatively low NDVI values. Specifically, if the estimated peak NDVI value was less than 0.2 (an empirically statistical value in this study), the corresponding fake-like peak will be removed form peak records. When y ′ (DOY − 1) < 0 and y ′ (DOY) ⩾ 0, we can obtain the y (DOY) values and corresponding DOY in valleys.

Crop-specific type mapping
Winter wheat and summer maize, two major crop types, were considered in this study. Referring to the NCP farming calendar, winter wheat is usually harvested in May and June and subsequently famers grow summer maize. For ease of classification, this study is subjectively divided into two periods by using DOY 177. That is, data from DOY65-177 (a total of eight images) were used for classification of winter wheat and non-wheat, and data from DOY177-289 (a total of eight images) for classification of summer maize and non-maize. This study collected 1767 winter wheat samples and 244 non-wheat samples in growing period of winter wheat, and 1396 summer maize samples and 655 non-maize samples in growing period of summer maize, from the field survey and high-resolution images (unmanned aerial vehicle and GF-2 images) in 2019. The NN algorithm was separately used for classifications of winter wheat and summer maize in view of the best performance in three algorithms (see the following section 3.1). This study took the classification of winter wheat in 2019 as a case to state the specific steps. This study first used the latitude and longitude information of samples to extract the corresponding time series of NDVI anomaly data in 2019. Then, this study employed 2/3 winter wheat and non-wheat samples to train the NN classifier, and the remaining 1/3 data were used for validating the performance of NN classier. The NN classier included eight input nodes corresponding to eight time-series based on NDVI anomaly images from DOY 65-177. In the hidden layer, this study set four nodes. The output layer included two nodes, namely winter wheat and non-wheat. Third, the trained NN classifier as well as 2019 cropland map and peak information was employed to map winter wheat and non-wheat across the entire study area in 2019. Fourth, based on the derived winter wheat map in 2019 as well as cropland map and peak information in 2018, we updated the winter wheat mask map in 2018. Similar to updating cropland map, this study calculated the correlation and significance of NDVI of DOY 65-177 and the bias absolute value of NDVImax between 2018 and 2019. For a specific pixel, if it has a significant correlation (p < 0.05) and the bias absolute value of NDVImax less than 0.1, its land-use type will not be updated. Otherwise, the pixel will be reclassified by using the NN classifier as well as corresponding NDVI anomaly data in 2018. Repeating the fourth step, we finally produced all winter wheat and non-wheat maps in 2000-2019. Similar to winter wheat mapping, this study also produced all summer maize and non-maize maps from 2000 to 2019.

Statistical analyses
The linear relationship between estimated data and statistical data, determination coefficient (R 2 ) and statistical significance were used to assess the overall performance of satellite-based algorithm. Of them, the student's t-test was used to calculate the significance at the 0.05 or 0.01 level. The user's accuracy (UA), producer's accuracy (PA), overall accuracy (OA) and Kappa coefficient (k) were used to assess the algorithm performances.
where n is the total number of all classifications; X i,j is an element, locating at ith row and jth column of the confusion matrix (or error matrix). This study aggregated 250 m maps of croplands, double cropping area, and winter wheat and summer maize rotation area into corresponding fractional (%) maps within a 40 × 40 pixel window (∼10 × 10 km 2 grid cell) and then computed a map of change rates of croplands, double cropping area, and winter wheat and summer maize rotation area from 2000 to 2019 by using the linear regression method.

Analyses of cropland classifications, crop classifications and cropping patterns
Three classifiers used in this study all showed good performances in the training phase. The overall accuracies were separately 90.9%, 88.0% 88.9% for NN-based, DT-based and SVM-based classifiers, while corresponding Kappa coefficients were 0.81, 0.75, 0.78, respectively. The NN-based croplands gave the best user's accuracy (92.6%) and producer's accuracy (86.0%) compared to those from DTbased and SVM-based classifiers. Given the best performance of NN-based classifier, it was thus used for updating cropland mask and producing crop-specific maps referring to the described steps of method section 2.2.3. Regarding NN-based winter wheat, the overall accuracy, user's accuracy and producer's accuracy were 94.5%, 97.3% and 96.4%, respectively, and the Kappa coefficient was 0.76. Regarding NN-based summer maize, the overall accuracy, user's accuracy and producer's accuracy were 92.2%, 93.1% and 95.7%, respectively, and the Kappa coefficient was 0.82. In addition, this study also validated the robustness of satellite-based predicted areas of winter wheat and summer maize against 32 city-level statistical planting area of winter wheat and summer maize from the NCP and compared with previous winter wheat and summer maize maps. Satellite-based predicted area of winter wheat and summer maize both performed good agreements against statistical data, as R 2 showed high values with values of 0.89 and 0.81 (p < 0.01), and the linear regression-based slopes being close to line 1:1 with 0.84 and 0.69 (figure S2). Figure 3(a) shows the NCP cropland mask map in 2017 according to the predetermined decision rule that at least two croplands out of three predictions from trained NN-based, DT-based and SVM-based classifiers. Figure 3(b) shows that double-season cropping patterns were more located in the western and southern NCP, while the single-season cropping patterns were mainly located in the mountainous regions of mid-eastern NCP and coastal regions. Regions with winter wheat and summer maize rotation were mostly concentrated upon western and southwestern NCP, namely southern Hebei province, western Shandong province, eastern Henan province, and northern Anhui and Jiangsu provinces.

Annual areal changes of croplands, single-/double-season cropping, and winter wheat and summer maize rotation
Net cropland area in the NCP presented a significant decreasing trend with slope of −783.8 km 2 a −1 (p < 0.01) during 2004-2019, although it displayed a slightly increasing trend in the entire study period (p = 0.211) because of rapid cropland expansion before the year of 2004 (figure 4(a)). Figure 4(b) shows that single-season cropping area displayed a slightly decreasing trend, and that trend was much steeper after the year of 2004 (slope = −938.3 km 2 a −1 , p = 0.111) compared to that of the entire study period (slope = −61.0 km 2 a −1 , p = 0.893). In contrast, the cropland loss did not directly decrease double-season cropping area. Doubleseason cropping area always remained a slightly increasing trend during the entire study period (slope = 446.9 km 2 a −1 , p = 0.335). Particularly, when considering double-season cropping types, this study, however, observed that the area of the winter wheat and summer maize rotation presented a significantly increasing trend (slope = 3423.3 km 2 a −1 , p < 0.01). The area of other three double-season cropping types, including winter wheat and nonmaize rotation, non-wheat and summer maize rotation and non-wheat and non-maize rotation, were yet all significantly decreasing trends (figure 4(c)). In  other words, double-season cropping area remained relatively stable, but double-season cropping types had respective patterns and trends. This study also found that overall area of winter wheat and overall area of summer maize both exhibited significantly increasing trends with slope of 2953.8 km 2 a −1 and 2002.9 km 2 a −1 (p < 0.01) in the entire study period.
To further clarify the spatial changes of croplands (figures 5(a)-(d), double-season cropping (figures 5(b)-(e)), and winter wheat and summer maize rotation (figures 5(c)-(f)), this study also investigated their corresponding patterns and trends at a 10 km × 10 km grid cell scale. This study observed that unlike patterns of cropland and double-season cropping with areas being more than 40 km 2 , winter wheat and summer maize rotation area more focused on the southern and western NCP. In 2000-2019, areal trends of cropland and double-season cropping showed similar spatial distributions. Their increased regions were mainly concentrated on northeastern and southern NCP. Peri-urban regions however presented significantly decreased cropland and double cropping areal trends. In contrast, significantly increasing regions of winter wheat and summer maize rotation mostly covered plain regions with elevation being less than 50 m.

Reliabilities and uncertainties of cropland and crop-specific type mapping
In this study, 250 m MODIS data and the best performed machine learning approach (i.e. NN algorithm) were used to derive annual cropland and crop-specific type maps. We further proved the reliabilities of satellite-derived cropland and crop type mapping. The good performances could be more contributed to the following two reasons. On the one hand, different from some previous cropland mapping studies that only relied on one machine learning approach, this study employed three typical approaches (i.e. NN, DT and SVM algorithms), time series MODIS NDVI anomaly data and a new decision rule. Taking cropland mapping as an example, that is, cropland mask is defined as two croplands out of three predictions, at least. This kind of treatment declined the uncertainty of derived result that only relied on single machine learning approach to some extent. The methods used by previous studies that employed time series vegetation index values to map croplands or crop types, could be more effective in relatively short term (Wardlow et al 2007, Skakun et al 2017. Because they could overlook the greening effects of CO 2 fertilization and land management for improving the vegetation biomass (Chen et al 2019, Liu et al 2020. The annual vegetation index anomaly could largely minimize impacts of CO 2 fertilization and land management, and thus likely help to improve accuracies of cropland and specific-crop classifications. At the same time, the best performance of machine learning approach was selected to apply for updating cropland and crop-specific type mapping. On the other hand, the update approaches of annual cropland and crop-specific type mapping could play an important role. For each pixel, the correlation and significance of NDVI time series in two neighboring years as well as the bias absolute value of 2 year NDVImax were calculated. The derived information is subsequently used for change detection analyses, that is, determining whether a specific pixel of land-use or crop type will be updated. If a specific pixel has a significant correlation (p < 0.05) and the bias absolute value of NDVImax less than 0.1, its corresponding land-use or crop type will not be updated, or vice versa. This simple method not only improves the efficiency of classifications but it also makes sure the consistency of areal trends of land-use or crop type. Particularly, the latter point is largely very important for investigating the areal trends of croplands and crop-specific types. In view of the open access 250 m MODIS data and robust algorithms, this study successfully obtained good performance cropland and crop-specific type maps with high classification accuracies (user's accuracy, producer's accuracy and overall accuracy being greater than 86%).
Some previous studies had completed winter wheat mapping of the NCP in several specific years (Wang et al 2015, Qiu et al 2017, Dong et al 2020. For example, a study from Wang et al (2015) used the traditional widely-used method of Twi-difference algorithm to derive three-year winter wheat areas of the NCP (2001, 2006 and 2011). They found that winter wheat areas of the NCP increased from 10.43 km 2 in 2001 to 11.10 km 2 in 2006 and to 11.99 km 2 in 2011. In view of using a similar study area in this study and Wang's study, we compared the results of the two studies. The comparison of winter wheat areas showed that both results were relatively consistent in 2011, but our derived areas in 2001 (9.69 km 2 ) and 2006 (11.70 km 2 ) were slightly lower or greater than Wang's results. Besides of the difference in methods used, one of the reasons for differences is probably that they used MODIS Collection 5 (C5) EVI time series data, which could be inevitably affected by relatively lower values and data quality of C5 EVI compared to C6 NDVI (Didan et al 2015, Zhang et al 2017. In comparison to C5 vegetation index, C6 products have gone through four efforts and improvements, including use of pre-composited data, modifying the constrained view angle maximum value composite compositing method, updating the EVI backup algorithm from SAVI to a new robust two-band EVI, and adjusting the relative azimuth angle and using a new dynamic range from −180 • to 180 • (Didan et al 2015). Areal amounts had allowable differences (figures S3-S6), but winter wheat patterns in the NCP showed good spatial consistency with those of previous studies (Wang et al 2015, Qiu et al 2017, Dong et al 2020, Luo et al 2020. Yet different from above-mentioned studies, this study could provide annual ongoing winter wheat products (2000-2019) with the help of relationships between time series vegetation index information in two neighboring years and NN-based machine learning approach. Importantly, long-term maps take us to have the ability to estimate areal trends of the winter wheat. Besides, to investigate changes of double cropping, annual summer maize maps were produced that previous studies rarely reported (Luo et al 2020). This study also observed that the developed method had a relatively stronger ability to derived winter wheat area compared to summer maize area, such as, R 2 being 0.89 for winter wheat estimation and 0.81 for summer maize estimation at the validation stage, respectively. This is largely because that there are more various crop types in summer maize growing phrase compared to winter wheat growing phrase. In collected field training data, autumn-harvest crops cover summer maize, paddy rice, peanut, cotton and soybean. Complex crop types greatly increased the classification challenge of summer crops (Duveiller and Defourny 2010, Chen et al 2016, Bégué et al 2018. Besides, to improve the economic income, farmers usually use intercropping systems in the NCP, e.g. the summer maize and soybean intercropping system. Yet, we found that our summer maize maps well matched spatial distributions and inter-annual trends of Luo et al (2020), although total areas of summer maize were slightly higher than their results. This is likely because that 250 m MODIS NDVI data could more decrease the impacts of mixed pixels compared to 1 km LAI data used in their study. Although there are still some limitations and uncertainties, this study is of significance for understanding the changes of regional croplands, single-/double-cropping systems and specific crop types.

Reasons for areal changes of double cropping, and winter wheat and summer maize rotation
The entire NCP showed a relatively stable status of double cropping area during the years from 2000 to 2019. This stability appears as one of the important reasons that the decreased double cropping area in some regions were totally offset by increased double cropping area in some heavily agricultural counties ( figure 5(e)). This study observed that the significant decreased cropland and double cropping areas were mainly located in regions around cities, in which cropland and double cropping areal decreasing trends also gave similar spatial distributions (figures 5(d) and (e). Rapid urbanization, decreasing rural active labor force and changing regional cropping systems are likely to be several important reasons of areal decrease. Previous studies have documented that urban expansion occupied a mass of high-quality croplands around cities (Liu et al 2014a, Ning et al 2018. The1 km resolution NCP land-use data in 2000 and 2015 showed that cropland loss was mainly caused by urban expansion, accounting for 86.4% of cropland total loss (Liu et al 2014a, Ning et al 2018, Cui et al 2019. Another result based on 1992-2020 CCI data displays that 91.6% of decreased croplands are converted into urban and built-up land in the NCP (figure 6). Total area of converting from croplands to urban and built-up land accounts for 74.3% of all urban expansion area, and note that the rate has overpassed 80% since 2010. On the other hand, rural active labor forces, especially people of 18-50 years old, are inclined to work in cities in view of finer wage income, and therefore have little time to return form cities for being engaged in agriculture (Liu et al 2014b). Reports from China's National Bureau of Statistics present that China's urbanization rate has nearly doubled from 30.9% in the year of 1999 to 59.6% in the year of 2018 (Du and Lin 2019). This means that there are more than 10 million agricultural population entering cities every year in China. These farmers' cultivated land could suffer from the single-season or double-season abandonment risk if not cultivated by other family members (Liu et al 2014b, Liu 2018. Besides, plastic greenhouse usage and cash crops planting expansion in some special regions likely affect the magnitude of double cropping area. For example, plastic greenhouse is widely used for vegetables or flowers planting in Shouguang county and Qingzhou county of Shandong province (Yang et al 2017), given that the economic income from vegetables or flowers can better support rural livelihoods compared to those from grain crops.
Excitingly, this study however finds that although cropland area has presented a significantly decreasing trend since 2004, double cropping area of the entire NCP performs a relatively stable status. Particularly, the area of winter wheat and summer maize rotation presents a significantly increasing trend. This phenomenon is largely contributed to the respective  Interannual changes of averaged economic benefits (a), crop purchase prices (c), the shaded area represents the corresponding standard deviation, pre-unit yield (d) and the relationships between averaged economic benefits and area (b). The symbol * * represents the significance at the 0.01 level. area increases of winter wheat and summer maize, but is also greatly associated with relatively higher economic benefits of winter wheat and summer maize in the NCP (figure 7), including advantage of per-unit yield, the increasing purchase price of grain crops, as well as transitions beginning to subsidize rather than tax agriculture in China in 2004. Rice also has a relatively higher economic benefits, but available water resource of the NCP limits its planting area. Besides, this study observes a steep increasing areal trend of winter wheat and summer maize rotation in the NCP after the year of 2010. A potential reason is that national and some local land-use and subsidy policies of supporting agriculture and giving favorable treatments to farmers, e.g. improving protective grain purchase price, increasing grain subsidies and enhancing subsidies of purchasing agricultural machinery, have greatly promoted farmers' enthusiasm since the beginning of the year of 2010. These findings suggest that positive local land-use policies and stabilized grain price are very important for promoting cropland intensive level, and maintaining the stability of the double cropping area and grain crop areas.

Potential applications and future development
The developed method in this study has a large potential to be used for continuous and large-scale cropping frequency mapping owing to its robust performances. In view of the purposes of this study, we only used the developed method and products to analyze patterns of croplands, double cropping, and winter wheat and summer maize rotation, and their areal changes. However, some potentially fruitful extensions of this developed method and products are being mentioned. For example, the developed method and products could have the large potential to monitor single cropping pattern and land abandonment by considering annual cropping frequency changes, which would be greatly helpful for land-use management and policy-making (Prishchepov et al 2013, Renwick et al 2013, Spera et al 2014, Estel et al 2015. Yet, it should be noted that land abandonment monitoring is challenging owing to its transient and fragmentized features. Although the developed method has a robust performance and potentially fruitful extensions as AICc with averaged value and standard deviation of 6.76 ± 1.20, RMSE with 0.0135 ± 0.0075 and R with 0.988 ± 0.015 (figure S7), its corresponding limitations also need to be noted. First, the developed method is currently appropriate for regions where cropping frequency is two crops per year or less. Because it uses a six-degree polynomial function to fit vegetation growth curve and the first derivative of fitting curve to obtain peaks and valleys in its core algorithm. If cropping frequency is more than two crops per year, this method could be overfitting crop growth curve. And then, it could estimate some incorrect peaks and valleys. Certainly, estimating cropping frequency of three crops or more per year could be, to our knowledge, a widespread challenge for current satellite-based cropping frequency extraction methods (Wardlow and Egbert 2008, Cohn et al 2016. Second, cropping frequency in this study relies on peak number and their corresponding NDVI values. In normal years, it is relatively easy for accurately extracting cropping frequency for good crop phenological profiles. However, when crops suffer from extreme climatic events, e.g. severe drought event Yang 2009, Wang et al 2016), crop phenological profiles are largely affected. Similar phenomena could lead to an underestimation of cropping frequency. Thus, our further efforts will focus on overcoming these above-mentioned limitations.

Conclusions
This study introduced an approach, i.e. employing 250 m time-series MODIS NDVI anomaly data, the correlation and significance of NDVI time series in two neighboring years, crop peak information and typical machine learning algorithms (e.g. a feedforward NN algorithm, a DT algorithm and a SVM algorithm), to derive and update annual croplands and specific crop type area maps. The introduced method gave good performances against validated data and statistical data. More importantly, it estimated long-term continuous dynamics of croplands and cropping frequency. We found that cropland area presented a significantly decreasing trend since 2004 in the context of rapid urbanization. However, the double cropping area of the NCP showed a balanced magnitude but divergent patterns. Particularly, areas of the winter wheat and summer maize rotation presented significant increasing trends. Areas of other three double-season cropping types, including winter wheat and non-maize rotation, non-wheat and summer maize rotation and non-wheat and non-maize rotation yet followed significantly decreasing trends. To our knowledge, this study is the first endeavor in mapping long-term and spatial distributions of winter wheat and summer maize rotation in the most important granaries of China. The satellite-based maps and findings in this study could potentially provide some valuable information for regional landuse management and agriculture policy-making.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// figshare.com/articles/dataset/Winter_wheat_and_ summer_maize_datasets_over_the_North_China_ Plain_in_2000-2019/19228224.