Multi-Source Data Modeling of the Spatial Distribution of Winter Wheat Yield in China from 2000 to 2015

Yield gridded datasets are essential for agricultural land management, food security and harmonious human–land relationships. Many studies have developed yield spatialization models that are based on cropland areas. However, crop planting areas, phenological dates, and net primary production (NPP) have received minimal attention. This study proposes a novel method to simulate winter wheat yields in China from 2000 to 2015 using crop phenological datasets, phenological observations, and NPP. The results showed that the NPP in the growing season and statistical yield showed a significant positive correlation (R2 = 0.93, p < 0.01). The mean prediction error of the gridded yield dataset was 12.01%. The relative errors of the gridded yield dataset for approximately half of the samples were between −10% and 10%. Furthermore, the yield distribution was high in the east and low in the west. The high yield was primarily concentrated in the North China Plain, while low yield was observed in eastern Gansu, central Shanxi, southern Hebei, and eastern Sichuan. From 2000 to 2015, the yield mainly showed an increasing trend in the study area, with the average rate of 0.17 t ha−1 yr−1, especially in the North China Plain. This study suggests that NPP is a key indicator to evaluate the yield of winter wheat. Furthermore, this method can be used to generate gridded yield maps along with providing credible and fundamental data for climate change and sustainable agricultural development.


