Spatiotemporal Changes of Winter Wheat Planted and Harvested Areas, Photosynthesis and Grain Production in the Contiguous United States from 2008–2018

Winter wheat is a main cereal crop grown in the United States of America (USA), and the USA is the third largest wheat exporter globally. Timely and reliable in-season forecast and year-end estimation of winter wheat grain production in the USA are needed for regional and global food security. In this study, we assessed the consistency between the agricultural statistical reports and satellite-based data for winter wheat over the contiguous US (CONUS) at both the county and national scales. First, we compared the planted area estimates from the National Agricultural Statistics Service (NASS) and the Cropland Data Layer (CDL) from 2008–2018. Second, we investigated the relationship between gross primary production (GPP) estimated by the vegetation photosynthesis model (VPM) and grain production from the NASS. Lastly, we explored the in-season utility of GPPVPM in monitoring seasonal production. Strong spatiotemporal consistency of planted areas was found between the NASS and CDL datasets. However, in the Southern Great Plains, both the CDL and NASS planted acreage were noticeable larger (>20%) than the NASS harvested area, where some winter wheat fields were used as forage for cattle grazing. County-level GPPVPM was linearly related with grain production of winter wheat, with an R2 value of 0.68 across the CONUS. The relationships between grain production and GPPVPM in those counties without a substantial difference (<20%) between planted and harvested area were much stronger and their harvest index (HIGPP) values ranged from 0.2–0.3. GPPVPM in May could explain about 70–90% of the variance of winter wheat grain production. Our findings highlight the potential of GPPVPM in winter wheat monitoring, especially for those high harvested/planted ratio, which could provide useful data to guide planning and marketing for decision makers, stakeholders, and the public.


Introduction
Wheat is one of the most widely cultivated grain crops in the world, covering approximately one sixth of the total arable land area, and wheat grain contributes~20% of total dietary calories worldwide [1]. Wheat is also widely traded in the world food market with approximately~23% of the world's annual wheat grain production being traded internationally. The United States of America (USA) ranks the fifth for wheat production (after the European Union, China, India, and Russia) and the third for wheat grain export (after Russia and the European Union) [2]. Winter wheat dominates the USA wheat grain production, accounting for approximately 80% of national wheat production [3]. The interannual fluctuation of winter wheat production in the USA could have significant impacts GPP and grain production in the CONUS and on the performance of GPP for in-season grain production forecasting for winter wheat.
In this study, as one of a series of papers to explore the relationship between GPP VPM and grain production, we followed the general framework of the previous paper [33]. However, this paper focused on winter wheat croplands. First, we investigated the spatiotemporal dynamics of winter wheat cropping areas in the CONUS at national and county scales from 2008 to 2018. We quantified the spatiotemporal consistencies between NASS and CDL datasets for planted areas and for planted/harvested areas. Differences between winter wheat planted area and harvested area were expected because winter wheat fields can be used as dual-purpose fields: beef cattle grazing and/or grain production. This is a typical phenomenon in the Southern Great Plains (SGP). Second, we hypothesize that the relationships between winter wheat GPP and grain production are strong and linear. We analyzed the spatiotemporal dynamics of winter wheat GPP and grain production over the CONUS at national and county scales from 2008-2018. The GPP data are derived from the VPM (GPP VPM ), which were evaluated with GPP data from the in situ cropland eddy flux tower sites (GPP EC ) [30]. Third, we explored the relationships between GPP VPM and NASS grain production and quantified the ratio between GPP and grain production, namely, GPP-derived harvest index (HI GPP ). Finally, we assessed the linear regression models between county-level cumulated GPP over time and annual grain production and explored the potential to monitor winter wheat grain production in the CONUS within the growing season of winter wheat.

Study Area
The CONUS includes 3233 counties and 48 states. The climate is diverse, ranging from a temperate climate in the northern region to a subtropical climate in the southern region (e.g., Florida). Winter wheat is primarily cultivated in the Great Plains ( Figure 1). The Southern Great Plains (Texas, Oklahoma, and Kansas) are the largest winter wheat producing region, accounting for over 40% of national wheat grain production in the USA [34]. Winter wheat cultivation is largely rain-fed. Note that in the Southern Great Plains, varying areas of winter wheat fields are used as dual-purpose fields: beef cattle grazing and/or grain production, dependent upon weather, market conditions, and other factors [35].

