Mapping the Dynamics of Winter Wheat in the North China Plain from Dense Landsat Time Series (1999 to 2019)

: Monitoring spatio-temporal changes in winter wheat planting areas is of high importance for the evaluation of food security. This is particularly the case in China, having the world’s largest population and experiencing rapid urban expansion, concurrently, it puts high pressure on food demands and the availability of arable land. The relatively high spatial resolution of Landsat is required to resolve the historical mapping of smallholder wheat ﬁelds in China. However, accurate Landsat-based mapping of winter wheat planting dynamics over recent decades have not been conducted for China, or anywhere else globally. Based on all available Landsat TM/ETM+/OLI images (~28,826 tiles) using Google Earth Engine (GEE) cloud computing and a Random Forest machine-learning classiﬁer, we analyzed spatio-temporal dynamics in winter wheat planting areas during 1999–2019 in the North China Plain (NCP). We applied a median value of 30-day sliding windows to ﬁll in potential data gaps in the available Landsat images, and six EVI-based phenological features were then extracted to discriminate winter wheat from other land cover types. Reference data for training and validation were extracted from high-resolution imagery available via Google Earth™ online mapping service, Sentinel-2 and Landsat imagery. We ran a sensitivity analysis to derive the optimal training sample class ratio ( β = 1.8) accounting for the unbalanced distribution of land-cover types. We mapped winter wheat planting areas for 1999–2019 with overall accuracies ranging from 82% to 99% and the user’s/producer’s accuracies of winter wheat range between 90% and 99%. We observed an overall increase in winter wheat planting areas of 1.42 × 10 6 ha in the NCP as compared to the year 2000, with a signiﬁcant increase in the Shandong and Hebei provinces ( p < 0.05). This result contrasts the general discourse suggesting a decline in croplands (e.g., rapid urbanization) and climate change-induced unfavorable cropping conditions in the NCP. This suggests adjustments of the winter wheat planting area over time to satisfy wheat supply in relation to food security. This study highlights the application of Landsat images through GEE in documenting spatio-temporal dynamics of winter wheat planting areas for adequate management of cropping systems and assessing food security in China.


Introduction
Global population growth leads to increasing demand for agricultural crops, which is expected to continue in the coming decades [1]. Therefore, global-and regional-scale agricultural monitoring systems are needed to provide accurate information about agricultural such as setting thresholds based on vegetation indices [16,26,27]. Additionally, other methods, including supervised classification, transformation-based and metric-based approaches, have been used for crop classification [28,29]. Recently, a phenology-time weighted dynamic time warping approach has successfully been applied to map wheat [30][31][32], and deep learning has also proved to be effective to map crop types [33]. Among these approaches, the widely used Random Forest (RF) classifier was used for the NCP mapping of the past 20 years based on dense time series of Landsat medium resolution imagery. RF was chosen here due to its robustness, computational efficiency and the ability to make use of temporal information characterizing the growth cycle of winter wheat [34][35][36]. Moreover, it can efficiently handle and classify winter wheat fields with a non-normal distribution of spectral reflectance and RF is embedded in GEE that facilitates such classification endeavors across space and time based on Landsat data.
In this study, we aim to detect the dynamics of winter wheat and to fill the existing gap of long-term high spatial resolution observations of winter wheat in the NCP region. Specifically, we apply Landsat-5 TM/7 ETM+/8 OLI time series to produce EVI based temporal metrics capturing the specific growing season of winter wheat. The cloud platform of GEE provides access to the archived data from an array of Earth observation satellite systems (including the Landsat sensors). GEE allows working with images in an intrinsically parallel processing way using state-of-the-art cloud-computing and storage capability was used to process thousands of Landsat images [37]. The objectives of this study are to (1) present a phenology-based approach to map winter wheat using Landsat EVI time series; (2) implement the sensitivity analysis to account for the unbalanced distribution of land-cover types; (3) map winter wheat planting areas in NCP with 30-m Landsat time-series using the developed approach; and (4) assess the dynamics of Landsat-based winter wheat areas and detect the frequency of winter wheat planted in NCP from 1999 to 2019. Such information is a prerequisite for adequate management of cropping systems and assessing food security in China and other major wheat breadbaskets of the world.

Study Area
The study area is located between latitude 29.36 • N-42.67 • N and longitude 110.25 • E-122.70 • E and includes seven provincial-level administrative units (fine provinces and two direct-controlled municipalities) ( Figure 1a): Anhui, Jiangsu, Henan, Shandong, Hebei provinces, Beijing and Tianjin municipalities. The area includes the largest agricultural area in China, with 45,696,000 ha of croplands, of which wheat covers approximately 90%. The production of winter wheat in this area corresponds to 71% of the country's total wheat production [26]. The study area is separated by the Qin Mountains and the Huai River and thus has a topography which includes mountains in the west, hills in the south, plains in the east and middle areas, and an elevation ranging from around −100 to 3500 m (Figure 1a). The northern part of the study area is dominated by a temperate continental monsoon and the southern part is characterized by a subtropical continental monsoon.
In general, winter wheat is sown during October-November and harvested during June-July in the following year; pheno-phases includes seeding, tillering and heading during the growing season [16]. Winter wheat is mostly planted in flat areas with a relatively lower altitude that is easy to be cultivated and managed. Rapeseed is another crop type partly planted in the winter but constitutes only a small proportion of winter crops in this area. A double cropping rotation comprised of winter wheat-summer maize is the dominant cropping system, but there are also cropping rotations including winter wheat-corn, winter wheat-peanut and fallow-spring maize/peanut. The study area also has the highest population density in China with a rapid urban expansion at the expense of cropland [38,39]. area. A double cropping rotation comprised of winter wheat-summer maize is the dominant cropping system, but there are also cropping rotations including winter wheat-corn, winter wheat-peanut and fallow-spring maize/peanut. The study area also has the highest population density in China with a rapid urban expansion at the expense of cropland [38,39].

