Spatiotemporal Dynamics of the Northern Limit of Winter Wheat in China Using MODIS Time Series Images

: Studying the spatiotemporal changes of the northern limit of winter wheat (NLWW) in China is important to ensure regional food security and deal with the effects of climate change. Previous studies mainly used climate indicators to analyze the variation of the potential NLWW in different historical periods, while little attention has been paid to the actual migrations and changes of the NLWW. The objectives of the present study were three-fold: (i) to map the spatial distribution of winter wheat in northern China in 2001, 2007, 2014 and 2019; (ii) to extract the actual NLWW; and (iii) to quantitatively explore the dynamics of the NLWW. First, we adopted the “combining variations before and after estimated heading dates” method to map the winter wheat in northern China based on time series MODIS EVI2 data. Second, we used the kernel density estimation algorithm to extract the actual NLWW in four historical periods. Finally, the fishnet method was utilized to quantitatively analyze the direction and distance of the spatiotemporal changes of the NLWW. The results demonstrated that the NLWW has exhibited a marked fluctuating trend of migration southward, with a 37-km shift in latitude over the past 20 years. The elevation limit of winter wheat planting was around 1600 m; however, the centroid of winter wheat planting has shifted slowly to lower elevations. There was a gap between the actual NLWW and the potential NLWW. The reason for this gap was that the actual NLWW moved southward under the interacting effects of human activities and climate change, while the potential NLWW moved northward due to climate change. The results of this study are of great scientific value in the formulation of winter wheat planting strategies in climate-sensitive areas to respond to climate change and ensure food security.


