Skip to content
BY-NC-ND 3.0 license Open Access Published by De Gruyter Open Access December 30, 2016

Comparison of MODIS 250 m products for early corn yield predictions: a case study in Vojvodina, Serbia

  • Miro Govedarica , Dušan Jovanović , Filip Sabo EMAIL logo , Mirko Borisov , Milan Vrtunski and Ivan Alargić
From the journal Open Geosciences


The aim of the paper is to compare Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) products with different compositing periods (8-day, 16-day and 8-day-dual) at 250 m spatial resolution for early corn yield estimation. In order to achieve this objective, several regression models were used where the average yield was the dependent variable and the different values of NDVI were independent variables. The various inputs in the regression models included: (i) maximum NDVI value (peak value) during the season or heading date value, (ii) NDVI values from the first date after heading date, (iii), NDVI values from the second date after heading date, (iv) Seasonally integrated NDVI values. Results showed that the 16-day composite was better yield predictor than the 8-day composite when using maximum NDVI value during the season, which is the value from the most significant earliest period for yield estimation, which is called the heading date. The 8-day composites were more useful than 16-day composites later in the season for yield estimation when NDVI values from first date after heading date and values from second date after heading were used. However, the 8-day-dual was not useful for yield prediction. In order to validate the results, the authors used the leave-one-year-out approach, which trains the remaining years for the left out year and is used for yield prediction for missing year. It was found that the inverse regression model produced the best yield estimates. After excluding the anomalous 2012 year, the R2 values for the regression model were > 0.5 for all remaining years and products, with statistical significant at 0.05. The smallest difference between predicted and actual corn yield when using 8-day composite was 0.05% while the largest difference was 34.47%, whilst in the case of 16-day composite the smallest difference between predicted and actual yield was 1.67% and the largest difference was 44.12%.

1 Introduction

Corn (maize) crop production plays an important role in the economy of Serbia. Serbia was the fourth largest corn exporter in Europe in 2011 with US$0.5 billion [1], which is notable because of Serbia stumbling industry and because of the fact that irrigation is available only in total area of 84858 hectares (ha) [2, 3] which is significantly lower compared to other corn export leading countries [2]. The Serbian economy is greatly affected by crop production, yet also by climatic conditions. Corn is one of the main crops in whole South East Europe. Production quantity of corn was highest in 2010 and 2011 compared to other significant crops, while the production quantity in CIS Europe (Belarus, Ukraine, Russia and Republic of Moldova) placed corn in fourth place [1]. According to the Agricultural Census in 2012 [3], total agricultural land in Serbia was 3,437,423 ha, from which 73.1% was arable land. The main crops grown on arable land are corn, 39% of total arable land and wheat, 24% of total arable land.

The importance of arable land in Vojvodina, the northern province of Serbia, should be stated. From total arable land of 2,513,154 ha in Serbia, Vojvodina comprises 1,466,176 ha of arable land which is considered the more fertile land than any other arable land in Serbia. The total harvested area of corn in Serbia was 976,612 ha (610.304 ha in Vojvodina) in 2012 [3].

Crop yields in Serbia are not estimated officially. There are some unofficial predictions of crop yields prior to harvest by visually inspection of the status of crops on field and comparison to previous similar years. This can be very time consuming for large sown areas. Remote Sensing in agriculture can help in early prediction of crops and in near real-time monitoring because of wide coverage in area. Using this information, the policy makers can decide whether there are enough crops for export or whether there is a need for import, much earlier in the season. Timely and accurate predictions of corn yields or corn anomalies would certainly have a great benefit for a country like Serbia.

Remotely sensed data have been widely used for crop identification, monitoring and yield estimation. The most important remote sensing tools for monitoring crops are time series Vegetation Indices (VI) with different spatial resolutions. More recent studies [4-7] preferably chose Vegetation Indices from Moderate Resolution Imaging Spectroradiometer (MODIS) rather than from Advanced Very High Resolution radiometer (AVHRR), due to better spatial resolution and higher radiometric resolution (12-bit), improved geometric registration, atmospheric correction, and cloud screening [8-10].

The most common index used for crop yields estimates is the Normalized Difference Vegetation Index (NDVI) from both MODIS and AVHRR sensors [4-7, 11-13, 15]. Few studies compared results when using NDVI time series for crop yields estimate with results obtained when using other VI such are Enhanced Vegetation Index (EVI), EVI2 and Normalized Difference Water Index (NDWI) [4, 5]. Bolton and Friedl [4] used MODIS collection 5 nadir bidirectional reflectance distribution function adjusted surface reflectance data, with 500 m spatial resolution and eight day intervals, to compute NDVI, EVI2 and NDWI for estimating corn and soybean yield. Their results found that the best times to predict crop yields are 65-75 days after greenup for corn and 80 days after greenup for soybeans. Linear regression models were used for crop predictions. EVI2 proved to be better than NDVI for corn yield predictions, while they show similar results for soybean yield predictions. Son et al. [5] used MODIS Surface Reflectance 8-day composites with 500 m spatial resolution (MOD09A1 product) for rice yield estimation. Unlike previous studies, they tested quadratic regression model for crop yield estimate. They concluded that EVI is more sensitive to rice crop production than NDVI and that EVI better estimates yields of rice crop using the heading date of rice plant growth. Mkhabela et al. [6] also used MODIS data, 10-day composite NDVI for the Canadian Prairies. This study compared results from several models: linear, exponential, power and logarithmic in order to select the best model and found that the exponential model gives the highest determination coefficients. The study conclusion is that MODIS-NDVI can be used one to two months before harvest for good yield predictions. The authors emphasize is that regression models where the VI are used as inputs will differ depending on many factors including crop, soil type and environment. However, when other factors are used in regression, not only VI, the model will change and probably will not follow linear relationship with average yield as reported in [11]. Huang et al. [7] used 1 km NDVI data for yield estimation in the region where the crops were grown together. The study used several regression models and found that the S and growth regression models give the best prediction results. Ren et al. [12] evaluated the usefulness of 10-day MODIS NDVI for winter wheat yield estimation in China. Results show that NDVI follows significant linear relationship with production with R2 values greater than 0.6 in all models.

