Time Tracking of Different Cropping Patterns Using Landsat Images under Different Agricultural Systems during 1990–2050 in Cold China

Abstract

Rapid cropland reclamation is underway in Cold China in response to increases in food demand, while the lack analyses of time series cropping pattern mappings limits our understanding of the acute transformation process of cropland structure and associated environmental effects. The Cold China contains different agricultural systems (state and private farming), and such systems could lead to different cropping patterns. So far, such changes have not been revealed yet. Based on the Landsat images, this study tracked cropping information in five-year increments (1990–1995, 1995–2000, 2000–2005, 2005–2010, and 2010–2015) and predicted future patterns for the period of 2020–2050 under different agricultural systems using developed method for determining cropland patterns. The following results were obtained: The available time series of Landsat images in Cold China met the requirements for long-term cropping pattern studies, and the developed method exhibited high accuracy (over 91%) and obtained precise spatial information. A new satellite evidence was observed that cropping patterns significantly differed between the two farm types, with paddy field in state farming expanding at a faster rate (from 2.66 to 68.56%) than those in private farming (from 10.12 to 34.98%). More than 70% of paddy expansion was attributed to the transformation of upland crop in each period at the pixel level, which led to a greater loss of upland crop in state farming than private farming (9505.66 km² vs. 2840.29 km²) during 1990–2015. Rapid cropland reclamation is projected to stagnate in 2020, while paddy expansion will continue until 2040 primarily in private farming in Cold China. This study provides new evidence for different land use change pattern mechanisms between different agricultural systems, and the results have significant implications for understanding and guiding agricultural system development.

1. Introduction

Crops are basic requirements for feeding people, and cropping pattern information is important for human health and grain ration security. In recent decades, although cropland expansion has occurred in many places [1,2,3,4], the maintenance of cropland area is still a challenge worldwide [5]. Such issues are particularly prevalent in China, which feeds approximately twenty percent of the global population with no more than seven percent of the world’s cropland [6,7,8]. Urbanization, industrialization and rural expansion have reduced the cropland area in South China [9,10], and substantial cropland reclamation has been undertaken in North China due to a large amount of adaptive land resources [1]. The rapid expansion of agricultural land in North China has greatly affected the water cycles [11], climate change [12,13], and greenhouse gas emissions [14,15]. Crops in agricultural areas have varying effects on these environmental issues. Notably, upland crop (e.g., soybean) influences the carbon dioxide cycle [16], and paddy field contributes to methane emissions [17]. Therefore, long time series and high precision cropping pattern dataset is needed to study the differential contribution of crops on environment.

The global datasets such as the 1 km spatial resolution MODIS product from NASA (MOD12Q1) [18,19] are extensively employed by the international scientific community. While the cropland information regarding the crop types is still limited for these products. These datasets have been developed in specific phases and lack multitemporal cropping pattern mapping descriptions. The free accessibility of Landsat imagery (1972 to present [20]) provides opportunities to obtain large-scale and historical land change information. Then, the Chinese Academy of Sciences (CAS) launched a national land-use/cover classification project using Landsat images and auxiliary images (CBERS and HJ-1A). The visual interpretation and digitalization of individual images technology was employed for the whole of China to identify land cover types through interpretation symbols, such as the sizes, colors, textures, and shapes in the remote sensing images [21]. The result of this project, the National Land Cover Dataset of China (NLCD), was used to continuously monitor land-use/cover and the associated changes in China. At present, this dataset has been updated for the 1980s, 1990, 1995, 2000, 2005, 2008, 2010, 2013, and 2015, and it has become key data for land change research at the local, regional, and national scales in China [22,23,24]. A hierarchical classification system was applied in the NLCD [25], in which cropland belongs to the first class. Additionally, cropland is further classified as paddy field and upland crop in the second class. Thus, the dataset has the capability to reveal cropping patterns, while this pattern (mainly for paddy field) has been underestimated or overestimated in some places. Therefore, new cropping pattern data is needed.

In China, agricultural systems are based on two different management models. One model is a special organization that belongs to the central government and is distributed in Heilongjiang, Hainan and Xinjiang Provinces. This organization has extended to a social system with many business enterprises, more than 1500 state-owned farms, approximately 13.22 million people and 6.61 million hectares of cropland area [26,27]. In this region, the land ownership shows that all cropland is “state-owned” with annually detailed land planning. The farm governors can devolve land tenure rights to the workers via a paid contract. Then, farm workers plant cropland crops under the documents of the state and farm governors. The other model is called “household responsibility system”, which involves a 30-year free-paid contract and is implemented in other places under the jurisdiction of local governments. In this region, the farmers freely decide the crop types, management, and benefit rights in the agricultural market. Thus, the state and private farming are independent and based on separate land property right systems. Cold China contains many state and private farming, and the impact of different agricultural systems on cropping pattern needs to be revealed through the new cropping pattern dataset.

In North China, Cold China is a hotspot for land-use/cover changes because of the high-speed extension of agricultural land in conjunction with the large-scale loss of wetlands and unused land. Over the past three centuries, cropland in Cold China has expanded from the south to the north and from the west to the east [28]. In the past century, reclaimed cropland has migrated from the Songnen Plain to the Sanjiang Plain [29]. In the last half century, changes in evolved agricultural policy [30], the increased food demand, and economic returns have led to the northward expansion of cropland in Sanjiang Plain along with changes in the cropland planting structure. Moreover, a large number of state-owned farms are located in this region. Thus, the Sanjiang Plain is characterized by a concentration of cropland expansion and a variety of land ownership systems. This region is highly valuable for analyzing time series of cropland patterns under different agricultural systems.

