Long-term (2013–2022) mapping of winter wheat in the North China Plain using Landsat data: classification with optimal zoning strategy

ABSTRACT

Long-term mapping of winter wheat is vital for assessing food security and formulating agricultural policies. Landsat data are the only available source for long-term winter wheat mapping in the North China Plain due to the fragmented landscape in this area. Although various methods, such as index-based methods, curve similarity-based methods and machine learning-based methods, have been developed for winter wheat mapping based on remote sensing, the former two often require satellite data with high temporal resolution, which are unsuitable for Landsat data with sparse time-series. Machine learning is an effective method for crop classification using Landsat data. Yet, applying machine learning for winter wheat mapping in the North China Plain encounters two main issues: 1) the lack of adequate and accurate samples for classifier training; and 2) the difficulty of training a single classifier to accomplish the large-scale crop mapping due to the high spatial heterogeneity in this area. To address these two issues, we first designed a sample selection rule to build a large sample set based on several existing crop maps derived from recent Sentinel data, with specific consideration of the confusion error between winter wheat and winter rapeseed in the available crop maps. Then, we developed an optimal zoning method based on the quadtree region splitting algorithm with classification feature consistency criterion, which divided the study area into six subzones with uniform classification features. For each subzone, a specific random forest classifier was trained and used to generate annual winter wheat maps from 2013 to 2022 using Landsat 8 OLI data. Field sample validation confirmed the high accuracy of the produced maps, with an average overall accuracy of 91.1% and an average kappa coefficient of 0.810 across different years. The derived winter wheat area also has a good correlation (R² = 0.949) with census area at the provincial level. The results underscore the reliability of the produced annual winter wheat maps. Additional experiments demonstrate that our proposed optimal zoning method outperforms other zoning methods, including Köppen climate zoning, wheat planting zoning and non-zoning methods, in enhancing wheat mapping accuracy. It indicates that the proposed zoning is capable of generating more reasonable subzones for large-scale crop mapping.

KEYWORDS:

• Winter wheat mapping
• Landsat
• machine learning
• North China Plain
• optimal zoning

1. Introduction

Wheat is the world’s third-largest food crop, serving as a staple food for nearly 40% of the global population (FAO, 2014). China is the world’s largest producer and consumer of wheat, with an annual production of about 124 million tons and consumption of 120 million tons (Zhao et al., 2018). The North China Plain is the predominant wheat-producing region in China. Long-term mapping of winter wheat in this area is crucial as it reflects the primary spatiotemporal dynamic of winter wheat cultivation in China. Such maps are essential for identifying the driving force behind the winter wheat changes and for guiding the formulation of agricultural policy (He et al., 2019; Qiu et al., 2017; Shen et al., 2022; Zhang et al., 2019).

Remote sensing is the most effective tool for large-scale crop mapping. It has been extensively applied in mapping various crops in different regions (Li et al., 2021; Massey et al., 2017; Nabil et al., 2022; Pan et al., 2012; Song & Wang, 2019; Wang et al., 2019; Wardlow et al., 2007; Yang et al., 2023). Generally, classification methods for winter wheat mapping using remote sensing can be categorized into three types: index-based methods, curve similarity-based methods and machine learning-based methods (Dong et al., 2020).

The first type of method designs empirical indices for winter wheat based on its unique phenological features. For instance, Pan et al. (2012) proposed a crop proportion phenology index (CPPI) by capturing the unique two growth peaks at tillering and heading stages of winter wheat in Enhanced Vegetation Index (EVI) time-series. Tao et al. (2017) developed a peak before winter feature (PBWF) index to map winter wheat in the North China Plain by emphasizing the unique greenness of PBW phase. Qiu et al. (2017) proposed a method combining variations before and after estimated heading dates (CBAH) to identify winter wheat with a consideration of phenological variation across a large area. Qu and Zhang (2021) proposed a winter wheat index (WWI) based on two distinctive contrasts (i.e. the contrast between the over-wintering and sowing stages, and the contrast between the heading and harvesting stages) from Normalized Difference Vegetation Index (NDVI) time series. Mapping winter wheat based on the empirical indices is highly efficient, but heavily dependent on clear images acquired in key phenological stages. Thus, such methods are often inapplicable when the key observations are missing. The second type of method identifies winter wheat by measuring the similarity between the observed vegetation index (VI) curve and a standard VI curve of winter wheat (Belgiu & Csillik, 2018; Dong et al., 2020; Li et al., 2021; Li et al., 2021; Maus et al., 2016; Sun et al., 2012; Yang et al., 2019; Zhang et al., 2019; Zhang et al., 2019). Various similarity measurements, such as Euclidean distance, Kullback–Leibler divergence (KLD), and Jeffries-Matusita (JM) distance, have been used for winter wheat classification based on VI time-series (Sun et al., 2012). To address the issue of the shifts and distortions of VI curves induced by phenological variation across a large area, similarity based on dynamic time matching (DTW) was also applied to winter wheat mapping (Dong et al., 2020; Dong et al., 2020). The curve similarity-based method utilizes all available data in the time-series, thus has a relatively slight reliance on the effective observations of key phenological periods. However, it necessitates the dense time-series observations to match the observed and the reference curves. In summary, both the index-based and curve similarity-based methods require frequent observations to ensure the cloud-free observations in the key phenological periods or a time-series curve throughout the growing season. Thus, most of these studies utilize moderate resolution imaging spectroradiometer (MODIS) or Sentinel-2A/B data, which have a relatively high temporal resolution. However, the spatial resolution of MODIS is too coarse to support the accurate mapping of winter wheat in the North China Plain due to the fragmented landscape in this area. And the Sentinel-2A/B data, although having a higher spatial resolution, were not available before 2017. Thus, the Landsat data is the only available source that can be used for long-term mapping of winter wheat in this region. Unfortunately, the sparse time-series of Landsat poses a challenge in effectively applying the aforementioned methods.

