Using MODIS Data to Predict Regional Corn Yields

Abstract

A simple approach was developed to predict corn yields using the MoDerate Resolution Imaging Spectroradiometer (MODIS) data product from two geographically separate major corn crop production regions: Illinois, USA and Heilongjiang, China. The MOD09A1 data, which are eight-day interval surface reflectance data, were obtained from day of the year (DOY) 89 to 337 to calculate the leaf area index (LAI). The sum of the LAI from early in the season to a given date in the season (end of DOY (EOD)) was well fitted to a logistic function and represented seasonal changes in leaf area duration (LAD). A simple phenology model was derived to estimate the dates of emergence and maturity using the LAD-logistic function parameters b₁ and b₂, which represented the rate of increase in LAI and the date of maximum LAI, respectively. The phenology model predicted emergence and maturity dates fairly well, with root mean square error (RMSE) values of 6.3 and 4.9 days for the validation dataset, respectively. Two simple linear regression models (Y_P and Y_F) were established using LAD as the variable to predict corn yield. The yield model Y_P used LAD from predicted emergence to maturity, and the yield model Y_F used LAD for a predetermined period from DOY 89 to a particular EOD. When state/province corn yields for the validation dataset were predicted at DOY 321, near completion of the corn harvest, the Y_P model, including the predicted phenology, performed much better than the Y_F model, with RMSE values of 0.68 t/ha and 0.66 t/ha for Illinois and Heilongjiang, respectively. The Y_P model showed similar or better performance, even for the much earlier pre-harvest yield prediction at DOY 257. In addition, the model performance showed no difference between the two study regions with very different climates and cultivation methods, including cultivar and irrigation management. These results suggested that the approach described in this paper has potential for application to relatively wide agroclimatic regions with different cultivation methods and for extension to the other crops. However, it needs to be examined further in tropical and sub-tropical regions, which are very different from the two study regions with respect to agroclimatic constraints and agrotechnologies.

1. Introduction

Global warming is projected to accompany more frequent extreme weather events, such as heavy rainfall, high temperature, and drought [1,2]. Although these changes will positively affect crop production in some regions, crop production for food, feed, and fodder will be negatively affected in other regions [3,4], aggravating all dimensions of food security. Timely and reliable information on crop growth and yield at the regional, national and global scales is essential for food security and trade policies [5,6,7].

Notable advances in remote sensing have enabled reliable and timely prediction of crop yields [8,9]. Vegetation indices (VIs), which are calculated by combining the reflectance values of several spectral bands from multispectral satellite systems, have been used directly as crop yield estimators [8,10,11] and/or to estimate intermediate crop growth variables, such as LAI and biomass, for yield prediction [12,13,14,15]. In addition to VIs and crop growth variables, crop phenology information is required for reliable crop yield prediction because the effects of environmental conditions on crop yield differ by growth stage [16]. Crop phenology can be estimated using remote sensing data products [17,18,19,20]. For example, Islam and Bala [17] used NDVI and LAI derived from remote sensing data to identify the planting and ending dates of potato, Jönsson and Eklundh [18] developed TIMESAT for extracting phenological information from remote sensing time series estimates, Kandasamy et al. [19] compared the different methods for the performance to fill gaps and how this affects phenology estimations, and Kandasamy et al. [20] provides an approach to evaluate phenology algorithms and the uncertainty associated with its estimates.

Predicting crop yield using remote sensing data products often depends on an empirical approach that relates VIs alone or in combination with remote-sensing–derived meteorological variables to the reported crop yields [11,21,22,23]. Bolton and Friedl [21] used a simple linear regression model to predict corn and soybean yield using VIs (i.e., NDVI, EVI2, and NDWI) from MODIS in the Central USA and concluded that the EVI2 index exhibited better ability to predict maize yield than the NDVI and that the use of crop phenology information from MODIS improved the model predictability. Johnson [22] developed a regression tree model using the linear and/or exponential relationship of MODIS–derived NDVI and daytime land surface temperature with county-level yield statistics to predict corn and soybean yield for 12 states in central and northern USA. Shao et al. [23] developed simple linear regression model using multi-temporal NDVI from MODIS to predict county–level corn yields for the entire Midwestern USA and confirmed that using a digital elevation model climate data as additional model inputs slightly improved the performance of the regional corn yield model. Although simple models to predict crop yield can be developed [9,24], it is likely that the spatial portability of the model would be limited because the parameters of the empirical equation that were estimated in the study region would not be applicable to the other regions with different agroclimate and agrotechnologies. Additionally, VIs such as NDVI can detect inter-annual fluctuation of crop yield due to weather conditions while they cannot detect human-induced factors that result in increased crop yield [9], making the NDVI–crop yield regression model difficult to extend to other regions [10,25]. Sophisticated approaches requiring data in addition to remote sensing data have also been developed to improve predictions of crop growth and yield [9,26,27]. For instance, Prasad et al. [27] used a piecewise linear regression model with a break point to predict corn and soybean yield using monthly NDVI, soil moisture, surface temperature, and total rainfall. Nevertheless, it is preferable to develop a crop yield prediction model with wide spatial portability using a small dataset [28].

