Estimation of SOS and EOS for Midwestern US Corn and Soybean Crops

Understanding crop phenology is fundamental to agricultural production, management, planning, and decision-making. This study used 250 m 16-day Moderate Resolution Imaging Spectroradiometer (MODIS) Enhanced Vegetation Index (EVI) time-series data to detect crop phenology across the Midwestern United States, 2007–2015. Key crop phenology metrics, start of season (SOS) and end of season (EOS), were estimated for corn and soybean. For such a large study region, we found that MODIS-estimated SOS and EOS values were highly dependent on the nature of input time-series data, analytical methods, and threshold values chosen for crop phenology detection. With the entire sequence of MODIS EVI time-series data as input, SOS values were inconsistent compared to crop emergent dates from the United States Department of Agriculture (USDA) Crop Progress Reports (CPR). However, when we removed winter EVI images from the time-series data to reduce impacts of snow cover, we obtained much more consistent SOS estimation. Various threshold values (10 to 50% of seasonal EVI amplitude) were applied to derive SOS values. For both corn’s and soybean’s SOS estimation, a threshold value of 25% generated the best overall agreement with the CPR crop emergent dates. Root-mean-square error (RMSE) values were 4.81 and 5.30 days for corn and soybean, respectively. For corn’s EOS estimation, a threshold value of 40% led to a high R2 value of 0.82 and RMSE value of 5.16 days. We further examined spatial patterns of SOS and EOS for both crops—SOS for corn displayed a clear south-north gradient; the southern portion of the Midwest US has earlier SOS and EOS dates.