Introduction
The yield spatial distribution refers to the distribution of the yield in a certain region within a certain time. According to differences in data sources, the methods of estimating the spatial distribution of crop yield can be classified into two categories. The first category can estimate the spatial distribution of yields that are based on yield spatialization [1,2]. Yield statistical data, which are compiled according to political or administrative units, are the reliable data source of spatialization model. The second category usually estimates yield based on remote sensing images. Generally, the relationship between remote sensing products and harvest yield can be established for estimating yield [3,4]. These remote sensing products, including net primary production (NPP), normalized difference vegetation index (NDVI), and enhanced vegetation index (EVI), have been used for estimating yield [3,[5][6][7]. They have the characteristics of real time, rapidity, and wide cover range, which provide the fundamental data for yield estimation [4,8]. Furthermore, statistical and remote sensing data are instrumental in yield gridded data estimation, and can be auxiliary to one another. However, yield spatialization using

Study Area
Wheat, one of the major food crops in China, is planted in most regions of the country [25]. The country could be divided into three production areas of wheat: the spring wheat region, the northern winter wheat region, and the southern winter wheat region. The cultivation areas of wheat are mainly distributed in northern China, the northern region of Northeast China and the north-eastern area of Southwest China [26]. According to the Chinese National Bureau of Statistics, the planting area of wheat in China in 2018 was 24.51 million hectares, of which winter wheat accounted for 93% [27]. This study was performed in the main planting area of winter wheat in China (Figure 1). This region has an area of 2.86 million km 2 and is located at 92 • 44 ~122 • 41 E and 21 • 07 ~42 • 48 N. It is mainly located in temperate and subtropical monsoon climate zones of China, characterized by a synchronization of high temperature and ample precipitation. The average annual temperature and precipitation are approximately 10.4 • C and 800 mm, respectively.

Data Sources
A 1 km-grid crop phenological dataset, phenological observations from argo-meteorological stations, a 1 km-grid NPP dataset, statistical yields, and administrative data were used in this study, as shown in Table 1. The 1 km-grid crop phenological dataset, which was derived from figshare, was used to determine the planting area of winter wheat [28]. This included the key phenological date for three main crops (rice, wheat, and maize) of China from 2000 to 2015. The Day of Year (DOY) of greenup/emergence date, heading date, and the maturity date was recorded for wheat. Accuracy verification shows that the root mean square error (RMSE) averages of wheat phenological dates were 5.5 days, which met the requirements of planting area extraction [18].

Data Sources
A 1 km-grid crop phenological dataset, phenological observations from argo-meteorological stations, a 1 km-grid NPP dataset, statistical yields, and administrative data were used in this study, as shown in Table 1. The 1 km-grid crop phenological dataset, which was derived from figshare, was used to determine the planting area of winter wheat [28]. This included the key phenological date for three main crops (rice, wheat, and maize) of China from 2000 to 2015. The Day of Year (DOY) of green-up/emergence date, heading date, and the maturity date was recorded for wheat. Accuracy verification shows that the root mean square error (RMSE) averages of wheat phenological dates were 5.5 days, which met the requirements of planting area extraction [18].
The months of sowing and maturity of winter wheat were determined using phenological observations from argo-meteorological stations from 1994 to 2013, which were derived from the China Integrated Meteorological Information Service System (CIMISS). These data included crop phenological information, such as crop growth period and dates of phenological stages. This information was observed and recorded by well-trained agricultural technicians in the experimental field. This study extracted the average annual DOYs of sowing and maturity stages. A total of 58 stations were used in this study.
The 1 km-grid NPP dataset was derived from the Global Change Research Data Publishing and Repository. This dataset was collected to calculate the cumulative NPP during the growing season of winter wheat [29]. This dataset was calculated using the Carnegie Ames Stanford approach (CASA) model, which has been widely used to estimate NPP, with a monthly temporal resolution [30].
Statistical yield data included winter wheat yield at the provincial and city scale from 2000 to 2015. The provincial yield data were derived from the Nation Bureau of Statistics of China. The data accuracy was verified using yield data collected at the city scale. These data were derived from the locality statistical yearbook.

Methods
The flowchart of this study is shown in Figure 2. First, the planting areas were extracted based on a 1 km-grid crop phenological dataset along with the months of sowing and maturity that were determined using phenological observations. Second, the cumulative NPP during the growing season of winter wheat was calculated based on the planting area, the months of sowing and maturity, and the 1 km-grid NPP dataset. Third, the yield spatialization model was established using the NPP of the growing season, along with the yield gridded dataset generated from 2000 to 2015. Finally, the accuracy of the yield gridded dataset was verified. The months of sowing and maturity of winter wheat were determined using phenological observations from argo-meteorological stations from 1994 to 2013, which were derived from the China Integrated Meteorological Information Service System (CIMISS). These data included crop phenological information, such as crop growth period and dates of phenological stages. This information was observed and recorded by well-trained agricultural technicians in the experimental field. This study extracted the average annual DOYs of sowing and maturity stages. A total of 58 stations were used in this study.
The 1 km-grid NPP dataset was derived from the Global Change Research Data Publishing and Repository. This dataset was collected to calculate the cumulative NPP during the growing season of winter wheat [29]. This dataset was calculated using the Carnegie Ames Stanford approach (CASA) model, which has been widely used to estimate NPP, with a monthly temporal resolution [30].
Statistical yield data included winter wheat yield at the provincial and city scale from 2000 to 2015. The provincial yield data were derived from the Nation Bureau of Statistics of China. The data accuracy was verified using yield data collected at the city scale. These data were derived from the locality statistical yearbook.

Methods
The flowchart of this study is shown in Figure 2. First, the planting areas were extracted based on a 1 km-grid crop phenological dataset along with the months of sowing and maturity that were determined using phenological observations. Second, the cumulative NPP during the growing season of winter wheat was calculated based on the planting area, the months of sowing and maturity, and the 1 km-grid NPP dataset. Third, the yield spatialization model was established using the NPP of the growing season, along with the yield gridded dataset generated from 2000 to 2015. Finally, the accuracy of the yield gridded dataset was verified.

Extraction of Winter Wheat Planting Area
The 1 km-grid crop phenological dataset can be used to determine the planting area owing to different phenological dates, such as spring wheat and winter wheat. Luo et al. (2020) reported that green-up dates of winter wheat are earlier than the emergence dates of spring wheat [18]. Therefore,

Extraction of Winter Wheat Planting Area
The 1 km-grid crop phenological dataset can be used to determine the planting area owing to different phenological dates, such as spring wheat and winter wheat. Luo et al. (2020) reported that green-up dates of winter wheat are earlier than the emergence dates of spring wheat [18]. Therefore, the winter planting area could be extracted by setting a threshold of green-up dates. The green-up dates of winter wheat in China usually start in February to late March. Additionally, the annual average of green-up dates in the study area is recorded using phenological observations. The green-up dates ranged from 44 (DOY) in Jianhu station to 89 (DOY) in Qinxian station. The histogram of green-up dates was established using the 1 km-grid crop phenological dataset. Approximately 90% of the values were lower than 90. As a result, 90 was selected as the threshold of green-up dates to extract the planting area of winter wheat in this study.

Determination of the Months of Sowing and Maturity
The months of sowing and maturity were determined using phenological observations from Argo-meteorological stations from 1994 to 2013. First, the DOYs of sowing and maturity stages in each station were extracted every year. This was followed by calculating the annual averages of sowing and maturity stages that are shown in Figure 3. Third, the DOYs of sowing and maturity were converted to months. Finally, the Thiessen polygon method was used to obtain the spatial distribution of the sowing and maturity months. the winter planting area could be extracted by setting a threshold of green-up dates. The green-up dates of winter wheat in China usually start in February to late March. Additionally, the annual average of green-up dates in the study area is recorded using phenological observations. The greenup dates ranged from 44 (DOY) in Jianhu station to 89 (DOY) in Qinxian station. The histogram of green-up dates was established using the 1 km-grid crop phenological dataset. Approximately 90% of the values were lower than 90. As a result, 90 was selected as the threshold of green-up dates to extract the planting area of winter wheat in this study.

Determination of the Months of Sowing and Maturity
The months of sowing and maturity were determined using phenological observations from Argo-meteorological stations from 1994 to 2013. First, the DOYs of sowing and maturity stages in each station were extracted every year. This was followed by calculating the annual averages of sowing and maturity stages that are shown in Figure 3. Third, the DOYs of sowing and maturity were converted to months. Finally, the Thiessen polygon method was used to obtain the spatial distribution of the sowing and maturity months.

Estimation of NPP during the Growing Season
The sowing and maturity months were identified in 2.3.2. Subsequently, NPP during the growing season was estimated using Equation 1: where NPPG is the NPP during the growing season, ms and mm are the months of sowing and maturity, and NPPi is the ith pixel of the NPP.

Yield Spatialization Model
A linear regression model was used to simulate the spatial distribution of the yields. The model used the NPP to evaluate the yield of winter wheat. The model was calculated using Equation 2: where Yi is the statistical yield in the ith province, (NPPG)i,j is the jth pixel of NPPG in the ith province, is the yield distribution coefficient, and b is a constant. According to the principle of "no land, no yield", b equals to 0.

Estimation of NPP during the Growing Season
The sowing and maturity months were identified in Section 2.3.2. Subsequently, NPP during the growing season was estimated using Equation (1): where NPP G is the NPP during the growing season, m s and m m are the months of sowing and maturity, and NPP i is the ith pixel of the NPP.

Yield Spatialization Model
A linear regression model was used to simulate the spatial distribution of the yields. The model used the NPP to evaluate the yield of winter wheat. The model was calculated using Equation (2): where Y i is the statistical yield in the ith province, (NPP G ) i,j is the jth pixel of NPP G in the ith province, a is the yield distribution coefficient, and b is a constant. According to the principle of "no land, no yield", b equals to 0.

Methods of Accuracy Verification
The gridded yield dataset was simulated at the provincial scale. The most credible method to verify the accuracy of the results was to choose the lower level of the province. Therefore, statistical yields at the city scale were selected. This study randomly selected 10 cities each year to verify the accuracy because of data availability. Mean prediction error (MPE) in Equation (3) and relative error (RE) in Equation (4) were used for accuracy verification: where Y s and Y c are the simulated yield and statistical yield, respectively, i represents the ith city, and m is the number of cities.

Trend Analysis
Trend analysis was used to analyze the annual change rate for winter yields from 2000 to 2015: where slope is the change rate and yield i represents the yield in the ith year. n is the number of years, which equals 16 in this case. An F test was used to test the significance.

Spatial Distribution of NPP during the Growing Season of Winter Wheat from 2000 to 2015
The spatial distribution characteristics of the NPP during the growing season of winter wheat from 2000 to 2015 were demonstrated by calculating the annual averages of NPP ( Figure 4). The results showed significant differences in NPP for different regions. The NPP during the growing season of winter wheat from 2000 to 2015 was concentrated between 0.06 and 640 gC·m −2 , with an average of 165.39 gC·m −2 . The high NPP was mainly concentrated in the North China Plain, especially in southwestern Shandong, southeastern Henan, northern Anhui, and northern Jiangsu, where the NPP was more than 200 gC·m −2 . Additionally, the distribution of areas with high NPP (more than 200 gC·m −2 ) was relatively scattered over some areas, such as Yunnan, Guizhou and Sichuan. The low NPP was mainly distributed in eastern Gansu, central Shanxi, southern Hebei, and eastern Sichuan, where the NPP was below 100 gC·m −2 . Consequently, the NPP during the growing season could reflect the spatial heterogeneity of winter wheat.

Relationship between NPP and Yield
The indicator for the winter wheat yield was identified by selecting the cumulative NPP and statistical yield at the provincial scale to analyze the relationship ( Figure 5). The cumulative NPP and statistical yield were positively related with an R 2 of 0.926, and passed the 0.01 significance test. The Cook's distance ranged from 0 to 0.283 (less than 1), with an average value of 0.10, indicating that there are no outliers in the samples. Regression standardized residuals showed that the sample values corresponded naturally with the normal distribution and the sample observations are approximately in line with the assumption. Consequently, NPP is a beneficial indicator when evaluating the winter wheat yield. Furthermore, these results are also consistent with previous research [4,31,32]. NPP was selected to simulate the spatial distribution of the yield in this study.

Relationship between NPP and Yield
The indicator for the winter wheat yield was identified by selecting the cumulative NPP and statistical yield at the provincial scale to analyze the relationship ( Figure 5). The cumulative NPP and statistical yield were positively related with an R 2 of 0.926, and passed the 0.01 significance test. The Cook's distance ranged from 0 to 0.283 (less than 1), with an average value of 0.10, indicating that there are no outliers in the samples. Regression standardized residuals showed that the sample values corresponded naturally with the normal distribution and the sample observations are approximately in line with the assumption. Consequently, NPP is a beneficial indicator when evaluating the winter wheat yield. Furthermore, these results are also consistent with previous research [4,31,32]. NPP was selected to simulate the spatial distribution of the yield in this study.

Relationship between NPP and Yield
The indicator for the winter wheat yield was identified by selecting the cumulative NPP and statistical yield at the provincial scale to analyze the relationship ( Figure 5). The cumulative NPP and statistical yield were positively related with an R 2 of 0.926, and passed the 0.01 significance test. The Cook's distance ranged from 0 to 0.283 (less than 1), with an average value of 0.10, indicating that there are no outliers in the samples. Regression standardized residuals showed that the sample values corresponded naturally with the normal distribution and the sample observations are approximately in line with the assumption. Consequently, NPP is a beneficial indicator when evaluating the winter wheat yield. Furthermore, these results are also consistent with previous research [4,31,32]. NPP was selected to simulate the spatial distribution of the yield in this study.

Comparison of Regression Parameters from 2000 to 2015
The regression parameters were calculated using cumulative NPP and statistical yields at the provincial scale from 2000 to 2015 ( Table 2). The results showed that the coefficient of the model ranged from 2.610 to 3.429. The R 2 of the model was higher than 0.92 from 2000 to 2015, which indicated that the fitting effects of the model in each year were good. Additionally, both models passed the significance test (p < 0.001), which suggested that the equations were reliable. Therefore, the fitting parameters were used to simulate winter wheat yields from 2000 to 2015.

Accuracy Verification
Due to data availability, it is difficult to obtain verification data at the pixel scale. In this study, statistical yield at the city scale were collected. To verify the accuracy of the gridded yield dataset, a total of 160 cities (10 cities per year) were randomly selected as samples. Total production of each city was calculated based on the gridded yield dataset. Figure 6 shows comparisons between the gridded yield dataset and statistical yields at the city scale from 2000 to 2015.
The results showed that the simulated and statistical yields are linearly correlated ( Figure 6). R 2 values ranged from 0.90 to 0.99, indicating that simulated and statistical yields had high consistency at the city scale. In addition, MPE values ranged from 4.23% in 2014 to 24.45% in 2007, with an average MPE of 12.01%, indicating that the gridded yield dataset simulated by this method had high precision at the city scale. Furthermore, the REs of the gridded yield dataset for approximately half of the samples were between −10% and 10%. The number of cities with absolute RE values below 10%, from 10 to 20%, and above 20% were 78, 53, and 29, respectively. It is also indicated that the accuracy satisfied the requirements for simulating regional winter wheat yield. Thus, these results can be assumed as credible. The results showed that the simulated and statistical yields are linearly correlated ( Figure 6). R 2 values ranged from 0.90 to 0.99, indicating that simulated and statistical yields had high consistency at the city scale. In addition, MPE values ranged from 4.23% in 2014 to 24.45% in 2007, with an average MPE of 12.01%, indicating that the gridded yield dataset simulated by this method had high precision at the city scale. Furthermore, the REs of the gridded yield dataset for approximately half of the samples were between −10% and 10%. The number of cities with absolute RE values below 10%, from 10 to 20%, and above 20% were 78, 53, and 29, respectively. It is also indicated that the accuracy satisfied the requirements for simulating regional winter wheat yield. Thus, these results can be assumed as credible. Figure 7 shows the spatial distribution of winter wheat yields from 2000 to 2015. The results showed that the distribution of yields has apparent spatial differences. The yield had a pattern that indicated a high yield in the east and a low yield in the west. Similar to NPP, the high yield was mainly concentrated in the North China Plain, especially in southwestern Shandong, southeastern Henan, northern Anhui, and northern Jiangsu, where the yield was more than 6 t/ha. This region is the main winter wheat production area of China. Abundant light, heat, and water resources are the favorable factors for crop production [33,34]. The low yield was mainly distributed in eastern Gansu, central Shanxi, southern Hebei, and eastern Sichuan, where the yield was below 4 t/ha. Complicated topography and light deficiency are the major factors restricting yields from increasing in these regions.  The results showed that the distribution of yields has apparent spatial differences. The yield had a pattern that indicated a high yield in the east and a low yield in the west. Similar to NPP, the high yield was mainly concentrated in the North China Plain, especially in southwestern Shandong, southeastern Henan, northern Anhui, and northern Jiangsu, where the yield was more than 6 t/ha. This region is the main winter wheat production area of China. Abundant light, heat, and water resources are the favorable factors for crop production [33,34]. The low yield was mainly distributed in eastern Gansu, central Shanxi, southern Hebei, and eastern Sichuan, where the yield was below 4 t/ha. Complicated topography and light deficiency are the major factors restricting yields from increasing in these regions. The spatial distribution of the annual change rate of winter yields from 2000 to 2015 is shown in Figure 8. Notably, areas that passed the significance test (p < 0.05) participated in the analysis. The trends mainly showed an increase of an average rate of 0.17 t ha -1 yr -1 from 2000 to 2015 in the study area. There was also significant spatial heterogeneity in the trends of winter wheat yield in different regions. The region that increased was mainly concentrated in the North China Plain, while the distribution of the region that decreased was relatively dispersed. The difference in the North China Plain between the north and the south is apparent. Additionally, the increase rate is larger in the south (mainly >0.2 t ha -1 yr -1 ) than that in the north (mainly < 0.2 t ha -1 yr -1 ). The yield benefited from many factors, such as climate warming, cultivars shift, and fertilization [33,35,36]. The spatial distribution of the annual change rate of winter yields from 2000 to 2015 is shown in Figure 8. Notably, areas that passed the significance test (p < 0.05) participated in the analysis. The trends mainly showed an increase of an average rate of 0.17 t ha −1 yr −1 from 2000 to 2015 in the study area. There was also significant spatial heterogeneity in the trends of winter wheat yield in different regions. The region that increased was mainly concentrated in the North China Plain, while the distribution of the region that decreased was relatively dispersed. The difference in the North China Plain between the north and the south is apparent. Additionally, the increase rate is larger in the south (mainly > 0.2 t ha −1 yr −1 ) than that in the north (mainly < 0.2 t ha −1 yr −1 ). The yield benefited from many factors, such as climate warming, cultivars shift, and fertilization [33,35,36].

Comparisons with Other Studies
This study proposed a novel yield spatialization method based on a crop phenological dataset, phenological observations, and NPP. Generally, yield estimation methods can be divided into two categories that include yield spatialization and remote sensing estimation.
The yield spatialization method uses statistical yield, which is an important data source obtained by governments.  16,17]. Although the study object of these studies was not winter wheat, it still provides a spatialization method to this study. Additionally, this novel method has obvious advantages. This method extensively uses remote sensing images that include 1 km-grid crop phenological and NPP datasets. The yield gridded dataset can reflect the spatial heterogeneity of winter wheat very well and has a high temporal resolution.
Remote sensing estimation is another yield estimation solution. The relationship between NPP and actual crop yield is established, and the harvest index (HI) is critical to estimate yield accurately.  [4]. HIs were assigned to constant values in these studies. However, this parameter varies considerably because of regional differences, a variety of crops, environmental conditions, and management [7,37]. It could affect the accuracy of yield estimation on a regional scale. Although there are a few uncertainties in HI, NPP is a critical indicator to estimate yield on a regional scale [4,32]. The statistical yield was selected to establish the relationship between NPP and yield in this study, which can reduce a few of the uncertainties.