Most of the current dynamic changes in cropland rarely contain time series descriptions of the cropping patterns (paddy field and upland crop), and this issue limits their value for research on cropland changes and the associated effects on grain ration security, hydrologic processes, climate change, and the ecological environment. A high-precision cropping database covering multiple temporal periods is valuable for identifying cropping patterns, especially patterns occur under different agricultural systems in high latitudes regions. Thus, the objectives of this study are as follows: (1) developing a methodology for creating time series of high precision cropping pattern information through poor-quality observation detection, image co-registration, atmospheric correction, data extraction, information correction, small surface feature identification, and data subsetting and stacking; (2) mapping and further tracking the changes in cropping patterns and trajectory transformations in five-year increments (1990–1995, 1995–2000, 2000–2005, 2005–2010, and 2010–2015) under different agricultural systems in Cold China; and (3) predicting the potential cropping development patterns over ten-year intervals beginning in 2020 and ending in 2050 based on potential demand and patterns.

2. Materials and Methods

2.1. Study Area

Cold China was defined in the northernmost province of China (43°26′–53°33′N and 121°11′–135°05′E). This region is situated in cold temperate zone, with an annual average temperature below 1 °C. Our study area covers the Sanjiang Plain, which spans the border areas of China and Russia between 45°01′05″–48°27′56″N and 130°13′10″–135°05′26″E (Figure 1a,b) [31], and it covers 23 counties (Figure 1c), with a total area of 10.87 × 10⁴ km². Of this area, 68.55% (i.e., 7.45 × 10⁴ km²) is assigned to local governments (here, we define it as “private farming”), and the other 31.45% (i.e., 3.42 × 10⁴ km²) is managed by the central government (the definition of “state farming”) (Figure 1d) based on the administrative extent of the state and private farming (the administrative boundaries of the state and private farming were provided from Northeast Agricultural University, China). The Sanjiang Plain (a total of 52 first-level farms) comprises 51.89% (=3.42 × 10⁴ km²) of the state farming area in China (6.61 × 10⁴ km²); thus, this region is critical in the central government’s efforts to adjust the balance between the grain supply and demand in the Chinese agricultural market.

This study area is built upon sediment deposited by the Amur, Ussuri and Songhua rivers (Figure 1d), which provides convenient water conditions for paddy rice planting. The climate is a moderate temperate humid monsoon zone with an annual average temperature of approximately 1 °C–4 °C and an annual average precipitation of 500–650 mm [32]. The black soil is fit for crops with better organic matter. The landform is flat, with elevations ranging from 17 m to 968 m (Figure 1c).

2.2. Data Collection and Analysis

2.2.1. Landsat Image Collection and Preprocessing

The United States Geological Survey (USGS) allows for free data accessibility, which enables the use of all Landsat imagery. Because good-observations for Landsat images in each month were not always available due to the clouds, bad bands, and shadows, it is difficult or impossible to create an annual cropping pattern map through one year’s Landsat images. Thus, it was perhaps to collect the images from many years for crop mapping. We designed the images in an effort to produce an annual map and measured changes in cropping patterns at the time extension ranges of every three years (periods) (1989–1991, 1994–1996, 1999–2001, 2004–2006, 2009–2011, and 2014–2016, Figure 2). All the Landsat TM/ETM + /OLI images (Total number: 4357 scenes,

Table 1. Statistical information regarding the Landsat images for time series cropping pattern extraction.

Periods	1989–1991	1994–1996	1999–2001	2004–2006	2009–2011	2014–2016	Total
Total Landsat images	794	755	747	690	558	813	4357
Good-observations images	135	154	145	113	100	124	771
Epochs	1990	1995	2000	2005	2010	2015	Total
Extracted images	25	33	54	21	32	36	201
Auxiliary extracted images	5	6	3	13	5	6	38
Verified images	10	11	--	--	--	--	21
Auxiliary verified images	3	3	--	--	--	--	6

) from the USGS (https://glovis.usgs.gov/) during these six periods were first visually examined to evaluate the image quality (e.g., clouds, cloud shadows, snow cover, bad strips, etc.). Then, high-quality image observations from the months between April and September were determined based on phenology and crop calendar, which provided good temporal coverage of growing seasons of paddy field and upland crop [33]. A total number of 771 good-observations images (Table 1) for the scenes (paths/rows: 113/27, 114/26, 114/27, 114/28, 114/29, 115/26, 115/27, 115/28, 115/29, 116/27, 116/28, 116/29,

Table 2. Statistical information regarding the good-observations Landsat images for the path/row in each period.

	1989–1991	1994–1996	1999–2001	2004–2006	2009–2011	2014–2016	Total
113/027	12	11	14	13	9	10	69
114/026	13	16	14	10	5	10	68
114/027	14	14	10	10	8	11	67
114/028	15	13	13	7	6	10	64
114/029	13	14	11	9	6	12	65
115/026	7	12	14	8	9	13	63
115/027	9	12	12	7	12	8	60
115/028	11	9	11	7	8	7	53
115/029	7	11	12	8	9	9	56
116/027	10	16	11	12	8	11	68
116/028	10	14	13	10	11	13	71
116/029	14	12	10	12	9	10	67
Total	135	154	145	113	100	124	771

) were downloaded. The Universal Transverse Mercator coordinate system of the images was re-projected to the Albers Conical Equal Area projection to accurately calculate the dynamic area of cropland in this cold region. Meanwhile, radiometric calibration and atmospheric correction were also performed for all the available Landsat images. In order to obtain accurate annual cropping pattern information, we first focused on Landsat images in the special epochs of 1990, 1995, 2000, 2005, 2010, and 2015. If these extracted images failed to cover cropping pattern map for each epoch, the auxiliary extracted images in other epochs were used to fill the blank area where good-observation images were not available in these special epochs. Based on our statistics, 201 extracted images and 38 auxiliary extracted images (Table 1) were used to extract time series cropping pattern information.

2.2.2. Tracking Cropping Information