Factors that affect crop status and crop yield are varied and thus some studies combined NDVI with other factors. Prasad et al. [11] used monthly composite continental NDVI datasets from AVHRR together with surface parameters to develop a regression model which can be used for corn and soybean estimation in Iowa. They used nonlinear regression with breakpoint for prediction. The R2 value between predicted and observed corn yield is 0.78. The study also indicated that a smaller dataset (7-10 years) provides better results than larger dataset (19 years’ data) due to fewer variations in the input parameters. Doraiswamy et al. [13] also used MODIS NDVI time series together with MODIS Land Surface Temperature in order to predict corn and soybean yields at 250 m spatial resolution.

2 Methods

2.1 Study area

The study region is located in Srem, south west province of Vojvodina (figure 1]). The methodology was tested on corn parcels which are in possession of private agricultural organization “Nova Budućnost” (NB), located in the central part of Srem. This organization provided the authors with ground truth data for monitored years. Due to MODIS spatial resolution, only larger parcels (>10 ha)cound be monitored. This type of agricultural land is not common in Vojvodina, as most land parcels are usually smaller than 2 ha and could not be detected using MODIS. However, private agricultural organizations in Vojvodina usually have larger land parcels (> 50 ha) so the pixels will not be mixed with other crops or other land cover and are therefore suitable for this kind of analyses. Other MODIS products, with spatial resolutions of 500 and 1000 m would not be appropriate for monitoring crops in a small country such is Serbia. The total area of agricultural cropland which is in possession of NB is 1400 ha. The total area of monitored parcels where the corn was sown was: 785 ha, 824 ha, 329 ha, 400 ha, 188 ha, 894 ha, 168 ha and 270 ha for the monitored years 2007, 2008, 2009, 2010, 2011, 2012, 2013 and 2014 respectively. The size of each parcel where the corn was sown was greater than 40 ha for each year.

Figure 1 Location map of Vojvodina and Srem with RapidEye image, true color display. Rapid Eye image shows the location of study area in Srem.
Figure 1

Location map of Vojvodina and Srem with RapidEye image, true color display. Rapid Eye image shows the location of study area in Srem.

Corn is sown at the beginning of April and the harvest begins approximately in October, depending on the weather conditions. The time period from sowing to heading date is approximately 75 days, with the heading date is in the middle of June or at the end of June. The duration from heading to harvesting dates is approximately 120 days.

2.2 Satellite data

The complete 8-year time series (2007-2014) of the MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m product and MODIS/Terra Surface Reflectance 8-Day L3 Global 250m SIN Grid version 5, abbreviated as MOD13Q1 and MOD09Q1 respectively, were downloaded through the online Data Pool at the NASA Land Processes Distributed Active Archive Center (LP DAAC) [15] free of cost.

The data was processed online prior to download. The projection was set to Universal Transverse Mercator, zone 34 and image format was GeoTIFF. The M0D13Q1 product provides both EVI and NDVI 16-day composite images at 250-m spatial resolution, which is the highest possible MODIS spatial resolution. The pixels for VI generation are the pixels with the best quality from 16-day composite period. The MOD09Q1 product provides infrared band (band 2) and red band (band 1) 8-day composite images at 250 m spatial resolution. Bands 1 and 2 were used to generate 8-day NDVI composites. A modified 8-day composite called 8-day-dual composite was also used. This product was created from 8-day composite in order to match the 16-day composite and was previously used by Lee [14]. The 8-day-dual composite is generated by applying maximum value composite (MVC) algorithm to two consecutive images of 8-day composite in order to match the date from 16-day composite. For example, two 8-day composite images with dates: 2007.01.01 and 2007.01.09 were composited with MVC and matched with the first 16-day composite which is 2007.01.01 [14]. The equation for NDVI is as follows:


Landsat satellite images from Landsat 8, Landsat 7 and Landsat 5 were used monitor and interpret locations of sown crops for all of the mentioned years. Landsat data was used only as auxiliary data to MODIS, as it is impossible to visually discriminate crops from MODIS images or from MODIS Vis, due to the coarse spatial resolution. It was vital for us that pixels used in study were not mixed, which means homogenous corn pixels for all the monitored years. In order to avoid mixed pixels, a grid with cell size 250 χ 250 m (which corresponds to MODIS spatial resolution) was used and overlayed with Landsat and RapidEye multispectral images for all the monitored years and for all the crop parcels. If the grid pixel was entirely within the specific crop parcel, this pixel was sampled and was taken into further consideration. If the grid pixel was not entirely within the crop parcel, the pixel was ruled out of the statistical analysis. This kind of approach should provide the most reliable results for the study. Except for Landsat images, RapidEye multispectral images with 5 m spatial resolution (5 bands) were also used as an auxiliary data only for 2013 and 2014 years.