Winter Wheat Planted and Harvested Areas, and Grain Production Data from 2008-2018 from the USDA NASS Statistical Dataset
The USDA National Agricultural Statistic Service (NASS) Quick Stats database (https://quickstats.nass.usda.gov/, accessed on 28 April 2021) provides annual statistics on crop acreages (planted and harvested area) and grain production at national and county scales. The NASS crop statistical data are derived from the surveys of farm operators, grain processors, commercial storage firms, in situ yield samples, etc. As an example, the NASS winter wheat planted area and harvested area estimates come mostly from the Agricultural Survey in December and June. Farmers are asked by NASS survey enumerators about the acreage they planted and the acreage they plan to harvest [3]. The grain yield statistical data were derived from both the Objective Yield Survey (OYS) and Agricultural Yield Survey (AYS). The AYS is based on farmers' reported yield and the OYS is based on crop biophysical measurements of selected samples. Grain production is calculated by yield estimates at survey dates and to-be-harvested area, and predicted by assuming normal weather conditions for the remaining part of the growing season.

Winter Wheat Planted and Harvested Areas, and Grain Production Data from 2008-2018 from the USDA NASS Statistical Dataset
The USDA National Agricultural Statistic Service (NASS) Quick Stats database (https://quickstats.nass.usda.gov/, accessed on 28 April 2021) provides annual statistics on crop acreages (planted and harvested area) and grain production at national and county scales. The NASS crop statistical data are derived from the surveys of farm operators, grain processors, commercial storage firms, in situ yield samples, etc. As an example, the NASS winter wheat planted area and harvested area estimates come mostly from the Agricultural Survey in December and June. Farmers are asked by NASS survey enumerators about the acreage they planted and the acreage they plan to harvest [3]. The grain yield statistical data were derived from both the Objective Yield Survey (OYS) and Agricultural Yield Survey (AYS). The AYS is based on farmers' reported yield and the OYS is based on crop biophysical measurements of selected samples. Grain production is calculated by yield estimates at survey dates and to-be-harvested area, and predicted by assuming normal weather conditions for the remaining part of the growing season.

Winter Wheat Planted Area Data from the USDA NASS Cropland Data Layer Dataset (CDL)
The annual CDL dataset at a 30 m spatial resolution includes more than one hundred crop types in the CONUS [15]. The CDL dataset is produced by using a machine learning method, in situ reference data, and multi-date satellite images to identify crop types. The training data used in the crop classification include the National Land Cover Dataset (NLCD) and the USDA Farm Service Agency (FSA) Common Land Unit (CLU) data. The CLU attributes include the crop types and acreages reported by producers to the FSA county offices in the early growing season. The CDL classifier utilizes multiple temporal  The annual CDL dataset at a 30 m spatial resolution includes more than one hundred crop types in the CONUS [15]. The CDL dataset is produced by using a machine learning method, in situ reference data, and multi-date satellite images to identify crop types. The training data used in the crop classification include the National Land Cover Dataset (NLCD) and the USDA Farm Service Agency (FSA) Common Land Unit (CLU) data. The CLU attributes include the crop types and acreages reported by producers to the FSA county offices in the early growing season. The CDL classifier utilizes multiple temporal remote sensing images from multiple sensors, including Landsat TM/ETM+, AWiFS, Deimos-1, UK-DMC-2, and MODIS images [15]. The spatial resolution of the CDL dataset is 30 m from 2010-2018 and 56 m from 2008-2009. In 2018, the CDL data for 2008 and 2009 were reprocessed to 30 m. The CDL data from 2008 to 2018 at a 30 m spatial resolution were used, which ensures the spatial resolution is consistent during the study period.
The CDL dataset provides detailed spatial information for individual crop types, and high classification accuracies (>90%) for major crop types (e.g., winter wheat, soybean, and maize) were reported [15]. Numerous agriculture-related studies and land cover change studies have used the CDL dataset [36][37][38][39]. CropScape, which is a web-based data portal, provides interactive tools for people to download, visualize, and analyze the CDL data in a more effective and efficient way [40]. The annual total planted area of winter wheat in each county was calculated from CropScape. The national total areas in each year were then calculated as the total of annual planted crop areas over all the counties in the nation.

Gross Primary Production Estimates for Winter Wheat from the Vegetation Photosynthesis Model (GPP VPM )
The VPM estimates the daily GPP from the amount of PAR absorbed by chlorophyll (APAR chl ) and the light use efficiency [41,42]. GPP derived from the eddy-covariance tower sites (GPP EC ) was used to assess daily GPP VPM in different croplands, including winter wheat [43,44], soybean [45,46], paddy rice [47], maize [45,48], and sugarcane [49]. Strong temporal consistency between GPP EC and GPP VPM was reported in the aforementioned publications, with R 2 values ranging from 0.70 to 0.98. The driving data of the VPM include vegetation indices (land surface water index (LSWI) and EVI) and meteorological data (air temperature, shortwave radiation). The theory, model structure, and parameters of the VPM were documented in [30,31]. We followed the same procedure to run VPM simulations as reported in an earlier publication [30], which used different maximum light use efficiency (LUE 0 ) parameter values for C 3 and C 4 crops to improve GPP estimation in croplands, as follows: where FPARchl is the fraction of PAR absorbed by chlorophyll in the canopy; LUE 0 is the maximum light use efficiency without environmental stress; Tscalar and Wscalar are the temperature and water limitation scalars; FA C3 and FA C4 are the area fractions of C 3 and C 4 plants within each 500 m MODIS pixel (range of 0-1.0), LUE 0-C3 and LUE 0-C4 are the maximum LUE values for C3 and C4 plants, respectively. The growing season total GPP of winter wheat at each pixel was calculated by summarizing daily GPP between the USDA planted and harvested dates. The mean GPP values in each county during the growing season were area weighted according to the area fraction of winter wheat in a county derived from the CDL datasets within 500 m pixels. More details about the VPM and GPP calculation can be found in [33]. The growing season total GPP of winter wheat in each county was calculated by the mean GPP multiplied by the total area of winter wheat within all pixels in the county (Figure 1f).

Statistical Analysis
To explore the interannual changes of winter wheat cropping areas and grain production from 2008-2018, we calculated the anomaly of each variable as the difference between the value in a specific year and the multi-year average from 2008-2018, and then normalized it by the multi-year average. We used linear regression models to quantify the relationships between grain production and cropping areas, and between GPP VPM and grain production at national and county scales. The performance of models was assessed using the coefficient of determination (R 2 ), bias, and root mean square error (RMSE) between the modeled grain production estimates and the NASS grain production statistics. The regression slope between grain production and GPP VPM , representing the conversion coefficient from GPP VPM to grain production, was also termed as the GPP-based harvest index (HI GPP ).
Remote Sens. 2021, 13, 1735 6 of 20 where GP sim is the simulated grain production at a county scale based on linear models, GP sim is the mean GP sim across counties in the CONUS, i is the county number, and n is the total number of counties. As there are missing NASS statistics data for some counties over the years, we only kept those counties with continuous data from 2008-2018. The total number is 795. GP NASS is the NASS grain production, GP N ASS is the mean GP NASS across counties in the CONUS.

In-Season Forecasting of Winter Wheat Grain Production Using Cumulated GPP VPM Data
Accurate and timely prediction of grain production and yield is one requirement in crop production monitoring programs. To explore the in-season prediction skill from GPP VPM , linear regression models were applied between cumulated GPP (GPP cum ) over time and grain production at the county scale (see Equation (1)).
where t is the time intervals in a year (ranging from 1 to 46), as GPP VPM data are at an 8-day temporal resolution and have 46 estimates within each year; k is the number of days within a time interval, being 8 days for t = 1 to 45, and 5 or 6 days (non-leap year, leap year) for t = 46. In this study, we used a calendar year schedule to run the statistical models for all counties, though winter wheat is usually planted in late fall of the previous year (e.g., late September to October). The GPP cum in the previous fall season is usually low and was not included in Equation (1), as we tried to make the task simple, without delineating the exact planting dates of winter wheat in the fall season of the previous year. We also assumed that it has a negligible effect on the modeling skills.  (Figure 1a). The disagreements between the NASS planted area and CDL planted area were relatively small for most of the years (<10%). We calculated the anomaly for the mean values of 2008-2018, and the planted area from both NASS and CDL datasets showed a significant decrease from 2013 ( Figure 2b). The NASS harvested area had a similar change pattern to the planted area but with different magnitudes over the years (Figure 2b).
At the county scale, for all the counties with winter wheat planted acreages from both NASS and CDL datasets in the period of 2008-2018, winter wheat planted areas calculated from the CDL were highly consistent with those from the NASS's officially reported dataset (R 2 = 0.98, p < 0.001) (Figure 3a). There are relatively small interannual variations in the slope values (<10%) and R 2 values among the individual years (Table 1). In many counties, winter wheat harvested areas were much smaller than winter wheat planted areas (Figure 3b,c). Most of those counties with a large discrepancy between winter wheat planted area and harvested area were distributed in the Southern Great Plains and California, where some winter wheat fields were used as cool-season forage for beef cattle production (grazing or haying).
At the county scale, the interannual trends of winter wheat planted areas from 2008-2018 were calculated ( Figure 4). As the NASS dataset has missing data in each year for various counties, only 626 counties were included (Figure 4b many counties in Oklahoma, Texas, and several other midwest states had moderate to large gains of winter wheat planted area ( Figure 4a). For most states, there is a reasonable spatial agreement of loss and gain of winter wheat planted area between the NASS and CDL datasets, however, in Oklahoma and Texas, there are noticeable differences between these two datasets (Figure 4a,b). A number of counties in Oklahoma and Texas had losses of winter wheat planted area according to the NASS dataset but gains of winter wheat area according to the CDL dataset. This discrepancy could be related to the imprecise estimates of winter wheat area from the CDL dataset because of image spatial resolution, misclassifications, and misalignment of field boundaries [36]. At the county scale, for all the counties with winter wheat planted acreages from both NASS and CDL datasets in the period of 2008-2018, winter wheat planted areas calculated from the CDL were highly consistent with those from the NASS's officially reported dataset (R 2 = 0.98, p<0.001) (Figure 3a). There are relatively small interannual variations in the slope values (<10%) and R 2 values among the individual years (Table 1). In many counties, winter wheat harvested areas were much smaller than winter wheat planted areas (Figure 3b,c). Most of those counties with a large discrepancy between winter wheat planted area and harvested area were distributed in the Southern Great Plains and California, where some winter wheat fields were used as cool-season forage for beef cattle production (grazing or haying).    At the county scale, the interannual trends of winter wheat planted areas from 2008-2018 were calculated ( Figure 4). As the NASS dataset has missing data in each year for various counties, only 626 counties were included (Figure 4b,d). According to the CDL dataset, 2168 counties had a loss of winter wheat planted area and 998 counties had a gain of winter wheat planted area (Figure 4a). There was a geographic shift of winter wheat plant area, as shown by the loss and gain of winter wheat planted area at the county scale For all the counties with winter wheat grain production, from planted and harvested area data from 2008-2018, grain production had noticeably stronger relationships with harvested areas than with planted area ( Figure 5) due to part of the planted acreages being used for grazing and hay. Table 2 lists the statistics of the linear regression models (slope, R 2 , and RMSE) between grain production and cropping areas in the period of 2008-2018, and the poorest relationship occurred in the spring drought year of 2014. tial agreement of loss and gain of winter wheat planted area between the NASS and CDL datasets, however, in Oklahoma and Texas, there are noticeable differences between these two datasets (Figure 4a,b). A number of counties in Oklahoma and Texas had losses of winter wheat planted area according to the NASS dataset but gains of winter wheat area according to the CDL dataset. This discrepancy could be related to the imprecise estimates of winter wheat area from the CDL dataset because of image spatial resolution, misclassifications, and misalignment of field boundaries [36].  For all the counties with winter wheat grain production, from planted and harvested area data from 2008-2018, grain production had noticeably stronger relationships with harvested areas than with planted area ( Figure 5) due to part of the planted acreages being used for grazing and hay. Table 2 lists the statistics of the linear regression models (slope, R 2 , and RMSE) between grain production and cropping areas in the period of 2008-2018, and the poorest relationship occurred in the spring drought year of 2014.

Spatiotemporal Dynamics of GPP VPM and Grain Production from NASS Dataset from 2008-2018
At the national scale, winter wheat grain production showed a decreasing trend from 2008-2018 ( Figure 6), which is largely determined by the decrease in planted and harvested area (Figure 1). GPP VPM showed a similar decreasing trend as grain production from 2008 to 2018, but with strong response in some drought years, such as 2011 and 2014.

Spatiotemporal Dynamics of GPPVPM and Grain Production from NASS Dataset from 2008-2018
At the national scale, winter wheat grain production showed a decreasing trend from 2008-2018 (Figure 6), which is largely determined by the decrease in planted and harvested area (Figure 1). GPPVPM showed a similar decreasing trend as grain production from 2008 to 2018, but with strong response in some drought years, such as 2011 and 2014.
At the county scale, we quantified the interannual trends of GPPVPM and grain production from 2008-2018 (Figure 7). Out of 291 counties with continuous NASS grain production data during the 11 years, winter wheat grain production showed increasing trends in 90 counties and decreasing trends in 201 counties (Figure 7a,c). The interannual trend of NASS grain production had a strong linear relationship with the that of GPPVPM (Figure 7e). Out of 2609 counties in the CONUS, GPPVPM data show 821 counties with a substantial increasing trend and 1788 counties with a decreasing trend (Figure 7b,d).
Many counties with a decrease in GPPVPM were located in the Great Plains (South Dakota, Kansas, Oklahoma) and the states along the Mississippi River (Figure 7b,d). Many counties in the northwestern states and the western edge of the Southern Great Plains showed an increasing trend of GPPVPM (Figure 7b,d). The probability distributions of the interannual trends of NASS grain production and GPPVPM were similar, with narrow ranges of normal distribution (Figure 7f).  At the county scale, we quantified the interannual trends of GPP VPM and grain production from 2008-2018 (Figure 7). Out of 291 counties with continuous NASS grain production data during the 11 years, winter wheat grain production showed increasing trends in 90 counties and decreasing trends in 201 counties (Figure 7a,c). The interannual trend of NASS grain production had a strong linear relationship with the that of GPP VPM (Figure 7e). Out of 2609 counties in the CONUS, GPP VPM data show 821 counties with a substantial increas-ing trend and 1788 counties with a decreasing trend (Figure 7b,d). Many counties with a decrease in GPP VPM were located in the Great Plains (South Dakota, Kansas, Oklahoma) and the states along the Mississippi River (Figure 7b,d). Many counties in the northwestern states and the western edge of the Southern Great Plains showed an increasing trend of GPP VPM (Figure 7b,d). The probability distributions of the interannual trends of NASS grain production and GPP VPM were similar, with narrow ranges of normal distribution (Figure 7f).

The Relationships between County-Level GPP VPM and Winter Wheat Grain Production from 2008 to 2018
At the county scale, the relationships between GPP VPM and winter wheat grain production from 2008 to 2018 varied from one year to another (Figure 8a), with R 2 values ranging from 0.59 (2008) to 0.75 (2016) ( Table 3). The slope values in the models with no intercept (GP = HI GPP × GPP VPM ) varied from 0.25 (2012) to 0.33 (2011) ( Table 3). When all the county-year data in the CONUS were used in the linear regression model for GPP VPM and grain production (GP), GPP VPM could explain 68% of the variance in winter wheat grain production (Figure 8a). The relationships (Figure 8b) were also affected by the differences between the planted area calculated from CDL and harvested area from NASS statistics, which varied according to county and year (Figure 3). When including only those counties with differences less than 20% between the planted area from CDL and the harvested area from NASS, the relationship between GPP VPM and grain production was much stronger (R 2 > 0.8, Table 4).
Remote Sens. 2021, 13, 1735 13 of 21 ranging from 0.59 (2008) to 0.75 (2016) ( Table 3). The slope values in the models with no intercept (GP = HIGPP × GPPVPM) varied from 0.25 (2012) to 0.33 (2011) ( Table 3). When all the county-year data in the CONUS were used in the linear regression model for GPPVPM and grain production (GP), GPPVPM could explain 68% of the variance in winter wheat grain production (Figure 8a). The relationships (Figure 8b) were also affected by the differences between the planted area calculated from CDL and harvested area from NASS statistics, which varied according to county and year ( Figure 3). When including only those counties with differences less than 20% between the planted area from CDL and the harvested area from NASS, the relationship between GPPVPM and grain production was much stronger (R 2 > 0.8, Table 4). HIGPP values had a moderate range of variation across most of the county-years (Figure 8c), which was not affected much by the difference between the two cropping areas (Figure 8c). For those counties with small differences (<20%) between the two cropping areas, HIGPP values ranged from 0.2-0.3 (Figure 8d). Geographically, high HIGPP values occurred in the northern part of the CONUS, and the winter wheat belt region in the Southern Great Plains generally had a lower HIGPP, except for part of western Kansas (Figure 9).    Table 4. Statistics of linear regression results between county-level GPP VPM and winter wheat grain production from NASS over CONUS from 2008-2018 by considering the relative differences between CDL-derived planted area and NASS harvested area. HI GPP values had a moderate range of variation across most of the county-years (Figure 8c), which was not affected much by the difference between the two cropping areas (Figure 8c). For those counties with small differences (<20%) between the two cropping areas, HI GPP values ranged from 0.2-0.3 (Figure 8d). Geographically, high HI GPP values occurred in the northern part of the CONUS, and the winter wheat belt region in the Southern Great Plains generally had a lower HI GPP , except for part of western Kansas (Figure 9).   year slope R 2 bias (10 3 ton) RMSE (10 3 ton) Figure 9. Spatial distribution of harvest index derived from GPP VPM and NASS grain production (HI GPP ) in 2010 for (a) all the counties and (b) counties with small differences (<20%) between CDL-derived planted area and NASS harvested area.

In-Season Forecasting of Winter Wheat Grain Production Using Cumulative GPP Data
In the CONUS, winter wheat is usually planted from September-November, and harvested from June-August of the following year. We assessed the potential of the simple linear regression model in forecasting grain production based on cumulative GPP. When using all the county-year data in the CONUS, the model prediction skill increased over time, reaching 60% to 80% by the end of June (Figure 10a). When using all the county-year data in a state, the model prediction skills for those states located in the cold northern part of the CONUS, where there were few differences between CDL-derived planted area and NASS harvested area, reached 90% by the end of May (Figure 10b,c). For those states located in the Southern Great Plains with big differences between the two cropping areas, that is, CDL-derived planted area and harvested area from NASS statistics, the model prediction skill varies over the years, ranging from 70% to 90% (Figure 10d,e). After excluding those counties with a difference larger than 20% between the two cropping areas, the model prediction skill for the CONUS increases and varies between 80% and 90% over the years (Figure 10f).

In-Season Forecasting of Winter Wheat Grain Production using Cumulative GPP Data
In the CONUS, winter wheat is usually planted from September-November, and harvested from June-August of the following year. We assessed the potential of the simple linear regression model in forecasting grain production based on cumulative GPP. When using all the county-year data in the CONUS, the model prediction skill increased over time, reaching 60% to 80% by the end of June (Figure 10a). When using all the county-year data in a state, the model prediction skills for those states located in the cold northern part of the CONUS, where there were few differences between CDL-derived planted area and NASS harvested area, reached 90% by the end of May (Figure 10b,c). For those states located in the Southern Great Plains with big differences between the two cropping areas, that is, CDL-derived planted area and harvested area from NASS statistics, the model prediction skill varies over the years, ranging from 70% to 90% (Figure 10d,e). After excluding those counties with a difference larger than 20% between the two cropping areas, the model prediction skill for the CONUS increases and varies between 80% and 90% over the years (Figure 10f).

Spatiotemporal Dynamics of Winter Wheat Planted Area, GPP, and Grain Production
Many factors contributed to the spatiotemporal changes of winter wheat planting, GPP, and grain production, including the grain market volatility (e.g., international trade market, fluctuating market price), and climate-induced variation. For international trade, the US share of global exports for wheat was~20% before 2013, which then shrank tõ 15% in recent years, driven by rising exports from Russia, the European Union, and South America [50]. The global wheat price also experienced an increase from~5 dollars/bushel in 2008 to~9 dollars/bushel in 2012, then decreased to~USD 5/bushel in 2018 [51]. The interannual variation of winter wheat planted area, GPP, and grain production is highly consistent with the wheat price dynamics (Figures 2 and 6). Domestically, as winter wheat is usually planted from September to October, and is harvested in June and July after maize and soybean planting, this makes it compete directly with maize and soybean. With less profit created by winter wheat over the years, maize and soybeans are preferred by US farmers. The recurring droughts throughout the Great Plains also exacerbated the shrinkage of winter wheat acreage. The 2011 drought in the Southern Great Plains and the 2017 drought in the Northern Great Plains impacted wheat growing, resulting in a large production reduction in 2012 ( Figure 6). After the 2017 drought, wheat was replaced by soybean as the predominant crop in North Dakota, and more US states are increasingly planting maize and soybeans rather than wheat. For example, in Kansas, the largest wheat growing state in the US, wheat planted areas become smaller than maize and soybean areas in 2017 [52].

Spatiotemporal Consistency of Winter Wheat Cropping Areas from the CDL and NASS Datasets
Crop planted area and its spatial distribution are the first data layer for crop production monitoring and forecasting. A few studies compared the planted area derived from the CDL dataset and NASS statistics over a few years or over a few states [15,36]. The CDL pixel-based planted area estimates generally showed a satisfying accuracy for major crops (85-95% for maize, soybean, and wheat), but had a tendency of underestimation compared with NASS planted area. The pixel counting estimation of crop acreage from the CDL may underestimate crop acreages for certain counties and states but may yield overestimation for other states and counties. It should be indicated that the R 2 value does not indicate that the estimation is biased up or down. In this study, the winter wheat planted area estimates in the CONUS for the CDL dataset showed a strong spatiotemporal consistency (R 2 = 0.98) with those from the NASS dataset from 2008-2018 at the county scale ( Figure 3). This high spatiotemporal consistency can be largely attributed to the accuracy of the CDL and NASS statistical data, especially for major crop types (winter wheat, soybean, and maize) [15] in the major crop production areas. The high accuracy and robustness of the CDL dataset for winter wheat planted area estimation makes it reliable for crop monitoring.
Crop harvested area and its spatial distribution were directly related to crop production because some crop fields could not be harvested due to damage, failure, and other factors. However, there is still a lack of pixel-based harvested area at a regional scale. Thus, it is important to calculate the difference between planted and harvested area for monitoring crop production. In our study, we found that the spatiotemporal consistency of winter wheat planted and harvested area in the CONUS varied substantially, especially in the Southern Great Plains and in some spring drought years (2011, 2013-2014), with a difference of more than 20% (Figures 2 and 3). The Southern Great Plains are one of the key regions for the nation's wheat and beef production. In many counties of the Southern Great Plains region, winter wheat is used as part of a dual purpose graze-grain system, serving as both a grain crop and a forage crop to supply beef cattle [35,53]. Drought is an important factor that can lead to the loss of grain yield and abandonment of winter wheat acreage in the Southern Great Plains. Winter wheat crops in this area normally break from dormancy in late February or early to mid-March. A lack of precipitation or snow before the dormancy period could significantly reduce the growth of winter wheat crops [54]. These factors could have led to the significant differences between planted and harvested area, especially over the Southern Great Plains.

Harvest Index-The Relationship between Winter Wheat Grain Production and GPP VPM
The "harvest index" (HI) is used in crop production research in different ways: the ratio of crop grain production over (1) aboveground biomass (AGB), namely HI AGB [28], (2) net primary production (NPP), namely HI NPP [27,29], and (3) gross primary production (GPP), namely HI GPP [26,33]. These three His are different but related to each other, as GPP, NPP, and AGB are closely related. GPP represents the total carbon fixed by photosynthesis during the growing season. NPP is part of GPP, deducted by autotrophic respiration, and is often estimated as a sum of aboveground biomass and belowground biomass at harvest time. AGB is the total biomass allocated to the leaves, stems, branches, and seeds. The grain production (yield) is part of AGB at harvest. HI AGB derived from field-level studies showed a large variability in previous studies [28,55], which is contributed to the seed variety and growing condition. HI NPP is usually calculated from HI AGB and the root:shoot ratio at the field scale. The root:shoot ratio for wheat ranges from 0.19-0.21 [28,56]. For HI GPP , Done et al. [57] and He et al. [26] reported a value of 0.29-0.37 based on modeled GPP and yield data at a regional scale. Our results showed a median HI GPP of 0.30 from 2008-2018, and it fluctuated over counties and years ( Figure 5). This variation of the GPP-yield relationship (HI GPP ) could result from environmental factors, crop management, and crop varieties (genomics). For example, winter wheat grain production in the Southern Great Plains is relatively low and more variable than for other regions, because some winter wheat fields are used as pasture for cattle grazing. The rate of abandonment could vary in some extreme climate years, as adverse weather reduces wheat yield. In addition, fertilizer application, tillage, and row spacing could also affect the yield both spatially and temporally [58]. Additional studies are needed to quantify the responses of GPP and grain production under extreme climate events.

In-Season Forecasting for Winter Wheat Grain Production
The LUE-based GPP models provide a simple but efficient way to estimate biomass or yield [59]. Unlike empirical models, they are semi-empirical and can be applied at both site or regional scales with calibration and validation [59,60]. Our study combined the GPP estimated from the VPM and CDL crop area and to predict winter wheat grain production. At the county and state scales, the predicted grain production estimates correlated well with the NASS grain production. This in-season forecasting pilot study was run at county and state scales because the NASS grain production data are readily available at those scales. Lack of access to the field-scale grain production data limits the validation of GPP VPM -based grain production data at the field scale. GPP data at a moderate spatial resolution (500 m) were used in this study. For the major production regions of winter wheat, which produced more than 40% of the total winter wheat grain, the field sizes are usually large and there are few mixed pixels in these counties and states. However, for other regions, the mixed pixels problem can cause significant discrepancies in calculating wheat GPP and grain production. GPP products with high spatial resolutions (e.g., 30 m Landsat, and 10 m Sentinel-2) and suitable temporal frequency (weekly) may overcome the above limitations for field-scale agricultural applications.
This study highlighted the importance of accurate and in-season or annual maps of crop types for the prediction of crop production, including both planted and harvested areas. Annual or in-season maps of crop planted areas provided the basic information at the beginning of the growing season. Annual to-be-harvested area maps in the growing season would be helpful in predicting grain production as grain production is calculated as a function of "to-be-harvested" crop areas. For some major crops like maize and soybean, there could be no significant differences between planted area and harvested area in the CONUS [33]. Thus, early mapping of the planted area in the growing season could be used for grain production estimation. However, for winter wheat, our study showed that there were significant differences between planted and harvested areas over the CONUS, especially the Southern Great Plains, where the relative difference could be greater than 20%. Though numerous studies have tried to provide in-season crop type maps as early as possible based on remote sensing data [61,62], there is not much attention paid to the research of planted and harvested area mapping. A separate classification of planted area and harvested area can help to get a more accurate prediction of grain production. Another issue is to have updated (in-season) and high spatial resolution (e.g., 30 m or 10 m) crop classification maps. The CDL dataset usually has a six-month lag time after harvest, before being released to the public, and a delay of six months after harvest would preclude in-season grain production prediction. Time series satellite observations at high spatial resolutions, for example, from Sentinel-2 (5-day to 10-day revisit) and Sentinel-1 (6-day or 12-day revisit), offer rich data to map different types of crops at a field scale in a timely fashion, especially in those regions with complicated crop landscapes and frequent cloud cover.

Conclusions
We evaluated the spatiotemporal consistency among NASS winter wheat statistics (planted area, harvested area, grain production), CDL-derived planted area, and GPP VPM for 2008-2018 at national and county scales. High spatiotemporal consistencies of planted area from the CDL and NASS datasets were found. There were large disagreements between the CDL planted area and NASS harvested area in many counties, especially in the SGP. The correlation between annual total GPP VPM and grain production varied with individual counties, depending upon the differences between CDL planted area and NASS harvested area. It was found that the cumulative GPP VPM at an 8-day interval and HI GPP were able to accurately predict winter wheat grain production at the county scale from early May to late June. HI GPP , calculated as the ratio of grain production over GPP, varied over a small range over different counties and years, when excluding those counties with large differences of CDL planted and NASS harvested areas, with a median value of 0.27. The results from this study highlight the importance of differentiating planted area and harvested area in crop mapping. Our results also demonstrate the potential of using GPP VPM products in estimating and monitoring grain production of winter wheat in the CONUS.
Author Contributions: Conceptualization, X.X. and X.W.; Methodology, X.W. and X.X.; data analysis. X.W. and X.X.; interpretation and discussion, X.W., X.X., J.S., Z.Y., Y.Q. and J.W.; writing, X.W. and X.X. led the early drafts, and all authors contributed to editing and revision of various versions of the manuscript; project administration, X.X., and funding acquisition, XX. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: The gross primary production (GPP) data for the CONUS from the VPM model, driven by MODIS imagery, CDL cropland area data, and NCEP re-analysis climate data, are available upon request to the corresponding author.