Paddy field information acquisition. Paddy field has unique property in that it is transplanted in compound surface environments of water and soil. When rice is transplanted, a period of “flood/open-canopy” regularly occurs and is sustained for a few weeks. Other crops, such as wheat, soybean, and corn, do not require transplantation and are not grown in environments with flooded soil. The rice “flood/open-canopy” stage in compound surface environments can be observed with a spectral sensor. After this period, the optical sensors are unable to obtain the water body information below paddy field canopy. Paddy field enters a period of rapid growth and quickly reaches the maturity and harvest phases, from which the spectral characteristics of the paddy field are consistent with those of other crops. Therefore, the best time to remotely extract paddy field data is during the “flood/open-canopy” period, and pixels of from this stage can be identified if the observations meet one of the following conditions: the land surface water index (LSWI, Equation (1) [34]) > the normalized difference vegetation index (NDVI, Equation (2) [35]) or the LSWI > the enhanced vegetation index (EVI, Equation (3) [36]) [37,38], otherwise, the pixels are not identifiable. Multitemporal paddy planting information was acquired through supervised classification in this stage, because the training sites of this stage in Landsat images displayed obvious symbols of paddy field in terms of the hue (dark blue) and regular size. In addition, flooded natural wetland and water body were removed to improve the accuracy. If this step provided insufficient data to meet the annual cartographic needs, auxiliary Landsat images adjacent to the “flood/open-canopy” stage were used through extraction and artificial visual correction.

$$\mathrm{LSWI}=({\rho }_{green}-{\rho }_{\mathrm{SWIR}})/({\rho }_{green}+{\rho }_{\mathrm{SWIR}})$$

(1)

$$\mathrm{NDVI}=({\rho }_{NIR}-{\rho }_{red})/({\rho }_{NIR}+{\rho }_{red})$$

(2)

$$\mathrm{EVI}=2.5\times \frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR}+6\times {\rho }_{red}-7.5\times {\rho }_{blue}+1}$$

(3)

$$\mathrm{NDWI}=({\rho }_{\mathrm{green}}-{\rho }_{NIR})/({\rho }_{green}+{\rho }_{NIR})$$

(4)

Excluding flooded natural wetland. Paddy field is also called artificial wetland because it is characterized as both wetland and cropland. Thus, differentiating these land cover types after the rice “flood/open-transplanting” stage using an optical sensor is difficult. Natural wetland is always flooded in late spring due to land surface snowmelt, and the natural wetland grows quickly and reaches NDVI values greater than 0.3 before the period of crop planting [39]. We set a threshold of 0.32 to exclude natural wetland in this “flood/open-transplanting” stage using the matching Landsat images.

Multitemporal water body masks. Multiple masks of water bodies were generated to reduce the error associated with cropping pattern maps, particularly for paddy field, due to the similar spectral signatures in the “flood/open-canopy” stage. Remotely extracted water bodies were chosen from August to September based on multiple Landsat images between 1990 and 2015 in 5-year intervals (1990, 1995, 2000, 2005, 2010, and 2015). During these months, crops entered rapid growth and mature stages, and the sensors cannot monitor water information under the canopy in crops. Normalized Difference Water Index (NDWI, Equation (4) [40]) was performed to achieve spatial water information and the threshold was set to be greater than zero.

Upland crop information acquisition. Upland crop information was obtained from the NLCD dataset in five-year increments from 1990 to 2015. The cropland classification was highly accurate [25], while it does not necessarily support the accuracy of upland crop classification. In this research, we assessed the quality of the multitemporal upland crop data by overlaying it with a corresponding 30-m false-color map (R/G/B = NIR/Red/Green) derived from Landsat TM/ETM + /OLI images. A visual inspection was conducted on all images to correct the misinformation using professional knowledge. Because the classification accuracy of visual interpretation and digitalization images was dependent on the professional knowledge of these individuals. Thus, quality control was implemented via repeated interpretations. In this step, 20% of the dynamic vector polygons in each epoch were randomly selected and reinterpreted to check the quality.

Time series cropping pattern excluding tiny surface features. The new time series cropping pattern (paddy field and upland crop) dataset was created through the NLCD dataset and 266 Landsat images over time. The images were processed based on poor-quality observation detection, image co-registration, atmospheric correction, information extraction, information correction, and data subsetting and stacking (Figure 2). The final product was still a 30-m spatial resolution mixed dataset that embedded the tiny surface features of the cropland areas. The multitemporal and large-scale tiny surface features of the cropland, such as canals, ridges, country roads, tree belts, and construction land for hydraulic engineering, were very difficult to achieve due to the large associated workload and restricted historical high-precision images. In this research, tiny surface features of the cropland obtained through the platform of the Chinese Academy of Sciences and local scientific institutions were used to precisely quantify the cropland map for 2015. The resulting map was combined with the extracted dynamic vector polygons of cropland from the new dataset in each period to refine the status maps of cropland in 2010, 2005, 2000, 1995, and 1990. The new cropping pattern maps were then used to analyze the historical cropland changes under different agricultural systems.

2.2.3. Accuracy Assessment