Introduction
Winter wheat is one of the most important food crops for ensuring food security, and is planted extensively worldwide, providing approximately 20% of human calorie demand [1,2]. Climate change is significantly altering the spatial distribution of winter wheat [3,4]. Globally, climate change in northern China is more significant, and the pattern of winter wheat planting is prone to spatiotemporal changes [5,6]. The northern limit of winter wheat (NLWW) is the northern boundary of winter wheat distribution in China where humans can engage in safe winter wheat planting. Understanding the spatiotemporal dynamics of the NLWW is an important foundation for an effective climate change response and winter wheat production strategy formulation [7,8].
The meteorological data offered insight into the response of the NLWW to climate change [9]. The potential NLWW was determined by agricultural climate indicators calculated from meteorological data. Previous studies have demonstrated that the potential NLWW in China has shown a clear northward shift in the historical period, and this shift will continue in the future [10]. This indicates that the area suitable for winter wheat planting has increased [11]. However, the planting area of winter wheat in northern China decreased significantly on the provincial (municipal) statistical level (http://www.stats.gov.cn/). Therefore, a gap may exist between the potential NLWW and the actual NLWW in China.
At present, the studies on the NLWW focus on the potential NLWW shifts, while studies on actual NLWW dynamics are rare. The distribution of winter wheat is essential for mapping the actual NLWW. Remote sensing techniques have gradually emerged as an important scientific and technical method for monitoring the winter wheat planting area [12,13],and they offer the possibility of obtaining the actual NLWW through spatial statistical methods [9]. The remote sensing mapping of crops can use image statistics-based approaches in certain stages, and use time series remote sensing data through different classification algorithms [14]. However, there are many challenges to these methods, including a large demand for training samples, the high cost of sample collection and the time-consuming and laborious processing of vast image datasets.
Phenology-based methods can effectively avoid these problems based on the unique characteristic during the growing season [15]. Crop mapping was achieved by using a time series remote sensing vegetation index to calculate the similarities of standard vegetation index time profiles for different crops [16]. However, the time profile of the crop vegetation index usually shows considerable intraclass variability [17], which poses great challenges to crop classification [18]. Due to the long growth and development period of winter wheat, the characteristics of phenological changes reflected in winter wheat are clearly unique, and the key phenological period, with multiple characteristic peaks and troughs, can be used to identify winter wheat [19]. The use of the time series remote sensing data to better reflect the change law of winter wheat phenology has become the mainstream technical method [18].
To improve the mapping accuracy of winter wheat, phenological information of high temporal resolution was formed from the image data of moderate-high spatial resolution. For example, Landsat, RapidEye and Sentinel-2 images were integrated to fully use the information derived from the high-spatiotemporal-resolution data of winter wheat to identify winter wheat [20,21]. Additionally, combining the characteristics of synthetic aperture radar and optical data can improve the classification of winter wheat [22,23]. The mapping accuracy of winter wheat can be slightly improved through the integration of various remote sensing data, but these data are limited by a short time period and the difficulties associated with obtaining high-quality data for a long time period.
MODIS can provide sufficient phenology information due to its high temporal resolution (1-day temporal resolution) [24,25]. Thus, the characteristics of the growth phenology of crops can be highly reflected, and this approach is widely used in crop mapping [26,27]. From small to large spatial research scales, studies have employed MODIS time series data for the mapping of maize, paddy rice, winter wheat and soybean planting areas [28][29][30]. This study used MODIS time series data combined with a phenology-based method to identify the spatial distribution of winter wheat, and further extracted the actual NLWW, then analyzed the reasons for the gap between the actual NLWW and the potential NLWW. The objectives of this study were the following: (i) to determine the distribution of winter wheat in northern China during the period 2001-2019 using MODIS EVI2 time series images; (ii) to detect the actual NLWW using kernel density estimation (KDE) based on the generated winter wheat maps; and (iii) to investigate the spatiotemporal dynamics of the actual NLWW during the period 2001-2019.

Study Area
The study area encompasses Hebei, Shanxi, Shaanxi and Gansu Provinces, in addition to the Beijing and Tianjin municipalities, covering approximately 380,000 km 2 between 34.33-41.00°N and 106.40-119.82°E (Figure 1). The area's characteristics include the climatic sensitivity and strong spatiotemporal fluctuation of the ecological environment, with severe cold in winter and limited precipitation. The accumulated temperature greater than 0 °C during the winter wheat growth stage is approximately 2200 °C, and the absolute minimum temperature is around −24 °C. Annual precipitation is 440-660 mm. Generally, the winter wheat growth period extends from late September to early July of the following year. The elevation of the study area is within the range of 0-3059 m above sea level, and the terrain is reasonably complex. The research area is located in the Chinese farming-grazing ecotone, and the agricultural system belongs to a mixed agriculture system [31]. The cropping system is dominated by single-cropping, and the area of double-cropping is increasing [32]. The wheat planting area (including winter wheat and spring wheat) is approximately 22,600 km 2 , accounting for 8% of China's wheat area. The crops in the study area are varied, but consist predominantly of wheat and miscellaneous grains. The planting area of winter wheat accounts for 30-40% of the planting area of grain crops, and other crops include maize, sorghum and potato [33].

Datasets and Preprocessing
The data used in this study include MODIS EVI2 datasets, Landsat images and sampling point data. To extract actual NLWW, the planting areas of winter wheat in the study areas during the four periods are mapped based on the data of the MODIS time series EVI2. Landsat data are used to map winter wheat in the study area in order to evaluate the mapping accuracy of MODIS-derived winter wheat. Sampling point data are used to determine constants for MODIS-derived mapping, as well as perform accuracy evaluations of the MODIS-derived and Landsat-derived winter wheat maps.  Figure 1). Four tiles (h26v04, h26v05, h27v04 and h27v05), which covered the entire study area and spanned the winter wheat growing season, were acquired from the NASA website (https://ladsweb.modaps.eosdis.nasa.gov/). The MODIS Reprojection Tool was used for the projection, conversion and assembly of the images. The data were processed uniformly into the TIF format, and Albers projection WGS84 coordinates were adopted. The MOD09Q1 product was utilized to calculate the two-band Enhanced Vegetation Index (EVI2) (Equation (1)) time series datasets [34].
where EVI2 is the two-band enhanced vegetation index, NIR1 is the reflectance of the near-infrared band, and RED is the reflectance of the red band. Compared with the normalized difference vegetation index (NDVI), the EVI2 is more suitable for mapping crop areas. Additionally, EVI2 is less susceptible than NDVI to biases resulting from cloud and haze contamination, and can better reflect vegetation information [34]. The MOD09Q1 EVI2 time series images underwent preprocessing, which comprised calibration, correction, cloud removal and filtering, prior to use for winter wheat mapping. Specifically, the pre-processing procedure for MODIS data included three components: First, the EVI2 was calculated through red and near-infrared wavelengths. Second, we used the quality control flag layer in the MOD09Q1 products to extract the clouds and cloud shadows from each image. Third, the harmonic analysis of the time series was used to reconstruct the essential curve of the EVI2 [35].

Landsat Images
Landsat data were used to map winter wheat in the study area, so as to evaluate the mapping accuracy of MODIS-derived winter wheat. Landsat images for the Path 123/Row 033 (123033) tile area were used for winter wheat extraction ( Figure 1 The imaging times broadly covered the period of the winter wheat growing season, which increased the ease of distinguishing the winter wheat crop. The specific acquisition dates of the images and the concurrent extent of cloud cover are listed in Table 1. The data were downloaded from the U.S. Geological Survey website (https://glovis.usgs.gov/). Prior to information extraction, all aforementioned remote sensing data were preprocessed (e.g., radiation correction, atmospheric correction and geometric correction) using ENVI TM 5.3 to eliminate possible errors. Specifically, the pre-processing procedure for Landsat images included three components: First, the digital number values of the Landsat images were converted to radiance values. Second, the land surface reflectance was derived using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) model. Third, the Landsat TM and ETM+ images were geometrically corrected using Landsat OLI as the reference standard. For the accurate extraction of winter wheat data, Landsat images from three periods were chosen to characterize the winter wheat phenology changes. These were realized through the following three steps: First, in the tillering-stage images, the near-infrared band is the clearest band for winter wheat with regards to distinguishing other non-vegetative features, and the reflectivity of this band can distinguish vegetation and non-vegetation. Second, according to the phenological changes of the vegetation index, shifting "up-down" before and after the heading period of winter wheat, the characteristics that make NDVI increase (from the tillering to the heading period) and NDVI decrease (from the heading to the harvesting period) can be used to mask winter crops, forests and nonvegetation types. Finally, after the identified types above were eliminated, winter wheat was extracted by the Random Forest classification method according to the remote sensing images of the winter wheat heading period [36].

Sampling Point Data
We used UniStrong 50A handheld GPS receivers to conduct the ground surveys of field sites during the period 16-24 April 2019. The field surveys for sample information collection covered 78 towns in the Hebei, Shanxi, Shaanxi and Gansu Provinces. Overall, 1709 field survey sites were used in this study after careful consideration and processing. We randomly selected 371 field survey sites of winter wheat, which were distributed across the study area, for use as a training set to determine constants. A further 137 winter wheat field survey sites were used to evaluate the results of winter wheat extraction from the Landsat images. Finally, 1201 field survey sites were used to evaluate the MODIS-estimated map of winter wheat, including 578 winter wheat sites and 623 non-winter wheat sites.
We referred to the geographic location and type information of the field survey sites in 2019. Considering the quality of the Google Earth TM images and the land cover type around the reference points, 441, 435 and 413 sites were obtained by visual interpretation from Google Earth TM

Other Data Sources
Agrometeorological observation records of winter wheat during the period 2000-2013 were obtained from agricultural meteorological stations governed by the National Meteorological Information Center (www.data.cma.cn/). The specific content included the crop names recorded on a 10-d timescale, the name of the development period, the date of the development period and the degree of development. For this study, we selected 44 agrometeorological sites with records of key phenological dates extending over a period exceeding 10 years, which included the seedling dates, heading dates and harvesting dates of winter wheat ( Figure 1). Agricultural census data at the county level for 2000-2014 were obtained from the Ministry of Agriculture and Rural Affairs. The elevation data were derived from a 90 × 90-m spatial resolution digital elevation model, provided by the geospatial data cloud platform of China (www.gscloud.cn/).

Winter Wheat Extraction
On the large spatial scale, the "combining variations before and after estimated heading dates" (CBAH) method is most equipped to deal with the intraclass variability of the temporal profiles of vegetation indices, and thus produces winter wheat mapping with a high accuracy [18]. We considered the range of critical stages and dates of winter wheat planting, and the characteristics of the sensitive climate of the study, area to adjust the CBAH method. According to the records of key phenological dates of winter wheat over the required 10-year (or greater) period, the ranges of three critical stages and dates for the study area were obtained by boxplot analysis.
There were considerable dissimilarities between winter wheat and others crops with regard to the changes of phenological growth curve shapes. These changes are likely due to several factors, including the following: winter wheat grows longer than other crops; the vegetation index value of winter wheat is the minimum at the seedling stage; the vegetation index value of winter wheat reaches the maximum at the heading stage; and the vegetation index of winter wheat at the harvesting stage decreases to valley value. Therefore, the growth stage of winter wheat can be divided into early growth stage and late growth stage, with the heading date as the boundary. The winter wheat was mapped by combining variations in the vegetation index before and after the estimated heading stage.
In this study, this was achieved through the following steps ( Figure 3). First, the seedling stages (SDS), heading stages (HDS), harvesting stages (HTS), and the dates of winter wheat growth, were determined based on boxplot analysis (represented by bars and dots in the upper subplot in Figure  3). Second, we calculated the minimum, maximum and minimum EVI2 values in the seedling stage, heading stage and harvesting stage, respectively. Then, the amplitude variation between the heading and seedling stages, and between the heading and harvesting stages, was identified ( Figure 3). Third, the EVI2 changes in amplitude were extracted according to the derived seedling, heading and harvesting dates for winter wheat ( Figure 3). Finally, the EVI2 variations during the early growth stage (EVE) and EVI2 variations during the late growth stage (EVL) were derived based on the EVI2 variations during the early and late growth stages. The EVE was the summation of two items: the EVI2 change amplitude in the maximum value of EVI2 at the heading stage (EVI2headingmax) and the minimum value of EVI2 at the seedling stage (EVI2seedlingmin), and the EVI2 difference between the estimated heading (EVI2heading) and seedling dates (EVI2seedling) during the early growth stage (Equation (2)). Similarly, the EVL was derived based on two other items: the EVI2 change amplitude in the maximum value of EVI2 at the heading stage (EVI2headingmax) and the minimum value of EVI2 at the harvesting stage (EVI2harvestingmin), and the EVI2 difference between the estimated heading (EVI2heading) and harvesting dates (EVI2harvesting) during the late growth stage (Equation (3)).
where EVI2seedlingmin, EVI2headingmax and EVI2harvestingmin represent the minimum value of EVI2 at the seedling stage of winter wheat, the maximum value of EVI2 at the heading stage, and the minimum value of EVI2 at the harvesting stage, respectively; and EVI2seedling, EVI2heading and EVI2harvesting represent the value of EVI2 at the seedling date (median of seedling stage) of winter wheat, the value of EVI2 at the heading date (median of heading stage), and the value of EVI2 at the harvesting date (median of harvesting stage), respectively. For a winter wheat pixel, the EVE and EVL can simultaneously gain large values, but not in nonwinter wheat pixels. This is primarily due to three reasons: First, the amplitude variation between the heading and seedling stages, and between the heading and harvesting stages, was large. Second, the EVI2 values at the heading stages would be much higher than those of the non-winter wheat crops. Finally, due to the length of the growing period of winter wheat, the increment range of the EVI2 values from three critical periods can be distinguished from non-winter wheat. In addition, the CBAH method solved the intraclass variability of winter wheat in different regions. The ranges and dates of three critical stages of winter wheat growing were determined through the monitoring information of the agricultural meteorological sites, eliminating intraclass variation caused by phenological changes. Next, the EVI2 curve shape's change characteristics, with regard to winter wheat growth, were used to greatly reduce the intraclass variability of other vegetation. In different regions, the EVE and EVL values calculated in winter wheat pixels remain the largest, although there are differences in the curve shapes of winter wheat growth phenology [18]. Therefore, it is feasible to use EVE and EVL incremental change values to identify winter wheat.
In the corresponding unknown pixels, almost no consistency exists with the curve change shown in the growth process of winter wheat. For the phenological changes of other crop pixels, the EVE and EVL cannot obtain large values at the same time. Hence, it is feasible to identify winter wheat using the EVE and EVL (Equation (4)): If (EVE > v1 and EVL >v2), it is winter wheat (4) where v1 and v2 are constants. This algorithm determines whether each pixel in the entire study area represents winter wheat or not. We

Accuracy Assessment
The accuracy assessment of the MODIS-estimated map of winter wheat was undertaken using field survey data, the interpretation of Landsat images of the 123033 tile area, and agricultural census data. First, the degree of spatial consistency of the MODIS-estimated map of winter wheat in 2019 was evaluated using the results from the 1201 field survey sites. Second, a total of 508 field survey sites and sites interpreted from Google Earth TM images were utilized to validate the results of winter wheat extraction from the Landsat images [37], of which 487 sites were classified correctly. To evaluate the MODIS-estimated map of winter wheat, the Landsat-interpreted maps of winter wheat were aggregated to 250 × 250-m pixels using the majority resampling method [18]. The Landsatinterpreted maps of the 123033 tile area in the four periods were used to evaluate the MODISestimated maps of winter wheat in 2001, 2007, 2014 and 2019. Third, the overall areas of winter wheat estimated from the MODIS images were compared with those derived from agricultural census data at the county level in 2001, 2007 and 2014. For the MODIS-derived maps, the numbers of pixels of winter wheat were calculated at the county level, which allowed direct comparison with the agricultural census data. The accuracy assessment was presented in the form of a confusion matrix or as a coefficient of determination (R 2 ) [17].

Dynamic Degree Model of Winter Wheat Planting Area Change
To quantitatively describe the range, speed and difference of winter wheat planting area change in a specific period, the dynamic degree model of winter wheat planting area change was introduced by referring to the dynamic degree model of land use [38]. The dynamic degree model of the change of winter wheat planting area represents the change of winter wheat planting area in a certain period and within a certain research range. The specific calculation formula is as follows (Equation (5)): where D d represents the dynamic degree of the spatiotemporal change of winter wheat planting area in t period. W b is the early planting area of winter wheat, W a is the late planting area of winter wheat and t is the length of the period in units of years.

Extraction of the NLWW Using Kernel Density Estimation
The aim of the kernel density estimation (KDE) is to produce a smoother visualization of discontinuous points, building upon the principle of a heat map distribution between core areas and surrounding neighborhoods [39]. KDE has been increasingly used in the fields of agriculture, geography, archaeology and epidemiology [40]. Here, we used KDE (as expressed in Equation (6)) to extract the NLWW: where ( ) is the estimator of the probability density, is a coordinate vector of estimated points, N is the number of winter wheat planting sites, ℎ (ℎ > 0) is a user-defined smoothing parameter or bandwidth, K is a user-defined nonnegative kernel function that is taken as the quadratic Epanechnikov kernel, and is a coordinate vector of sample points. On the scale of an agricultural landscape, the bandwidth is typically set to 2.5 km [39]. In this study, we used the percentile method to determine the threshold value of the kernel density estimator with the 95% percentile, and then extracted the NLWW [40].

Spatiotemporal Dynamics of the NLWW
To calculate the direction and distance of migration of the NLWW, a 5 × 5-km fishnet was generated for the study area in the longitudinal direction [41]. When the former period's limits intersected with the fishnet line, the latitudinal coordinates of the intersection were calculated and recorded. Similarly, the latitudinal coordinates where the latter period's limits intersected the fishnet line were also calculated and recorded. There were three scenarios in terms of direction. First, when the difference between the intersection points of the latitudinal coordinates in the former and latter periods was positive, the direction was coded as 1 (i.e., the NLWW moved northward). Second, the direction was coded as 0 when the intersection difference of the latitudinal coordinates was 0 (i.e., the NLWW did not move). Third, the direction was set as -1 if the difference between the latitudinal coordinates of the former and latter periods was negative (i.e., the NLWW moved southward). The distance signed with the direction of movement was computed using Equation (7): where fp-lp represents the distance northward or southward, fp is the former period, lp is the latter period, L lp represents the coordinate value of the latitudinal intersection in the latter period, and L fp represents the coordinate value of the latitudinal intersection in the former period.

Accuracy Evaluation of MODIS-Derived Winter Wheat Maps
The MODIS-derived map of winter wheat for the study area in 2019 was evaluated using 1201 ground survey sites. Of the 578 sites of winter wheat cultivation, 549 were classified correctly. The producers' accuracy (PA) and users' accuracy (UA) were 94.98% and 94.17%, respectively. Of the 623 non-winter wheat sites, 589 were identified correctly (PA: 94.54%). An overall accuracy (OA) of 94.75% was obtained when compared with in-situ observation data. The kappa index was 0.895 (Table 2).
For the MODIS-derived maps of winter wheat in 2001, 2007, 2014 and 2019, accuracy assessments were conducted using Landsat-interpreted maps of the 123033 tile area in the corresponding years [42]. Visually, the MODIS-estimated maps appeared to be consistent with the corresponding Landsat-interpreted maps (Figure 4).
Landsat-interpreted images were used to evaluate the accuracy of the MODIS-derived winter wheat maps from 2001 to 2019. The OA ranged between 92.93% and 96.06%. The PA ranged between 69.73% and 75.90%, and the UA ranged between 73.93% and 82.08%. The kappa index values ranged from 0.689 to 0.750. The OA, PA and UA average values were 94.87%, 78.30% and 72.70%, respectively, and the kappa index was 0.726 (Table 3).
The planting area of MODIS-derived winter wheat in 2001, 2007 and 2014 was aggregated to the county level. The area of winter wheat cultivation estimated from the MODIS images correlated with the agricultural census data. The coefficient of determination (R 2 ) between the MODIS-estimated dataset and the agricultural census data at the county level was in the range of 0.76-0.83, suggesting reasonable agreement between the two datasets ( Figure 5).

Spatiotemporal Dynamics of Winter Wheat Planting Area
The winter wheat planting area in northern China displayed a decreasing trend from 2001 to 2019. The planting area of winter wheat decreased significantly, from 16.0 × 10 3 km 2 to 12.4 × 10 3 km 2 during the period 2001-2014, and increased in 2019. In terms of spatial variation, the planting area in the transition zone of winter wheat planting decreased most obviously, particularly in the Shanxi and Shaanxi Provinces ( Figure 6). The Jenks optimization method was used to visually reflect the dynamic degree of spatial and temporal changes of the winter wheat planting area. Over the past 20 years, in general there has been a significant reduction in the dynamic degree of the change of winter wheat planting area (71% reduction) (Figure 7a). The trend was clearest in the northernmost planting transition zone of winter wheat (Figure 7b). The change area and proportion of winter wheat dynamic degree in a 5 × 5-km grid. (Dark orange is < −5% dynamic degree, orange is −5~0% dynamic degree, light green is 0~5% dynamic degree, green is > 5% dynamic degree.).

Latitudinal Dynamics of the NLWW during the Period 2001-2019
During the period 2001-2019, the NLWW showed a clear trend of southward movement, i.e., as an average for the entire study area, the NLWW moved 37 km southward. Specifically, it moved southward by 83 km in the area of transition between the Shaanxi Basin and Shaanxi Plateau (represented by S(a) in Figure 8a  In the past 20 years, the shift of the NLWW has exhibited a marked fluctuating trend of southward movement. During this migration, the amplitude of the shift reflected in each period has varied. The period with the most substantial southward movement was 2007-2014.
During the period 2001-2007, the average movement of the NLWW was 11 km southward. On the Loess Plateau of eastern Gansu, there was modest southward movement of 29 km (represented by S(a) in Figure 9a,b). The areas in which the NLWW moved southward more obviously were the Shanxi Central Plain (93 km; represented by S(b) in Figure 9a,b), and the border area of the Hebei Plain and Hebei mountain ranges (35 km; represented by S(c) in Figure 9a,b). However, northward movements of 44 and 30 km occurred at the Shaanxi-Gansu border and the Shaanxi-Shanxi border, respectively (represented by N(a) and N(b), respectively, in Figure 9a,b).  Figure 10a,b). The area where the NLWW moved more obviously southward was the Shanxi Central Plain (163 km; represented by S(b) in Figure 10a,b). The NLWW shifted 53 km southward on the border area of the Hebei Plain and Hebei mountain (represented by S(c) in Figure 10a,b). In addition, there was northward movement of 29 km on the Loess Plateau of eastern Gansu (represented by N(a) in Figure 10a,b). In comparison with the previous period, the NLWW generally showed only a slight trend of southward movement (4 km) during the period 2014-2019. The areas with southward movement were primarily the Shanxi Central Basin (47 km; represented by S(a) in Figure 11) and in the transition area between the Shaanxi Basin and Shaanxi Plateau (20 km). However, northward movement occurred at the northern end of the North China Plain (34 km; represented by N(a) in Figure 11).

Elevational Gradient Distributions of the NLWW during the Period 2001-2019
Over the past 20 years, a fluctuating trend of increase from east to west has been demonstrated in the elevation of the NLWW for each period (Figure 12a). In comparison with western areas, eastern parts exhibited clear changes in the elevation of the NLWW (i.e., winter wheat cultivation in highelevation areas in the east gradually decreased), particularly in the border area of the Hebei Plain and Hebei mountain ranges. The elevation limit of winter wheat planting was approximately 1600 m. Overall, the amount of winter wheat planting in low-elevation areas gradually increased, and its planting centroid shifted slowly toward low-elevation areas (Figure 12b).

Discussion
We analyzed the latitudinal dynamics of the NLWW during the period 2001-2019 in northern China. This analysis indicated that the actual NLWW showed a clear trend of southward movement. In contrast, the potential NLWW showed considerable northward movement and westward expansion [10,43]. There was indeed a gap between the potential NLWW and the actual NLWW in China. Therefore, we analyzed the reasons for this gap.

Analysis of Driving Factors for the Gap Between the Potential NLWW and the Actual NLWW
Climate change has led to a northward shifting of the potential NLWW in China [44]. In general, the limit moves to the north during warm and wet periods, and it expands to the south during the cold and dry period [45]. Over the past five decades, the record of temperature in the growing season has demonstrated a significant increase in temperature [46]. Simultaneously, extremely arid areas on land have become increasingly humid, and increasingly favorable climatic conditions have made it possible to plant winter wheat at higher latitudes [47]. This has provided favorable climatic conditions for the northward migration of the potential NLWW [48]. The northward shift of the potential NLWW shows that the climatic resource conditions suitable for winter wheat planting are gradually expanding northward, but it does not mean that the actual winter wheat planting is moving to the north.
Through the interacting effects of human activities and climate change, the actual NLWW was formed under suitable climatic conditions, and moved southward. Human activities play a decisive role in the formation and change of the actual NLWW. Urbanization development, the migration of agricultural populations and comparative incomes, among other factors, affect the shift of the actual NLWW. In the areas where groups of cities are located, farmland resources have been repurposed to provide additional space for city development [49]. For example, the Development of Western China further accelerated urbanization after the year 2000 [50]. Rural labor has undergone a large-scale transfer to urban areas, resulting in the loss of famers engaged in winter wheat cultivation, which further led to a decline, or even non-planting, of winter wheat [51]. Large areas of farmland have been abandoned in some regions of northern China [52]. In comparison with other economic crops, winter wheat often succumbed to drought, freezing and insects [9], which elevated the risks and costs associated with winter wheat planting [53]. Hence, the large input-output gap associated with winter wheat planting often reduces the willingness of farmers to plant winter wheat [54]. Many farmers in transition areas of winter wheat cultivation prefer to plant less, or even no, winter wheat. Government policy has also been a key driving force for actual NLWW changes. A series of ecological engineering projects has been applied in northern China [55]. The Grain for Green project encouraged farmland shrinkage by ecological restoration, especially in northern China [50]. This policy dictated that farmland that is unsuitable for cropping should be converted to woodland or grassland for ecological restoration. The winter wheat planting area significantly decreased during the period 2000-2010, mainly owning to the deployment of protection policies [41]. Agricultural incentive policies guided farmers to change cropping systems by subsidizing the reduction of winter wheat planting [54,56]. In addition, local governments have undertaken the supply-side structural reform of agriculture to adjust the planting structure [57]. For example, based on the natural resource conditions and production characteristics, Hebei Province's central and southern regions have been regarded as the key development areas for strong gluten winter wheat. The winter wheat planting area in the region has increased in recent years. Further, winter wheat production was also affected by expanding irrigation infrastructure, the update of winter wheat varieties, purchase prices, and so forth [58].

Uncertainties in MODIS-Derived Winter Wheat Maps
Considering the impacts of unstable climate change, weather-related agricultural disasters and the quality of the MODIS and Landsat remote sensing data, the research period of 2001-2019 was divided into four years: 2001, 2007, 2014 and 2019. The MODIS-derived map of winter wheat for each of these years can represent the stable state of the winter wheat planting area. Compared with the curve shape-based algorithm and the winter bump method, the CBAH method incorporated the intraclass variability over large regions, and thus achieved higher accuracy for winter wheat mapping [18]. For the former two algorithms, clear commission errors occurred in winter wheat mapping (i.e., overestimations were confirmed in provincial census statistics and ground survey data). The characteristics for vegetation indices were relatively stable in smaller regions, but the two algorithms had difficulties in obtaining satisfactory results in larger regions [16]. Therefore, in our larger region of 380,000 km², we used the CBAH method to attain the most appropriate results for winter wheat mapping.
It is important to note that our actual NLWW was the NLWW of the core area of winter wheat planting extracted using the KDE method on the spatial scale of the MODIS data (i.e., 250 × 250 m). Some sections of the study area contain complex terrain, fragmented land and mixed crops, which would have an impact on the accuracy of the MODIS-derived mapping. For example, a small proportion of winter wheat planting was scattered and distributed within mountainous areas, where the natural growth environment is poor.
We applied the CBAH method to map winter wheat in the study area. There were uncertainties in the verification methods adopted. First, owing to the low spatial resolution of the MODIS EVI2 data, it was not possible to fully capture the scattered areas of winter wheat planting. Second, the extraction results of the winter wheat planting area were verified using agricultural census data; however, the agricultural census data have certain inherent errors. Moreover, an assessment of the degree of consistency in the spatial distribution could not be performed. For example, at the county level, in comparison with the agricultural census data, the MODIS-estimated area of winter wheat was underestimated or overestimated in certain counties [18]. Finally, Landsat validation was only performed for an area in the middle of the Hebei Plain where winter wheat planting is reasonably concentrated. In Shanxi, Shaanxi and eastern Gansu, no validation was possible due to a lack of data or data quality problems.

Conclusions
In the context of climate change, it is crucial to understand the spatiotemporal shifts of the NLWW in order to develop adaptation strategies and ensure food security in China. In this study, we extracted the actual NLWW for 2001, 2007, 2014 and 2019 in China, and quantitatively analyzed the directions and distances of the temporal and spatial changes. The CBAH method achieved good results for winter wheat mapping in this study area based on the MODIS-derived maps of the spatial distribution of winter wheat. In the past 20 years, there has been a significant decrease of winter wheat in the transition zone of winter wheat planting. It was found that the NLWW has exhibited a marked fluctuating trend of southward migration, and the average distance of movement was 37 km in latitude. Our results revealed that the elevation limit of winter wheat planting was approximately 1,600 m, and the centroid of the area of winter wheat planting has slowly shifted toward lowelevation areas. Human activities were the decisive factor for the formation and southward movement of the actual NLWW, while the northward movement of the potential NLWW was caused by climate warming and humidification, and there was a large gap between these two limits. Our study provides an important scientific reference for adaptation of the winter wheat cropping system to future climate change.
This study represents a pioneering investigation of NLWW extraction for a large area using MODIS time series data. The actual NLWW in northern China has exhibited a fluctuating trend of southward movement as a result of anthropogenic activities under suitable climatic conditions. Further research is needed to ascertain the mechanism whereby climatic factors and human activities might influence the formation and change of the actual NLWW.