Tracking the spatio-temporal change of cropping intensity in China during 2000–2015

Improvement in the efficiency of farmland utilization and multiple cropping systems are of prime importance for achieving food security in China. Therefore, spatially-explicit analysis detecting trends of cropping intensity are important preconditions for sustainable agricultural development. However, knowledge about the spatiotemporal dynamics of cropping intensity in China remains limited. In this study, we generated annual cropping intensity maps in China during 2000–2015 using a rule-based algorithm and MOD09A1 time series imagery. We then analyzed the spatio-temporal changes of cropping intensity. The results showed single-cropping and double-cropping areas were about 1.28 ± 0.027 × 106 km2 and 0.52 ± 0.027 × 106 km2 in China in 2015 and their areas were relatively stable from 2000–2015. However, cropping intensity had substantial spatial changes during 2000–2015. About 0.164 ± 0.026 × 106 km2 of single-cropping area was converted to double-cropping area, which mainly occurred in the Huang-Huai-Hai Region. About 0.193 ± 0.028 × 106 km2 of double-cropping area was converted to single-cropping area, which mainly occurred in the southern part of China. About 85% of croplands with decreases in cropping intensity were located in the southern part of China, and about 80% of croplands with increases in cropping intensity was distributed in the Huang-Huai-Hai Region and the northern part of the Middle and Lower Reaches of the Yangtze River region (p < 0.05). The landscapes of different cropping systems tended to be homogenized in major agricultural production regions.