Accuracy evaluation of new dataset was based on the randomly select samples (pixels), and a number of 400 samples were randomly allocated in the whole study area in each epoch (1990, 1995, 2000, 2005, 2010, and 2015). First, because the areas among different land-use types changed in different epochs, it was better to use the stratified random sampling to select samples in each land use type (e.g., paddy field, upland crop, and other types) in each epoch. Second, for epochs of 1990 and 1995, the samples were referred by the color composite images from Landsat TM due to the difficulty in obtaining good quality and high-resolution images before the year of 2000. The color composites (R/G/B = Bands 5/4/3) in Landsat images displayed obvious symbols of paddy field in terms of the hue (dark blue) and regular size, and the upland crop areas could be identified by their uniform texture (Figure 3a–h). The verified Landsat images were separated from the images that was used to obtain time series cropping pattern information, that means, we used some images (the extracted images and auxiliary extracted images, total number: 239 scenes, Table 1) to obtain cropping pattern information and employed other images (the verified images and auxiliary verified images, total number: 27 scenes, Table 1) to conduct accuracy verification. When an auxiliary extracted image was used, the auxiliary verified image in the same location and same year was also used. The verified images and auxiliary verified images could completely cover the study area and provided all the available samples. Thus, we collected 800 samples for the epochs of 1990 and 1995. Third, for the epochs of 2000, 2005, 2010, and 2015, we obtained the historical imagery archive of Google Earth for these four epochs through the professional authorized payment software (software name: 91 bitmap assistant, website: http://www.91weitu.com/). In order to guarantee the temporal consistency between Landsat source data and validation imagery, we filtered the same years of Google Earth imagery with that of Landsat images used for classification. Specifically, the time of Google images focused on the growing seasons of paddy field and upland crop in the special epochs of 2000, 2005, 2010, and 2015. Based on the boundary in research area, all Google images were downloaded in these four special epochs (resolution: 2 m, date: from April 1st to September 30th in each epoch, sensor: QuickBird, vendor: DigitalGlobe & Google). Based on the stratified random sampling for 400 samples in each epoch, we assigned true/faults for each sample and obtained 217 samples, 266 samples, 347 samples, and 382 samples, because other random samples were unavailable due to the spatial inconsistency in Google Earth images. Meanwhile, we conducted field surveys and marked survey information on the maps, accompanied by confirmed photographs (such as: Figure 3e–h) in the spring and autumn of 2015. There were 3139 field photos collected during field investigation, and these photos provided auxiliary information for accuracy evaluation. Finally, for the land system in the accuracy evaluation, the land classification system of the created dataset in the study was the same with that of the NLCD. In this classification system, all land use types were divided into six first class types, which included cropland, woodland, grassland, water body, construction land, and unused land, and cropland was divided into two categories representing paddy field and upland crop. Then, these verified land types were merged into three types (paddy field, upland crop, and noncropland) based on the needs of this study, and the non-cropland contained five first class types (Figure 3i–l).

2.2.4. Trajectory Transformations

A clustering analysis of the land-use types was conducted for the trajectory transformations. These trajectory transformations include continuous paddy field, continuous upland crop, reclaimed paddy field, reclaimed upland crop, paddy field to upland crop, upland crop to paddy field, paddy field to noncropland, and upland crop to noncropland. The discrimination abilities of the trajectory transformations in the five-year increments (1990–1995, 1995–2000, 2000–2005, 2005–2010, and 2010–2015) were quantified through the land-use confusion matrix approach method. Information on the detailed trajectory transformations for different land use categories can be described in

Table 3. Detailed information regarding the trajectory transformations for different categories.

Types	Epoch 1	Epoch 2	Epoch 3	Epoch 4	Epoch 5	Types	Epoch 1	Epoch 2	Epoch 3	Epoch 4	Epoch 5
c p	p0→p1	p1→p2	p2→p3	p3→p4	p4→p5	c u	u0→u1	u1→u2	u2→u3	u3→u4	u4→u5
r p	nc0→p1	nc1→p2	nc2→p3	nc3→p4	nc4→p5	r u	nc0→u1	nc1→u2	nc2→u3	nc3→u4	nc4→u5
p→nc	p0→ nc1	p1→nc2	p2→nc3	p3→nc4	p4→nc5	u→nc	u0→nc1	u1→nc2	u2→nc3	u3→nc4	u4→nc5
p→u	p0→u1	p1→u2	p2→u3	p3→u4	p4→u5	u→p	u0→p1	u1→p2	u2→p3	u3→p4	u4→p5

Notes: p: paddy field; u: upland crop; nc: noncropland; c p: continuous paddy field; r p: reclaimed paddy field; p→nc: paddy field to noncropland; p→u: paddy field to upland crop; c u: continuous upland crop; r p: reclaimed upland crop; u→nc: upland crop to noncropland; u→p: upland crop to paddy field; period i (i = 1, 2, 3, 4, 5) represents the time period stages of 1990–1995, 1995–2000, 2000–2005, 2005–2010, and 2010–2015, respectively.

.

2.2.5. Cropping Pattern Prediction

The Cellular Automata/Markov Change prediction model (CA_MARKOV) is a combined Cellular Automata, Markov Chain, Multi-Criteria, and Multi-Objective Land Allocation land cover prediction system [41]. CA_MARKOV allocates land based on the suitability of the land for end covers along with a cellular automaton rule to promote spatial contiguity. It works well when historical land cover data is not available or is not a good predictor of future land cover [42,43]. Thus, we used CA_MARKOV to extrapolate the cropping patterns for the epochs of 2020, 2030, 2040, and 2050, and further provided some information that CA_MARKOV was well used in the satellite remote sensing analysis domain. In CA_MARKOV, the Markov chain-predicted transitions are usually generated based on only two training land use maps. This approach may provide an inadequate number of training samples. We applied the LUCC budget technique [44,45] to obtain the land changes and the associated component changes, including the total change, gain, loss, and transition, according to the historical land use maps. These training maps were used to determine the Markov chain transition matrices through equal weight coefficients. Based on the Markov chain transition matrices, the Markovian Transition Estimator was used to obtain suitability images for each of the land cover types. In this process, the multiple influencing factor maps in the study area, including the natural factors (streams, slopes, and dem) and socio-economic factors (urban areas, towns, roads, and railways), were added to rectify suitability images, and the Multi-criteria Evaluation (MCE) was conducted to further regulate spatial distribution of suitability images through the thematic maps, including water distribution, ecological protection areas, natural reserves, and land use planning schemes in the study area. Based on the suitability images, the issues with concurrent transitions were solved through the combination of multi-objective land allocation (MOLA) and cellular automata (CA). Then, the predicted maps can be created. In order to test the accuracy of this method in the study area, we utilized the 1990–2010 dataset to predict the map of 2015 (comparison data) and compared it against the real map of 2015 (reference data) in our new dataset (Figure 4). The verification was conducted in the terrset geospatial monitoring and modeling system, and the accuracy verification indicators included quantitative consistency (Kappa location (Klocation)) and spatial consistency (Kappa standard (Kstandard)). The result exhibited high accuracy (Klocation: 0.9714; Kstandard: 0.9302) in the test involving the whole study area. Thus, the method was suitable for land use prediction in this region. We used the method to obtain the suitability images of each land cover type in each epoch (Figure 5) and extrapolated the land prediction maps for the epochs of 2020, 2030, 2040, and 2050.

3. Results

3.1. Validation Accuracy of the New Cropping Pattern Maps

Accuracy assessments of time series of cropping pattern maps were conducted based on the randomly select samples (pixels). The results for the cropland, including paddy field and upland crop, and noncropland maps were highly accurate (

Table 4. Confusion matrix of land-use validation based on Landsat TM and Google Earth imagery in the special epochs of 1990, 1995, 2000, 2005, 2010, and 2015. The overall accuracy (OA), user accuracy, producer accuracy, and kappa coefficients are also provided.

		Ground Truth (GT) Samples (Pixels)			Total Classified Pixels	User’s Accuracy
Year	Land-Use Type	Paddy Field	Upland Crop	Non-Cropland
1990	Paddy field	89	4	7	100	89.00%
	Upland crop	3	151	10	164	92.07%
	Non-cropland	5	13	118	136	86.76%
	Total GT pixels	97	168	135	400	OA = 89.50%
	Producer’s Accuracy	91.75%	89.88%	87.41%		Kappa = 0.87
1995	Paddy field	126	4	7	137	91.97%
	Upland crop	1	112	9	122	91.80%
	Non-cropland	9	5	127	141	90.07%
	Total GT pixels	136	121	143	400	OA = 91.25%
	Producer’s Accuracy	92.65%	92.56%	88.81%		Kappa = 0.88
2000	Paddy field	54	2	4	60	90.00%
	Upland crop	1	81	4	86	94.19%
	Non-cropland	3	6	62	72	86.11%
	Total GT pixels	58	89	69	217	OA = 90.78%
	Producer’s Accuracy	93.10%	91.01%	89.86%		Kappa = 0.86
2005	Paddy field	78	2	4	84	92.86%
	Upland crop	1	69	5	75	92.00%
	Non-cropland	6	4	97	107	90.65%
	Total GT pixels	85	75	107	266	OA = 91.78%
	Producer’s Accuracy	91.76%	92.00%	91.59%		Kappa = 0.88
2010	Paddy field	134	5	5	144	93.06%
	Upland crop	1	93	6	100	93.00%
	Non-cropland	7	3	93	103	90.29%
	Total GT pixels	142	101	104	347	OA = 92.22%
	Producer’s Accuracy	94.37%	92.08%	89.42%		Kappa = 0.88
2015	Paddy field	153	3	5	161	95.03%
	Upland crop	3	106	4	113	93.81%
	Non-cropland	6	2	100	108	92.59%
	Total GT pixels	162	111	109	382	OA = 93.98%
	Producer’s Accuracy	94.44%	95.50%	91.74%		Kappa = 0.89

). The overall accuracies were 89.50%, 91.25%, 90.78%, 91.78%, 92.22%, and 93.98% in the different epochs, with corresponding kappa coefficients of 0.87, 0.88, 0.86, 0.88, 0.88, and 0.89, respectively. This finding indicated that the final epoch had the highest overall accuracy and kappa coefficient because of the improved Landsat image intensity and the availability of high-quality Google Earth. The paddy field displayed higher producer accuracy (PA) and user accuracy (UA) in 2015 (PA, 94.44%; UA, 95.03%). The year 1990 exhibited the lowest PA (91.75%) and UA (89.00%), and this result was attributable to the limited intensity of Landsat data. The precision of upland crop classification was similar with that of the paddy field. The time series maps of cropland had high accuracies and could be used for mapping cropping patterns.

3.2. Comparison of the Cropping Patterns between NLCD-Based and New-Based Datasets

Time series maps of cropland in six epochs were generated using a new dataset (Figure 6(a2)–(f2)). The cropland planting area constantly increased from 1990 to 2015, with a total cropland area of 41,167.26 km² in the first year and 52,507.57 km² in the final year. During this period, both the NLCD maps (Figure 6(a1)–(g1)) and new dataset (Figure 6(a2)–(g2)) exhibited a continuous and consistent increasing trend for cropland area, despite the differences in the cropping pattern in each epoch. Based on the new dataset, most of cropland area in this region was located on flat plains, with low-lying landscapes and high proportions of swamps, rivers and ditches. Suitable cropland area was generally already cultivated, and further land reclamation was and will continue to be limited (Figure 1d and Figure 6).

The paddy field expanded persistently from the study period. Based on the new dataset, in the first year, paddy field covered very limited areas (2850.87 km² in total). In these areas, the spatial distribution of paddy field showed in a scattered manner to Songhua River and valley plain along the Khanka. After 1990, the paddy planting area acutely increased and reached 26,537.14 km² by the final year. The area covered by paddy field in 2015 was 9.3 times greater than that in 1990. During this period, paddy field area reached 4691.24 km², 9539.90 km², 12,586.86 km², and 18,866.99 km² in epochs of 1995, 2000, 2005, and 2010, respectively (Figure 6(g2)). Paddy expansion first occurred in the western part of the plain along the Songhua river (Figure 6(a2)) and then rapidly increased along the Songhua, Amur, Ussuri and Khanka rivers, where the water resources were abundant near the tributaries and streams of these rivers (Figure 6(b2)–(f2)). Therefore, paddy expansion was mainly based on the accessibility of water resources. Both the NLCD maps and new dataset displayed a continuous increasing trend for the paddy field area in six epochs, while the paddy expansion patterns exhibited apparent inconsistency in each period in both datasets (Figure 6(a1) vs. Figure 6(a2); Figure 6(b1) vs. Figure 6(b2); Figure 6(c1) vs. Figure 6(c2); Figure 6(d1) vs. Figure 6(d2); Figure 6(e1) vs. Figure 6(e2); and Figure 6(f1) vs. Figure 6(f2)). The new dataset had better spatial pattern information for upland crop and paddy field.

3.3. Different Cropping Patterns in Different Agricultural Systems

The area and the percentage of cropland were used to compare different changes under different agricultural systems. The area changes were calculated based on the area in latter epoch subtracted the area in previous epoch, and the percentage were obtained based on the total croplands occupied by paddy field or by upland crop in each epoch in state farming or private farming. The net cropland expansion area on state farming (+6709.64 km²) was faster than that in private farming (+4630.68 km²) during the period of 1990–2015, (Figure 7a,b). Meanwhile, the cropping pattern differed significantly in both regions over time. In private farming (Figure 7a), the percentage of upland crop decreased from 89.88 to 65.02%, with a loss of area of 2840.29 km², and the percentage of paddy field increased from 10.12 to 34.98%. Thus, the area covered by paddy field increased by 7470.96 km². While in state farming (Figure 7b), the changes of cropping patterns (including upland crop and paddy field) in state farming changed more rapidly than that in private farming, with the percentage of upland crop obviously decreased from 97.34 to 31.44%, which led to a large loss of area of 9505.66 km²; moreover, the paddy area, as well as the corresponding percentage of total cropland, linearly increased from 467.99 km² (2.66%) to 16,683.30 km² (68.56%), thus reflecting an increase of 16,215.30 km².

3.4. Trajectory Transformations of Cropping Pattern at the Pixel Level in Different Agricultural Systems

Maps of trajectory transformations of cropping patterns, as well as the cropping pattern and noncropland times, were generated based on the new data from all epochs in the study area (Figure 8a,b). For paddy field, the paddy field maps displayed a successive and steadily increasing pattern at the pixel level (Figure 8a). Increase in paddy area was mainly attributable to the transformation of upland crop, which accounted for 83.47%, 90.31%, 85.22%, 76.77% and 73.38% of the reclaimed paddy field in five periods (

Table 5. Trajectory transformations of land use for the periods of 1990–1995, 1995–2000, 2000–2005, 2005–2010, and 2010–2015 (km²).

		1990–1995	1995–2000	2000–2005	2005–2010	2010–2015
Cropland	c p	2785	4643	9522	12,544	18852
	p→u	55	34	4	33	13
	c u	35,748	35,045	34,309	30,563	25,842
	u→p	1590	4421	2612	4853	5638
Cropland and Noncropland	nc→p	314	474	453	1468	2045
	nc→u	4238	2493	1371	932	114
	p→nc	9	13	14	10	5
	u→nc	977	575	651	266	48

Notes: p: paddy field; u: upland crop; nc: noncropland; c p: continuous paddy field; p→nc: paddy field to noncropland; p→u: paddy field to upland crop; c u: continuous upland crop; u→nc: upland crop to noncropland; u→p: upland crop to paddy field.

). The areas of abandoned paddy field, including the transformation of paddy field to upland crop and paddy field to noncropland, were limited. For upland crop, the trajectory transformation of continuous upland crop was variable at the pixel level (Figure 8b). Increase in upland crop area was mainly due to the transformation of noncropland. This transformation accounted for 98.70%, 98.66%, 99.68%, 96.57% and 89.60% of reclaimed upland crop in five periods. Moreover, areas of abandoned upland crop were converted to paddy field on a large scale. For noncropland, the areas of abandoned cropland converted to paddy field had an increasing trend, with areas of 314 km², 474 km², 453 km², 1468 km², and 2045 km² in each period (Table 5). While the opposite trend can be observed in the conversion of abandoned cropland converted to upland crop, with areas of 4238 km², 2493 km², 1371 km², 932 km², and 114 km² in each period.

Trajectory transformations of continuous paddy field was inconsistent under different agricultural systems (Figure 8c). In 1990, most paddy field displayed an aggregate distribution in private farming close to the Songhua river and the valley plain (Figure 8(a2)). Additionally, there was very little distribution of crop fields that had paddy field in state farming. Then, on state-run farms, paddy field have acutely expanded to neighboring crops since the year of 2000 (Figure 6 and Figure 8a). So, most paddy field had cropland ages of no more than twenty years following continuous cultivation during 1990–2015 (Figure 8(c3)–(c5)). This situation was also demonstrated in time series cropland maps among six epochs (Figure 6(a2)–(f2)). While in private farming, the reclaimed paddy field lagged behind than that of state-run farming (Figure 8c), although most of the paddy fields were first cultivated more than 25 years before and have been continuous cultivated. A common feature of trajectory transformations of continuous paddy field under different agricultural systems is that paddy transformations presented a northward migration, such as in the directions of c1→c2 for private farming and c3→c4→c5 for state farming, due to the climate warming, increased food demand, and better economic profit for paddy field.

3.5. Cropping Patterns for the Period of 2020–2050

The cropping pattern maps of 2020, 2030, 2040, and 2050 were created based on the prediction (Figure 9a–d). Cropland planting area in this region was predicted to increase by 2602.80 km² over the thirty-five years, with net changes of +1224.68 km², +839.57 km², +565.10 km², −26.55 km² in 2015–2020, 2020–2030, 2030–2040 and 2040–2050, respectively. The cropland reclamation area in the future period is much smaller than that in the past period (1990–2015, with a total increase of 11,340.31 km²). It can be noted that cropland growth from 2015 to 2020 accounted for nearly half (47.05%) of the predicted net change. Thus, rapid cropland reclamation is expected to stagnate after the epoch of 2020 due to terrain and spatial land use planning restrictions, such as areas of ecological protection and natural reserves. The paddy planting area is predicted to expand continuously by 7361.25 km², 9852.57 km², 3932.54 km², and 781.15 km² in each period (Figure 9e–h). This result indicates that rapid paddy expansion become slowness in 2040. Future paddy expansion pattern is first observed in the flat areas neighboring the continuous paddy field along the Songhua, Amur, Ussuri and Khanka Rivers and then extend to several tributaries of these rivers and other areas.

From the perspective of agricultural systems, incremental areas covered by paddy field are still faster in state farming in the period of 2015–2020, while the opposite phenomenon is observed in private farming after that period. In contrast, the areas covered by upland crop will continuously decrease by 6136.57 km², 9013.00 km², 3367.44 km², and 807.70 km² in each period, with a similar trend in both land ownership regions. Moreover, greater losses of upland crop will occur in private farming that in state farming beginning in 2020.

4. Discussion

4.1. First Region for Continuous and Large-Scale Land Use Transformations from Upland Crop to Paddy Field in China

This study first reported the dynamic tracing of cropping patterns and showed the reasons for paddy expansion in a region of northern Asia. We found that cropping patterns have changed significantly, and continuous rice paddy expansion has occurred, and the land source for paddy expansion was mainly attributed to the large-scale transformations of upland crop. Although a considerable number of paddy fields are located in northern Asia such as northern Japan, northwestern Korea, and China, while dynamic tracking of paddy expansion has not been performed and the source has not been identified in northern Asia. Therefore, we first provide new satellite evidence for these issues in the study area of Cold China. This study showed that continuous and large-scale land use transformations from upland crop to paddy field occurred throughout Cold China, and it represent a new and acute process of land use change in China based on our national survey. At present, China has a large population and limited cropland area [46]. Given the rapid speed of industrialization, urbanization and rural development in the 21st century, the total cropland area in China has decreased [47,48]. This cropland evolution pattern has been characterized by a substantial loss of cropland areas and a reduction in the cropping intensity (such as the cropping intensity from triple cropping to double cropping; and from double cropping to single cropping) in South China, or a return to forest and grassland via ecological restoration projects in Central China. Conversely, the cropping area of single crops has continuously increased in North China, especially in Cold China. Thus, cropland expansion (net change area of +11,340.31 km² in 1990–2015, Figure 6) in the study area has played an essential role in the maintenance of cropland area in China.

This study showed that cropping patterns in Cold China have varied over time. That is, the area of upland crop shrunk following a large increase, and the area of paddy field continuously expanded during the period of 1990–2015 (Figure 6). This region has been characterized by the rapid expansion of agricultural land, particularly paddy field (net change of +23,686.27 km² from 1990 to 2015). Therefore, the Cold China has gradually become an important planting base for paddy field in China. Currently, suitable cropping areas are generally already cultivated in areas of low-lying terrain. Further land reclamation will be very limited. Therefore, the cropping pattern displays a new land use phenomenon in which large areas of upland crop are converted to paddy field (Figure 8). This phenomenon first occurred in Cold China due to the large number of state-run farms and appropriate natural resources, such as dense rivers, black soil, and flat terrain. Paddy development in Cold China will continue to guide cropping pattern reformation (Figure 9) in China because this region includes more than half of the state-owned farm areas, which is conducive to the role of the state in adjusting the crop structure in the agricultural market. Meanwhile, paddy field provide a major dietary staple [49] for over 1.3 billion Chinese people and for more than half of the global population.

4.2. Reasons for Different Paddy Patterns under Different Agricultural Systems in China

Agricultural systems are based on two different models in China, and the results presented here indicate that paddy field have expanded at a faster rate and more regularly in state farming than in private farming (Figure 7 and Figure 9). The considerable differences in cropping patterns in both regions can be summarized based on the following factors. (1) Land management regimes. All cropland on state farming is state-owned, and the central government generates annual detailed crop planting plans to adjust the balance in the national agricultural market. According to state regulations, farm governors own land rights and devolve the rights to workers. These workers manage cropland planting, such as paddy field planting in a specific region, under the guidance of the central government and farm governors. While the cropland in private farming is collective-owned and implemented by farmers. Farmers can freely decide the crop planting types and control the land use, management, and benefit right strategies; (2) Agricultural policy. The purpose of agricultural policy has changed over time in China. The goal of the first stage was to provide the basic food for over 1.3 billion Chinese people, and these policies were to promote cropland areas. The “Food First” policy in the First Five-Year Plan (1953–1957) [50] increased the land reclamation area. The “Agricultural Modernization” policy (1978–1985) [51,52] promoted advanced agricultural machinery, and many modernized farms were formed. In the second stage, the goal was to improve people’s agricultural living standard. Although subsequent agricultural policies, such as “Layout Planning for the Advantageous Regions of Rice (2008–2015),” “Sanjiang Plain Land Renovation Project (since 2008)” and “Building High-Standard Basic Farmland”, also generated large-scale cropland reclamation in Cold China; however, these policies focused on cropland quality and cropping patterns (mainly for paddy field). A total of 80–90% of these projects were located in state farming, which effectively improved the agricultural infrastructure and transformed scattered paddy field into an aggregated pattern. Whereas the paddy pattern is scattered in private farming due to the lack of necessary engineering facilities; (3) Restricted water resources and cropland area. Water resources are an important limiting factor for paddy field planting. In state farming, advanced irrigation technology and excellent water facilities provide sufficient water environment for rice planting. Paddy expansion exhibits a “push-and-pull” pattern under surface canal irrigation. However, due to insufficient surface irrigation facilities, paddy expansion in private farming occurs primarily due to the spontaneous drilling of growers, leading to a small-scale and sporadic pattern. Another restricting factor is the cropland area. In state farming, which is characterized by abundant cropland areas and scarce populations, workers obtain cropland via paid annual contracts from farm governors, and they can purchase more land as desired. Notably, paddy field can easily be planted on a large scale (approximately 30–40 hectares/household). However, in private farming, which is characterized by a limited cropland area and a large population, farmers obtain cropland via a 30-year free paid contract, and obtaining more cropland is difficult or impossible. Moreover, a paddy field in private farming can only be planted on a small scale (approximately 2–6 hectares/household) based on our survey.

4.3. Driving of Physical Conditions and Human Factors on New Cropping Pattern Changes in Cold and High Latitudes of China

The driving forces underlying the new land use changes can be summarized based on two aspects. The first aspect is the physical conditions, such as water resources, soil resources, and topography. Paddy expansion is mainly based on the accessibility of water resources. The region is built upon sediment deposited by the Amur, Ussuri and Songhua rivers, which provides convenient water conditions for paddy rice planting. This region is one of the three black soil regions on Earth, and the black soil is more fit for rice planting because of its good organic matter content and high fertility. The terrain is flat, which is propitious to water storage in paddy patches. The second aspect is human factors, such as agricultural policy, agricultural technological progress, and crop economic returns. This region covers more than half of the cropland area of state-run farms in China, and many agricultural policies and water conservancy projects are focused on state farming regions, which improves the infrastructure in paddy field. Agricultural technological progress increases the adaptability of seeds and the efficiency of fertilization in cold regions; and thus, promotes rice paddy expansion. The economic returns that the net income of paddy field (average more than 8500 yuan/hm²) is higher than that of upland crop (average less than 6500 yuan/hm²) is another human factor underlying paddy expansion in cold regions.

4.4. Effect of Cropping Pattern Changes on the Environment in Cold Region

Cold China is a commodity grain production base in China and was covered by a large area of upland crop in 1990. Continuous and acute land use change from upland crop (mainly soybean) to paddy field will inevitably affect the terrestrial carbon cycle and lead to the large-scale loss of carbon dioxide in this cold region [53,54,55,56]. Additionally, paddy expansion influences the surface energy balance and greenhouse gas emissions in addition to virus transmission due to its characteristics as both wetland and cropland. Paddy field planting obviously reduces the surface temperature and changes the Bowen ratio during the growing season; and it is also a significant source of methane flux into the atmosphere and contributes more than ten percent of worldwide methane emissions. As a vital habitat for wild animals, such as free-ranging ducks and waterfowl species, paddy field is related to the spread of avian influenza viruses and other diseases. By growing in a flooded soil environment, paddy expansion may have the most apparent influence on water resource utilization. Rapid and acute paddy expansion in Cold China accelerated river water irrigation and groundwater overexploitation, and such change may increase the risk of groundwater depletion to a certain extent. Based on the study, the paddy expansion will continue in Cold China after 2020, especially in private farming region. In these regions, surface water resources are not abundant compared to those in state farming. A comprehensive evaluation of the carrying capacity of agricultural water resources is needed to provide an early scientific warning and coordinate the water resources and paddy planting area in cold regions. Meanwhile, with the exception of environmental issues involving the changes in paddy field expansion and upland crop shrinkage, the issues caused by the loss of noncropland should not be ignored. The expansion in cropland was mainly at the loss of woodland, grassland and unused land in noncropland based on our survey. Thus, these issues were multiple fields, such as the loss in carbon source/sink in woodland and the loss of biodiversity in unused land in ecology, and the changes in air temperature and evapotranspiration among different land use types in climatology.

4.5. Uncertainty of Crop Information Extraction and Future Prediction

The main uncertainties in this research was related to two factors. (1) The availability of good-observations in Landsat images. The accuracy of land use information acquisition depends on the quality of remote sensing images. Good-observations in Landsat images can comprehensively and accurately track dynamic information. Therefore, the cropping pattern based on remote sensing information during 1989–2016 is mainly limited by image quality in this research. Specifically, the rice “flood/open-canopy” stage is short. With restricted number of images per year (accompanied by 16-day revisit cycle with Landsat data) and uncertain image quality (clouds, bad strips, etc.), establishing an annual cropping pattern map via the Landsat images for one year is difficult or impossible. Thus, this research designed images to create an annual map and measured the changes in cropping pattern in the time resolution over three years. If the data in stage of rice” flood/open-canopy” failed to meet the cartographic requirements in the special epochs of 1990, 1995, 2000, 2005, 2010, and 2015, we extended the data time throughout the growing season in these special epochs. If the maps still failed to cover the whole study in each epoch, the auxiliary Landsat images adjacent to special epochs were used to fill the blank area where good-observation images were not available in these special epochs. Finally, 222 images in special epochs and 44 images in other epochs (Table 1) were used to finish time series cropping pattern mapping; (2) The uncertainty for cropping pattern prediction. Based on the past changes of cropping pattern areas, a potential and reliable future development scenario was given in this study. To reduce the uncertainty of land prediction, we used the LUCC budget technique to obtain the relations between the nonstationary processes and predictability needs based on Markov chain transition matrices through various training maps in each period. Moreover, the natural factors, social and economic factors, and thematic maps were added to accurately predict the suitability images of each land cover type in each epoch. In this scenario, the information regarding paddy field planning and water consumption was not included due to the sufficient water resources (total storage: 202.31 × 10⁸ m³; annual average precipitation: 567.50 × 10⁸ m³; annual average flow: 2634.26 × 10⁸ m³) in the study area. Given the apparent influence of paddy expansion on water resource utilization, the assessments of paddy field planting areas and available water resources are needed in future research. Meanwhile, the shortcoming in future prediction is that this study didn’t take into account the impacts of possible changes in climate and food consumption preferences on cropping patterns in the future, such as the new crop that will fit better for this region climate and will increase the growers’ profit in Cold China. These new crops (the exclusion of paddy field) also have the possibility to expand in the future. Because the changes in cropland are often nonlinear (e.g., progressive or regressive, slow or fast), and cropland often displays a high degree of space-time complexity and complicated feedback mechanisms [57,58,59]. However, due to the complexity of changes in land use, climate, and human activity, it is difficult to make reliable multi-scenario prediction at present. These issues will be our important research directions in the future. We’ll report them immediately as soon as we get any signals of the new crops.

5. Conclusions

Large-scale applications of time series Landsat images in cropland change studies beginning in 1990 face various challenges, such as variable image quality (e.g., clouds and bad strips) and inhomogeneous data availability at spatial and temporal scales. Based on time series of Landsat images, this study attempts to track cropping information (1990–2015) and predicts the future pattern (2015–2050) under different agricultural systems in Cold China by developing an updated method for determining cropland patterns. The results demonstrated that the available Landsat images in Cold China met the requirements for long-term cropping pattern studies; and the developed method exhibited high accuracy and good spatial pattern results. Based on the new cropping pattern dataset, the cropping pattern significantly differed in both agricultural systems, with paddy field in state farming expanding at a faster rate (from 2.66 to 68.56%) than those in private farming (from 10.12 to 34.98%). Paddy expansion was mainly attributable to the transformation of upland crop at the pixel level, which led to a greater loss of upland crops in state farming than private farming (9505.66 km² vs. 2840.29 km²) from 1990 to 2015. Meanwhile, the loss in non-cropland was first mainly converted to upland crop and then primarily converted to paddy field. Based on the prediction, rapid cropland reclamation will be projected to stagnate in 2020, while paddy expansion will continue until 2040 primarily in private farming areas in Cold China. This study provides new satellite evidence showing that different agricultural systems significantly affected cropping patterns through rice paddy expansion, and the results have significant implications for understanding and guiding agroecosystem development.