2.3 Ground truth data

All of the previously mentioned work used crop statistical data such are average crop yields for previous years and also ground truth data, which can be used to identify locations of planted crops. The ground truth data consisted of sowing structures for the years 2011, 2012, 2013 and 2014. Crop statistic data (average yield, sown corn area) was available for all years for the study area for the years 2007-2014. Sowing structures for missing years were reimbursed with MODIS classification and with auxiliary data, which consisted of high resolution images from Landsat mission and RapidEye. Sowing structures were critical in the study as they represented location maps of sown crops for mentioned years and were used for crop masking. The locations of crops change every year and thus it was important to locate planted corn in an image for each year. This data was enhanced with high resolution images. All of the data information used in the study can be seen in table 1.

Table 1

Satellite and crop data used in the study for available years.

Satellite dataCrop data
MODIS productsAuxiliary dataArea ofAverageSowing
MOD09Q1MOD12Q1Landsat 5,7,8RapidEyesown cropsyieldstructures
2007-20142007-20142007-20142013, 20142007-20142007-20142011-2014

2.4 Time series filtering

Before applying any kind of classification or analyses to NDVI time series data, the products were subjected to filtering. When using data for a longer data set (2007-2014), it is highly likely that the composites (16-day and 8-day) will be affected mainly by cloud contaminations, which can seriously influence final results. Figure 2]. indicates that one corn pixel from 8-day composites suffers from many variations around heading date and during the growth up period. Even though the 16-day composites already underwent the MVC method in order to minimize these effects, a lot of errors still exist. The 8-day composites were not subjected to any kind of filtering techniques.

Figure 2 Unfiltered Corn pixel and filtered corn pixel.
Figure 2

Unfiltered Corn pixel and filtered corn pixel.

Many filtering techniques have been developed to smooth the NDVI time-series. In this study, the Whittaker smoother [16] technique was used for smoothing NDVI composites. Recently, this smoothing technique received broader attention in remote sensing analyses of timeseries images [17-19]. It is based on penalized least squares and the smoother the results, the more it will deviate from actual input values [18]. Two parameters are controlled by the user, the number of iterations and the parameter k, which is the smoothing parameter [20]. This smoothing technique was used since the Whittaker smoother performed very well [18, 21] and even outperformed other techniques [18]. The software used for smoothing the NDVI time-series was SPIRITS [20, 22]. Smoothing was applied to the 16-day, 8-day and 8-day-dual composites.

Figure 3 Comparison of NDVI profiles for a) Corn and b) Soybean in 2013.
Figure 3

Comparison of NDVI profiles for a) Corn and b) Soybean in 2013.

2.5 Crop type separability

Lee [14] discussed crop type separability when comparing MODIS 16 day, 8 day and dual-8-day NDVI products and concluded there are significant differences in NDVI profiles for these products at 250 m spatial resolution. Similarly, this study showed that indeed there are significant differences between NDVI profiles when comparing mentioned composites. As it can be seen from Figure 2], the 16-day composites have a steady NDVI profile both for corn and soybean, without greater perturbations. However, 8-day composites show dissimilarities during the entire season. NDVI values were lower during the growth up period and during the senescence period. The maximum NDVI value occurred after the maximum NDVI value from 16-day composite and is evidently higher. Dual 8-day composites showed a similar profile during the growth up period, where the maximum value occurred before maximum values compared to 16-day and 8-day NDVI. As well, the peak greenness has the highest value. The senescence period is similar to 16-day composites, which it almost overlaps. The 8-day and 8-day-dual NDVI profiles show much greater perturbations than the 16 day composites for all crops, which can be the result of true NDVI values or affected by cloud contaminations or other atmospheric effect. Higher NDVI values for 8-day-dual composites may be the consequence of MVC algorithm applied to the raw dataset, i.e. bands 1 and 2. The authors sought to examine whether these dissimilarities can affect the final corn yield and in what measure.

The profiles are very similar for corn and soybean, with the only difference being that NDVI values for soybean are the highest for all products and for all crops that were analyzed in this study. This was also noticed by Lee [14] and was later used in classification purposes. Figure 4] shows that all crops extracted from 8-day-dual composite have higher NDVI values than the ones from 16-day composite. Also, clear separability can be noticed between sugar beet and corn, soybean is very notable after middle of August.

Figure 4 Comparisons of different NDVI profiles for 16-day and 8-day dual composites.
Figure 4

Comparisons of different NDVI profiles for 16-day and 8-day dual composites.

Prior to examining the VI time series for corn, MODIS NDVI phenology based classification was developed [23] and was used in this study. Phenology based classification was used for the years where the crop map data was not available (2007, 2008, 2009 and 2010). First, possible parcels of corn were identified by visually analyzing multispectral, multitemporal Landsat images from the same period as corresponding MODIS NDVI products. Parcels which are visually identified as corn are called potential corn locations and are used for masking. Then, based on the years for which crop maps are available, we defined threshold values for specific day of the year (DOY) when the discrimination between corn, soybean, sugar beet and sunflower is possible, since these crops are often mixed, especially corn and soybean. If the potential corn locations overlap with MODIS identified corn parcels and if the size of an identified parcel corresponds to that reported in ground truth data (crop statistic), then the pixels of that parcel belong to corn. In all other cases, the pixels are rejected as corn. Certain upper or lower threshold values for crops were noted for the years: 2011, 2013 and 2014. These threshold values were consistent throughout all the training years. Year 2012 was not used in this classification due to bad weather conditions and thus the NDVI resulted in low values. NDVI phenology based classification resulted in identification of corn pixels for 2007, 2008, 2009 and 2010 year. For the purposes of this classification, the MODIS 13Q1 product was used due to greater consistency during all monitored years.