The third type of method, machine learning-based methods, such as random forest (RF), support vector machine (SVM), and deep learning methods, have demonstrated effectiveness in crop classification using time-series data in recent years (Asgarian et al., 2016; Liang & Wang, 2019; Liu et al., 2018; Nasrallah et al., 2018; Skakun et al., 2017; Weiss et al., 2020; Xu et al., 2020). Zhang et al. (2021) utilized random forest to map winter wheat in the North China Plain based on the EVI time-series features extracted from Landsat data. Wang et al. (2019) applied random forest for winter wheat classification based on the harmonic fitting results of Landsat time-series and achieved more than 80% classification accuracy across different regions and years in the United States. These studies highlight the effectiveness of machine learning-based methods in exploring features from sparse Landsat time-series data for winter wheat mapping. Despite the promising performance of machine learning in crop classification, winter wheat mapping by machine learning encounters two major challenges: 1) obtaining sufficient labeled samples for training classifiers; and 2) zoning an entire large region into different subzones to reduce the spatial heterogeneity. Selecting high-confidence samples from existing crop maps is a practical and effective approach for the first challenge (Li et al., 2021). For instance, regions with consistent labels among different land cover maps are often considered as high-confidence and could be used for training classifier (Jin et al., 2019; Liu & Chen, 2018). Conversely, the edge pixels are often considered as unreliable and should be filtered out (Tran et al., 2022). Such a strategy is feasible for winter wheat mapping in the North China Plain because several crop maps in this area have been produced based on Sentinel-2 data (Dong et al., 2020; Dong et al., 2020; Yang et al., 2023; Zang et al., 2023). These maps could be used to generate reliable samples for training Landsat data. Although such strategies guarantee the reliability of the selected training samples, the sample representativeness becomes a concern because the difficult or confused samples may be filtered out (Li et al., 2021). Thus, it is important to carefully design the sample generation rules that is able to balance the tradeoff between the reliability and the representativeness, with consideration of the regional characteristics and the error types of the available crop maps. Regarding the second challenge, some studies divided the study area into a group of square or hexagonal grids (Dong et al., 2020; Huang et al., 2022) and trained different classifiers for different grids. However, such a strategy is somewhat arbitrary as it does not consider the varying spatial heterogeneity across different grids. Zoning with climatic or agro-ecological zones is another commonly used strategy for large-scale crop mapping (Li et al., 2021; Liu et al., 2022; Yang et al., 2021). The crop growth situation and the crop varieties within a single subzone tend to be more similar, thus a single classifier is adequate for crop classification in a single subzone. However, the classification process is influenced not only by the crop characteristics but also by the background land cover features and the data availability. Therefore, the climatic or agro-ecological zoning is not an optimal zoning strategy for crop classification either. To date, rare attention was paid to the development of an optimal zoning method specifically tailored for crop classification.

Due to the different challenges associated with different methods mentioned above, there are currently no publicly available data of long-term winter wheat maps in the North China Plain up to now. To fill this gap, this study proposes a novel machine learning-based approach for winter wheat mapping. Firstly, we generated a training sample set (comprising wheat and non-wheat samples) from several existing crop maps in the study area by considering both the sample reliability and representativeness. Secondly, we developed an optimal zoning method for crop classification by combining a measure of inter-regional classification feature consistency with the quadtree region splitting algorithm. Based on the generated sample set and the derived optimal subzones, we trained specific classifiers for Landsat data in each subzone. Then the classifiers were used to generate annual winter wheat maps from 2013 to 2022 in the North China Plain.

2. Materials

2.1. Study area

This research focuses on the North China Plain, the largest wheat-producing region in China, contributing 60% of the national wheat production (Yu et al., 2002). The region encompasses 10 provinces and municipalities, including Beijing, Tianjin, Hebei, Shandong, Henan, Anhui, Jiangsu, Hubei, Shaanxi, and Shanxi (Figure 1a). The study area exhibits significant climatic heterogeneity, including five subzones based on the Köppen climate zoning system (Peel et al., 2007) (Figure 1b): BSk (Arid Steppe Cold), Dwa (Cold Dry Winter Hot Summer), Cwa (Temperate Dry Winter Hot Summer), Dfb (Cold Without dry season Warm Summer), and Cfa (Temperate Without dry season Hot Summer). There are also large differences in the agricultural characteristics including farming systems, variety types, and production levels. It is commonly divided into five subzones (Figure 1c) according to the previous study on wheat planting zoning in China (Cai, 2010): NWS (Northwest spring wheat area), NW (Northern winter wheat area), HW (Huang-Huai winter wheat area), SWW (Southwest winter wheat area), YW (Yangtze River winter wheat area).

Figure 1. Study area (a), Köppen climate zoning (b) and wheat planting zoning (c).

2.2. Data

2.2.1. Landsat images and preprocessing

In this study, Landsat-8 surface reflectance imagery (Collection 2 Level 2 Tier 1) from 2013 to 2023 was acquired from Google Earth Engine. Given the growing season of winter wheat spans across 2 years, we selected Landsat time series from September 1 of the sowing year to June 30 of the subsequent year for winter wheat classification of each sowing year. To alleviate the cloud contamination effects, the reflectance value of the pixel with a cloud flag was replaced with an interpolated value from the two nearest valid observation. Other than the multi-spectral bands (Band 2–7, corresponding to blue, green, red, near infrared, shortwave infrared 1 and shortwave infrared 2 respectively), a suite of spectral indices were also calculated to enhance the crop spectral features for classification based on the previous studies on winter wheat mapping and winter rapeseed mapping (Table 1) (Chen et al., 2019; Dong et al., 2020; Qiu et al., 2015; Tian et al., 2021; Wang et al., 2017; Xiao et al., 2005; Yang et al., 2021; Zang et al., 2020, 2023; Zhang et al., 2021). The cloud probability band was also input as a classification feature to represent the quality of the valid observation.