Landsat Data Processing
The study area covers 64 tiles of Landsat images (Figure 1a). Landsat images used in this study including, Landsat-5 TM which spans the long period between 1984 and late 2011. Landsat-7 ETM+ is also available from the year of 1999. Landsat-8 OLI provide imagery starting from April 2013. Landsat TM/ETM+/OLI products provide the "Surface Reflectance Tier 1" products that have the highest quality imagery meeting certain radiometric and geometric requirements. They all contain 3 visible bands (blue, green and red), 1 near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands, as well as 1-2 thermal infrared (TIR) bands. This dataset is the atmospherically corrected surface reflectance and is all available in GEE cloud platform. We used all available atmospherically corrected surface reflectance tier 1 data from Landsat-5 TM, Landsat-7 ETM+ and Landsat-8 OLI between October 1st and June 30th from 1999 to 2019. Cloud, cloud shadows, cirrus and snow/ice were masked with the CFMASK algorithm [40]. In total, 28,826 Landsat scenes were processed via the GEE platform (Figure 1b), and the spatial distribution of the number of available Landsat images for each growing season of winter wheat is shown in Figure S1 and by satellites (Landsat-5/7/8) (Figure 1c). To circumvent the impacts from potential data gaps caused by the occasional scarcity of Landsat images available, we excluded any analysis of the growing seasons spanning 2011-2012, 2012-2013, when only Landsat-5 TM and Landsat-7 ETM+ Scan Line Corrector (SLC)-off were available ( Figure  1c and Figure S1).

Landsat Data Processing
The study area covers 64 tiles of Landsat images (Figure 1a). Landsat images used in this study including, Landsat-5 TM which spans the long period between 1984 and late 2011. Landsat-7 ETM+ is also available from the year of 1999. Landsat-8 OLI provide imagery starting from April 2013. Landsat TM/ETM+/OLI products provide the "Surface Reflectance Tier 1" products that have the highest quality imagery meeting certain radiometric and geometric requirements. They all contain 3 visible bands (blue, green and red), 1 near-infrared (VNIR) bands and 2 short-wave infrared (SWIR) bands, as well as 1-2 thermal infrared (TIR) bands. This dataset is the atmospherically corrected surface reflectance and is all available in GEE cloud platform. We used all available atmospherically corrected surface reflectance tier 1 data from Landsat-5 TM, Landsat-7 ETM+ and Landsat-8 OLI between October 1st and June 30th from 1999 to 2019. Cloud, cloud shadows, cirrus and snow/ice were masked with the CFMASK algorithm [40]. In total, 28,826 Landsat scenes were processed via the GEE platform (Figure 1b), and the spatial distribution of the number of available Landsat images for each growing season of winter wheat is shown in Figure S1 and by satellites (Landsat-5/7/8) (Figure 1c). To circumvent the impacts from potential data gaps caused by the occasional scarcity of Landsat images available, we excluded any analysis of the growing seasons spanning 2011-2012, 2012-2013, when only Landsat-5 TM and Landsat-7 ETM+ Scan Line Corrector (SLC)-off were available (Figures 1c and S1).