Comparisons with Other Studies
This study proposed a novel yield spatialization method based on a crop phenological dataset, phenological observations, and NPP. Generally, yield estimation methods can be divided into two categories that include yield spatialization and remote sensing estimation.
The yield spatialization method uses statistical yield, which is an important data source obtained by governments.  16,17]. Although the study object of these studies was not winter wheat, it still provides a spatialization method to this study. Additionally, this novel method has obvious advantages. This method extensively uses remote sensing images that include 1 km-grid crop phenological and NPP datasets. The yield gridded dataset can reflect the spatial heterogeneity of winter wheat very well and has a high temporal resolution.
Remote sensing estimation is another yield estimation solution. The relationship between NPP and actual crop yield is established, and the harvest index (HI) is critical to estimate yield accurately.  [4]. HIs were assigned to constant values in these studies. However, this parameter varies considerably because of regional differences, a variety of crops, environmental conditions, and management [7,37]. It could affect the accuracy of yield estimation on a regional scale. Although there are a few uncertainties in HI, NPP is a critical indicator to estimate yield on a regional scale [4,32]. The statistical yield was selected to establish the relationship between NPP and yield in this study, which can reduce a few of the uncertainties. Furthermore, the yield gridded dataset with statistical yields at the city level shows that the accuracy satisfies the requirements for simulating regional winter wheat yield. Therefore, regardless of the current method, this is a flexible approach for future research.