Simple classification of Landsat images could not be used dueto the vegetation dynamics during the season. That is, it is very hard to discriminate corn from other “green” crops by using only one to two Landsat images. The reflectance values of corn, soy bean, sunflower and sugar beet were very similar for the study area. However, MODIS 16-day composite VI’s could be used together with high and medium resolution images for this kind of classification.

2.6 Model formulation and validation

MODIS images were first subset over the study region. The images were then stacked for each year, and for each product. As previously mentioned, in order to avoid mixed, heterogeneous pixels, a grid that corresponds to MODIS spatial resolution (250 χ 250 m) was overlaid with Landsat or RapidEye images and the polygons that contained entire corn MODIS pixels for each year were drawn and used for sampling.

Based on the all identified parcels of corn for all the monitored years, several regression models were tested. These included linear, exponential, inverse, growth, S and quadratic regression models. Previously mentioned studies also used several regression models and the linear model consistently gave the best prediction results.

Some studies used seasonal maximum NDVI as the main input parameter [24] whereas other studies combined different approaches and used accumulated NDVI (EVI) for a specific period. Son et al. [5] used integrated VI values after the heading date while Ren et al. [12] used several methods that included accumulation of NDVI for different periods. The authors used several methods, seasonally integrated NDVI [13, 25, 26], the seasonal NDVI peak value [24], and single NDVI values after the heading (peak) date. That is, NDVI values were used from dates n, n+1 and n+2, where n is the heading date. The first methodology (entire season) captured the effect of climatic conditions before and after the flowering stage of corn. The second approach involved identifying maximum NDVI value for all composites and is very useful for early prediction of yield; however, what happens after the maximum value is not included in the model. The peak values for the study area and for years 2007-2014 occurred most frequently at the beginning of July; however, in the case of 2012 year the peak value was in the middle of June. Values after the heading date were used to predict yield at later stages and were very useful at this stage for comparison of products. Nevertheless, since the heading dates were at the beginning of July, n+1 and n+2 dates can also be characterized as early predictions. In order to validate regression models, the authors used leave one year out approach [5, 6]. This methodology consists of removing one year at a time and generating new regression models. Models were then used for predicting yield for the missing year. Figure 5] displays an overall workflow of this paper.

Figure 5 Flowchart of methodologies used in the paper.
Figure 5

Flowchart of methodologies used in the paper.

3 Results

Amongst all the tested regression models, inverse and linear model proved to be most accurate for predicting yield for each year. Inverse was better than linear model, so no results are shown for linear model. The equation for inverse regression models is:


Where Y denotes corn yield (dependent variable), a and b are the coefficients of the model which vary with years and i are the prediction dates, which were heading date (n), first date after heading date (n+1), second date after heading (n+2). NDVI (independent variable) is the value from corresponding date i.

Other used models were ruled out as inappropriate for the study area. Also, seasonally integrated NDVI values gave poor results for all used models thus were not considered further. As mentioned earlier, the authors strongly emphasized the prediction of yield for 2012 anomalous year and to early corn predictions. Early corn yield predictions include period with maximum NDVI values for each of products, but also two dates after the peak value (n+1 and n+2).

8-day-dual composite, did not produce good results in terms of final yield prediction. Coefficients of determination between actual and predicted yield were lower than 0.3 in all used models, which meant that there were high dissimilarities between modeled yield and true yield.

Table 2 presents a comparative analysis between two products with different composite period for all years and for 3 values used in inverse regression model. The authors also provided statistical indicators which can determine whether the model is adequate, coefficients of determination (R2) and statistical significance (p-values). According to R squared values, it was concluded that the model is adequate for all years except for the prediction 2012 anomalous year, where the values were much lower than for other years. Also the p-values for 2012 year are unusually high. The 2012 year resulted in very poor corn yield which was 1.748 t/ha. Therefore, removing the year 2012 had a great impact on all statistical values in the model. It can be said that when 2012 year was included in the model, the R2 values were >0.5 for all years and were significant at 0.05. Table 3 shows coefficients which multiply 1/NDVI as well as constants for used prediction dates for the inverse regression model. Also, summary of t-statistics (t), standard error and corresponding significance values is provided in the table. It can be seen that almost all of the coefficients are significant at 0.05 and the values of the coefficients can be accepted. Year 2012 evidently disrupts the inverse model where the p-values are again higher than 0.05 so the values of the coefficients are nor well represented.

Table 2

Comparison of actual and predicted yields for two products (16-day and 8-day) for all years (2007-2014) and for all values used in inverse regression model, where n represents heading date (period when NDVI is the highest), n+1, n+2 represent values from the next two following dates. R squared and p-values are also shown for inverse model as a measure of prediction quality.

NDVI Product

Table 3