Training/Validation Sample from Google Earth™ (GE), Sentinel-2 and Landsat
Field observations were collected in the study area from field visits during March and April 2019. This involved recording of GPS coordinates for large and homogeneous winter wheat fields ( Figure S2) to be used for analysis of corresponding Landsat image spectral signatures. We then collected training/validation sample of winter wheat fields from multiple reference image sources based on GE, Sentinel-2 and Landsat data ( Figure 2). Specifically, 18,432 points were generated (keeping the minimum distance of 200 m) from the 364 polygons that were randomly manually digitized from these reference data sources for 9 epochs  Table 1). The Sentinel-2 derived true-color composites (bands 4/3/2) and false-color composites (bands 8/4/3) from 2015 to 2018 were mainly used as a reference for sample collections. GE images were used as a reference prior to the year of 2005 and Landsat true and false-color composites were used as well for reference data collection. The Areas Of Interest (AOIs) were generated according to the field size/shape and reference GE images, Sentinel-2 and Landsat color composite images. The AOIs of the 9 epochs were equally distributed across the study area and the distribution of AOIs overlaid with Landsat tiles (path/row) fully covered all the land cover types analyzed. and reference GE images, Sentinel-2 and Landsat color composite images. The AOIs of the 9 epochs were equally distributed across the study area and the distribution of AOIs overlaid with Landsat tiles (path/row) fully covered all the land cover types analyzed.
Apart from making use of the information from the true/false-color images of Sentinel-2, GE and Landsat to identify winter wheat and the other land-cover types, we also depended on the visual interpretation of image time-series of EVI profiles to collect referring sample. Winter wheat has two distinct seasonal peaks during the cycle from October 1st to June 30th in the following year that are easily used to separate winter wheat from other land cover types, supporting collection of winter wheat sample correctly using satellite/GE imageries. Spring/summer/other crop types mainly including maize, bean and peanut, are expected to cause the primary confusion with winter wheat. We, therefore, collected additional sample of EVI in the land-cover class of croplands ( Figure S3). We also collected information on additional land cover classes that may affect the classification, including deciduous forest, evergreen forest, built-up land and water during the process of sample collections ( Table 1). As this study focuses on winter wheat, we did not evaluate the other land cover types in the confusion matrixes and merged them into a class termed non-winter wheat.  A field campaign has been carried out to visually interpret the winter wheat in the images ( Figure S2). Apart from making use of the information from the true/false-color images of Sentinel-2, GE and Landsat to identify winter wheat and the other land-cover types, we also depended on the visual interpretation of image time-series of EVI profiles to collect referring sample. Winter wheat has two distinct seasonal peaks during the cycle from October 1st to June 30th in the following year that are easily used to separate winter wheat from other land cover types, supporting collection of winter wheat sample correctly using satellite/GE imageries. Spring/summer/other crop types mainly including maize, bean and peanut, are expected to cause the primary confusion with winter wheat. We, therefore, collected additional sample of EVI in the land-cover class of croplands ( Figure S3). We also collected information on additional land cover classes that may affect the classification, including deciduous forest, evergreen forest, built-up land and water during the process of sample collections ( Table 1). As this study focuses on winter wheat, we did not evaluate the other land cover types in the confusion matrixes and merged them into a class termed non-winter wheat.

Method
The workflow of mapping winter wheat planting areas includes four steps ( Figure 3). Firstly, we extracted time series of EVI from Landsat images and used the cloud-free pixels to derive EVI based metrics that discriminate winter wheat against other land cover types. Secondly, we created training and validation sample sets using high-/medium-resolution imagery from GE online mapping service, Sentinel-2 and Landsat ( Figure 2) and implemented a sensitivity analysis to derive the best ratio of training sample (β) for an optimal balance between user's (overestimation) accuracy and producer's (underestimation) accuracy. Thirdly, we used the β-adjusted training sample to train the RF model with the input of EVI based metrics to generate winter wheat planting areas for 1999-2019. Finally, we used sample from step 2 for validation and analyzed spatio-temporal changes in winter wheat planting areas.

Metrics Derived for Classification
All the non-masked Landsat images were used to calculate EVI and gaps were filled in with the median value of a 30-day sliding window. At the same time, the median window was used to smooth the time series of EVI over each winter wheat growing season (1 October-30 June) ( Figure 4). The original time series of EVI are shown in Figure S4.
We averaged the smoothed EVI for the Phase I (20 October-31 December; covering the period of seedling and tillering), Phase II (10 January-25 February; over-wintering),

Metrics Derived for Classification
All the non-masked Landsat images were used to calculate EVI and gaps were filled in with the median value of a 30-day sliding window. At the same time, the median Temporally smoothed EVI of winter wheat and other land cover types for northern to southern zones (the corresponding original EVI spectral are shown in Figure S4). Light and darker grey shadings represent the respective periods of seeding and heading of winter wheat. (c) EVI of different land cover classes represented by the six phenological metrics.

Sensitivity Analysis of Sampled Training Data
The class distribution of the collected sample for winter wheat and non-winter wheat was unbalanced, which could potentially introduce a bias of classification with overestimation of the dominant (winter wheat) and underestimation of the rare class (non-winter wheat). To evaluate the impact of unbalanced class distribution on the accuracy, we applied a sensitivity analysis approach proposed by [35] that tests different class ratios (β) to define the best ratio of training sample to find an optimal balance between user's (overestimation) accuracy and producer's (underestimation) accuracy. For the analysis of class ratio optimization and the final classification and accuracy assessment, a 10-fold iterative procedure was performed using the selected metrics for the RF. Initially, 60% (11,059) of the collected pixels were randomly selected as the training pixels (Train 1) and the remaining 40% (7373) were used as validation pixels (Test 1). For each iteration, a fixed number of 1500 training points (Train 2) were randomly generated from the pool of training pixels (Train 1) and the rest of the training pixels (Train 1) were used as a test set (Test 2). T he parameter β refers to the ratio of winter wheat and non-winter wheat pixels in the training set and was changed iteratively to approach an optimal value βn balancing producer's and user's accuracies. The procedure started from an equal class distribution (βi = 1), which defines the number of pixels for winter wheat (1500/(1 + βi)) and non-winter wheat (1500 * βi/(1 + βi)), and in each step βi was increased by 0.1 until reaching a non-winter wheat-fivefold class distribution (βi = 5) [36]. The optimal ratio of βi was subsequently used to adjust the subsets of Train 1 and the remaining sample pixels were used as test pixels. Temporally smoothed EVI of winter wheat and other land cover types for northern to southern zones (the corresponding original EVI spectral are shown in Figure S4). Light and darker grey shadings represent the respective periods of seeding and heading of winter wheat. (c) EVI of different land cover classes represented by the six phenological metrics.