Limitations
The planting area of winter wheat might affect the accuracy of the simulation of the yield. This study extracted the planting area based on a 1 km-grid crop phenological dataset. However, it is complicated for the cultivated system in China [24,38]. There are often multiple objects in a 1 km-grid, such as crops, trees, and roads. It is inevitable to extract the planting area of winter wheat with mixed pixels. This dataset has a sufficient quality for China because of data availability. Therefore, this dataset can be used to identify the planting area of winter wheat with high spatial and temporal resolution in China. Although this may have mixed pixel effects, a feasible and flexible approach for planting area extraction at the regional scale is provided.
NPP is a key parameter of this novel method, and the accuracy of the method depends on the precision of NPP during the growing season, which was calculated using planting area, phenological data, and an NPP dataset. This study considers the monthly temporal resolution of the NPP dataset. The sowing and maturity months were determined using the average annual DOYs. This can lead to a few uncertainties in this study. Studies have shown that the phenology of winter wheat has changed under climate change [39,40]. Xiao et al. (2013) found that the date of sowing was delayed (by 1.5 days per decade), whereas the date of maturity advanced (by 1.4 days per decade) in the North China Plain from 1981 to 2009 [41]. Liu et al. (2018) reported that the dates of sowing and maturity for wheat showed a similar trend in China from 1981 to 2010 [42]. Additionally, human management, such as sowing adjustment, irrigation, and fertilization, can affect the phenology of crops [42][43][44]. Earlier sowing for crops could extend the growing season, and result in a higher production [44]. Irrigation could increase soil water content, and a suitable irrigation method could promote crop growth and yield [45,46]. Li et al. (2018) pointed out that warming could decreased the wheat yield in dry years, while irrigation nullified the warming effect [47]. Reasonable fertilization could improve soil fertility and increase crop yield [48][49][50]. Although phenological dates of winter wheat may change at different times and locations, there is minimal change in NPP products. Therefore, the average annual DOYs were selected. Furthermore, the method was only tested in the winter wheat of China in this study. Therefore, the application of this method to other crops should be further investigated. Due to data availability, accuracy verification was carried out at the city scale. Fine-scale gridded yield or field data should be employed in accuracy verification in the future. The model accuracy can be improved using fine-scale statistical yields (i.e., city and county level data) in this model. Additionally, the distribution of phenological dates and NPP with high spatial and temporal resolution should be considered in the future.