homepage.do). State-level CPRs serve as an important data source for government agencies and commodity investors, as they track crop development during the growing season. CPR data have also been used to assess long-term trends over the past several decades [11][12][13]. For example, Kucharik [12] found that initiation of corn planting in 2005 was approximately two weeks earlier compared to the early 1980s for Corn Belt states.
Crop phenological parameters such as SOS and EOS can be estimated through remote sensing using methods similar to those developed for general vegetation phenology analysis. Time-series remote sensing data are now routinely used for characterizing vegetation phenology [14][15][16][17][18][19]. Using time-series Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) data, researchers are now developing operational global land surface phenological data for near real-time monitoring [19,20]. Using sequential remote sensing imagery to characterize crop phenology, however, presents a significant challenge due to the high spatial-temporal dynamics of agricultural landscapes [21]. Following Zhang et al.'s [19] phenological detection method, Wardlow et al. [21] derived green-up onset dates for corn, soybean, and sorghum using 16-day MODIS Normalized Difference Vegetation Index (NDVI) data. They compared MODIS-derived green-up onset dates with the USDA CPR data (e.g., 50% crop emerged) and found large inconsistencies across agricultural statistics districts (ASD). One of the main confounding factors was pre-crop or volunteer vegetation in the field. Sakamoto et al. [22] developed a Two-Step Filtering (TSF) method to detect phenological stages of maize and soybean using 6 years' time-series Wide Dynamic Range Vegetation Index (WDRVI) data derived from MODIS. Their phenological results were highly consistent with ground-based observations for two irrigated sites and one rainfed site in Nebraska. Their crop phenology detection method, however, has yet to be examined for large geographical areas and for longer-term (e.g., >8~10 years) monitoring.
At state and country scales, it is still unclear whether there is a robust relationship between remote sensing-derived SOS (EOS) dates and crop progress estimates provided from CPRs, especially for major crops such as corn and soybean. Although CPRs have their own uncertainties with respect to the nature of subjective observation and sample size, they are widely considered as the official statistics. Remote sensing phenological algorithms can be calibrated to generate very accurate site-specific crop phenology estimates. However, it is always a significant challenge to generalize such algorithms to a broad-scale, complex, agricultural landscape (e.g., major crop producing region such as the Midwestern US) and for relatively long time periods with varying weather conditions. Effective remote sensing analytical methods are needed to generate reliable estimates of crop phenology at state-country scales. Good agreement between remote sensing-derived crop phenology estimates and official CPRs would serve as a starting point for future operational use of remote sensing to potentially reduce the cost and effort in producing CPRs.
Accurate detection of key crop phenological phases such as SOS and EOS depend on remote sensing input data and selection of detection algorithms. A large number of phenological models or algorithms have been developed for vegetation phenology, and for land surface phenology in general. Most of such methods involve a two-step procedure of data smoothing and phenological parameter estimation. Because cloud and snow cover can significantly affect NDVI or Enhanced Vegetation Index (EVI) values, data smoothing is important to reduce signal noise caused by clouds and snow in time-series remote sensing data [2,16,23]. Phenological parameters are then estimated by using derivatives of the time-series data or applying user-defined thresholds (e.g., 50% of seasonal amplitude) [20,[24][25][26]. Among the various methods or tools, TIMESAT is one of the most commonly used packages because it provides an integrated framework for data smoothing and phenological parameter estimation [27]. Using TIMESAT, Tan et al. [20] suggested that pre-processing of time-series data, especially for pixels with snow cover, plays a very important role in estimating phenological parameters. In addition to snow/ice flags embedded in the MODIS surface reflectance products, Tan et al. [20] integrated MODIS land surface temperature to define winter season before vegetation phenology detection. Following Zhang et al. [19] and Ganguly et al. [28], the MODIS science team developed ready-to-use Land Cover Dynamics products (MCD12Q2, formerly Global Vegetation Phenology products). They gap-filled snow/ice pixels before characterizing vegetation phenology. None of the above-mentioned phenology products have been compared with CPRs for a broad-scale agricultural region, because they often focus upon general vegetation dynamics and comparison with CPRs would need crop-specific map products and additional data processing effort.
Aside from efforts to reduce snow/ice impacts upon the effectiveness of phenological monitoring, few researchers have attempted to integrate knowledge such as the general growing seasons of major crops to improve crop phenology detection. In the Midwestern US, farmers typically start corn planting in mid-April. The earliest planting dates recommended by crop insurance policies are also in early April. Therefore, a crop phenology algorithm can be designed to follow such local knowledge to reduce the complexity in dealing with the impacts of snow and ice.
Performance of crop phenology detection methods may also vary depending on the choice of time-series vegetation index products. Commonly used NDVI data are affected by soil background and surface soil wetness which may lead to uncertainties in detecting SOS. Under this circumstance, EVI is appealing because it can reduce the canopy background signal [29]. Finally, remote sensing-based phenology estimates can vary substantially because of uses of different user-defined thresholds and model parameters. It is always preferred if good crop phenology estimates can be obtained without extensive model calibration and threshold adjusting.
The overall goal of this study was to characterize key crop phenological parameters (SOS and EOS) for corn and soybean for the Midwest US using 250 m MODIS 16-day vegetation index products. SOS and EOS estimates were generated for a relatively long time period (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and compared to the official CPRs for a state-country scale assessment. Specific objectives were to: (1) examine whether a reduction in input time-series data, based on knowledge of general growing season (early-April to mid-November for corn and soybean) and snow/ice cover, can be used to improve SOS and EOS detection compared to the use of full time series; (2) examine how SOS and EOS estimates vary when TIMESAT threshold values (e.g., percent of seasonal amplitude) are varied; and (3) analyze spatial patterns of SOS and EOS based on MODIS-derived results.

Data and Data Preprocessing
Our study focused upon the Midwest US, defined here by the 12 states shown in Figure 1. This region is characterized by extensive crop cultivation. It is the major producer of US corn and soybean. In 2011, the study area produced 85% of the 10.4 billion bushels of corn binned and 81% of the 3 billion bushels of soybeans in the United States. It also produces half of the nation's wheat [30]. The Corn Belt located in this region is characterized by climates favoring crop production.
Remote Sens. 2017, 9,722 3 of 14 phenology. None of the above-mentioned phenology products have been compared with CPRs for a broad-scale agricultural region, because they often focus upon general vegetation dynamics and comparison with CPRs would need crop-specific map products and additional data processing effort. Aside from efforts to reduce snow/ice impacts upon the effectiveness of phenological monitoring, few researchers have attempted to integrate knowledge such as the general growing seasons of major crops to improve crop phenology detection. In the Midwestern US, farmers typically start corn planting in mid-April. The earliest planting dates recommended by crop insurance policies are also in early April. Therefore, a crop phenology algorithm can be designed to follow such local knowledge to reduce the complexity in dealing with the impacts of snow and ice.
Performance of crop phenology detection methods may also vary depending on the choice of time-series vegetation index products. Commonly used NDVI data are affected by soil background and surface soil wetness which may lead to uncertainties in detecting SOS. Under this circumstance, EVI is appealing because it can reduce the canopy background signal [29]. Finally, remote sensingbased phenology estimates can vary substantially because of uses of different user-defined thresholds and model parameters. It is always preferred if good crop phenology estimates can be obtained without extensive model calibration and threshold adjusting.
The overall goal of this study was to characterize key crop phenological parameters (SOS and EOS) for corn and soybean for the Midwest US using 250 m MODIS 16-day vegetation index products. SOS and EOS estimates were generated for a relatively long time period (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and compared to the official CPRs for a state-country scale assessment. Specific objectives were to: (1) examine whether a reduction in input time-series data, based on knowledge of general growing season (early-April to mid-November for corn and soybean) and snow/ice cover, can be used to improve SOS and EOS detection compared to the use of full time series; (2) examine how SOS and EOS estimates vary when TIMESAT threshold values (e.g., percent of seasonal amplitude) are varied; and (3) analyze spatial patterns of SOS and EOS based on MODIS-derived results.

Data and Data Preprocessing
Our study focused upon the Midwest US, defined here by the 12 states shown in Figure 1. This region is characterized by extensive crop cultivation. It is the major producer of US corn and soybean. In 2011, the study area produced 85% of the 10.4 billion bushels of corn binned and 81% of the 3 billion bushels of soybeans in the United States. It also produces half of the nation's wheat [30]. The Corn Belt located in this region is characterized by climates favoring crop production.  These high-quality vegetation index products include EVI, NDVI, surface reflectance values for red and NIR bands, and associated reliability index layers. Six MODIS vegetation index tiles were needed to cover the 12-state study area. For a study period of 2007-2015, a total of 1242 MODIS scenes were downloaded. For each composite interval, we used the MODIS Reprojection Tool (MRT) to extract, mosaic, and then project individual vegetation index and reliability index layers to an Albers Equal Area Conic projection. These image mosaics were then clipped to the 12-state study area boundary and stacked to build time-series EVI data. For each MODIS compositing image, it should be noted that EVI values may come from any day within the 16-day compositing window. For simplicity, we assumed that such values are obtained from day 8 (close to midpoint) of the compositing period.
Corn and soybean cropland masks were created with the Cropland Data Layer (CDL) produced by the USDA National Agricultural Statistics Service (NASS) [31]. Complete CDL coverages for all 12 states are available since 2007. Longer-term complete CDL data collected since 2001 are available for only three states-Illinois, Iowa, and North Dakota. NASS reports that CDLs have high levels of user's and producer's accuracies for major crops such as corn and soybeans, usually above 85% for Midwestern states. Thus, CDLs provide a reasonably reliable source of data for crop-specific phenological analysis. We use CDLs to extract pixels for corn and soybeans from 2007 through 2015 for the entire study area. Note that CDLs have different spatial resolutions of 30-56 m, depending on the year of production. For each year, we computed corn and soybean proportions within each 250 MODIS grid cell. Only MODIS pixels with corn or soybean proportions greater than 90% were considered for crop phenology detection and analysis.
CPRs provided by USDA NASS were used to evaluate the accuracy of MODIS-derived SOS and EOS estimates for corn and soybean. These progress reports of crop developmental stages are expressed as percentages of completed phases (e.g., 50% corn planted, emerged, and matured).

Estimate Crop Phenology Using Time-Series MODIS Data
The TIMESAT software package [26,27] was used to smooth MODIS EVI time-series data and detect crop phenology. Three data smoothing algorithms are available in TIMESAT: adaptive Savitzky-Golay (SG), asymmetric Gaussian, and double-logistic function. The adaptive Savitzky-Golay algorithm was used as our primary smoothing algorithm, since a recent study of smoothing algorithm comparisons suggested that the SG algorithm better characterizes temporal signals for both corn and soybean [23]. The SG algorithm basically applies a moving-window quadratic polynomial function to the original time-series data and estimates new values for the center point of each moving-window. TIMESAT allows users to provide a separate time-series input to identify pixels with cloud/shadow or other quality issues. Such pixels can be assigned with weights of zero, so they therefore will not affect data smoothing. We used the MODIS Reliability Index (RI) to identify such pixels (RI > 1). Because EVI signals with cloud contamination are often negatively biased, the TIMESAT package has an additional option to adapt fitted values to the upper envelope of the time-series data. Such adaptions were repeated twice to reduce potential impacts from cloud and shadow issues.
We examined two sets of time-series data for data smoothing and crop phenology detection. In the first set, full-year EVI time-series data were used as input to TIMESAT package. This approach served as a standard practice to benchmark phenology detection performance. For the second set, we simply defined a partial year EVI analytical window based on local knowledge of general crop growing season and information on crop insurance for the earliest crop planting dates. Farmers in the Midwestern US typically plant corn around mid-April and soybeans are planted slightly later than corn. The earliest planting dates recommended by crop insurance policies are also in early April. Most states complete harvest by mid-late November. Therefore, corn/soybean SOS and EOS detection can be focused within the analytical window of early-April to mid-November. This interval corresponds to 16-day MODIS time series composite periods 7-20 for each year. Use of this temporal analytical window reduced the percent of snow/ice impacted pixels (based on RI) to an average of less than 4% of total corn/soybean pixels for any given composite period, compared to an average of 7-62% of snow/ice impacted pixels for compositing periods defined from mid-November to late-March ( Figure 2). Remote Sens. 2017, 9,722 5 of 14 analytical window reduced the percent of snow/ice impacted pixels (based on RI) to an average of less than 4% of total corn/soybean pixels for any given composite period, compared to an average of 7-62% of snow/ice impacted pixels for compositing periods defined from mid-November to late-March ( Figure 2). On the basis of the smoothed EVI time series with full and partial image sequences as input, we first derived MODIS pixel-level (250 m resolution) SOS and EOS values using a default threshold (50% of the seasonal amplitude), measured from the left minimum and right minimum, respectively [26]. It should be noted that such threshold values were user-defined and can be easily adjusted to evaluate their impacts on crop phenology detection. Wu conducted a sensitivity analysis for SOS and EOS estimation based on various threshold values (0.1-0.5 with 0.05 incremental step).
With 9 years' EVI time-series data for a 12-state study area, TIMESAT data smoothing and phenology parameter estimation were computationally expensive. We divided our image processing tasks into 10 small segments and performed parallel processing at the Virginia Tech Advanced Research Computing's Newriver cluster. All data sub-setting and subsequent processing were streamlined using TIMESAT's setting files and Python scripts.

Comparison of MODIS-Derived Crop Phenology Parameters with CPRs
MODIS-derived phenology metrics of SOS and EOS were compared to crop progress statistics from CPRs. For each year and each state, 50% of corn or soybean emerged dates were extracted from CPRs and serve as references for SOS comparison. Similarly, 50% corn maturity dates were extracted from CPRs and serve as reference for EOS comparison. There was no data available for soybean maturity dates in CPRs, so the EOS comparison was conducted for corn only. It should be noted that CPRs are survey-based and only provide weekly snapshots of crop progress. In many cases, the exact dates for 50% crop emergence or crop maturity were not available. We applied linear interpolation to estimate 50% events using information from the two CPRs immediately straddling the crossing of the 50% threshold for each state and each crop type. Since CPRs are reported at state level, we averaged MODIS-derived SOS and EOS values for corn and soybean pixels in each state for each study year of 2007 to 2015. The root-mean-square error (RMSE) and coefficient of determination (R 2 ) were used to compare MODIS-derived crop phenology estimates and CPR values. On the basis of the smoothed EVI time series with full and partial image sequences as input, we first derived MODIS pixel-level (250 m resolution) SOS and EOS values using a default threshold (50% of the seasonal amplitude), measured from the left minimum and right minimum, respectively [26]. It should be noted that such threshold values were user-defined and can be easily adjusted to evaluate their impacts on crop phenology detection. Wu conducted a sensitivity analysis for SOS and EOS estimation based on various threshold values (0.1-0.5 with 0.05 incremental step).
With 9 years' EVI time-series data for a 12-state study area, TIMESAT data smoothing and phenology parameter estimation were computationally expensive. We divided our image processing tasks into 10 small segments and performed parallel processing at the Virginia Tech Advanced Research Computing's Newriver cluster. All data sub-setting and subsequent processing were streamlined using TIMESAT's setting files and Python scripts.

Comparison of MODIS-Derived Crop Phenology Parameters with CPRs
MODIS-derived phenology metrics of SOS and EOS were compared to crop progress statistics from CPRs. For each year and each state, 50% of corn or soybean emerged dates were extracted from CPRs and serve as references for SOS comparison. Similarly, 50% corn maturity dates were extracted from CPRs and serve as reference for EOS comparison. There was no data available for soybean maturity dates in CPRs, so the EOS comparison was conducted for corn only. It should be noted that CPRs are survey-based and only provide weekly snapshots of crop progress. In many cases, the exact dates for 50% crop emergence or crop maturity were not available. We applied linear interpolation to estimate 50% events using information from the two CPRs immediately straddling the crossing of the 50% threshold for each state and each crop type. Since CPRs are reported at state level, we averaged MODIS-derived SOS and EOS values for corn and soybean pixels in each state for each study year of 2007 to 2015. The root-mean-square error (RMSE) and coefficient of determination (R 2 ) were used to compare MODIS-derived crop phenology estimates and CPR values.

Analysis of Phenological Spatial Pattern
Our MODIS crop phenological parameter estimation was conducted at 250 m spatial resolution. Because of immense logistical, scale, and temporal obstacles, it will always be difficult to match field-based crop phenology data in order to obtain true validation, especially for such a large and complex agricultural region. Without such field-level validation, detailed analysis on spatial and temporal patterns of SOS and EOS estimates will always be uncertain. For this study, we only included a county-scale analysis based upon average SOS and EOS values for 9 years (2007-2015) for visual assessment of spatial patterns. Although validation data at the county-scale were not available, the spatial patterns of SOS and EOS can be loosely compared with broad-scale temperature and terrain variations to verify overall trends.

Analysis of Phenological Spatial Pattern
Our MODIS crop phenological parameter estimation was conducted at 250 m spatial resolution. Because of immense logistical, scale, and temporal obstacles, it will always be difficult to match fieldbased crop phenology data in order to obtain true validation, especially for such a large and complex agricultural region. Without such field-level validation, detailed analysis on spatial and temporal patterns of SOS and EOS estimates will always be uncertain. For this study, we only included a county-scale analysis based upon average SOS and EOS values for 9 years (2007-2015) for visual assessment of spatial patterns. Although validation data at the county-scale were not available, the spatial patterns of SOS and EOS can be loosely compared with broad-scale temperature and terrain variations to verify overall trends.

MODIS-Derived SOS and EOS Metrics Using a Default Threshold Value
Using a default 0.5 (50% of the EVI seasonal amplitude) threshold value, MODIS-derived corn and soybean SOS values were compared with 50 percent corn and soybean emerged dates from CPRs (2007-2015) ( Figure 3). With full year EVI time series as input, MODIS-derived SOS estimates for corn and soybean showed large inconsistencies compared to CPR data (Figure 3a,c). R 2 was 0.46 and 0.26 for corn and soybean, respectively. Corresponding RMSE values were 14.65 and 15.95 days. R 2 increased to 0.77 and 0.63 when partial EVI time series (early-April to mid-November) were used for TIMESAT data smoothing and SOS metric estimation (Figure 3b,d). However, RMSE values increased to 19.31 and 21.38 days for corn and soybean SOS estimation. Results appeared to be systematically shifted away from the 1:1 line.  With partial MODIS EVI time-series data as input, on average, the MODIS-predicted SOS value for corn was DOY (Day of Year) 162 in the Midwestern region, compared to DOY 143 for 50 percent corn emerged dates from CPR data. The average MODIS-predicted SOS value for soybean was DOY 177, about 21 days later than the CPR average value (DOY 156). Figure 4 compares MODIS EOS estimates with 50% corn maturity dates from CPRs, 2007-2015. With full-year EVI images as TIMESAT input, the R 2 value was 0.79 and the RMSE value was 5.69 days. By using partial year EVI images, the R 2 value increased to 0.84 and the corresponding RMSE value was 8.24 days. There were a couple outliers where late corn mature dates were reported in CPRs. For instance, 50% corn maturity was recorded on DOY 298 (October 25) in CPR, while MODIS predicted a much earlier EOS date at DOY 278 or 272 for full-year and partial-year input, respectively. There is no clear explanation for such a large discrepancy. In general, however, MODIS-predicted EOS values show good agreement with CPR data. With partial MODIS EVI time-series data as input, on average, the MODIS-predicted SOS value for corn was DOY (Day of Year) 162 in the Midwestern region, compared to DOY 143 for 50 percent corn emerged dates from CPR data. The average MODIS-predicted SOS value for soybean was DOY 177, about 21 days later than the CPR average value (DOY 156). Figure 4 compares MODIS EOS estimates with 50% corn maturity dates from CPRs, 2007-2015. With full-year EVI images as TIMESAT input, the R 2 value was 0.79 and the RMSE value was 5.69 days. By using partial year EVI images, the R 2 value increased to 0.84 and the corresponding RMSE value was 8.24 days. There were a couple outliers where late corn mature dates were reported in CPRs. For instance, 50% corn maturity was recorded on DOY 298 (October 25) in CPR, while MODIS predicted a much earlier EOS date at DOY 278 or 272 for full-year and partial-year input, respectively. There is no clear explanation for such a large discrepancy. In general, however, MODIS-predicted EOS values show good agreement with CPR data.

MODIS-Derived SOS Metrics with Various Threshold Values
Using full year images as input, for Corn's SOS estimation with various threshold values from 0.1 to 0.5, RMSE values were in the range of 8.37-51.87 days ( Figure 5). A threshold value of 0.35 resulted in the smallest RMSE value of 8.37 days. However, the corresponding R 2 value was 0.23, which suggested weak-moderate linear relationships between MODIS-derived SOS dates and CPR dates. Using partial year images as input, RMSE values were in the range of 4.56-19.31 days. The best results were obtained when a threshold value of 0.25 was applied for SOS detection. At this point, the RMSE value (4.81 days) was close to the smallest RMSE (4.56 days, threshold = 0.2) and the R 2 value (0.72) was slightly higher than the value (0.70) obtained using a threshold value of 0.2. We note that the R 2 value continued to increase from 0.72 to 0.77 when threshold values were increased from 0.25 to 0.5. However, MODIS-derived SOS values were pushed further away from the 1:1 line and RMSE values increased from 4.81 days to 19.31 days.
For soybean's SOS detection with partial year image as input, a threshold value of 0.25 also resulted in the best overall agreement with CPR data. Figure 6 shows the scatter plots comparing the results from the MODIS-derived SOS values and CPR data using a 0.25 threshold value. Overall, SOS detection for both crops had reasonable agreement with CPR data. RMSE values were less than one week and R 2 values were higher than 0.55. Soybean had lower consistency compared to corn when CPRs were used as references.

MODIS-Derived SOS Metrics with Various Threshold Values
Using full year images as input, for Corn's SOS estimation with various threshold values from 0.1 to 0.5, RMSE values were in the range of 8.37-51.87 days ( Figure 5). A threshold value of 0.35 resulted in the smallest RMSE value of 8.37 days. However, the corresponding R 2 value was 0.23, which suggested weak-moderate linear relationships between MODIS-derived SOS dates and CPR dates. Using partial year images as input, RMSE values were in the range of 4.56-19.31 days. The best results were obtained when a threshold value of 0.25 was applied for SOS detection. At this point, the RMSE value (4.81 days) was close to the smallest RMSE (4.56 days, threshold = 0.2) and the R 2 value (0.72) was slightly higher than the value (0.70) obtained using a threshold value of 0.2. We note that the R 2 value continued to increase from 0.72 to 0.77 when threshold values were increased from 0.25 to 0.5. However, MODIS-derived SOS values were pushed further away from the 1:1 line and RMSE values increased from 4.81 days to 19.31 days.
For soybean's SOS detection with partial year image as input, a threshold value of 0.25 also resulted in the best overall agreement with CPR data. Figure 6 shows the scatter plots comparing the results from the MODIS-derived SOS values and CPR data using a 0.25 threshold value. Overall, SOS detection for both crops had reasonable agreement with CPR data. RMSE values were less than one week and R 2 values were higher than 0.55. Soybean had lower consistency compared to corn when CPRs were used as references.  We further compared the MODIS-derived SOS values and CPR data for individual states. At the individual state level, the sample size was reduced to a small number of 9 (i.e., [2007][2008][2009][2010][2011][2012][2013][2014][2015]. For each state, the threshold value of 0.25 also resulted in the best (or close to optimal) agreement with the CPR data. Table 1 shows RMSE and R 2 values for each state using a single threshold value of 0.25. At the individual state level, RMSE values were in the ranges of 2.84-6.83 and 3.11-7.12 days for corn's and soybean's SOS estimation, respectively. We note that both RMSE and R2 statistics could be noisy at the individual state level, because the sample size (n = 9) is small.  We further compared the MODIS-derived SOS values and CPR data for individual states. At the individual state level, the sample size was reduced to a small number of 9 (i.e., [2007][2008][2009][2010][2011][2012][2013][2014][2015]. For each state, the threshold value of 0.25 also resulted in the best (or close to optimal) agreement with the CPR data. Table 1 shows RMSE and R 2 values for each state using a single threshold value of 0.25. At the individual state level, RMSE values were in the ranges of 2.84-6.83 and 3.11-7.12 days for corn's and soybean's SOS estimation, respectively. We note that both RMSE and R2 statistics could be noisy at the individual state level, because the sample size (n = 9) is small. We further compared the MODIS-derived SOS values and CPR data for individual states. At the individual state level, the sample size was reduced to a small number of 9 (i.e., [2007][2008][2009][2010][2011][2012][2013][2014][2015]. For each state, the threshold value of 0.25 also resulted in the best (or close to optimal) agreement with the CPR data. Table 1 shows RMSE and R 2 values for each state using a single threshold value of 0.25. At the individual state level, RMSE values were in the ranges of 2.84-6.83 and 3.11-7.12 days for corn's and soybean's SOS estimation, respectively. We note that both RMSE and R2 statistics could be noisy at the individual state level, because the sample size (n = 9) is small.

MODIS-Derived EOS Metrics with Various Threshold Values
Similar to the SOS detection, the choice of threshold values clearly had impacts on MODIS EOS estimates (Figure 7). For corn's EOS estimation, a threshold value of 0.4 combined with partial year images as input led to a very high R 2 value (0.82) and the smallest RMSE value (5.16 days). At the individual state level, the same threshold value of 0.4 led to the smallest RMSE values for eight out of 12 states. For the remaining four states, the RMSE values were close to optimal, with less than 2 days difference compared to the smallest RMSE values. MODIS-derived EOS values for soybean were not compared to the CPRs, because there was no percent soybean maturity data in the CPR database.

MODIS-Derived EOS Metrics with Various Threshold Values
Similar to the SOS detection, the choice of threshold values clearly had impacts on MODIS EOS estimates ( Figure 7). For corn's EOS estimation, a threshold value of 0.4 combined with partial year images as input led to a very high R 2 value (0.82) and the smallest RMSE value (5.16 days). At the individual state level, the same threshold value of 0.4 led to the smallest RMSE values for eight out of 12 states. For the remaining four states, the RMSE values were close to optimal, with less than 2 days difference compared to the smallest RMSE values. MODIS-derived EOS values for soybean were not compared to the CPRs, because there was no percent soybean maturity data in the CPR database.

Spatial Patterns of SOS and EOS
After evaluating all combinations of MODIS input data and threshold value choices, it was determined that a combination of partial year MODIS time-series and a threshold value of 0.25 generated the most consistent SOS estimates for both corn and soybean. A threshold value of 0.4 generated the most consistent EOS estimates for corn. The annual SOS and EOS maps derived via this approach were used to analyze spatial patterns. For ease of illustration, we focused on countyscale analysis and used average SOS and EOS values for 9 years (2007-2015) for visual assessment

Spatial Patterns of SOS and EOS
After evaluating all combinations of MODIS input data and threshold value choices, it was determined that a combination of partial year MODIS time-series and a threshold value of 0.25 generated the most consistent SOS estimates for both corn and soybean. A threshold value of 0.4 generated the most consistent EOS estimates for corn. The annual SOS and EOS maps derived via this approach were used to analyze spatial patterns. For ease of illustration, we focused on county-scale analysis and used average SOS and EOS values for 9 years (2007-2015) for visual assessment (Figure 8). For corn, SOS and EOS values were observed to increase following the south-north gradient (Figure 8a,c). Earlier SOS in the southern portion of the Midwest US was expected because of the relatively warmer winter and longer growing season. SOS values for soybean vary substantially (i.e., late April-late June) across the study region (Figure 8e). There was no clear south-north gradient and counties located at the edges of the Corn Belt appeared to have late SOSs. For both corn and soybean, the standard deviation values of SOS (and EOS) were higher in the southern portion of the Midwest US, especially for Illinois, Indiana, and Missouri (Figure 8b,d,f). Such high inter-annual variability can be attributed to favorable climates in these states, which provide farmers in this sub-region with more flexibility in managing crop planting and harvesting dates.
Remote Sens. 2017, 9,722 10 of 14 ( Figure 8). For corn, SOS and EOS values were observed to increase following the south-north gradient (Figure 8a,c). Earlier SOS in the southern portion of the Midwest US was expected because of the relatively warmer winter and longer growing season. SOS values for soybean vary substantially (i.e., late April-late June) across the study region (Figure 8e). There was no clear southnorth gradient and counties located at the edges of the Corn Belt appeared to have late SOSs. For both corn and soybean, the standard deviation values of SOS (and EOS) were higher in the southern portion of the Midwest US, especially for Illinois, Indiana, and Missouri (Figure 8b,d,f). Such high inter-annual variability can be attributed to favorable climates in these states, which provide farmers in this sub-region with more flexibility in managing crop planting and harvesting dates.

Discussion
Consistent determination of SOS and EOS for corn and soybean depends on many factors such as input time-series data, the smoothing algorithms chosen, and specific analytical approaches and parameters used to pinpoint phenological metrics. Time-series EVI can be noisy due to clouds,

Discussion
Consistent determination of SOS and EOS for corn and soybean depends on many factors such as input time-series data, the smoothing algorithms chosen, and specific analytical approaches and parameters used to pinpoint phenological metrics. Time-series EVI can be noisy due to clouds, shadows, and variations in viewing and illumination geometry. Using MODIS RI as ancillary quality data may significantly improve SOS (and EOS) estimation. However, we suspected that snow pixels and pixels with partial snow cover still affect seasonal EVI amplitude estimates, as well as SOS and EOS estimates defined by thresholding of seasonal amplitude.
In contrast to previous studies in which snow pixels were gap-filled and winter seasons were specifically defined, pixel-by-pixel, our approach was to simply remove all winter images from mid-November to late-March to reduce snow/ice impacts. This simplified method was designed for major summer crops such as corn and soybean, because their growing seasons are well-defined (e.g., early-April to mid-November) and crop insurance also provides guides for identifying the earliest planting dates. Local knowledge of crop calendar and insurance policies thus serve as an alternate approach in determining valid data processing windows. The partial year time-series input data provided much improved performance in detecting SOS for both crops with respect to demonstrating improved statistical relationships with CPR information (Figures 3 and 5). We emphasize that snow pixels and winter images are particularly challenging for deriving vegetation dynamics and phenology [14]. For example, for a high-latitude environment where the growing season is short and snow cover is persistent, Beck et al. [14] developed a new method to estimate winter NDVI and then applied the double logistic function to model vegetation phenology. Such a method outperformed several existing algorithms using Fourier series and asymmetric Gaussian functions. Other advantages of their method also include automated growing season detection and effectiveness in dealing with outliers. For our future study, such a method needs to be explored and compared to our approach.
Threshold values used to detect SOS and EOS were also important. For a large study area such as the Midwest US, we were searching for a threshold applicable to the entire study area, rather than variable thresholds across the region. Compared to a 50% seasonal amplitude as the default threshold, a 25% value appeared to be better associated with CPR 50% crop emerged dates for both crops. The overall less consistent SOS estimates for soybean might be attributed to its relatively late planting dates compared to corn. Late planting dates potentially lead to higher variations in pre-crop vegetation conditions, which is a main confounding factor for SOS estimation [21]. Other potential confounding factors include surface soil wetness and local weather conditions. These factors may act interactively to affect vegetation/crop condition and increase the uncertainty of SOS estimation.
We also examined the ready-to-use MODIS Land Cover Dynamics product (MCD12Q2, 500 m resolution) to help benchmark performance of our crop phenology estimates. This MCD12Q2 product employed feature point or transition date detection algorithms developed by Zhang et al. (2003Zhang et al. ( , 2006. The main inputs for this MODIS product include an 8-day EVI derived from Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR) and MODIS land surface temperature data. We compared the onset of greenness increase dates derived from the MCD12Q2 and CPR 50% crop emerged dates for both crops. The R 2 value was 0.43 and 0.19 for corn and soybean, respectively. Corresponding RMSE values were 7.59 and 13.89 days. MCD12Q2 tends to show earlier SOS estimates compared to CPRs. Such RMSE and R 2 statistics suggested that MCD12Q2 products may not be ideal dataset for crop phenology monitoring, because they are designed for monitoring vegetation dynamics at a global scale.
We used 16-day MODIS time-series data to derive SOS and EOS estimates. MODIS EVI values were assumed to be obtained from the eighth day (or midpoint) of the compositing period for simplicity. The actual Date-of-Acquisition (DOA) information is available to further improve the temporal precision. For example, Brown et al. [32] found that integration of DOA and temporal interpolation of the 16-day composite vegetation index improved their multi-temporal crop classification results. The only concern with including DOA information is that the pixel-by-pixel temporal interpolation can be time-consuming, especially for a large study area and for long-term analysis. We also conducted additional analysis to examine how the use of different vegetation indices (EVI and NDVI) affects SOS detection. For each pair of EVI-derived and NDVI-derived SOS values, we conducted a correlation analysis and found high correlation for corn (r~0.96) and soybean (r~0.92). For both EVI and NDVI products, the same threshold value of 0.25 resulted in the best overall agreement with the CRP data. Furthermore, there was no significant difference of accuracy statistics (RMSE and R 2 ) using two input datasets.
Spatial patterns of SOS (EOS) estimates were analyzed at the county scale. Such county-level crop phenological metrics may augment the USDA CPRs by providing improved spatial details. We note that county-level metrics have not been validated due to the unavailability and doubtful quality of reference data. Spatial pattern analysis of SOS and EOS at finer resolutions (e.g., 250 m) is possible but not encouraged, because we did not conduct pixel-scale validation.

Conclusions
We used 250 m MODIS EVI time-series data to estimate annual SOS and EOS for corn and soybean within the 12-state Midwest US, 2007-2015. MODIS-derived SOS and EOS values were compared with the USDA CPR 50% crop emerged dates and 50% crop mature dates, 2007-2015. Inconsistent SOS and EOS values were derived when full year EVI images were used as input to TIMESAT phenological parameter estimation. Using corn/soybean's general crop calendar or growing season (early-April to mid-November) as a guide, we removed winter images from MODIS time-series data to reduce the impacts from snow and ice. Use of partial year time-series data substantially improved consistency between MODIS-derived SOS (and EOS) dates and CPR data. We also examined various threshold values (10 to 50% of seasonal EVI amplitude) for SOS and EOS detection. For both corn and soybean, a threshold value of 25% generated much improved SOS estimates compared to the default 50% for the Midwest US, at both entire study region and individual state scales. For corn's EOS estimation, a threshold value of 40% led to a very high R 2 value of 0.82 and RMSE value of 5.16 days. Spatial analyses of SOS and EOS values revealed a clear south-north gradient for corn-earlier SOS (and EOS) in the south and later SOS (and EOS) in the north portion of the study region.