Introduction
Increasing human population has made it difficult for agricultural systems to meet demand for crop and livestock production, thereby making food security an especially serious challenge for China (Foley et al 2005). Increasing cropland area and yield have been the primary and direct solutions for the rising agricultural production needs over the past decades (Godfray et al 2010, Ray et al 2012. Cropland expansion is mainly converted from natural ecosystems and agricultural intensification is mainly based on excessive fertilizer applications, which have resulted in substantial threats to ecosystems (Lambin and Meyfroidt 2011). Rising concerns are focusing on exploring other means for achieving the balance between food security and environmental security (Foley et al 2011, Ray andFoley 2013), and increasing cropping intensity is regarded as one of critical and effective approaches (Turner and Doolittle 1978, Ray and Foley 2013, Stephan et al 2016. Cropping intensity refers to the cropping frequency in a given cropland area per year, which is often measured using multiple cropping indices (MCI). Mean cropping intensity of a region represents the average intensity of areas cropped across the whole region during one year.
Cropping systems in China are diverse due to complex, varied climate and socio-economic conditions. Single-cropping rice, maize, or soybean systems dominate Northeast China . Wheat/ maize and wheat/soybean rotations are the major cropping systems in most of North China Plain . Rice-based cropping systems are prevalent in the Middle and Lower Reaches of the Yangtze River Region (Leeming 1979) with winter wheat/rice in the northern Yangtze River region and rice/rice/green manure cropping rotation in the southern Yangtze River region. Double season rice and rotations of rice, potato, and sugarcane are the major cropping systems in South China.
Several previous studies were carried out to identify cropping practices using various data sources. The National Agricultural Census datasets were used to calculate MCI for the administrative provinces (Xie and Liu 2015), which reflects the general cropping intensity but does not capture detailed spatial information on cropping intensity (Yan et al 2014). Satellite remote sensing can repeatedly observe the land surface, which makes it possible to monitor crop phenology and the spatial pattern and temporal changes of cropping intensity at large spatial scales. The spectral reflectance bands from optical satellites can be used to calculate vegetation indices (VIs), such as Landsat (Jain et al 2013), SPOT Normalized Difference Vegetation Index (NDVI) (Ding et al 2015), GIMMS-NDVI , and MODIS (Yan et al 2014). Although Landsat has a relatively high spatial resolution of 30 m, it has a temporally coarse 16 d revisit cycle and observations are easily affected by frequent cloud contamination. MODIS provides daily observations at a spatial resolution of 500 m but has a higher chance to capture cloud-free observations over the same time period, which are used to composite data products that are less affected by cloud cover and shadows. Considering the trade-off between spatial and temporal resolution, MODIS data is the better choice for identifying and tracking the spatio-temporal changes of cropping intensity across broad extents due to its higher frequency and higher number of cloud-free observations.
Additionally, different techniques and algorithms have been developed to quantify the cropping frequency (Jain et al 2013). The MODIS peak value method has been widely used, which defines the number of peaks on MODIS time series VIs as the number of cropping cycles (e.g. Biradar and Xiao 2011). Correspondingly, various techniques have been used to extract the peak frequency, such as the wavelet transform method (Qiu et al 2017) and the two-difference algorithm . However, using only these statistical methods to interpret the remote sensing data will likely introduce uncertainty because they do not consider in situ knowledge about crop phenology and intensity.
Existing annual MCI products of China cover various spatial resolutions (500 m, 1 km) and time spans (the latest year is 2013) (Ding et al 2015, Qiu et al 2017. Previous studies did not achieve consensus on the overall trend of cropping intensity change in China since 2000. Ding et al (2015) found that cropping intensity increased significantly during 1999-2013 in China based on 1 km SPOT NDVI imagery, while Qiu et al (2017) found that the cropping intensity slightly declined during 2001-2013 based on 500 m MODIS imagery. Therefore, knowledge about the spatio-temporal changes in cropping intensity in China remains limited and more accurate cropping intensity maps are needed.
The objectives of this study were two-fold: (1) to develop annual cropping intensity maps with high accuracy in China from 2000-2015 using a rule-based algorithm on time series MODIS enhanced vegetation index (EVI) and in situ crop phenology data from field observations; and (2) to assess the spatio-temporal changes of cropping intensity in China from 2000-2015. The results from this study can provide informational support for making cropland use policy, which is important for food security in China.

Materials and methods
2.1. Data 2.1.1. Cropland data Cropland data was extracted from China's National Land-Use/Cover Change Dataset (NLCD-China) at a mapping scale of 1:100 000 by Chinese Academy of Sciences. The NLCD-China dataset was generated based on visual interpretation of the 30 m Landsat TM/ETM+images (Liu et al 2005. The cropland layer used in our study had high accuracy (∼95%) at the national scale based on substantial ground references. The interpreted NLCD-China was converted to gridded raster datasets in the following steps (Liu et al 2005). First, a fishnet of 1 km×1 km grid was created. Second, the vector map was overlain with the fishnet and the area fraction for each land cover and land use type was calculated for each grid cell. Third, a series of raster layers for six land cover types and 25 five land use types were developed. Cropland area decreased about 0.7% of the total cropland in China from 2000 to 2010 . In this study, to reduce the influence from cropland change, we used the spatial distribution map of cropland in 2000 generated from those pixels where cropland was the largest area proportion. Then the cropland map was resampled to match the 500 m MODIS-derived cropping intensity data. where, ρ nir , ρ red and ρ blue are the values of reflectance of near infrared (841-876 nm), red (620-670 nm), and blue (459-479 nm) bands, respectively.

Crop calendar data from agro-meteorological observations
Crop phenology information derived from agrometeorological observations in 2002 was incorporated to guide decision rules in the cropping intensity mapping algorithm because of its data availability (Yan et al 2014). A total of 394 in situ observations from the national agro-meteorological stations across China were used (figure S1 is available online at stacks.iop. org/ERL/14/035008/mmedia), which is freely available from the China Meteorological Data Sharing Service System (https://data.cma.cn/en). The agrometeorological observations included records of critical phenological phases (i.e. planting date, emergence date, and maturity date) covering 11 major crop types (early rice, late rice, single rice, winter wheat, spring wheat, spring maize, summer maize, cotton, soybean, sorghum and rapeseed). The crop phenology information provided the ground references we used to build decision rules and thresholds for mapping cropping intensity.

Statistic data
We collected provincial agricultural and socio-economic datasets in China. We used these datasets to compare our MODIS-based cropping intensity product and analyze the potential factors influencing the spatio-temporal changes of cropping intensity. Table 1 is a brief summary for the provincial agricultural and socio-economic statistic datasets. Cost-profit ratio for major crops refers to the ratio of net profit to total cost for the agricultural production of major crops. The cost-profit ratio not only considered the farmer profit but also the total cost during the agricultural production including labor cost, land cost, and material cost such as fertilizer input. The major crop types in our study included Indica rice, Japonica rice, corn, and wheat.

A rule-based algorithm for mapping cropping intensity
Annual time-series EVI data can clearly reflect crop phenological stages including planting, emergence, heading, maturity, and harvest (Sakamoto et al 2005). Correspondingly, the dynamic process of 'rise-peakfall' in EVI time series represents the cropland growth cycle. EVI reached its peak on the heading date. We also showed the capacity of MODIS EVI time series to catch the crop growth cycle at three cropland ecosystem field stations with multiple cropping areas in China (Yan et al 2014).
Based on the phenology data of in situ field observations, we developed a rule-based algorithm to identify and map the spatial distribution of cropping intensity in China in 2002 (Yan et al 2014). The algorithm consisted of three following procedures: (1) smoothing the EVI time series using the Harmonic Analysis of Time Series method to remove noise and construct gapless images; (2) detecting the potential cropping cycles by determining the number of peaks of the smoothed EVI time series profile using a moving window method; and (3) determining the reasonable peak values as the number of cropping cycles. Decision rules and thresholds were further set to determine reasonable peak in the EVI profiles based on crop calendar information from agro-meteorological observations. First, the reasonable peak value of EVI should be higher than 0.35, which represents the peak value of semi-arid grasslands. Second, the temporal interval between two probable peak dates of crop growth should be larger than 80 d. Third, the possible earliest and latest dates that reasonable peaks occur should be earlier than the 320th day of the year and later than the 56th day of the year. Figure 1 showed the smoothed EVI with the crop calendar information from Yugan station, which is in a typical double cropping region. For more detailed information about the rule-based algorithm of cropping intensity mapping, please refer to our previous study (Yan et al 2014). In this study, we applied this rule-based algorithm and generated annual cropping intensity maps in China from 2000-2015.

Accuracy assessment and comparison of MODISbased cropping intensity maps
To assess the accuracy of MODIS-based cropping intensity maps at both pixel and regional scales, we collected ground reference samples and the government agricultural census data. The ground reference samples were visually interpreted from time series multiple-source high spatial resolution images, including Landsat 4-5 TM, Landsat 7 ETM+, Landsat 8 OLI, and Sentinel-2 time series images in 2000 and 2015. The government agricultural census data at the provincial level in 2000, 2005, and 2010 were used in this study. The combined Landsat-like sensors provide relatively high frequency of land surface observations at the spatial resolutions from 10 to 30 m, which provides a chance to visually interpret cropping cycles at relatively fine scales. Considering the very small area proportion (∼2% of total cropland area) of triplecropped area in China, the accuracy assessment was conducted for single-and double-cropping practices. Firstly, 200 ground samples were generated using stratified random sampling in ArcGIS based on cropping intensity map in 2015 (figure S2). The sample numbers for single-(146 pixels) and double-cropping (54 pixels) areas were approximately consistent with the area proportion of these cropping systems. Secondly, visual interpretation of cropping intensity for these samples was conducted for 2000 and 2015 using a combination of all available Landsat-like images. Based on the phenology characteristics of single-and double-cropping systems (see figures S3 and S4), the visual interpretation was performed using the LandLook Viewer (https://landlook.usgs.gov/), which is a prototype tool for rapid online viewing and interactively explore Landsat and Sentinel-2 archives at up to full resolution. Thirdly, we generated a confusion matrix for cropping intensity maps in 2000 and 2015 and cropping intensity change from 2000 to 2015, and we then calculated the overall accuracy, user accuracy, and producer accuracy. Area estimates (uncertainties) were also calculated from the confusion matrix using the method described in Olofsson et al (2014). At the regional scale, agricultural sown areas estimated from this study and the government agricultural census were compared at the provincial level for 2000, 2005, and 2010. Pixel-level agricultural sown area was calculated by multiplying cropland area by cropping intensity for each pixel in 2000, 2005, and 2010. Provincial agricultural sown areas were the sum of agricultural sown area of each grid cell. Finally, linear relationships between agricultural sown areas estimated from this study and the government agricultural census were analyzed for 2000, 2005, and 2010.

Methods for characterizing the spatial pattern of cropping intensity trends
To characterize the changes in the spatial and temporal patterns of cropping intensity, we aggregated the 500 m cropping intensity maps into 5 km gridcells using two aggregation approaches. In the first approach, we generated annual 5 km cropping intensity maps using the average aggregation method in the period of 2000-2015. In the second approach, we generated for the same period annual 5 km area fractional maps of single-and double-cropping based on the sum aggregation method. Then, we performed trend analysis for annual 5 km cropping intensity maps using linear regression models with t test at the 5% significance level. Two landscape metrics were applied to quantify spatial fragmentation of cropping pattern, which were calculated using software FRAGSTATS 4.2 (McGarigal et al 2012). Patch density and Landscape Fragmentation Index were used to quantify the degree of fragmentation at the landscape and regional scales, respectively where N i represents the patch numbers of landscape i; S i represents the area of landscape i where N refers to the total number of patches; S refers to the total landscape area of the region.

Accuracy assessment and comparison of MODIS-based cropping intensity maps in China
According to the confusion matrix (  (table 3). The overall accuracy of change in cropping intensity was 91.5%. Both the user accuracy and producer accuracy of single-cropping converted to double-cropping were 75% and the user accuracy and producer accuracy of double-cropping to single-cropping were 73.7% and 82.4%. According to figure 2, the agricultural sown areas calculated from the MODIS cropping intensity maps and from the government agricultural census data had reasonably good and stable linear relationships in 2000 (R 2 =0.88, p<1e −5 ), 2005 (R 2 =0.86, p<1e −5 ), and 2010 (R 2 =0.89, p<1e −5 ), respectively (figure 2). The slopes ranged from 0.75 to 0.86 and the RMSE was 8.44×10 6 km 2 , 6.73×10 6 km 2 , and 7.95×10 6 km 2 for 2000, 2005, and 2010, respectively. The crop planting area estimated by MODIS data in our study was larger than the reported area by agricultural census data, which was also found in previous studies (Frolking et al 1999). Two factors may contribute to the uncertainty: (1) the underestimate of cropland area in China from official agricultural census (Frolking et al 1999, Seto et al 2000, Xiao et al 2003, Liu et al 2005, which may be due to political and policy factors (Seto et al 2000, Xiao et al 2003, and (2) the inherent overestimate of cropland area in remote-sensing estimates due to the coarse spatial resolution for small field patches (Frolking et al 1999, Seto et al 2000.

Spatial distribution of MODIS-based cropping intensity in China in 2015
The MODIS-based cropping intensity map showed clear spatial patterns of cropping intensity in China in 2015. Single-cropping had the largest area of approximately 1.28±0.027×10 6 km 2 , accounting for 70.14% of the total cropland area, which was concentrated in northern China and slatternly distributed in southern China ( figure 3). The double-cropping area was 0.52±0.027×10 6 km 2 , accounting for 27.74% of the total cropland area. Double-cropping was mainly distributed in the Huang-Huai-Hai Region (37.71%) and the Middle and Lower Reaches of the Yangtze River Region (30.93%). Figure 3 showed triple-cropping areas occupied a small proportion of cropland (∼2% of total cropland area), so we mainly focused on the single-and double-cropping areas in the following analyses.

Changes in the spatio-temporal distribution of cropping intensity
At the country scale, the total single-cropping (1.24± 0.027×10 6 km 2 in 2000 to 1.28±0.027×10 6 km 2 in 2015) and double-cropping (0.54±0.027×10 6 km 2 in 2000 to 0.52±0.027×10 6 km 2 in 2015) areas remain relatively stable. However, there was a noticeable spatial transition of cropping intensity between 2000 and 2015. About 0.164±0.026×10 6 km 2 of single-cropped land was converted to double-cropped land, mainly in the Huang-Huai-Hai Region. About 0.193±0.028×10 6 km 2 of double-cropped land was converted to single-cropped land, which mainly occurred in the southern part of China ( figure 4(a)). The spatial pattern of change in cropping intensity during 2000-2015 was shown in figure 4(b). Approximately 85.42% of cropland with decreased cropping intensity was in the southern part of China. About 79.96% of cropland with increased cropping intensity was distributed in the Huang-Huai-Hai Region and the northern part of the Middle and Lower Reaches of the Yangtze River region. In addition, we determined the spatial pattern of area changes of single-and double-cropping systems during 2000-2015. Significant expansion of singlecropping occurred in the southern part of China (the Middle and Lower Reaches of the Yangtze River Region and the Southwest Region), while singlecropping area significantly shrank in the Huang-Huai-Hai Region ( figure 4(c)). The reverse trend was observed for double-cropping areas ( figure 4(d)). Single-cropping was dominate in Northern China and was relatively stable.  At the regional scale, regional variations in cropping intensity changes were also detected (table 4, figures 4(e) and (f)). The Huang-Huai-Hai Region had a significant increased cropping intensity (p<0.05) through the expansion of double-cropping area and a slight reduction in single-cropping area. The Middle and Lower Reaches of the Yangtze River Region, the Southwest Region, and the South China Region experienced significant decreases in cropping intensity (p<0.05), with an increase in single-cropping area and a decrease in multiple-cropping area. In the other agricultural regions, there was a general absence of any significant trends both in average cropping intensity and area of different cropping systems.

Improvement to previous cropping intensity methods and products
The cropping intensity calculated based on the agricultural census data at the administrative level (Xie and Liu 2015) has large uncertainty, restricted by spatial resolution, cost, and the knowledge of census people. Our remote sensing approach provided a detailed spatial and spatiotemporal dynamics of cropping intensity and has high efficiency to track cropping intensity through time (Yan et al 2014). Compared with previous research on change in cropping intensity in China based on remote sensing since 2000 (Ding et al 2015, Qiu et al 2017, our study incorporated crop-phenology information from in situ field observations into the algorithm and quantified the area, area uncertainties, and spatiotemporal changes of cropping intensity. Our results showed improvements in two aspects: higher accuracy and more reasonable results for cropping intensity in China. First, our improved methodology incorporated crop-phenology information and led to higher accuracy. Ding et al (2015) performed the NDVI peak-based algorithm using SPOT NDVI at 1 km spatial resolution with an overall accuracy of 91.95%. Qiu et al (2017) applied the wavelet feature-based method to calculate cropping intensity using MODIS EVI with an overall accuracy of 91.63% (in 2013). We incorporated the crop-phenology information from in situ field stations to improve the EVI peak-based algorithm, which helped make our rule-based algorithm and helped determine reasonable peak numbers in the EVI profiles. We set the thresholds for possible minimum peak value, the minimum temporal interval between two probable peaks, and the possible earliest and latest dates that the peaks reasonably occur according to the field-observed crop calendar information to avoid misclassification. The overall accuracy of cropping intensity in 2000 and 2015 in our study were 95.5% for each year, which was higher than the previous two studies. Additionally, the previous studies did not have accuracy assessments for the spatiotemporal changes of cropping intensity in China, which is important to analyze the uncertainty of their cropping intensity change results. In our study, we not only evaluated the accuracy of cropping intensity in 2000 and 2015, but also the accuracy of change in cropping intensity from 2000-2015 based on the ground references. The accuracy of change in cropping intensity from 2000-2015 was 91.5%.
Second, a more accurate cropping intensity product leads to more reasonable results of the spatiotemporal change of cropping intensity in China. For the trend in cropping intensity in China, Ding et al (2015) found that cropping intensity increased significantly by 1.3% per year during 1999-2013 in China, while Qiu et al (2017) found that the cropping intensity slightly declined during 2001-2013. According to the accuracy assessment, we carried out uncertainty analysis for the cropping intensity change in China (Olofsson et al 2014). Our results showed that the total cropping intensity was relatively stable. The single-cropping area changed from 1.24±0.027× 10 6 km 2 in 2000 to 1.28±0.027×10 6 km 2 in 2015 and the double-cropping area changed from 0.54± 0.027×10 6 km 2 in 2000 to 0.52±0.027×10 6 km 2 in 2015. This meant the total cropping intensity was relatively stable in China from 2000-2015. For example, our product did not show many significant changes in cropping intensity in the Loess Plateau. However, both Ding et al (2016) and Qiu et al (2017) showed that significant increases in cropping intensity occurred in the revegetation area in Loess Plateau , Feng et al 2016, where large areas of farmland and barren land have been converted to woodland and grassland under the 'Grain for Green' project . Additionally, rising labor scarcity leads to decline in cropping intensity, the rising out-migrant population in southern part of China (figures 6(d) and (e)) will aggregate the labor shortage and then result in decreasing cropping intensity (Jiang et al 2013). Therefore, it is reasonable that our product and Qiu, et al (2017) showed significant decreasing trends of cropping intensity in Southern China.
Our previous work mainly focused on the algorithm and spatial distribution of cropping intensity in 2002 (Yan et al 2014). This study extended our previous efforts to generate the annual cropping intensity maps during 2000-2015 and address the trend of cropping intensity since 2000. The accuracy assessment of the algorithm was firstly provided in this paper not only in terms of cropping intensity mapping in 2000 and 2015, but also that of the change analysis, which have demonstrated the consistency and robustness of the algorithm for tracking spatiotemporal changes of cropping intensity in China. Based on the annual cropping intensity products, we analyzed the spatiotemporal changes of cropping intensity and improved the previous knowledge about the cropping intensity in China since 2000.

Change in landscape pattern of cropping systems
Spatiotemporal changes in cropping intensity not only include the land use changes of different cropping systems, but also on the change in the spatial configuration pattern of cropping systems. The annual fragmentation of cropland decreased in China but the change was statistically insignificant (p>0.05). However, the spatial fragmentation in the Huang-Huai-Hai Region and the Middle and Lower Reaches of the Yangtze River Region significantly declined by 61.51% and 34.16% (p<0.05) ( figure 5(a)), implying that cropping systems tended to be homogenized and cropland with the same cropping system tended to be more distributed in 2015 than in 2000 ( figure 6).
We observed a negative trend in patch density in single-cropped lands of China (p<0.05), declining from 2000 to 2015 by 32.49% ( figure 5(b)). Regionally, a significant decreasing trend was also detected in the Huang-Huai-Hai Region, the Middle and Lower Reaches of the Yangtze River Region, and the Southwest Region, where the patch densities have declined by 66.75%, 69.07%, and 34.55%, respectively ( figure 5(c)). This decrease implied that the degree of fragmentation of single-cropped land declined and that these areas tended to be more concentratedly distributed ( figure 5(b)). We detected no trend in patch density of double-cropping cropland across China ( figure 5(b)). Regionally, double-cropped lands tended to be more concentrated in the Huang-Huai-Hai Region and the Middle and Lower Reaches of the Yangtze River Region and more dispersed in the South China Region ( figure 5(c)).

Possible factors influencing the cropping intensity change
The potential factors contributing to the spatiotemporal changes of cropping intensity in China are complex. In the long run, climate changes have also been found to have effects on cropping intensity (Iizumi and Ramankutty 2015). Climate warming in China have increased the potential cropping intensity in the last 50 years (Liu et al 2013a), including the Tibetan Plateau (Zhang et al 2013). Meanwhile, crop mixes and crops distribution were also influenced by climate change (Yin et al 2016, Zhang et al 2018. For example, climate warming has made the planting boundary shift northward for maize and rice in Northeast China and rice in Southern China (Liu et al 2013b, Ye et al 2015. Socio-economic factors also effected the trend of cropping intensity change in two opposite ways. On one hand, a series of new agricultural policies implemented the early 21st century benefited farmers and increased farmers' enthusiasm for growing grain, including the elimination of the agricultural tax, increasing agricultural subsidies, and increasing grain price ( figure 7(a)). These policies were supposed to increase farmers' incomes and investment in agriculture, which would positively contribute to intensification of cropland (Gale et al 2005, Zuo et al 2014. Also, institutional innovation and policy support will increase farm operational scale in China through land transfer, land consolidation, and mechanization (Huang and Ding 2016), which would help the distribution of croplands become homogenized. On the other hand, urban expansion and growth in out-migrant agricultural workers would lead to lower cropping intensity (Jiang et al 2013, Xie andJiang 2016). The increasing out-migrant population and rising wages (figures 7(d), (e)) caused declines in available farm workers and a continuous growth in labor cost for major crops ( figure 7(b)). Although there were fewer farm workers, the land owner often preferred to retain their land for family farming (Xie and Jiang 2016), which would lead to decreased cropping frequency (Su et al 2016). Among agricultural regions, the out-migrant population trend was consistent with the distribution of hotspots of decline in cropping intensity. From 2000-2010, more provinces in the southern part of China, except for coastal provinces, became the main sources of migrant workers (figures 7(d) and (e)). Meanwhile, the cost-profit ratio of major crops, i.e. the rate of profit, soared to the highest level in 2004 and then exhibited a decreasing trend ( figure 7(c)). The decline in the rate of profit would reduce the enthusiasm of farming. This variation might explain the decrease in the area of multiplecropping starting in 2004.

Linkage with food security
In an effort to achieve food security, a series of agricultural programs were launched to provide agricultural subsidies, increase grain price, and improve soil quality and crop varieties to increase farmers' enthusiasm for growing grain . Although the agricultural areas with multiple-cropping systems decreased, grain production continued to increase according to the agricultural census, mainly due to increased grain yield and cropland reclamation. The average grain yield increased from about 5000 kg/ ha in 2000 to 6000 kg/ha in 2013, which can be largely attributed to increased investment in agriculture . However, a large gap remains between grain production and demand in China (Wei et al 2015) and about 1.2×10 8 t of grain (20% of grain production in China) was imported in 2015.
Cropping intensity is an important indicator of agricultural land use efficiency (Rudel et al 2009, Licker et al 2010, Neumann et al 2010. Although improving cropping intensity is an effective way to increase grain production, suitable cropping systems still depend on the local natural environment and socio-economic conditions. The spatio-temporal change of cropping intensity from our study clearly showed the regions that had significant decreases and increases in cropland area from 2000-2015. According to the changes in cropping intensity, increased cropping intensity in some regions would be a potential way to increase grain if the ecological and environmental conditions are suitable (Zuo et al 2014). For example, we should evaluate the potential for improving cropping intensity in Southwest China despite ecological and environmental limitations, such as intensive soil erosion (Zuo et al 2014). Some agricultural regions show increased cropping intensity, such as the Huang-Huai-Hai Region, a region with high grain production. In this region, the increasing grain production is based on the use of excessive fertilizer, underground water, and biocides, which reduce the sustainability of cropland production and food security (Guo et al 2010, Rodell et al 2018, as well as ecosystems (Sanderson et al 2002, Foley et al 2005, Kareiva et al 2007, Ellis and Ramankutty 2008. Therefore, a balance should be reached between cropland sustainable production and ecosystem services.

Conclusions
This study explored changes in cropping intensity and fragmentation of different cropping systems across China during 2000-2015. In 2015, single-cropping had the largest area of approximately 1.28±0.027×10 6 km 2 , accounting for 70.14% of the total cropland area and double-cropping area was 0.52±0.027×10 6 km 2 , accounting for 27.74% of the total cropland area. From 2000-2015, the total area of single-cropped and doublecropped lands was relatively stable at the country scale. However, the distribution of cropping intensity significantly changed between 2000-2015. About 0.164± 0.026 ×10 6 km 2 of single-cropped land was converted to double-cropped land, mainly occurring in the Huang-Huai-Hai Region. About 0.193±0.028×10 6 km 2 of double-cropped land was converted to single-cropped land, which occurred mainly in the southern part of China. Substantial changes occurred through spatial comparison: cropping intensity stabilized in Northern China, increased in the Huang-Huai-Hai Region, and decreased in Southern China. Single-cropping increased in Southern China, which was opposite to the decrease observed in the Huang-Huai-Hai Region. We also found that the trend in double-cropping reversed. The decline in the cropping intensity in Southern China was correlated with the increased out-migrant population trend in the region. Regional cropping system patterns tended to increase in spatial consistency and landscapes of the same cropping practices became less fragmented in the main agricultural regions (the Huang-Huai-Hai Region and the Middle and Lower Reaches of the Yangtze River Region).