Potential Applications
Regional yield estimation is critical for many applications, such as agricultural lands management, food security warning systems, agricultural trade policy, and carbon cycle research [22,51,52]. Generally, crop yield is influenced by the interaction of climate resources and human activities [53]. Climate resources that include temperature, precipitation, and light have changed significantly in China because of climate change. The impact of climate change on crop yield has been of considerable attention [32,[54][55][56]. Furthermore, human activities that include technological progress, irrigation, and fertilization, can promote crop yield [57,58]. However, crop production can cause a series of challenges, which include irrational utilization of land, resource consumption, and environmental deterioration [59][60][61][62]. Therefore, coordinating the relationship between population, resources, environment, and crop production is essential for ensuring sustainable agricultural development [59,62]. Gridded yield datasets can show the spatial characteristics of distributions and facilitate applications using spatial analysis when compared with statistical yield. This can provide credible and fundamental data that increase the understanding of the relationship between climate change and yield. Furthermore, human-land relationships can be coordinated to promote sustainable agricultural development at the pixel scale.

Conclusions
Many studies have developed yield spatialization models that are based on cropland areas. However, minimal attention has been received by planting areas, phenological dates, and the NPP of crops. Therefore, this study proposed a novel yield spatialization method that is based on planting area, phenological data, and NPP.
The average annual NPP during the growing season of winter wheat was 165.39 gC·m −2 . The NPP and yield were positively related, with an R 2 of 0.926 (p < 0.01). This indicates that NPP is an advantageous indicator to evaluate the yield of winter wheat. A gridded yield dataset was generated by the novel method that has high precision (MPE = 12.01%) at the city scale. The number of cities with absolute RE values below 10%, from 10 to 20%, and above 20% were 78, 53, and 29, respectively. The distribution of yields has apparent spatial differences, and the yield indicated a high yield in the east and a low yield in the west. The high yield was mainly concentrated in the North China Plain (more than 6 t/ha), whereas the low yield was mainly distributed in eastern Gansu, central Shanxi, southern Hebei, and eastern Sichuan (lower than 4 t/ha). The yield trend increased with an average rate of 0.17 t ha −1 yr −1 in the study area from 2000 to 2015. The region that increased was mainly concentrated in the North China Plain. Therefore, this study suggested that NPP is a key indicator of the yield spatialization model. Additionally, the method can generate gridded yield datasets with higher accuracy, which could show the spatial and temporal variation of yields.
Author Contributions: All of the authors have contributed to the manuscript. D.H. analyzed the results and wrote the manuscript. X.Y. conceived and designed the study. H.C. and X.Y. contributed to editing the manuscript and provided many suggestions. X.X. assisted with data processing. All authors have read and agreed to the published version of the manuscript. We also thank Zhi Qiao for his comments on this study.

Conflicts of Interest:
The authors declare no conflicts of interest.