Table 1. Spectral indices used for classification input.

Index		Expression
NDVI	Normalized Difference Vegetation Index	$\frac{{\mathit{\rho }}_{\mathit{NIR}}-{\rho }_{\mathit{RED}}}{{\rho }_{\mathit{NIR}}+{\rho }_{\mathit{RED}}}$
NDPI (Chen et al., 2019; Wang et al., 2017)	Normalized Difference Phenology Index	$\frac{{\mathit{\rho }}_{\mathit{NIR}}-\left(\mathit{\alpha }\times {\rho }_{\mathit{RED}}+\left(1-\mathit{\alpha }\right)\times {\rho }_{\mathit{SWIR}1}\right)}{{\rho }_{\mathit{NIR}}+\left(\mathit{\alpha }\times {\rho }_{\mathit{RED}}+\left(1-\mathit{\alpha }\right)\times {\rho }_{\mathit{SWIR}1}\right)}$ (α = 0.74)
EVI	Enhanced Vegetation Index	$2.5\times \frac{{\mathit{\rho }}_{\mathit{NIR}}-{\rho }_{\mathit{RED}}}{{\rho }_{\mathit{NIR}}+6{\rho }_{\mathit{RED}}-7.5{\rho }_{\mathit{BLUE}}+1}$
LSWI (Xiao et al., 2002, 2004)	Land Surface Water Index	$\frac{{\mathit{\rho }}_{\mathit{NIR}}-{\rho }_{\mathit{SWIR}1}}{{\rho }_{\mathit{NIR}}+{\rho }_{\mathit{SWIR}1}}$
DYI (Zang et al., 2020)	Differential Yellow Index	${\rho }_{\mathit{GREEN}}-{\rho }_{\mathit{BLUE}}$

2.2.2. Field survey data

To verify the accuracy of the mapping results, we collected survey samples from multiple sources (Table 2 and Figure 2): 1) Field surveys were conducted in the sowing years of 2017, and 2022, resulting in 1001, and 587 point samples, respectively; 2) Visual interpretation of ultra-high-resolution images acquired by Unmanned Aerial Vehicles (UAVs) during the 2018 and 2019 sowing years yielded 93 and 223 polygon samples; and 3) the open-access point samples (Qiu et al., 2022) for the sowing years of 2015, 2016, 2020, 2021 were downloaded.

Figure 2. Distribution of field survey samples.

Table 2. Information of field survey samples.

Sowing Year	Sampling Time	Sample size		Sample Type	Equivalized number of 30m pixels		Sample Sources
		Winter wheat	Non-winter wheat		Winter wheat	Non-winter wheat
2015	/	12	/	Point	12	/	Qiu et al. (2022)
2016	/	95	/	Point	95	/	Qiu et al. (2022)
2017	2018.06	633	368	Point	633	368	GPS field survey
2018	2018.12	45	48	Polygon	2457	1323	UAV Interpretation
2019	2020.03	102	121	Polygon	3394	7144	UAV Interpretation
2020	/	28	/	Point	28	/	Qiu et al. (2022)
2021	/	18	/	Point	18	/	Qiu et al. (2022)
2022	2023.03–05	343	244	Point	343	244	GPS field survey

2.2.3. Other data