Classification and Accuracy Assessment
We averaged the smoothed EVI for the Phase I (20 October-31 December; covering the period of seedling and tillering), Phase II (10 January-25 February; over-wintering), Phase III (20 March-15 May or 10 April-31 May; jointing/heading) and Phase IV (20 May-30 June or 1 June-30 June; maturing) to derive four EVI based metrics that are used to classify the winter wheat against other land cover types (Figure 4a,b). Two different dates for Phase III and Phase IV were applied for the northern and southern study area, defined by the transection of the latitude 34.57N degree and determined by the change in winter wheat phenology (Figure 4a) from northern to southern zones. The different dates are implemented using the grids (80 × 80 km) for reference samplings shown in Figure 2a. We further calculated the Feature 1 (F1) and Feature 2 (F2) as the difference in EVI average values between phase I-II and III-IV, respectively. These six image composite metrics represent the inputs for RF model training. A distinct difference between winter wheat and the other land cover is observed from several of the EVI metrics derived ( Figure 4) and particularly so for the F1. The use of average EVI over a given phase has the advantage of reducing the influence from gaps in Landsat image availability.

Sensitivity Analysis of Sampled Training Data
The class distribution of the collected sample for winter wheat and non-winter wheat was unbalanced, which could potentially introduce a bias of classification with overestimation of the dominant (winter wheat) and underestimation of the rare class (non-winter wheat). To evaluate the impact of unbalanced class distribution on the accuracy, we applied a sensitivity analysis approach proposed by [35] that tests different class ratios (β) to define the best ratio of training sample to find an optimal balance between user's (overestimation) accuracy and producer's (underestimation) accuracy. For the analysis of class ratio optimization and the final classification and accuracy assessment, a 10-fold iterative procedure was performed using the selected metrics for the RF. Initially, 60% (11,059) of the collected pixels were randomly selected as the training pixels (Train 1) and the remaining 40% (7373) were used as validation pixels (Test 1). For each iteration, a fixed number of 1500 training points (Train 2) were randomly generated from the pool of training pixels (Train 1) and the rest of the training pixels (Train 1) were used as a test set (Test 2). The parameter β refers to the ratio of winter wheat and non-winter wheat pixels in the training set and was changed iteratively to approach an optimal value βn balancing producer's and user's accuracies. The procedure started from an equal class distribution (βi = 1), which defines the number of pixels for winter wheat (1500/(1 + βi)) and non-winter wheat (1500 * βi/(1 + βi)), and in each step βi was increased by 0.1 until reaching a non-winter wheat-fivefold class distribution (βi = 5) [36]. The optimal ratio of βi was subsequently used to adjust the subsets of Train 1 and the remaining sample pixels were used as test pixels.