Standardized coefficients, t-statics, standard errors and corresponding p-values for two NDVI products in regards of prediction dates (n, n+1 and n+2). NDVI in the denominator corresponds to specific prediction date.

NDVI product
dateCoeficientsSth. Errortp-valueCoeficientsSth. Errortp-value
n1 / NDVI-46.9013.56-3.460.0181 / NDVI-69.5728.39-2.450.058
2007n+11 / NDVI-35.527.42-4.790.0051 / NDVI-75.6720.04-3.780.013
n+21 / NDVI-21.722.92-7.440.0011 / NDVI-67.4213.34-5.050.004
n1 / NDVI-49.2913.36-3.690.0141 / NDVI-68.6128.33-2.420.060
2008n+11 / NDVI-39.676.02-6.590.0011 / NDVI-75.3920.07-3.760.013
n+21 / NDVI-23.342.24-10.410.0001 / NDVI-71.3012.43-5.740.002
n1 / NDVI-44.3012.81-3.460.0181 / NDVI-69.3523.36-2.970.031
2009n+11 / NDVI-34.196.25-5.470.0031 / NDVI-71.9118.88-3.810.013
n+21 / NDVI-20.573.04-6.770.0011 / NDVI-64.8611.93-5.440.003
n1 / NDVI-46.2613.66-3.390.0201 / NDVI-67.6328.13-2.400.061
2010n+11 / NDVI-35.537.51-4.730.0051 / NDVI-74.8520.12-3.720.014
n+21 / NDVI-21.843.30-6.610.0011 / NDVI-69.9713.40-5.220.003
n1 / NDVI-56.499.21-6.130.0021 / NDVI-72.2927.24-2.650.045
2011n+11 / NDVI-37.286.21-6.000.0021 / NDVI-78.8618.40-4.290.008
n+21 / NDVI-21.313.30-6.460.0011 / NDVI-67.9512.89-5.270.003
n1 / NDVI-27.8419.40-1.440.2111 / NDVI-2.5851.47-0.050.962
2012n+11 / NDVI-27.9812.45-2.250.0751 / NDVI-57.7249.81-1.160.299
n+21 / NDVI-17.354.70-3.690.0141 / NDVI-78.3134.77-2.250.074
n1 / NDVI-45.8811.36-4.040.0101 / NDVI-83.2616.96-4.910.004
2013n+11 / NDVI-34.317.03-4.880.0051 / NDVI-81.149.35-8.680.000
n+21 / NDVI-21.553.56-6.050.0021 / NDVI-66.0310.48-6.300.001
n1 / NDVI-44.1814.95-2.960.0321 / NDVI-61.7230.13-2.050.096
2014n+11 / NDVI-34.518.20-4.210.0081 / NDVI-70.9421.29-3.330.021
n+21 / NDVI-20.763.61-5.740.0021 / NDVI-63.2512.60-5.020.004

However, even with poor R squared and p values the model used for the prediction of yield for 2012 year in case of 16-day composite gave good prediction results, when the input was maximum NDVI value in the season (n). This prediction result is much better than for 8-day composite, 5.432 t/ha predicted yield. Prediction results tend to be better when the input value was later in the season, n+1 and n+2, for both MODIS products, as indicated in figures 6a, 6b, 6c and 8a, 8b, 8c.

Figure 6 Relationship between final Yield (t/ha) and predicted Yield (t/ha) using 8-day composite and a) NDVI heading date value b) Value from first date after heading date and c) Value from second date after heading date.
Figure 6

Relationship between final Yield (t/ha) and predicted Yield (t/ha) using 8-day composite and a) NDVI heading date value b) Value from first date after heading date and c) Value from second date after heading date.

The study focused primarily on early prediction date, i.e. maximum NDVI values which most often occurred at the beginning of July. In this case, the NDVI 16-day composite gave slightly better results than NDVI 8-day composite, figure 6a) and figure 8a). However, due to the fact that 2013 year resulted in low average yield and none of the models successfully predicted this early in the season, the authors analyzed the removal of this year from comparison between predicted and actual yield for both products, figure 7a) and b). The removal influenced significantly on the increase of R squared value for the 16-day product and not so much for the 8-day product.

Figure 7 Relationship between final Yield (t/ha) and predicted Yield (t/ha) using a) Maximum NDVI value from 8-day composite and b) Maximum NDVI value from 16-day composite when the actual yield and predicted yield for 2013 year is removed.
Figure 7

Relationship between final Yield (t/ha) and predicted Yield (t/ha) using a) Maximum NDVI value from 8-day composite and b) Maximum NDVI value from 16-day composite when the actual yield and predicted yield for 2013 year is removed.

Figure 8 Relationship between final Yield (t/ha) and predicted Yield (t/ha) using 16-day NDVI composite and a) NDVI heading date value b) Value from first date after heading date and c) Value from second date after heading date.
Figure 8

Relationship between final Yield (t/ha) and predicted Yield (t/ha) using 16-day NDVI composite and a) NDVI heading date value b) Value from first date after heading date and c) Value from second date after heading date.

Unlike the 2012 anomalous year, weather conditions during 2014 year (from May) were very advantageous for corn in Serbia. Very large amounts of rain resulted in very good corn yield, 11.4 t/ha for NB. The 8-day product produced the best predictions for this year when using maximum VI value and, when predicting later in the season, the values of yield dropped but not significantly. The highest value of NDVI came early in the season and reflected the final yield very well (absolute difference of 2.39% between predicted and actual yield).