We obtained existing crop mapping products in this area (Table 3), including three winter wheat maps and a rapeseed map, for automatic sample generation. In addition, the provincial census data on the annual sowing area of winter wheat from 2013 to 2021 were acquired from the China Agriculture and Forestry Database (CAFD) (https://data.stats.gov.cn/easyquery.htm?cn=C01) for comparison.

Table 3. Crop mapping products used for sample generation.

Sowing Year	Product Type	Spatial Resolution (m)	Data source	Officially reported accuracy (OA)
2017	Winter wheat map	30	Mapped by the CBAH (Qiu et al., 2017) method with the Sentinel-2 data from Sep. 2017 to Jun. 2018	87.9%*
	Winter wheat map	30	30 m winter wheat distribution map of China for four years (2016–2019) (Dong et al., 2020)	89.9%
	Winter wheat map	20	Winter wheat map in Northern China (2017–2018) (Dong et al., 2020)	89.9%
	Rapeseed map	20	China rapeseed map (2017–2021) (Zang et al., 2023)	94.9%

*This accuracy is reported in Dong et al. (2020).

3. Method

The proposed winter wheat mapping methodology consists of three parts (Figure 3): (a) a large number of wheat and non-wheat samples are automatically generated using existing crop mapping products; (b) the study area is optimally zoned into different subzones based on the classification feature heterogeneity; (c) a specific classifier is trained in each subzone with corresponding training samples and is used for mapping winter wheat from 2013 to 2022 based on Landsat-8 OLI data.

Figure 3. The workflow of the proposed winter wheat mapping approach ((a) Sample Generation; (b) Optimal Zoning; (c) Mapping and Validation).

3.1. Sample generation

We used three winter wheat maps and a rapeseed map in the sowing year of 2017 to generate training samples for Landsat data. The sample generation consists of two main steps (Figure 3a).

In the first step, areas with consistent labels among the four mapping data were selected to guarantee the reliability of the generated samples. Specifically, pixels classified as winter wheat in all three winter wheat maps and as non-rapeseed in the rapeseed map were selected for the wheat sample set. Similarly, pixels classified as non-wheat in all three winter wheat maps were selected for the non-wheat sample set. Such a selection rule effectively minimizes the inclusion of potentially erroneous samples, but it also results in an under representation of the generated sample set. A major problem is that there are few rapeseed samples that can be selected into the non-wheat sample set because winter wheat and rapeseed have similar growth seasons and phenological features during growth, leading to overlap between spectral features (Figure 4), thus they are often misclassified as the winter wheat in several winter wheat maps (Yang et al., 2023). This lack of rapeseed samples hinders the effective training of the classifier to distinguish the winter wheat from the rapeseed.

Figure 4. Landsat-8 satellite-based normalized vegetation index (NDVI) time-series variation curves for winter wheat and rapeseed.

In the second step, thus, a part of rapeseed pixels was randomly selected from the rapeseed map and added into the non-wheat sample set to increase its representativeness. The number of added rapeseed pixels ($N_{\mathrm{rapeseed}}$) was determined based on census data:

(1) ${\mathit{N}}_{\mathit{rapeseed}}=\mathit{N}\times \frac{{\mathit{S}}_{\mathit{rapeseed}}}{\mathit{S}-{\mathit{S}}_{\mathit{winterwheat}}}$(1)

where $S_{\mathrm{rapeseed}}$ and $S_{\mathrm{winterwheat}}$ are the census areas of rapeseed and winter wheat in the study area in 2017; S is the total area of the study area; and N is the number of wheat samples generated. The inclusion of rapeseed samples ensured the diversity and representativeness of the non-wheat sample set.

3.2. An optimal zoning method for crop classification

Classification feature heterogeneity means that the classifier relies on classification features with different levels of importance in different zones. The underlying reason for this variation is that the classification process is influenced by a combination of crop spectral features, background feature type and remote sensing data availability. Therefore, the core principle of zoning is to maximize the internal similarity of classification features within one zone while minimizing the similarity across different zones. However, traditional climatic zoning and agro-ecological zoning primarily consider climate or crop suitability features, which are inconsistent with the remote sensing classification features. To address this issue, we introduced a novel classification feature consistency index (CCI) to measure the similarity of classification features between two zones. Then, we developed an optimal zoning method utilizing the classical quadtree region splitting algorithm with a splitting and merging criterion of the CCI.

3.2.1. Classification feature consistency index (CCI)

The CCI is used to measure the similarity of classification features between two zones. If two zones exhibit good consistency of classification features, they should be merged into a single zone where only one classifier is trained. In contrast, if a zone exhibits internal inconsistency in the classification features, it should be split into multiple subzones where different classifiers are trained. However, due to the high dimensionality of the classification features and their complex nonlinear relationship with the classification results, the consistency of the classification features between different zones is difficult to be expressed explicitly. To avoid this difficulty, we suggest that the similarity in classification features between two zones can be represented by the consistency degree of the outputs of two classifiers trained separately with the samples in these two zones.

(2) $\mathit{O}{A}_{\mathit{A},\mathit{B}}=\frac{{\sum }_{i=1}^n{c}_i}{\mathit{n}},\mathit{where}{c}_i=\left\{\begin{array}{c}1,\mathit{if}{f}_A\left(m_i\right)=f_B\left(m_i\right)\\ 0,\mathit{if}{f}_A\left(m_i\right)\ne f_B\left(m_i\right)\end{array}\right.$(2)

Assuming that Zone A and Zone B contain training sample sets T_A and T_B respectively, we use T_A and T_B to train two classifiers f_A and f_B. Here, the training sample set of each zone consists of 2000 wheat and 2000 non-wheat samples randomly selected from the sample set generated in Section 3.1. Then, a test sample set $\mathbf{M}=\left\{m_1,\dots ,m_n\right\}$ is randomly selected from all other samples in Zone A and Zone B, and are classified using f_A and f_B, respectively. The consistency between the two classification results can be quantified by the overall accuracy metric (OA_A,B):

A higher OA_{A, B} indicates a higher consistency between two classification results and thus a higher similarity of the classification features of the two zones. However, we must also consider the presence of Bayesian errors (Error_Bayes) of the classifiers, i.e. completely random classification errors. Even classifiers trained on sample sets with identical features will have random inconsistencies in their classification results on the same test samples. It may cause the OA_A,B index to underestimate the consistency of classification features between different zones. Therefore, in this study, CCI_A,B between two zones is defined as the sum of OA_A,B and Error_Bayes:

(3) $\mathit{CC}{I}_{\mathit{A},\mathit{B}}=\mathit{O}{A}_{\mathit{A},\mathit{B}}+\mathit{Erro}{r}_{\mathit{Bayes}}$(3)

where Error_Bayes can be approximated as the training error due to the fact that the training error of a model infinitely converges to the Bayesian error when there are adequate training samples (Zhou, 2016). In this study, we computed the average training error rates of the two classifiers as the Bayesian error.

3.2.2. Quadtree region splitting based on classification feature consistency

Building upon the quadtree region splitting algorithm (Horowitz & Pavlidis, 1976), we developed a zoning process based on the CCI criterion. It consists of two main steps (Figure 3b).

In the first step, a top-down splitting process was applied to the study area using a quadtree algorithm. The entire study area is firstly divided into four subzones labeled as A, B, C, and D. The CCIs between every pair of subzones (CCI_A,B, CCI_A,C, CCI_A,D, CCI_B,C, CCI_B,D, and CCI_C,D) were calculated using the method described in Section 3.2.1. If the minimum value of these indices (CCI_min) is below a predefined threshold, it indicates a significant difference in the classification features among the four subzones. In such case, the quadtree region splitting result is retained. Conversely, if CCImin is above the threshold, it suggests no difference in classification features among subzones, and the region splitting process is halted. Subsequently, the quadtree region splitting process is repeated for the newly obtained zone through splitting until either there are no zones that can be further split (CCI_min > Threshold). In addition, the splitting will also be stopped if the available samples in the subzone are less than 2000 or the maximum number of splits (Split_Num) has been reached.

The second step involves a bottom-up merging process. After the splitting process, a queue of zones to be grown is obtained. Initially, all subzones are considered as seed points, meaning they are potential zones to be merged. From the seed queue, a zone Z₀ is selected. Its neighboring zones, denoted as $\left\{{\mathrm{Z}}_1,\dots ,{\mathrm{Z}}_{\mathrm{p}}\right\}$, are then searched. The CCI between zone Z₀ and its neighboring zones is sequentially compared. The neighboring zone with the highest CCI, denoted as Z_max, is selected.

(4) $Z_{\mathit{max}}=\underset{Z_i}{\mathit{argmax}}\left\{\mathit{CC}{I}_{Z_0,Z_1},\mathit{CC}{I}_{Z_0,Z_2},\dots .,\mathit{CC}{I}_{Z_0,Z_p}\right\}$(4)

If the $\mathrm{CC}{\mathrm{I}}_{Z_0,Z_{\max }}$ exceeds the predefined threshold, indicating no significant differences between two neighboring zones, the two zones are merged. The resulting merged zones are subsequently appended to both the “to be grown” queue and the seed queue. Conversely, if $\mathrm{CC}{\mathrm{I}}_{Z_0,Z_{\max }}$ is below the threshold, it signifies substantial classification feature difference between zone Z₀ and all its neighboring zones. Consequently, zone Z₀ is eliminated from the seed queue, preventing its merge with other zones. This iterative process continues until the seed queue becomes empty, ultimately yielding the final merged zone results.

3.2.3. Details of the zoning method implementation

Random forest was chosen as the classifier to compute the CCI, owing to its robustness and efficiency in large-scale land cover and crop classification studies (Li et al., 2021; Liu et al., 2018; Tran et al., 2022; Xing et al., 2021). Two parameters should be set during the quadtree region splitting process: the maximum number of splits (Max_split_number) and the split-merge threshold (Threshold). For our study area, Max_split_numer was set to 5, which correspond to the smallest zoning unit in the crop classification. And the CCI threshold was set to 0.97 through empirical trials.

3.3. Annually mapping of winter wheat and accuracy validation

After generating a substantial number of training samples and obtaining a set of optimal subzones, we randomly selected 4,000 wheat samples and 4,000 non-wheat training samples in each subzone to train its own classifier. Although the subzones have different spatial sizes, we use the equal number of training samples in each subzone because the heterogeneity of the classification features is comparable among different subzones generated by the proposed splitting-merging process. It is also noted that we use a larger training sample size than that used in the zoning stage (2,000 wheat samples and 2,000 non-wheat samples) because the subzones become larger and contain more available samples after the merging process. These trained classifiers of the subzones were then employed for annually winter wheat mapping from 2013 to 2022 using Landsat 8 OLI data (Figure 3c). To assess the mapping accuracy, we calculated various metrics including overall accuracy (OA), producer accuracy (PA), and user accuracy (UA) using the available validation samples for each year. Furthermore, we compare the winter wheat area derived from our produced maps with the census area obtained from the China Agriculture and Forestry Database (CAFD) at province level.

4. Results

4.1. Optimal zoning result

Based on the proposed optimal zoning method, the entire study area was divided into six subzones (Figure 5), and it is evident that the zoning differs significantly from the commonly used climatic and wheat planting zones (Figure 1b,c). The intermediate steps of the optimal zoning process are shown in Figure 6a, where the study area was first split and then merged until the proposed CCI-based criterion was reached. For each splitting and merging step, the classification accuracies (OA) were also calculated based on the validation samples (Figure 6b). In the splitting process, the OA first increased and then decreased. This is reasonable because the classification difficulty commonly decreases when a large study area is divided into smaller subzones. However, when the divided subzones are too small, the generated samples probably become less representative because the sample diversity decreases, leading to a low classification accuracy. In the merging process, the accuracy continuously increases and reaches the highest value in the final merging step. It indicates that the proposed optimal zoning method is able to acquire the best zoning map for winter wheat classification.

Figure 5. Optimal zoning results and the distribution of training samples in each subzone.

Figure 6. (a) The splitting and the merging process of the study area. (b) The change of classification accuracy (OA) during the splitting and merging process.

4.2. Winter wheat mapping results

Based on the optimal zoning result and the generated samples (Figure 5), six different classifiers for each subzone were trained by the corresponding training samples. Then, these trained classifiers were respectively used for classifying winter wheat from 2013 to 2022 in their corresponding subzones. Finally, the 30 m winter wheat maps in the North China Plain from 2013 to 2023 were generated.

The winter wheat maps generated from Landsat data (Figure 7a) show similar spatial patterns with previous studies (Dong et al., 2020; Qiu et al., 2017; Zhang et al., 2021). However, our data additionally reflect the long-term changes of winter wheat thanks to the long time series of Landsat 8 OLI since 2013. We calculated winter wheat coverage in 5-km grid for each year and conducted linear trend analysis for each grid from 2013 to 2022. As shown in Figure 7b, there is a significant spatial heterogeneity in the long-term trend of winter wheat coverage. In general, obvious decreasing trends are observed in southern boundary of winter wheat planting region, including southern Hubei and central Anhui, which might be induced by the dryland-to-paddy policy in these two provinces (Yang, 2021; Zhu, Li, et al., 2021; Zhu, Zou, et al., 2021). In comparison, areas with increasing trends are more prominent in the northern part of the study area, such as Hebei, northern Shandong, and northeastern part of Guanzhong region (a plain region across Shaanxi and Shanxi provinces). This northward shift of winter wheat planting areas might be attributed to a combination of factors, such as expanded irrigation and climate warming (Fan et al., 2020). Furthermore, scattered areas with decreasing winter wheat could also be found around the urban areas due to the urban expansion.

Figure 7. (a) Produced 30 m winter wheat map in 2022; (b) change trend of winter wheat from 2013 to 2022 at 5-km grid. The urban boundary is accessed from He et al. (2019).

Figure 8 illustrates four detailed examples of the winter wheat changes. Overall, the extracted winter wheat match well with the Landsat images acquired in April or May when the winter wheat should show a dark green color. Figure 8a (located in Gu’an County, Hebei Province) shows a decrease of winter wheat recently, due to the cropland fallow since 2020; Figure 8b (located in Xiong’an New Area, Hebei Province) shows a significant reduction of winter wheat, because of the large-scale construction in Xiong’an New Area an; Figure 8c (located in Zhoushan County, Henan Province) illustrates a rapid increase of winter wheat, driven by the increased irrigation in the region; Figure 8d (located in Jiangling County, Hubei Province) shows a decrease of winter wheat, because of the conversion of winter wheat fields into paddy rice fields.

Figure 8. Spatial details of winter wheat maps in different years.

4.3. Quality evaluation of winter wheat map

4.3.1. Accuracy assessment based on field survey samples

The classification accuracy of winter wheat was calculated for each year from 2015 to 2022, based on the field survey sample. Since the samples accessed from Qiu et al. (2022) only contained positive samples without negative samples, only the PA was calculated. As shown in Table 4, the OA for each year ranged from 85.2% to 94.3%; the Kappa coefficient ranged from 0.692 to 0.873; and the PA ranged from 75.0% to 97.0%. The accuracy validation results indicate that our method can map the winter wheat with reasonable accuracies across different years.

Table 4. The accuracy indices of the winter wheat maps from 2015 to 2023.

Sowing Year	PA	UA	OA	Kappa
2015	75.0%	/	/	/
2016	87.4%	/	/	/
2017	93.5%	93.7%	92.4%	82.6%
2018	97.0%	94.3%	94.3%	87.3%
2019	96.4%	89.1%	92.3%	84.7%
2020	89.3%	/	/	/
2021	83.3%	/	/	/
2022	89.8%	85.6%	85.2%	69.2%
Average	89.0%	90.7%	91.1%	81.0%

4.3.2. Comparison with census data

The winter wheat areas obtained from the winter wheat maps were compared with the census area at the provincial level for the years of 2013–2020. Figure 9 shows that there is a good consistency between the mapped area and the census area with an R-squared value of 0.949.

Figure 9. The comparison of the mapped area of winter wheat with the census area at the province level during 2013–2020.

4.3.3. Comparison with existing mapping products derived from Sentinel-2 data

Since Sentinel-2 data have higher spatial and temporal resolution, it is reasonable to assume that the winter wheat maps derived from Sentinel-2 data have higher reliability and thus can also serve as references for evaluating our produced winter wheat maps based on Landsat data. As shown in Table 5, our produced winter wheat maps have good agreement with the maps derived from Sentinel-2 data. The OA calculated based on these Sentinel-2 derived maps range from 80.4% to 89.3%. Considering that these Sentinel-2 derived maps also contain classification errors (officially reported OA ranged from 88% to 94%), such a level of agreement indicates that our produced map based on Landsat data is not much worse than the Sentinel-derived maps.

Table 5. Accuracies assessment based on Sentinel-2 derived maps.

Reference map	Sowing year	OA
TWDTW (Dong et al., 2020)	2016–2018	80.4%
PTDTW (Dong et al., 2020)	2017	87.9%
CBAH (Qiu et al., 2017)	2017	89.3%
ATDG (Yang et al., 2023)	2021–2022	85.0%

Figure 10 further shows several detailed comparisons of our winter wheat maps and other Sentinel-2 derived maps for the years of 2017 and 2022. In general, our map shows similar spatial pattern to the Sentinel-2 derived maps in most examples. However, the Sentinel-2 derived maps, especially ATDG, exhibit more detailed parcel boundaries due to the higher spatial resolution.

Figure 10. Details of winter wheat maps obtained using the optimal zoning method and other methods for the years 2017(a) and 2022(b).

5. Discussion

5.1. Superiority of the proposed optimal zoning method in crop mapping

Previous studies have made many improvements in various aspects of crop mapping with remote sensing data, including fusing more data sources (Blickensdörfer et al., 2022; Gu et al., 2023; Werner et al., 2024), developing more advanced classifiers or matching methods (Belgiu & Csillik, 2018; Li et al., 2024; Wang et al., 2021); generating training samples with less labor (Belgiu et al., 2021; Xuan et al., 2023; Zang et al., 2023); and designing more effective classification features (Qu et al., 2021; Xu et al., 2023; Zang et al., 2020). However, little attention was paid to improving the zoning strategy commonly used in large-scale crop mapping. The zoning strategies used in previous studies were often originally designed for other purposes (e.g. climate classification, agricultural management), thus might not be suitable for crop mapping.

This study firstly developed a specific zoning method for crop mapping by maximizing the internal similarity of classification features within one zone while minimizing the similarity across different zones. To verify the superiority of the proposed optimal zoning method over traditional zoning methods in winter wheat mapping, we compared the classification accuracy based on the proposed zoning method with those based on Köppen climate zoning, the wheat planting zoning (Figure 1), and the non-zoning. For each zoning method, we repeated the 10 times of the classification experiments (4,000 positive and negative training samples of winter wheat each randomly selected in each subzone). The mean and standard deviation of OA are reported in Table 6. Accuracies were only evaluated for the years of 2017, 2018, 2019, and 2022, because the sample sets in these years include both the positive and negative samples. The superiority of the proposed optimal zoning methods over other zoning methods varies across different years due to the different spatial distribution of the validation samples. In 2017 and 2018, the validation samples are mainly located in northern areas (Figure 2) with relatively simple crop structures, resulting in high classification accuracies for all zoning methods. In 2019, the validation sample set includes samples from rapeseed planting areas in Hubei province. Due to the similar phenology of rapeseed and winter wheat in this region, traditional zoning methods show particularly low classification accuracies for this year. In contrast, the proposed zoning method achieves much better accuracy by generating a subzone D (Figure 5) in this area to train a specific classifier that can effectively discriminate between these two crops.

Table 6. Accuracy assessment of the classification results based on different zoning methods. Bold font and underlining indicate the highest and second highest accuracy, respectively.

Sowing year	Mean OA ± Standard derivation
	Non-zoning	Climate zoning	Wheat planting zoning	Optimal zoning
2017	93.0%±0.2%	92.4%±0.2%	92.1%±0.3%	92.4%±0.5%
2018	94.6%±0.2%	94.1%±0.1%	92.8%±0.1%	94.4%±0.2%
2019	77.9%±0.3%	86.1%±1.7%	85.8%±1.0%	92.5%±1.0%
2022	84.0%±0.2%	84.1%±0.4%	85.2%±0.3%	85.4%±0.3%
Average	87.4%±0.2%	89.2%±0.6%	89.0%±0.4%	91.2%±0.5%

We also investigate the accuracy improvements of the optimal zoning method relative to non-zoning method in different subzones (Table 7). Significant improvements can be found in subzones D and F due to their diversified planting types such as winter rapeseed, paddy rice, and other crops with growing seasons overlapping with winter wheat. Such complicated cropping structure greatly reduces the classification accuracy and thus provides improvement opportunities for the optimal zoning method. For subzones A and B, the crop diversity is low (there is few other winter crop), leading to high classification accuracies. Thus, the potential benefits of optimal zoning are limited. For subzones C and E, both methods achieve 100% accuracy; however, it is meaningless due to the very limited validation samples. In summary, the optimal zoning method mainly contributes to reducing the classification difficulty in the areas with complicated cropping types, which is consistent with its designing objective.

Table 7. Accuracy improvements of the optimal zoning compared to non-zoning in different subzones.

Zone	Optimal Zoning	Non-zoning	Improved OA	Number of samples
A	94.1%	94.4%	–0. 3%	3673
B	93.7%	93.9%	–0. 2%	1108
C	100%	100%	0	24
D	73.0%	67.5%	5.5%	5170
E	100%	100%	0	2
F	87.1%	83.0%	4.1%	5929
Average	91.3%	89.8%	1.5%

We further examine the similarity and the difference in the spatial pattern of different zoning results (Figure 11). Large mountain ranges often form the important boundaries among different zones due to the significant climate difference between two sides of the mountains. For example, all of the three zoning maps have similar boundaries near Qinling mountains range, a famous boundary between north and south China. And the proposed zoning and Köppen climate zoning maps also show consistent boundaries near Taihang mountain range. But Taihang mountain range is not included in the wheat planting zoning because of no wheat planting in its west side. Huaihe river, also a famous landmark for separating north and south China due to the 0°C isotherm in January, is also considered to be a zoning boundary in both the Köppen climate zoning and the wheat planting zoning. However, there is actually no sharp climate transition across this river. Thus, it is reasonable to be excluded in our proposed zoning result. Instead, our method generates a north-south boundary between zones C and F, which does not exist in the other two traditional zoning results. It might be induced by the elevation rise from east to west in this region. In the northeast part of the study area, the proposed zoning boundaries are somehow coincided with both the Köppen climate zoning boundaries and the administrative boundary of Shandong and Hebei province. Both the climate and the agricultural policy difference could contribute to separate subzones in this area. In summary, the proposed zoning method directly considers the classification features that could vary across different natural and artificial conditions, thus can better benefit the crop mapping compared to the traditional zoning methods that mainly consider the natural factors.

Figure 11. Different spatial patterns of different zoning results: optimal zoning (a), Köppen climate zoning (b) and wheat planting zoning (c).

5.2. Sensitivity analysis of CCI threshold

The CCI threshold (Threshold) is the most important parameter in the splitting and merging process, deciding whether to split down or merge up. In order to explore the effect of threshold setting on zoning and crop mapping accuracy, we compared the zoning and classification results at different CCI thresholds, i.e. 0.99, 0.97, and 0.95 (Figure 12 and Table 8). It shows that the classification accuracy achieves the best when the CCI threshold is 0.97.

Figure 12. Optimal zoning results based on different CCI thresholds: 0.95 (a), 0.97 (b) and 0.99 (c).

Table 8. Comparison of classification result accuracy of different CCI thresholds. Bold font indicates that accuracy ranks first, and underlining indicates that accuracy ranks second.

Sowing year	Mean OA ± Standard derivation
	0.99	0.97	0.95
2017	92.2% ± 1.0%	92.4% ± 0.5%	91.8% ± 1.0%
2018	92.4% ± 0.8%	94.4% ± 0.2%	93.9% ± 0.3%
2019	90.8% ± 1.0%	92.5% ± 1.0%	86.7% ± 0.7%
2022	85.0% ± 1.0%	85.4% ± 0.3%	84.7% ± 1.0%
average	90.1% ± 1.0%	91.2% ± 0.5%	89.3% ± 0.8%

Overall, if the CCI threshold is too high (0.99), the generated subzones become relatively small, leading to inadequate diversity of generated samples for classifier training. And if the threshold is too low (0.95), the subzones become too large, resulting in high heterogeneity of classification features. However, no matter how the CCI threshold is set, 0.97, or 0.99, or 0.95, the corresponding mapping accuracy is better than that of other zoning methods.

5.3. The potential and challenges of winter wheat mapping in a longer term

In this study, we mapped the winter wheat only after 2013 because Landsat-8 OLI was launched in 2013. However, the Landsat 5 TM and Landsat 7 ETM+, which have similar spatial, spectral and temporal resolutions with Landsat 8 OLI, have been available since 1984. They thus can support for a much longer term of winter wheat mapping in the study area. Unfortunately, it remains several challenges in transferring the proposed mapping approach to the elder Landsat data. Firstly, while the USGS assumes that all Level 1 Landsat data are radiometrically calibrated and consistently geolocated, spectral inconsistency is still inevitable due to the different wavelength sensitivities (Sulla-Menashe et al., 2016; Teillet et al., 2001). Secondly, the proposed optimal zoning for winter wheat mapping was determined based on the CCI calculated in 2017. While it is reasonable to assume that the optimal zoning is stable in a ten-year term, it is likely to shift in the longer term due to the climate change, crop technique advancement and economic development. Thirdly, we trained the classifiers by the training samples in 2017 and used them to classify winter wheat in the other years. It is acceptable if the classification features keep stable across different years. However, the wheat growth condition and phenology could vary across different years, thus also alter the classification features and lead to a decrease in classification accuracy (Gadiraju & Vatsavai, 2020). Therefore, it is crucial to develop transferring technology with the consideration of the sensor inconsistence, zoning shifting and classification feature varying, for a longer-term mapping of winter wheat in the study area.

5.4. Remaining issues

There are some remaining issues to be resolved in this study.

First, although the 30 m resolution is the highest available for long-term mapping of winter wheat, mixed pixels are still unavoidable. Mixed pixel effect not only introduces about classification uncertainty but also brings about area bias in the fragmented agricultural landscapes (Dong et al., 2022; Ozdogan & Woodcock, 2006). Thus, it should be noticed that the comparison with census data provides only a reference rather than a rigorous validation (Dong et al., 2020). Spectral unmixing or sub-pixel mapping techniques offer potential solutions to the issue. However, obtaining the required training samples with coverage information presents significant challenges, thereby limiting their applications at a large scale.

Second, insufficient observation frequency is another serious challenge for crop mapping based on Landsat data. Previous studies have demonstrated that insufficient effective observations can significantly reduce the mapping accuracy of crops such as winter wheat and winter rapeseed (Dong et al., 2020; Zang et al., 2020). Spatio-temporal fusion of Landsat and MODIS data might offer a potential solution to this issue.

Third, although we particularly consider excluding winter rapeseed in the proposed approach, there are still other crops that could potentially be confused with winter wheat. For example, garlic is an important winter crop with a similar phenology to winter wheat in this study area, which is likely to introduce mapping errors in our produced winter wheat maps. However, such mapping errors should be marginal for two reasons. On the one hand, garlic accounts for a very small proportion (less than 1%) of winter crops in the study area. On the other hand, although garlic is not specifically considered in the mapping approach, the strict selection rule of the winter wheat samples (i.e. the intersection of three wheat maps) helps to exclude the misclassified garlic pixels in the sample set.

Finally, the validation samples were mainly collected along roadsides for saving cost, which may lead to biases in the representativeness of the samples. Rigorous accuracy assessment requires better sampling strategies (Olofsson et al., 2014; Stehman & Foody, 2019), but these are constrained by labor and time costs.

6. Conclusion

This study established a mapping approach for annual winter wheat in the North China Plain based on Landsat data. Firstly, a sample generation rule is specifically designed for winter wheat mapping in this area by considering the confusion between winter wheat and rapeseed in the existing crop mapping products. Secondly, an optimal zoning method for crop classification is proposed to reduce the classification feature heterogeneity over a large area. Comparative experiments showed that the proposed zoning method effectively helps to improve the mapping accuracy compared to the traditional zoning methods. It indicates that the proposed optimal zoning method can reduce the classification difficulty on the sparse Landsat time-series data and thus has potential to support long-term crop mapping.

With the proposed approach, we produced annual 30 m winter wheat map from 2013–2022 in the North China Plain, which achieves a good classification accuracy and a good agreement with census data. The produced map is opened for public to support further analysis of spatiotemporal dynamics of winter wheat in the North China Plain.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The Landsat-8 data used in this paper are available in GEE. The produced wheat maps are shared in GEE (https://code.earthengine.google.com/0bedafeeb6899c8a92e23a4e64782206?noload=1) and Zenodo (https://doi.org/10.5281/zenodo.10522027).

The validation data in this study are available from the corresponding author, X.H.Chen, upon reasonable request.

Additional information

Funding

This study is supported by the National Key Research and Development Program of China [No. 2022YFD2001100 and No. 2017YFD0300201].

Notes on contributors

Yifei Liu

Yifei Liu received the B.E. degree in remote sensing science and technology from China University of Geosciences, Wuhan, China, in 2021. She is pursuing the M.E. degree with the school of Beijing Normal University, Beijing, China. Her research interest is crop mapping with remote sensing.

Xuehong Chen

Xuehong Chen received the B.S. degree in physics and the M.S. degree in civil engineering from Beijing Normal University, Beijing, China, in 2006 and 2009, respectively, and the Ph.D. degree in earth and environmental science from Nagoya University, Nagoya, Japan, in 2012. He is currently an Associate Professor with the Institute of Remote Sensing Science and Engineering, Faculty of Geographical Science, Beijing Normal University. His research interests include data fusion, change detection and crop mapping with remote sensing.

Jin Chen

Jin Chen received the B.A. and M.S. degrees in geography from Beijing Normal University, Beijing, China, in 1989 and 1992, respectively, and the Ph.D. degree in civil engineering from Kyushu University, Fukuoka, Japan, in 2000. He was a Postdoctoral Researcher with the University of California, Berkeley, CA, USA, from 2000 to 2001, and with the National Institute of Environmental Studies, Tsukuba, Japan, from 2001 to 2004. He is currently a Professor with the Institute of Remote Sensing Science and Engineering, Faculty of Geographical Science, Beijing Normal University. His research interests include remote sensing modeling and vegetation parameter retrieval through the inversion of remote sensing model.