In addition to VIs, remote-sensing–derived LAI has also been used to predict crop biomass and yield [29,30,31]. LAI, the one side total leaf area per unit of ground area [32], is a key biophysical variable used to determine crop growth and yield [33,34]. Son et al. [31] developed a regression model using MODIS–LAI and EVI as input variables to predict rice yield in the Mekong delta, Vietnam. LAD is the integrated value of LAI over time. Although LAD is an important crop parameter that has a strong correlation with dry matter production and grain yield [35,36], it has not been used in a yield prediction model using remote sensing data. LAD has been reported to have positive correlation with corn yield under water and nitrogen stress conditions imposed at different growth stages [37] and under varying planting densities of three maize hybrids [38], and genetic differences in photosynthetic duration (longer LAD) were reported to be associated with a longer grain filling duration and higher yield [39]. These findings suggest that LAD would have greater potential to represent corn yield variability in regions with diverse agroclimate and agrotechnologies compared to VIs and LAI at a particular crop growth stage.

The objectives of this study were to develop a simple approach to predict crop phenology and corn yield using LAI and LAD derived from a minimum set of remote sensing data products and to examine the spatial portability of the simple method. This study focused on predicting corn yield in major production areas of the USA and China with different agroclimates and agrotechnologies.

2. Materials and Methods

2.1. Study Area

A simple approach using remote sensing data products was developed to predict crop yields in major crop production areas in the USA and China, where varieties and cultivation methods differ considerably. Illinois, USA (40°0′0″N, 89°0′0″W) (Figure 1a) and Heilongjiang Province, China (47°50′0″N, 127°40′0″E) (Figure 2a) were selected as the regions of interest. Illinois and Heilongjiang are located on opposite sides of the Earth. The latitude of Heilongjiang is higher than that of Illinois. Accordingly, the environmental conditions (e.g., temperature, precipitation, and solar radiation) differ between sites. The annual mean temperature in Illinois is approximately 10 °C, while in Heilongjiang it is −4 °C to 4 °C. These different environmental conditions explain why different crop varieties and cultivation methods are used at the two sites [40].

The corn production area and the quantity of corn produced as a proportion of the national total in Illinois were approximately 32% and 15% in 2013, respectively. The corresponding values in Heilongjiang were approximately 9% and 12%.

2.2. Data and Processing

2.2.1. Crop Yield and Phenology Data