The R2 values during all monitored years and for both products were quite positive (> 0.5 in almost all the years) and increased later in the season (n+1, n+2). This can be seen on figure 9a) and b). The anomalous year (2012) resulted in significantly lower R2 values; however, they also increased when using dates for prediction later in the season. Considering only 2014, which is considered an extremely good, year, the R2 values were also lower than other but fairly good. Overall, the 8-day composite gave better coefficients of determination for all the monitored years.

Figure 9 Comparison of R2 values for all years in regards of the NDVI values used in the regression model where n denotes heading date, n+1 first date after heading date and n+2 the second date after heading date for a) 8-day composite b) 16-day composite products.
Figure 9

Comparison of R2 values for all years in regards of the NDVI values used in the regression model where n denotes heading date, n+1 first date after heading date and n+2 the second date after heading date for a) 8-day composite b) 16-day composite products.

4 Discussion

This paper explored the potential of MOD13Q1 and MOD09Q1 products for early corn yield prediction in the small province of Vojvodina, Srem. The study used only pure corn pixels, with 250 m spatial resolution. This was achieved by using high resolution satellite images and ground truth data. The parcels which were monitored were very convenient for this application because of large size (> 40 ha) and thus it can be concluded that the results are reliable.

Regarding early corn yield estimation (2-3 months before harvest), Mkhabela et al. [27] used linear model and found that cumulative average NDVI during flowering and grain filling produced good results. However, this study did not find that accumulating VI values gave good predictions. Wang et al. [28] also used linear model for yield prediction. When using NDVI at the time of peak correlation, which varied with years, NDVI explained 38.5% of the yield in the linear regression [28]. This methodology of using values in time of peak correlation is not consistent and sometimes can be deep during the senescence period of corn. Research [13] found that the relationship between actual and predicted yield resulted in 0.78 value of R2 when using NDVI composites from NOAA AVHRR, whilst this study found similar value earlier in the season or higher value for 16-day composite when ruling out 2013 year, but for a much smaller area and much higher spatial resolution.

Son et al. [5] used several input values in the regression model, including values from heading date and accumulated values from heading date towards five extended dates. This temporal accumulation resulted in lowering the R2 values later in the season for each date after the heading date. Unlike their approach, the authors used only single NDVI values after the heading date (without accumulation), which resulted in increasing the R2 values later in the season, figure 9. The conclusion is that the senescence period of corn is crucial for yield prediction (n+2). However, the authors did not inspect the R2 values after the n+2 date. It can be concluded from all the products that the 8-day-composite with NDVI value from n+2 date gave the best results, Table 2, figure 6c).

Considering the comparisons of MODIS or, and AVHRR products in the previously related work, few papers dealt with the possibilities of comparative analysis of MODIS or other products with different compositing period for yield estimate. Research [29] previously discussed the comparisons of EVI and NDVI (MOD13Q1) indices for crop mapping. Gallo et al. [30] compared MODIS and AVHRR NDVI 16-day composites. Lee [14] used 8-day, 16-day and 8-day-dual 250 m composites in order to evaluate crop type separability. Chen et al. [31] compared MODIS daily and 8-day composite products 356 for floodplain and wetland mapping. Yi et al. [32] evaluated two versions (4 and 5) of MODIS surface reflectance values at 500 m spatial resolution for wheat leaf area index retrieval.

There is a gap in the literature regarding the comparative analyses of MODIS products with different temporal resolution (compositing period) at any spatial resolution for yield estimation and crop monitoring. Researchers tended to use MOD09Q1 derived NDVI or EVI2 with 8-day composite period since this product can capture variations of crops during the whole season better than MOD13Q1. However this product is more vulnerable to cloud or other atmospheric contaminations than MOD13Q1 and it is possible that it will not always produce reliable results. Future studies can focus more on comparing MODIS products composites (8-day and 16-day) not only with 250 m spatial resolution but with coarser resolution, i.e. 500 m and 1 km such are MOD09A1, M0D13A1, MOD13A2, etc.

5 Conclusion

All of the previously mentioned studies focused mainly on finding the most reliable regression model for crop yield prediction and also on identifying the period which is most correlated with final yield. However, the main aim of this study was to compare products with different compositing period, i.e. MODIS 250-meter products (M0D13Q1 and MOD09Q1) for early corn yield estimates during an 8-year period from 2007 to 2014 in Srem, a small region in Vojvodina. More specifically, the authors compared 250 m NDVI 16-day composite from M0D13Q1 product, 250 m NDVI 8-day composite derived from MOD09Q1 surface reflectance values and 250 m 8-day-dual NDVI also from MOD09Q1 [14]. The earliest prediction date for this study is the beginning of July, approximately two and a half months before the beginning of the harvest. Several regression models were used in order to achieve our objective. Linear and inverse models proved to be the most efficient for this study area. Results imply that MOD09Q1 did not outperform M0D13Q1 when considering the earliest period for corn yield estimate. However, when using values after the heading date, the MOD09Q1 proved to be more efficient. 8-day-dual composite was not a good predictor of corn yield in any of models The MODIS phenology based classification needs to be considered in order to be able to monitor all necessary years. The quality of predictors is strongly associated with the prediction of yield for 2012 year, as 2012 was anomalous year, which suffered from drought from the beginning of June until the end of August. Another conclusion is that corn yield was more related with NDVI later in the season, figure 9 and the values of R2 tended to increase for both products. Therefore, by selecting only single values of NDVI after the seasonal maximum, good predictions are achieved which on the other hand have the limitation that the yield is predicted later in the season.