Classification and Accuracy Assessment
We applied a RF model to map winter wheat planting areas with Landsat EVI time series data consisting of six EVI based metrics (October-June). The RF classifier is an ensemble classifier that generates multiple decision trees and is trained using bagging, thereby enabling the trees to determine the probability of the class membership [31]. RF performs multiple criteria classifications and has the advantage of fast processing speed and is less sensitive to overfitting. We used the default GEE parameters (the number of variables per split is set to the square root of the number of input features and the fraction of input to bag per tree is 0.5) for RF with 100 trees to train the time series of Landsat EVI with the optimal β-adjusted Train 1 sample pixels (a total of 11,059 (60%)) ( Figure 5). The remaining sample pixels (a total of 7373 (40%)) were subsequently used to validate the performance of the developed model for each epoch separately. We developed one model for the entire period of analysis with the 60% of the entire collected samples, as separate RF models for each epoch cannot be trained since the collected training data only include one land cover class (winter wheat) during the early period (e.g., the epoch for 1999-2000).
Remote Sens. 2021, 13, x FOR PEER R EV IEW 9 of 18 We applied a RF model to map winter wheat planting areas with Landsat EVI time series data consisting of six EVI based metrics (October-June). The RF classifier is an ensemble classifier that generates multiple decision trees and is trained using bagging, thereby enabling the trees to determine the probability of the class membership [31]. RF performs multiple criteria classifications and has the advantage of fast processing speed and is less sensitive to overfitting. We used the default GEE parameters (the number of variables per split is set to the square root of the number of input features and the fraction of input to bag per tree is 0.5) for RF with 100 trees to train the time series of Landsat EVI with the optimal β-adjusted Train 1 sample pixels (a total of 11,059 (60%)) ( Figure 5). The remaining sample pixels (a total of 7373 (40%)) were subsequently used to validate the performance of the developed model for each epoch separately. We developed one model for the entire period of analysis with the 60% of the entire collected samples, as separate RF models for each epoch cannot be trained since the collected training data only include one land cover class (winter wheat) during the early period (e.g., the epoch for 199 . The classification accuracy was estimated using contingency (confusion) matrices. We calculated the area-weighted overall accuracy, producer's and user's accuracies for the two categories of winter wheat and non-winter wheat for each year. We also calculated error-adjusted area estimates and the proportion of each class based on the results of accuracy assessments following state-of-the-art accuracy assessment recommendations [41]. Figure 5. Estimates of the optimal class ratio (β) that achieves a balance between user's accuracy (black line) and producer's accuracy (red line). The lines indicating the mean of the accuracies were produced from 10-fold subsampling runs for each β. The shadings show 95% confidence intervals. An optimal ratio of β = 1.8 was derived and applied for all years of classification.

Analysis of Spatio-Temporal Patterns of Winter Wheat Planting Areas
The years of 2011-2012 are excluded from further analyses because of limited Landsat data availability (Figure 1c and Figure S1), whereas for the years of 1999,2002,[2005][2006][2007]2009 and 2015, an uncertainty measure is estimated from the unobserved component of the maximum extent of wheat during the full period following the approach of [42]. As part of the wheat planting area in those years is potentially not taken into account (not observed), this represents a source of underestimation of the reported area. We combined the measured values of wheat planting area with an estimate of the area of unobserved Figure 5. Estimates of the optimal class ratio (β) that achieves a balance between user's accuracy (black line) and producer's accuracy (red line). The lines indicating the mean of the accuracies were produced from 10-fold subsampling runs for each β. The shadings show 95% confidence intervals. An optimal ratio of β = 1.8 was derived and applied for all years of classification.
The classification accuracy was estimated using contingency (confusion) matrices. We calculated the area-weighted overall accuracy, producer's and user's accuracies for the two categories of winter wheat and non-winter wheat for each year. We also calculated error-adjusted area estimates and the proportion of each class based on the results of accuracy assessments following state-of-the-art accuracy assessment recommendations [41].

Analysis of Spatio-Temporal Patterns of Winter Wheat Planting Areas
The years of 2011-2012 are excluded from further analyses because of limited Landsat data availability (Figures 1c and S1), whereas for the years of 1999,2002,[2005][2006][2007]2009 and 2015, an uncertainty measure is estimated from the unobserved component of the maximum extent of wheat during the full period following the approach of [42]. As part of the wheat planting area in those years is potentially not taken into account (not observed), this represents a source of underestimation of the reported area. We combined the measured values of wheat planting area with an estimate of the area of unobserved but potential wheat planted area to account for this by computing the maximum wheat area extend during the entire period. The true wheat planting area will lie somewhere within this unobserved range, but the actual limits cannot be established. To characterize the spatio-temporal patterns of winter wheat planting areas in the NCP, we examined their dynamics for 1999-2019 using a linear trend analysis for the entire area and (years with insufficient Landsat data coverage were excluded from the trend analysis) at the province level (Beijing and Tianjin are excluded from this analysis as they occupy a small proportion of the winter wheat planted area across the NCP (~0.62%). We then calculated the per-pixel frequency of Landsat derived winter wheat planting, which is the total occurrences of winter wheat planting from 1999 to 2019 divided by the number of years. Moreover, we assessed the changes in winter wheat planting areas by calculating the classes of stable, gain and loss of winter wheat areas for the period 2000-2019. The classes were defined accordingly: We split the entire period into two sub-periods, covering 1999-2009 and 2010-2019, respectively, and calculated the per-pixel frequency of winter wheat fields within each sub-period. Frequencies of winter wheat fields of the two sub-periods were compared and if the increase rate was above 20%, the pixel was denoted as gain and if the decrease rate was lower 20%, it was denoted as loss. Pixels of change frequencies in-between 20% decrease and 20% increase were denoted as stable.

Winter Wheat Area Mapping Using Landsat Time Series
The winter wheat planting areas were mapped for the growing season of 2017-2018 ( Figure 6), showing that winter wheat was mostly occurring in Henan, southern Hebei, northern Anhui, northern Jiangsu and Shandong (except the Taishan mountain areas), with an entire area coverage of 3.32 (±0.25) × 10 7 ha in NCP. Specifically, the Henan province has the largest area of winter wheat with 8.47 × 10 6 ha (25.5% of the winter wheat planting areas in NCP), followed by Shandong (22.11%), Jiangsu (20.79%), Anhui (19.93%), Hebei (11.09%), Tianjin (0.41%) and Beijing (0.15%). Figure 6b and c show four different subsets of the classification results representing the varying density of winter wheat and the corresponding Landsat-8 false-color composites (bands-5/4/3).

Accuracy Assessment of Landsat-Based Winter Wheat Area Maps
The overall accuracies of the Landsat-based classification of winter wheat planting areas ranged from 0.824 to 0.999 (Tables 2 and S1). The user's accuracies varied between 0.900 and 0.988 and the producer's accuracies between 0.940 and 0.999, indicating that the model can accurately map winter wheat planting areas throughout all epochs/years. It can be seen that the relatively low overall accuracy of 0.824 for the growing season of 2016-2017 was mainly attributed to the errors of commission for deciduous forest (38.2%), built-up land (19.5%) and other crops (18.8%). More detailed confusion matrices for each land cover type are shown in Table S2. northern Anhui, northern Jiangsu and Shandong (except the Taishan mountain areas), with an entire area coverage of 3.32 (±0.25)  10 7 ha in NCP. Specifically, the Henan province has the largest area of winter wheat with 8.47  10 6 ha (25.5% of the winter wheat planting areas in NCP), followed by Shandong (22.11%), Jiangsu (20.79%), Anhui (19.93%), Hebei (11.09%), Tianjin (0.41%) and Beijing (0.15%). Figure 6b and c show four different subsets of the classification results representing the varying density of winter wheat and the corresponding Landsat-8 false-color composites (bands-5/4/3).

Accuracy Assessment of Landsat-Based Winter Wheat Area Maps
The overall accuracies of the Landsat-based classification of winter wheat planting areas ranged from 0.824 to 0.999 (Table 2 and Table S1). The user's accuracies varied between 0.900 and 0.988 and the producer's accuracies between 0.940 and 0.999, indicating that the model can accurately map winter wheat planting areas throughout all    (Figure 7a). These periods, therefore, are not included in the calculation of the trend analysis. Further, we detect the frequency of winter wheat planting for the observed periods, showing a distinct higher frequency of winter wheat planting in eastern Henan, northern Anhui, southern Hebei, western Shandong and the most area in the northern of the Yangtse River in Jiangsu province (Figure 7c). Moreover, 1.02 × 10 7 ha (occupying~18% of the area where winter wheat is encountered during 1999-2019, i.e., frequency of winter wheat planted areas >0) of the study areas are observed always as winter wheat ( Figure S5). Finally, we analyzed changes in winter wheat planting areas (see Section 2.2.4.) during 1999-2019 ( Figure 7d) and found that 44% of the area is characterized by the class denoted stable, whereas 29% of the area showed a gain and 27% a loss in winter wheat planting areas.

Landsat Images for Phenology-Based Winter Wheat Planting Area Mapping
We developed an approach using freely available time series data of Landsat imagery to reconstruct the dynamics of winter wheat planting areas in China at a spatial resolution of 30 m dating back to 1999. The results thereby advance the monitoring of winter wheat from MODIS 250-1000 m to 30 m spatial resolution [17,27]. As wheat growing in China is characterized by a high cropland fragmentation from the small and patchy nature of fields, the use of coarser-resolution MODIS data (250 m) includes considerable uncertainties. Satellite data of even finer spatial resolution with the ability to resolve temporal dynamics in seasonal spectral signatures, such as Sentinel-2 MSI (10 m), is promising to further improve mapping of winter wheat over such fragmented croplands, which however is limited to the mapping of recent years (since the launch of the Sentinel-2 A in June 2015). Four features that are important for the high-quality results are: (I) The spatial resolution (30 m) benefits the mapping of winter wheat in areas characterized by smallholder agriculture and heterogeneous landscape structures. (II) GEE, as a cloud computing platform, enables the means for fast data preprocessing and computation of massive amounts of Landsat scenes (28,826 Landsat images used in this study). (III) The use of phenological features derived from Landsat time series during growing seasons (as compared to the use of a single scene) was helpful to distinguish winter wheat from other crop and land cover types. (IV) The sensitivity detection analysis to derive the optimal class ratio (β = 1.8) of the training sample pixels, which achieves a balance between user's accuracy and producer's accuracy, was used to reduce the bias that is brought about by the unbalance of the class distribution.

Landsat Images for Phenology-Based Winter Wheat Planting Area Mapping
We developed an approach using freely available time series data of Landsat imagery to reconstruct the dynamics of winter wheat planting areas in China at a spatial resolution of 30 m dating back to 1999. The results thereby advance the monitoring of winter wheat from MODIS 250-1000 m to 30 m spatial resolution [17,27]. As wheat growing in China is characterized by a high cropland fragmentation from the small and patchy nature of fields, the use of coarser-resolution MODIS data (250 m) includes considerable uncertainties. Satellite data of even finer spatial resolution with the ability to resolve temporal dynamics in seasonal spectral signatures, such as Sentinel-2 MSI (10 m), is promising to further improve mapping of winter wheat over such fragmented croplands, which however is limited to the mapping of recent years (since the launch of the Sentinel-2 A in June 2015). Four features that are important for the high-quality results are: (I) The spatial resolution (30 m) benefits the mapping of winter wheat in areas characterized by smallholder agriculture and heterogeneous landscape structures. (II) GEE, as a cloud computing platform, enables the means for fast data preprocessing and computation of massive amounts of Landsat scenes (28,826 Landsat images used in this study). (III) The use of phenological features derived from Landsat time series during growing seasons (as compared to the use of a single scene) was helpful to distinguish winter wheat from other crop and land cover types. (IV) The sensitivity detection analysis to derive the optimal class ratio (β = 1.8) of the training sample pixels, which achieves a balance between user's accuracy and producer's accuracy, was used to reduce the bias that is brought about by the unbalance of the class distribution.

Uncertainty in Winter Wheat Mapping
While the accuracies of the image classifications were generally high, potential uncertainties can be associated with the applied data and algorithms. Although the combined use of Landsat TM/ETM+/OLI images delivers a reasonably dense time series, the availability of Landsat images during the period of 2011-2012 is low (Figures 1c and S1), resulting in larger uncertainties of the results of the area classified as wheat during this period as the low frequency of observations cause a failure when computing all the phenological metrics used for the RF classifier. Therefore, 2011-2012 were omitted from the classification results reported, while for other years (1999,2002,(2005)(2006)(2007)2009 and 2015) uncertainty is estimated from the unobserved component of the maximum extent of wheat planting (Figure 7a,b). Reducing uncertainties in relation to production of historical reference data for the image classification is common challenge, but the use of historical VHR data in GE provides a valuable complementary data source to the use of Sentinel-2 and Landsat imagery. Uncertainties in the training/validation workflow due to possible Landsat coregistration errors is expected to be minimal due to the use of the Level-1TP (Precision and Terrain Correction) product in the wheat growing areas that are mostly flat areas. Furthermore, other uncertainties may exist in the design of the study: The phenology of winter wheat may have shifted due to global warming over the long time series. Applying different dates for Phase III and Phase IV calculations north and south of 34.57N latitude (due to differences in phenology) inevitably introduces the risk of creating artefacts in the classification output. Yet, no border effect could be observed. Finally, intra-class variability caused by field conditions and agricultural practices was not considered here [26].
We compared the winter wheat map (3.32 × 10 7 ha) developed in this study with the recently published winter wheat maps in the season of 2017-2018, with 3.73 × 10 7 [31] and 2.07 × 10 7 ha [32] (Figure 8). While the overall area estimate of our classification matches the area estimate of in Figure 8a relatively well, the spatially explicit mapping comparison show a 68% and 53% match of the winter wheat areas mapped in the two studies as compared to the winter wheat map developed here. Whereas our study generally maps less areas as winter wheat as compared to ref [31] (orange areas in Figure 8a), it is clear that our study maps considerably larger areas of winter wheat than ref [32] (pink areas in Figure 8b). This suggests that uncertainties of winter wheat maps are impacted by the differences in remote sensing data source and methods used.
While the accuracies of the image classifications were generally high, potential uncertainties can be associated with the applied data and algorithms. Although the combined use of Landsat TM/ETM+/OLI images delivers a reasonably dense time series, the availability of Landsat images during the period of 2011-2012 is low (Figure 1c and Figure  S1), resulting in larger uncertainties of the results of the area classified as wheat during this period as the low frequency of observations cause a failure when computing all the phenological metrics used for the RF classifier. Therefore, 2011-2012 were omitted from the classification results reported, while for other years (1999, 2002, 200 5-2007, 2009 and 2015) uncertainty is estimated from the unobserved component of the maximum extent of wheat planting (Figure 7a,b). Reducing uncertainties in relation to production of historical reference data for the image classification is common challenge, but the use of historical VHR data in GE provides a valuable complementary data source to the use of Sentinel-2 and Landsat imagery. Uncertainties in the training/validation workflow due to possible Landsat co-registration errors is expected to be minimal due to the use of the Level-1TP (Precision and Terrain Correction) product in the wheat growing areas that are mostly flat areas. Furthermore, other uncertainties may exist in the design of the study: The phenology of winter wheat may have shifted due to global warming over the long time series. Applying different dates for Phase III and Phase IV calculations north and south of 34.57N latitude (due to differences in phenology) inevitably introduces the risk of creating artefacts in the classification output. Yet, no border effect could be observed. Finally, intraclass variability caused by field conditions and agricultural practices was not considered here [26].
We compared the winter wheat map (3.32 × 10 7 ha) developed in this study with the recently published winter wheat maps in the season of 2017-2018, with 3.73 × 10 7 [31] and 2.07 × 10 7 ha [32] (Figure 8). While the overall area estimate of our classification matches the area estimate of in Figure 8a relatively well, the spatially explicit mapping comparison show a 68% and 53% match of the winter wheat areas mapped in the two studies as compared to the winter wheat map developed here. Whereas our study generally maps less areas as winter wheat as compared to ref [31] (orange areas in Figure 8a), it is clear that our study maps considerably larger areas of winter wheat than ref [32] (pink areas in Figure 8b). This suggests that uncertainties of winter wheat maps are impacted by the differences in remote sensing data source and methods used.  [32]. A+B+ denotes both studies agree on the winter wheat mapping, A+B-denotes winter wheat mapped in the cited study but not in our study, A-B+ denotes winter wheat mapped  [32]. A+B+ denotes both studies agree on the winter wheat mapping, A+B− denotes winter wheat mapped in the cited study but not in our study, A−B+ denotes winter wheat mapped in our study but not in the cited study, A−B− denotes that winter wheat are not mapped in any of the studies. The corresponding winter wheat maps are shown in Supplementary ( Figure S6).

Dynamics of Winter Wheat
An increase of winter wheat, on average, 1.42 × 10 6 ha in the entire NCP region, was observed in comparison to the year 2000 and significant increases were observed in the Hebei and Shandong provinces. Surprisingly, our findings contrast with other studies, which indicated that rapid economic development and urbanization in China resulted in the reduction of arable land in China available for agriculture/winter wheat production [38]. China experienced dramatic economic growth during recent decades, with a strong impact on livelihood strategies with many people moving from rural regions to cities, which has caused widespread urbanization and urban expansion [43]. This phenomenon has been particularly pronounced in NCP that includes the Jingjinji and Yangtze river delta urban agglomeration having the highest density of human population and the fastest rate of urban expansion in China [44]. Therefore, on the one hand, urbanization results in the loss of cultivated cropland due to urban sprawl [38]. On the other hand, urbanization may also have influenced cropping systems by increasing/decreasing winter wheat planting areas at the expense of other crop types. As urbanization can lead to a shift in cropping systems, for instance, increasingly more farmers plant market-driven crop species to optimize the household income [45]. Hence, in such circumstances, the winter wheat planting areas can be sustained and show an increasing tendency (3.32(±0.25) × 10 7 ha) characterized even by a significant increase in some provinces (Hebei and Shandong). This study thereby provides a quantification supporting the important role of winter wheat for sustaining food security in China and that the high-quality arable land of NCP is being protected for agricultural use.
Climate change leads to an increase in the frequency of extreme weather events [46], which would pose a severe challenge for cropping systems in NCP. As winter wheat has a long growing season that spans autumn, winter and spring, it is sensitive to a changing climate, in particular increased temperature [47]. This may potentially lead to a change in winter wheat planting due to an adaption to ongoing climate change. In fact, heat stress can severely affect the wheat's growing season and yield [48], so farmers would be forced to change the future cropping types to plant crop types having a higher resilience to warming in order to mitigate economic losses. The availability of water (rainfall pattern, river runoff, irrigation) could be another factor controlling changes in winter wheat planting and the cropping systems [49]. Although suffering from these different adverse impacts of climate change and water availability, this study suggests that winter wheat planting is steadily increasing in the NCP, mitigating the increased demand from a rapidly growing population and reducing the dependency on grain import.

Implication and Future Development of Winter Wheat Mapping
Our study demonstrated the possibility of mapping winter wheat planting areas and tracking its dynamics using Landsat data in China, showing the potential to be applied at a global scale. The study contributes to remote sensing-based monitoring of wheat production, forming an important improved basis for global food security programs, and helps to explain the agriculture dominated greening Earth [50] and may also contribute to the understanding of the decline in terrestrial water storage in NCP [51]. Moreover, our study can be used as inspiration for developing workflows for winter wheat mapping in other major wheat breadbaskets of the world, such as smallholder farming systems in India, and the more mechanized large-scale production systems of European countries and Russia. While a finer-resolution of, e.g., Sentinel-2-based winter wheat mapping is a logical next step, the extensively temporal span of Landsat based mapping dating back to the 1980s is unmatched. This work was faced by the challenge of the uneven temporal Landsat data coverage. Yet, from the approach of forming epochs, high accuracy of the wheat mapping was achieved throughout the period of analysis. The limited Landsat data availability prior to the turn of the new millennium did not allow for a robust assessment following the epoch methodology taken here. The increase in winter wheat planting areas during recent decades, seemingly associated with a decline in areas under agriculture [38], suggests changed cropping systems, which should be further studied.

Conclusions
An approach for mapping of winter wheat planting area at 30 m spatial resolution was developed for the breadbasket of China from 1999 to 2019, using Landsat imagery time series. We observed an increase in winter wheat of 1.42 × 10 6 ha in the entire area of the NCP in comparison to the year 2000, with a significant increase in the Shandong and Hebei provinces (p < 0.05). In the context of an overall decrease in cropland areas and adverse impacts from climate change and water availability, these results suggest that cropping systems in NCP sustain winter wheat planting to meet the increased food demand from an increasing population. The results underline the value and potential for Landsat images for long-term monitoring of winter wheat in areas characterized by small agricultural plot sizes and efficiently quantify the past and ongoing changes in wheat cropping systems in one of the most densely populated regions of China. The approach developed has the potential to be reproducible elsewhere and/or at the global scale, as winter wheat generally has distinct temporal phenological spectral features as compared to other crop-types, and the entire Landsat archive is accessible via cloud-based computing services such as the one provided by GEE.  Figure S4. Temporal EVI of (a) winter wheat and (b) other land cover types in the study area; Figure S5. Map of stable and changed winter wheat planting areas for the period 1999-2019; Figure S6. Maps of winter wheat planting areas for the period 2017-2018 based on the reference studies [1,2]; Table S1. Accuracy error matrix. Table S2. Confusion matrix of land cover validation based on the remaining sample pixels (Test 1) from GE, Sentinel-2 and Landsat.