Corn yields in Illinois and Heilongjiang (Figure 3) were obtained from official reports provided by statistical services in each country to evaluate the reliability of a model for predicting corn yield. In Illinois, crop yields from 2000 to 2013 were gathered from the National Agricultural Statistics Service (NASS) by agricultural district (AD) (available at https://quickstats.nass.usda.gov/). Crop area and total production were also obtained from NASS by county and state. The unit system for crop yield in Illinois was converted from “bushels per acre” to “kilograms per hectare” to fit the unit system used in Heilongjiang.

The yields and planted area for 2002 to 2012 in Heilongjiang Province were collected from the Heilongjiang Statistical Yearbook by prefecture and province. The corn yields in Heihe and Daxinganling prefectures in Heilongjiang Province were excluded from the calculation to minimize the computing resources because mean corn production was considerably lower than in other areas of the province.

The data of corn emergence and maturity date for Illinois were obtained from weekly crop progress reports by the NASS Illinois Field Office (IFO). These data, which were provided by AD, were only available for five ADs (northwest, northeast, central, western, and eastern districts) between 2003 and 2012 due to missing data in some ADs. The dates of emergence and maturity were obtained. Because of the limited availability of phenology data in Heilongjiang Province, a model was established and validated to predict the phenological stages using data only for the ADs in Illinois, USA and was applied to predict the phenological stages in Heilongjiang, China.

The phenology data were used to evaluate the reliability of a model to predict corn phenological dates using remote sensing data. The date on which a given phenological stage reached 50% in the study area, e.g., an AD, were used to represent when a phenological stage occurred. It was assumed that the proportion of area in a given stage would increase linearly over a two-week period. The weekly data for the proportion of area in which the crop was at the phenological stage of interest were collected from the NASS IFO. Those weekly data were used to compare the estimated dates on which a phenological stage occurred in the AD. For example, the day on which the corn crop reached a phenological stage was determined by linear interpolation between the weekly proportion of area for the given phenological stage.

2.2.2. Crop Cover Data

Crop cover data for Illinois (Figure 1b) were obtained from cropland data layers (CDL) provided by NASS to identify the area where a given crop was grown (available at https://www.nass.usda.gov/Research_and_Science/Cropland/SARS1a.php). Crop cover data were obtained from 2000 to 2013. The crop cover data from 2002 to 2012 in China (Figure 2b) were generated using remote sensing [30]. Crop areas were identified from the MODIS land cover type product (MCD12Q1). Major crops, including corn, within the crop areas were classified using a maximum likelihood classifier and time-series MODIS 16-day NDVI dataset (MOD13Q1). Crop cover data, which have a spatial resolution of 250 m, were subjected to reprojection using ENVI (ExelisVIS: Exelis Visual Information Solutions, Boulder, CO, USA). The projections of the crop cover maps in both regions were converted to a Universal Transverse Mercator (UTM) projection and WGS-84 coordinates at 1 km spatial resolution.

2.2.3. Remote Sensing Data

The MODIS surface reflectance data for 2000–2013 were obtained from the MOD09A1 data, which have a spatial resolution of 500 m. These data were downloaded from Reverb, which is a web-based remote sensing data portal operated by the National Aeronautics and Space Administration (available at http://modis.gsfc.nasa.gov/). The h10v04, h10v05, h11v04, and h11v05 tiled grid data for Illinois were collected from DOY 89 to 329. The h26v04 and h27v04 tiled grid data were obtained for the same period in Heilongjiang.

The near-infrared (NIR; band 2) and red (band 1) band of reflectance data were prepared to estimate LAI using a series of data tools (Figure 4). First, all tiled reflectance data were mosaicked into a single dataset using interface description language (IDL; ExelisVIS). The projection of mosaic data was converted to UTM projection at 1 km spatial resolution and WGS-84 geographic latitude and longitude coordinates using IDL, which applies the triangulation wrap and nearest neighbor resampling methods. Reflectance data were resized to fit the size and georeference of crop cover data using FWTools, which is a collection of open-source GIS applications. Gridded data for the area where a given crop was grown by year were prepared by extracting the reflectance data for the given extent to compare grids belonging to the crop of interest in the crop cover data using MATLAB (MathWorks Inc., Natick, MA, USA).

Cropland was identified using the proportion of crop cover in each pixel of the MODIS product using a mixed problem and low resolution [21,23]. Crop cover data were overlapped with the MODIS data to calculate the percentage of pixels where a given crop was grown in the research region. It was assumed that the crop of interest would be grown within pixels where the percentage of corn was >60% and >90%. Although identifying a crop using the percentage of crop cover could be used for pure growing regions of the crop, a percentage value was not set for the identification method, so the mixed problem differed by region and crop. This method can cause problems when integrating from a smaller to a larger scale due to the pixels below the threshold percentage being excluded. In this study, the pixels for corn growing in Illinois and Heilongjiang were identified according to the data standardization procedure described in Figure 4 and the number of identified pixels in Illinois and Heilongjiang were approximately 610,000 and 860,000, respectively.

2.3. Estimation of LAI

2.3.1. Estimation of LAI Using Remote Sensing Data

LAI was calculated using the NIR and red bands of the reflectance data. Nguy-Robertson et al. [41] suggested that the combined vegetation index (CVI) is more accurate than a single vegetation index set to estimate LAI. Instead of using the LAI products from the MODIS data, the LAI value was calculated as follows [31]:

$$\mathrm{LAI}=\{\begin{array}{c}\left(\mathrm{NDVI}-0.28\right)/0.28, \mathrm{NDVI} \le 0.7\\ \left(\mathrm{SR}-1.0\right)/3.5, \mathrm{NDVI} >0.7 \end{array},$$

(1)

where NDVI and SR are the normalized difference vegetation index [42] and simple ratio [43], respectively. NDVI and SR were determined as follows:

$$\mathrm{NDVI}=\left(\mathrm{NIR}-\mathrm{red}\right)/\left(\mathrm{NIR}+\mathrm{red}\right)$$

(2)

$$\mathrm{SR}=\mathrm{NIR}/\mathrm{red},$$

(3)

where NIR and red indicate the near-infrared and red band spectra, respectively.

2.3.2. Estimation of Daily LAI Using a Logistic Function

Various smoothing algorithms have been used to reduce the noise in remote sensing time-series data [44]. The Savitzky–Golay, asymmetric Gaussian, double-logistic, Whittaker smoother, and discrete Fourier transformation smoothing algorithms have been applied to the NDVI data in the MODIS product. In this study, MODIS–derived LAI was smoothed using a simple method. It was assumed that the total sum of the LAI over time would be shaped like a logistic function. For example, the LAI would be negligible early in the season until the emergence date. However, the daily LAI sum would increase rapidly during the vegetative growth stage, but the sum of the LAI would become negative after flowering. Finally, the sum of the LAI over a season would remain relatively constant until harvest after physiological maturity (Figure 5).

A logistic function was used to represent the temporal change in leaf area duration (LAD), which is the integral of LAI over time as follows (Equation (4)):

$$\mathrm{LAD}\left(\mathrm{t}\right)={\displaystyle \int }LAI\left(t\right) dt=\frac{b_3}{1.0+\exp \left[-b_1\left(t-b_2\right)\right]},$$

(4)

where b₁ represents the rate of LAI growth, b₂ represents the date of the maximum LAI, b₃ represents the cumulative LAI at physiological maturity, and t indicates days after planting.

Although it would be preferable to determine t as the date after the exact planting date, it was challenging to obtain this date in each grid cell. Because a logistic function was used in Equation (4), it was assumed that t could be determined using a date earlier than the actual planting date. In this study, DOY 89 was used as the beginning of the cropping season, which was a date prior to the actual planting for the regions of interest. The order of the remote sensing data products since DOY 89 was used to determine t and conveniently estimate the logistic function parameters. For example, t was 2 when data production on DOY 97 was used for Equation (4). The simplex method [45] was used to determine the parameter values of Equation (4) for each grid cell. The sum of the LAI derived from the MODIS product until a given date was compared with that obtained from Equation (4). The LAI values from 89 to 337 DOY were accumulated at eight-day intervals for each grid cell at a 1-km resolution. The sum of the square error between the observed and simulated values of cumulative LAI was minimized to obtain parameter values for b₁, b₂, and b₃ in the simplex algorithm.

The derivative of Equation (4) on a given date represented the LAI on that date:

$$\mathrm{LAI}=\frac{b_3·b_1·\exp \left[-b_1\left(d - b_2\right)\right]}{{\left\{1.0 + \exp \left[-b_1\left(d - b_2\right)\right]\right\}}^2}.$$

(5)

Daily changes in LAI (Figure 6) were determined for each grid cell.

2.4. Prediction of Phenological Dates

The dates of emergence and maturity were estimated using a simple empirical equation (Equation (6)). It was assumed that date D_P for a given phenology P could be determined using b₁ and b₂ from the LAD-logistic function (Equation (5)), which represent the rate of LAI increase and the date of maximum LAI at a given site, respectively, as follows:

$${\mathrm{D}}_{\mathrm{P}}={\mathrm{b}}_2+{\mathsf{\tau }}_{\mathrm{P}}+{\mathsf{\rho }}_{\mathrm{P}}/{\mathrm{b}}_1,$$

(6)

where ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ are empirical coefficients, ${\mathsf{\tau }}_{\mathrm{P}}$ represents the overall time difference between the date of the maximum LAI and a given phenological stage P, and ${\mathsf{\rho }}_{\mathrm{P}}$ represents the impact of the increase in LAI on the phenological change over time. ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ were estimated using the simplex method. The NASS–derived dates for a given phenology were compared with those obtained from Equation (5).

2.5. Prediction of Crop Yield

Daily LAD was also determined from the date of emergence to that of maturity for each grid cell, as follows [46]:

$$\mathrm{LAD}=\frac{\left(LA{I}_{d+1 }+ LA{I}_d\right)}{2},$$

(7)

where LAI_d is the LAI value on d. The sum of the LAD values was obtained for the growing periods, e.g., from emergence to maturity, for each grid cell. Then, the sum of those values was averaged by region, e.g., AD and prefecture in a season as follows:

$$\mathrm{ALAD}=\frac{1}{n}{{\displaystyle \sum }}_{\mathrm{c}=0}^{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{d}=0}^{\mathrm{m}}LA{D}_{dc},$$

(8)

where c and d represent the grid cell index in the region in which the corn was grown and the date index from emergence to maturity at c, respectively.

Huang et al. [24] reported that crop yield can be determined with the LAI using simple linear regression. In this study, ALAD was used as an independent variable in the simple linear regression to predict grain yield. The coefficients of the regression line were obtained for the reported yields and yield predictions for a given district D in a season as follows:

$${\mathrm{Yield}}_{\mathrm{D}}={\mathsf{\alpha }\mathrm{ALAD}}_{\mathrm{D}}+\mathsf{\beta },$$

(9)

where α and β are coefficients estimated by the least squares difference method.

2.5.1. Y_P Model Using LAD Accumulated from the Estimated Emergence Date

The Y_P model used the predicted phenological stages and smoothed LAI values to calculate LAD. LAI values were accumulated from the emergence to maturity dates predicted using Equations (4) and (6). The last date on which the remote sensing data products were used was denoted EOD. The EOD value was the date elapsed from the initial date of analysis when the approach was used in other regions, e.g., in the southern hemisphere. However, DOY was used to indicate the EOD for simplicity. The coefficients of the simple linear regression for the Y_P model were obtained for the relation between the reported yields and yield estimates. No remote sensing data was available from EOD to maturity if the EOD was earlier than the maturity date. In such cases, Equation (5) was used to estimate daily LAI until maturity.

2.5.2. Y_F Model Using LAD Accumulated from an Arbitrarily Fixed Date

To estimate the dates of emergence and maturity, phenology data in the region of interest must be available to determine ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ for Equation (6). Therefore, the Yp model would not be applicable to a region where few phenology data are available. Instead of using the smoothed LAD values from the predicted emergence to maturity date, the Y_F model used the LAD accumulated with LAI values calculated from Equation (1) for the period representing the entire growing season, e.g., from DOY 89 to 337. The ALAD value was calculated to determine yields using Equation (9) between DOY 89 and a particular EOD or between DOY 89 and 337. The coefficients of the simple linear regression for the Y_F model were obtained for the relations between the reported yields and yield estimates.

2.5.3. Comparison between the Y_P and Y_F Models

The Y_P and Y_F models were different with respect to the period and the LAI resource data used to calculate LAD. The Y_P model used LAD from the emergence to maturity dates predicted using Equations (4) and (6), and LAD was calculated from the LAI values smoothed using Equations (4) and (5). The Y_F model used LAD from DOY 89 to a particular EOD, and LAD was calculated with raw LAI values derived from Equation (1).

2.6. Classification of the Calibration and Validation Datasets

Yield data were classified into two groups for the calibration and validation of phenology and yield. Three years in which yield data were available for both Illinois and Heilongjiang were chosen randomly for validation of corn yield prediction at the state scale. Calibration datasets for predicting phenology and yield at the district scale were selected randomly from approximately 70% of the remaining datasets, including data for the 2000, 2001, and 2013 seasons, during which yield data were available in either China or the USA. In total, 69 and 62 sets of yield data were used as the calibration datasets for Illinois and Heilongjiang, respectively (Figure 7). Yield data at the district scale, e.g., AD and prefecture for the USA and China, respectively, were pooled to determine the empirical parameter values, including ${\mathsf{\tau }}_{\mathrm{P}}$, ${\mathsf{\rho }}_{\mathrm{P}}$, α, and β. The other datasets (i.e., 30 and 26 sets in Illinois and Heilongjiang, respectively) were used to validate the yield models.

${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ were determined for districts where phenology data were available in the calibration datasets. For example, all phenological dates were available in only five ADs in the USA for 2003 to 2012. The simplex method was used to determine ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ for emergence and maturity, respectively. The crop phenology prediction model was validated with districts where phenology data were available in the validation datasets. Because no phenological data were available in China, ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ obtained from the USA were used to determine whether Equation (5) could be used to predict crop yield.

2.7. Degree of Agreement Analysis

The degree of agreement statistics was determined by spatial scale, season, and variable. The RMSE value was determined to compare the observed and estimated phenology dates, as follows:

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(P_i-O_i\right)}^2 },$$

(10)

where n represents the number of comparisons, and P_i and O_i are the estimated and observed data.

Four statistical measures correlation (R²), normalized root mean square error (NRMSE), concordance correlation coefficient (CCC), and RMSE, were determined for crop yield. The yield of each grid cell was summarized by individual season and district, e.g., AD and prefecture, to compare with the reported yields at the regional scale. Yields by region were aggregated to compare the predicted and reported yields at the state scale, e.g., Illinois and Heilongjiang, by season. The NRMSE was determined as follows [47]:

$$\mathrm{NRMSE}=\mathrm{RMSE}\times \frac{100}{\mathrm{M}},$$

(11)

where M is the mean yield reported by the statistical agencies in China and the USA. The simulated results were considered to be either excellent (NRMSE < 10%), good (10% < NRMSE < 20%), fair (20% < NRMSE < 30%), or poor (NRMSE > 30%). The CCC value, which was used to represent precision and accuracy, was determined as follows [48]:

$$\mathrm{CCC}=\frac{2{\mathsf{\rho }\mathsf{\sigma }}_s{\mathsf{\sigma }}_y}{{\mathsf{\sigma }}_{\mathrm{x}}^2+{\mathsf{\sigma }}_{\mathrm{y}}^2+\left({\mathsf{\mu }}_{\mathrm{x}}-{\mathsf{\mu }}_{\mathrm{y}}\right)}$$

(12)

where $\mathsf{\rho }$ is the correlation coefficient between the estimated and reported data, and $\mathsf{\sigma }$ and $\mathsf{\mu }$ are the standard deviation and mean of the estimated and observed data, respectively. CCC ranges from −1 to 1, where −1 and 1 represent perfect disagreement and agreement, respectively, and 0 represents independence between the estimated and reported data [49].

3. Results

3.1. Crop Phenology

The estimated ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ in Equation (6) for the calibration datasets differed by the last date on which the remote sensing product was used (EOD) (Figure 8). ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ tended to group together by EOD. For example, ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ for the emergence and maturity dates were similar after EOD 257.

The RMSEs for the occurrence date of the phenological stage for the calibration datasets differed by phenological stage (Figure 9). For example, the RMSEs for maturity (<7 days) were relatively higher than those for emergence (<5 days). Temporal changes in RMSE differed by phenological stage. Although RMSE decreased with increasing EOD until EOD 257, the error values after EOD 257 were relatively similar.

The crop phenology prediction model was validated for EODs 209, 257, and 321. Corn flowering occurred on approximately DOY 209 in Illinois, and the corn harvest was completed near DOY 321. DOY 257 was selected to examine the reliability of the yield predictions in advance because the date was one of the earliest EODs on which maturity dates for corn were reliably predicted in the calibration dataset.

The occurrence dates of a given phenology stage were estimated within a reasonable range of error (Figure 10). For example, emergence and maturity dates were estimated within seven days for most districts, e.g., 67% and 75%, respectively, when these phenological stages were estimated on EOD 257.

The degree of agreement tended to be higher for the estimates of the timing of maturity than those of emergence (

Table 1. Statistical indices for the phenological stage estimate model in the validation dataset.

EOD	Emergence			Maturity
	R2	RMSE (Days)	NRMSE (%)	R2	RMSE (Days)	NRMSE (%)
209	0.51	5.31	3.98	0.27	7.61	2.95
257	0.41	6.30	4.72	0.70	4.91	1.90
321	0.35	6.29	4.72	0.78	4.43	1.71

). For example, the RMSE values for maturity were lower than those for emergence. Temporal changes in the degree of agreement also differed.

3.2. Crop Yield at the District Level

Estimated α and β values of Equation (9) for the calibration datasets changed with EOD in the Y_P model (Figure 11) and grouped together after EOD 257 because they depended on an estimated date for a phenological event, e.g., emergence or maturity. The α and β values changed over the EODs in the Y_F model using an arbitrarily fixed starting date of DOY 89.

During calibration, the Y_P model had a greater degree of agreement for predicting crop yield than the Y_F model (Figure 12). The R² values for all EODs of the Y_P model were higher than those for the Y_F model when predicting crop yield by district. Although both models predicted differences in crop yield between Illinois and Heilongjiang, the Y_P model was more precise than the Y_F model (Figure 13).

The crop yield prediction model was validated for EODs 209, 257, and 321. The corn flowering date occurred on approximately DOY 209 in Illinois, and the corn harvest was completed near DOY 321. DOY 257 was selected to examine the reliability of the yield predictions in advance because the date was one of the earliest EODs on which maturity dates for corn were reliably predicted in the calibration dataset.

The Y_P model based on the phenology dates identified from the logistic function had a lower error than that of the Y_F model based on fixed dates (

Table 2. Statistical indices for the corn yield prediction models at the agricultural district/prefecture level in Illinois and Heilongjiang.

EOD	YP Model			YF Model
	R2	RMSE (kg/ha)	NRMSE (%)	R2	RMSE (kg/ha)	NRMSE (%)
209	0.65	1158.82	13.20	0.57	1283.03	14.62
257	0.68	1083.74	12.35	0.68	1086.66	12.38
321	0.70	1042.43	11.88	0.66	1127.67	12.85

). The Y_P and Y_F models predicted differences in yield between the two regions, and corn yield was considerably lower in Heilongjiang than in Illinois. Although the Y_P and Y_F models for EOD 257 had similar R², the Y_P model had greater R² than the Y_F model for EODs 209 and 321. The Y_F model always had a higher RMSE than the Y_P model.

The errors in yield prediction were similar for EODs 257 and 321 in the Y_P model. The NRMSE of the crop yield prediction on EOD 257 was slightly greater than that on EOD 321. However, the Y_F model had the lowest error on EOD 257.

3.3. Crop Yield at the State/Province Level

The statistical indices for the corn yield prediction model over three randomly selected years (2003, 2009, and 2012) are shown in

Table 3. Statistical indices for the yield prediction models at the state level in Illinois and at the province level in Heilongjiang.

Region	EOD	YP Model				YF Model
		R2	RMSE (kg/ha)	NRMSE (%)	CCC	R2	RMSE (kg/ha)	NRMSE (%)	CCC
IL	209	0.43	1785.62	19.27	0.21	0.18	2137.85	23.07	−0.12
	257	0.87	687.68	7.42	0.93	0.99	1006.67	10.86	0.78
	321	0.95	684.72	7.39	0.91	0.94	1068.36	11.53	0.74
HE	209	0.99	1115.68	15.72	0.59	0.96	1008.13	14.20	0.68
	257	0.99	964.88	13.59	0.68	0.99	839.75	11.83	0.79
	321	0.99	664.07	9.36	0.87	0.99	816.82	11.51	0.79

IL: Illinois, HE: Heilongjiang.

. The Y_P model showed high statistical measures for the degree of agreement in corn yield prediction in Illinois and Heilongjiang. In particular, the Y_P model performed best in Illinois, while the Y_F model had a similar performance in both regions. The R² and CCC values of both prediction models and regions, except for EOD 209, were >0.87 and 0.68, respectively. The NRMSE values were 7.39–13.59. Both prediction models exhibited good corn yield prediction performance in the two regions.

The reported and predicted corn yields by state/province for the prediction models at EOD 257 in Illinois and Heilongjiang are compared in Figure 14.

4. Discussion

Simple models were developed to predict corn phenological stages and yield using only the red and NIR band surface reflectance data of the MODIS products. Sakamoto et al. [50] suggested that the temporal discontinuity of remote sensing data makes it difficult to estimate crop properties, e.g., phenological dates; therefore, different approaches have been developed. In this study, the combination of Equations (4) and (5) allowed us to overcome data discontinuity to estimate the phenological date. The remote sensing data revealed fluctuations in LAI during a season, as shown in Figure 6. However, the integrative approach using Equation (4) made it possible to identify the date of maximum LAI accurately when EOD was later than DOY 249.

The phenology prediction model performance for predicting emergence and maturity dates differed depending on EOD. The errors for predicting emergence and maturity were greater when an EOD before the maximum LAI was used. For example, a reliable and consistent value of b₂ in Equation (4) could be obtained using remote sensing data only after the flowering period until near harvest. Because ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$ depend on b₂, long-term data were used to represent the growth conditions and to obtain reliable estimates of ${\mathsf{\tau }}_{\mathrm{P}}$ and ${\mathsf{\rho }}_{\mathrm{P}}$. The phenology model accuracy for maturity date prediction was relatively higher than for the emergence date prediction. Although the phenology, including the maturity date, depends on the characteristics of the cultivar, such as the physiological responses to environmental conditions, like photoperiod and temperature [51], changes in phenology are closely related to changes in LAI [52]. For example, flowering occurs near the time of maximum LAI, which is b₂ in Equation (4). Once the time of maximum LAI is predicted accurately, the maturity date can be reliably identified. Estimating the maximum LAI using remote sensing products is more accurate than measuring the LAI early in the season when LAI is <1 [53]. Emergence dates are related to planting dates, and the planting date is influenced by the soil temperature and soil moisture and varies widely by district and season (see Appendix A). For example, 33.6 °C is the optimal soil temperature for corn germination [54]. Thus, the planting date is delayed until the temperature condition is met in a specific field [55,56], which results in large variability in the length of the effective growing season in a particular region. As a result, the errors in the emergence dates were relatively higher than the errors in the maturity dates. The present crop phenology prediction model based on MODIS–derived LAD–logistic function provided more reliable and accurate maturity date predictions than those of previous studies. Sakamoto et al. [57] reported that the RMSE of the maturity date for corn in Illinois was 5.9 days when using a two-step filtering method and WDRVI derived from MODIS data. The RMSE value of the maturity date on EOD 257 was approximately 17% less than in a previous study. In contrast, the model resulted in greater prediction error for the emergence date. For example, the RMSE of the emergence date predicted on EOD 257 was 6.3 days, whereas the value in a previous study [57] was 4.9 days. The specific shape of seasonal LAI pattern by crop were used to predict phenological dates in the previous study. The present phenology prediction model, which is based on MODIS–derived LAD-logistic function and does not need specific shape of seasonal LAI pattern by crop, can be applied to other crops.

The α and β values in the Y_P model using the predicted phenological dates grouped together after EOD 257 (Figure 11). The α and β values of the model depended on the reliability of b₃, which represents the maximum cumulative LAI, LAD. When an EOD before the maximum LAI was used, b₃ was less reliable because the maximum LAI had not been reached in the field; thus, data up to a specific period after the maximum LAI are essential to reliably predict b₃. Due to errors in remote sensing data, e.g., caused by clouds and rainfall, maximum LAI estimates can have considerable errors, even after the maximum LAI. Nevertheless, α and β were similar after a set of EOD, suggesting that a small range of α and β values could be obtained to predict corn yield. The α and β values for the Y_F model using an arbitrarily fixed starting date depended on the last DOY (i.e., EOD) used for the LAD calculation, and the error for the predicted corn yield also depended on EOD because LAI decreased after the maximum LAI. This model also required data up to a specific period (i.e., EOD 257) after maximum LAI to predict corn yield reliably.

The errors for predicting corn yield in the Y_P and Y_F models at the AD level were relatively small. The RMSE for the Y_P model was 1.08 t/ha and 1.09 t/ha at EOD 257 for 20 districts (Table 2). In previous studies, the RMSE of the predicted corn yield ranged from 1.2 t/ha to 1.7 t/ha, depending on the region of interest and season [22,58,59]. Doraiswamy et al. [59] reported an RMSE of 1.21 t/ha for corn yield in several Illinois counties. Because of differences in spatial extent and scale, e.g., district or county and extent, caution is needed when interpreting the differences in RMSE values between the current and previous studies. Nevertheless, the RMSE of both the Y_P and Y_F models tended to be smaller than those of the previous studies, which merits further validation in a variety of additional regions.

For the validation datasets for each region, the errors for predicting corn yield at the state/province scale were relatively small in both the Y_P and Y_F models. When yield prediction was performed after EOD 257, the Y_P model exhibited RMSE values of 0.68–0.69 t/ha in Illinois and 0.66–0.96 t/ha in Heilongjiang, while the Y_F model had RMSE values of 1.0–1.1 t/ha in Illinois and 0.82–0.84 t/ha in Heilongjiang (Table 3). Even for the validation dataset for both regions, the two models revealed little decay in corn yield prediction performance compared with the validation dataset for each region, with RMSEs of 0.84 and 0.93 t/ha across the two regions for Y_P and Y_F, respectively. In addition, the current model performances were comparable to previous studies: the RMSE values for predicting corn yield ranged from 0.62 to 1.45 t/ha in major USA production regions [23,27] and Sakamoto et al. [58] reported that the RMSE of the corn yield prediction was 0.83 t/ha in Illinois. The corn yield prediction errors using only LAD result from changes in specific leaf area (SLA) in a given season and region. SLA is the ratio of leaf area to leaf biomass [60] and is affected by environmental conditions [61], such as weather and disease [52]. As a result, SLA varies by season, even when the same crops are cultivated at a given site [30]. Because remote sensing products represent LAI instead of biomass, the change in SLA affects the accuracy of the biomass and yield estimations when using LAI. Therefore, the inclusion of other factors affecting SLA would improve the performance of yield prediction model using LAI and LAD as yield estimators.

Corn management practices were different between Illinois and Heilongjiang. Irrigation practices have been widely adopted in Illinois, while rainfed cultivation is a common practice in Heilongjiang. Corn cultivars were also different by latitude and temperature, and corn yield in Illinois was higher than corn yields in Heilongjiang. Nevertheless, the current yield prediction models using only LAD showed similar high performance for the two study regions, Illinois and Heilongjiang, which have different crop cultivar characteristics, crop management, and environments. This result is due to the nature of LAD used as predictor; the fluctuation in LAD accurately represents the yield variation caused by environmental stress, such as drought, cultivation management, including nitrogen fertilization and planting density, and the genetic improvement of corn hybrids [37,38,39].

5. Conclusions

Simple approaches to predicting corn phenological stages and yields were developed using a minimum MODIS product dataset. Only the red and NIR band surface reflectance data were used to estimate the LAI. Rather than using the reported techniques for filtering/smoothing the LAI data, the LAI data summed over a cropping season (LAD) were fitted to a logistic function for smoothing. A phenology prediction model was established using the MODIS–derived LAD-logistic function parameters, and it was used to predict emergence and maturity dates within a reasonable range of error. Simple linear regression models were developed to predict yield using LAD over the predicted period from emergence to maturity as a predictor variable and LAD for a predetermined period from DOY 89 to a particular EOD. Our results indicate that these simple models using LAD as a predictor variable could predict yields for the two regions of interest with considerable precision and accuracy. The model including information related to phenology exhibited slightly better performance and could be applied from a fairly early pre-harvest stage of EOD 257. In addition, the model performance showed no difference between the two regions with very different climates and cultivation methods, including cultivar and irrigation management. Irrigation practices have been widely adopted in Illinois, USA, while rainfed cultivation is a common practice in Heilongjiang, China. The approach described in this paper has potential to be applied to relatively wide agroclimatic regions with different cultivation methods and to be extended to additional crops. However, it needs to be examined further in tropical and sub-tropical regions, which are very different from the two study regions with respect to agroclimatic constraints and agrotechnologies.