The main contribution of this study is that it showed progress in terms of early crop yield forecasting in Serbia, also showing that the frequently used MOD09Q1 product is not always the best solution for crop yield estimation and using only single NDVI values from heading date will improve the estimates for both products.

The methodologies mentioned should be used with caution, since any occurrences to the crop after the estimate will not be included in the final yield. The prediction models will probably change and adapt depending on the climatic conditions, environment and soil type and it should be said that the quality of MODIS VI’s is crucial in these situations.

Future analysis will be focused on estimating yields for major crops in Serbia on a regional scale or for entire country. The authors will also focus on developing models that examine the plant senescence period with more inputs in regression. One parameter that can affect final yield is LST used by Doraiswamy et al. [13] but the main limitation is its spatial resolution of 1 km. Future studies could focus more on finding crucial date for forecasting various crop yields as well on the comparison of multiple MODIS NDVI/EVI products with different spatial resolutions (500 m, 1 km) for much larger areas.

Finally, this study represents the first ever attempt to estimate yields in Serbia by using remote sensing methods and probably any kind of early yield prediction. Subject to constraints in number of years monitored, crop maps and the satellite data, study produced promising results.

Tel.: +381638015959


Special thanks to the employees at the agricultural organization “Nova Budućnost” for providing ground truth and crop statistic data. Also, Thanks to the employees of Nova Budućnost and the local farmers for providing information regarding crop calendar. Many thanks to USGS for making MODIS and Landsat data free and open access.


[1] Fao statistical yearbook Europe and Central Asia Food and Agriculture. Budapest, 2014.Search in Google Scholar

[2] Census of Agriculture. Agriculture in the Republic Of Serbia. Statistical Office of the Republic of Serbia, 2012.Search in Google Scholar

[3] Fao statistical yearbook World Food and Agriculture. Rome, 2013.Search in Google Scholar

[4] Bolton, D. K. and Friedl, M. A. 2013. Forecasting crop yield using remotely sensed vegetation indices and crop phenology metrics. Agricultural and Forest Meteorology, 173, 74-84.10.1016/j.agrformet.2013.01.007Search in Google Scholar

[5] N.T. Son, Chen C.F., Chen C.R., Minh V.Q. and Trung N.H. 2014 A comparative analysis of multitemporal MODIS EVI and NDVI data for large-scale rice yield estimation. Agricultural and Forest Meteorology, 197, 52-64.10.1016/j.agrformet.2014.06.007Search in Google Scholar

[6] Mkhabela, M. S., Bullock, P., Raj, S., Wang, S. and Yang,Y. 2011. Crop yield forecasting on the Canadian Prairies using MODIS NDVI data. Agricultural and Forest Meteorology, 151(3), 385-393.10.1016/j.agrformet.2010.11.012Search in Google Scholar

[7] Huang, J., Wang, H., Dai, Q. and Han, D.2014. Analysis of NDVI Data for Crop Identification and Yield Estimation. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 7(11), 4374-4384.10.1109/JSTARS.2014.2334332Search in Google Scholar

[8] Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X. and Ferreira, L. G. 2002. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1), 195-213.10.1016/S0034-4257(02)00096-2Search in Google Scholar

[9] Vermote, E. F., El Saleous, N. Z. and Justice, C. O. 2002. Atmospheric correction of MODIS data in the visible to middle infrared: first results. Remote Sensing of Environment, , 83(1), 97-111.10.1016/S0034-4257(02)00089-5Search in Google Scholar

[10] Wolfe, R. E., Nishihama, M., Fleig, A. J., Kuyper, J. A., Roy, D. P., Storey,J.C.and Patt,F.S.2002. Achieving sub-pixel geolocation accuracy in support of MODIS land science. Remote Sensing of Environment, 83(1), 31-49.10.1016/S0034-4257(02)00085-8Search in Google Scholar

[11] Prasad, A. K., Chai, L., Singh, R. P. and Kafatos, M. 2006. Crop yield estimation model for Iowa using remote sensing and surface parameters. International Journal of Applied Earth Observation and Geoinformation, 8(1), 26-33.10.1016/j.jag.2005.06.002Search in Google Scholar

[12] Ren, J., Chen, Z., Zhou, Q. and Tang, H. 2008. Regional yield estimation for winter wheat with MODIS-NDVI data in Shandong, China. International Journal of Applied Earth Observation and Geoinformation, 10(4), 403-413.Search in Google Scholar

[13] Doraiswamy, P. C., Akhmedov, B., Beard, L., Stern, A. and Mueller, R. 2007. Operational prediction of crop yields using MODIS data and products. International archives of photogrammetry, remote sensing and spatial information, Sciences Special Publications, Commission Working Group VIII WG VIII/10. Ispra, Italy: European Commision DG JRC-Institute for the Protection and Security of the Citizen, 1-5.Search in Google Scholar

[14] Lee, E. 2014. Analysis of MODIS 250 m NDVI Using Different Time-Series Data for Crop Type Separability. PhD thesis, University of Kansas, Kansas.Search in Google Scholar

[15] United States Geological Survey (USGS) – Reverb EOS Clearing House (echo). NASA Land Processes Distributed Active Archive Center (LP DAAC). Online Data PoolSearch in Google Scholar

[16] Eilers, P. H. 2003. A perfect smoother. Analytical Chemistry, 75(14), 3631-363610.1021/ac034173tSearch in Google Scholar PubMed

[17] Atzberger, C. and Eilers, P. H. 2011. Evaluating the effectiveness of smoothing algorithms in the absence of ground reference measurements. International Journal of Remote Sensing, 32(13), 3689-370910.1080/01431161003762405Search in Google Scholar

[18] Atkinson, P. M., Jeganathan, C., Dash, J. and Atzberger, C. 2012. Inter-comparison of four models for smoothing satellite sensor time-series data to estimate vegetation phenology. Remote Sensing of Environment, 123, 400-417.10.1016/j.rse.2012.04.001Search in Google Scholar

[19] Wijesingha, J. S. J., Deshapriya, N. L. and Samarakoon, L. 2015. Rice crop monitoring and yield assessment with MODIS 250m gridded vegetation product: A case study in Sa Kaeo Province, Thailand. International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences, 40(7), 121-127.10.5194/isprsarchives-XL-7-W3-121-2015Search in Google Scholar

[20] Rembold, F., Meroni, M., Urbano, F., Royer, A., Atzberger, C., Lemoine, G. and Haesen, D. 2015. Remote sensing time series analysis for crop monitoring with the SPIRITS software: new functionalities and use examples. Frontiers in Environmental Science, in Google Scholar

[21] Geng, L., Ma, M., Wang, X., Yu, W., Jia, S. and Wang, H. 2014. Comparison of eight techniques for reconstructing multisatellite sensor time-series NDVI data sets in the Heihe river basin, China. Remote Sensing, 6(3), 2024-2049.10.3390/rs6032024Search in Google Scholar

[22] Eerens, H., Haesen, D., Rembold, F., Urbano, F., Tote, C. and Bydekerke, L. 2014. Image time series processing for agriculture monitoring. Environmental Modelling & Software, 53, 154-162.10.1016/j.envsoft.2013.10.021Search in Google Scholar

[23] Govedarica, M., Jovanović, D. and Sabo, F. 2015. Corn yield estimation in Serbia using MODIS 13Q1 product. Proc. SPIE 9535, Third International Conference on Remote Sensing and Geoinformation of the Environment (RSCy2015), 95351D.10.1117/12.2192331Search in Google Scholar

[24] Becker-Reshef, I., Vermote, E., Lindeman, M. and Justice, C. 2010. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data. Remote Sensing of Environment, 114(6), 1312-1323.10.1016/j.rse.2010.01.010Search in Google Scholar

[25] Labus, M. P., Nielsen, G. A., Lawrence, R. L., Engel, R. and Long, D. S. 2002. Wheat yield estimates using multi-temporal NDVI satellite imagery. International Journal of Remote Sensing, 23(20), 4169-4180.10.1080/01431160110107653Search in Google Scholar

[26] Rojas,O. 2007. Operational maize yield model development and validation based on remote sensing and agro-meteorological data in Kenya. International Journal of Remote Sensing, 28(17), 3775-3793.10.1080/01431160601075608Search in Google Scholar

[27] Mkhabela, M. S., Mkhabela, M. S. and Mashinini, N. N. 2005. Early maize yield forecasting in the four agro-ecological regions of Swaziland using NDVI data derived from NOAA’s-AVHRR. Agricultural and Forest Meteorology, 129(1), 1-9.10.1016/j.agrformet.2004.12.006Search in Google Scholar

[28] Meng, W. A. N. G., Tao, F. L., Shi, W. J., 2014. Corn Yield Forecasting in Northeast China Using Remotely Sensed Spectral Indices and Crop Phenology Metrics. Journal of Integrative Agriculture, 13(7), 1538-154510.1016/S2095-3119(14)60817-0Search in Google Scholar

[29] Wardlow, B. D. and Egbert, S. L. 2010. A comparison of MODIS 250-m EVI and NDVI data for crop mapping: a case study for southwest Kansas, International Journal of Remote Sensing, 31(3), 805-830.10.1080/01431160902897858Search in Google Scholar

[30] Gallo, K., Ji, L., Reed,B., Dwyer, J. and Eidenshink, J. 2004. Comparison of MODIS and AVHRR 16-day normalized difference vegetation index composite data. Geophysical Research Letters, 31(7) doi: 10.1029/2003GL019385.Search in Google Scholar

[31] Chen, Y., Huang, C., Ticehurst, C., Merrin, L. and Thew, P. 2013 An evaluation of MODIS daily and 8-day composite products for floodplain and wetland inundation mapping. Wetlands, 33(5), 823-835.10.1007/s13157-013-0439-4Search in Google Scholar

[32] Yi, Y., Yang, D., Huang, J. and Chen, D. 2008. Evaluation of MODIS surface reflectance products for wheat leaf area index (LAI) retrieval. ISPRS Journal of Photogrammetry and Remote Sensing, 63(6), 661-677.10.1016/j.isprsjprs.2008.04.004Search in Google Scholar

Received: 2015-9-29
Accepted: 2016-7-19
Published Online: 2016-12-30
Published in Print: 2016-1-1

© 2016 M. Govedarica et al.

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.

Downloaded on 2.3.2024 from
Scroll to top button