Received 7 April 2016; accepted 3 June 2016; published 6 June 2016
Sugarcane (Saccharum spp. L.) is a widely grown semi-perennial crop that plays a major role in global agriculture for the production of sugar and other by-products/co-products (bioethanol, rum, bagasse, fertilizer, filter muds etc.)  . More than twenty eight million hectares of agricultural land is devoted globally to growing sugarcane, producing approximately 1.7 trillion tonnes of raw sugar each year  . Accurate and timely prediction of yield offers the global sugar industry improved efficiency and profitability by supporting decision making processes such as crop harvesting scheduling, marketing, milling and forward selling strategies. Currently, the in-season estimation of yield is undertaken using visual or destructive sampling techniques by either growers or mill funded productivity officers. However, this method is labour intensive with accuracies influenced by varied seasonal climatic conditions, crop age due to an extended harvest period and human error. Commercially available harvester mounted yield monitoring devices operated in different sugarcane growing regions  also offer varied degrees of yield mapping accuracy, but data is only available post-harvest and therefore not available to support in-season decision making. As an alternative, remote sensing technologies have been evaluated in recent years as a more accurate and cost effective method of sugarcane yield prediction  .
Time series of satellite-based remote sensing images offer a unique opportunity to assess land cover dynamics and monitoring productivity of agricultural crops at varying spatial and temporal resolutions  . Numerous studies have reported the potential of time series observations of satellite images at different resolution such as AVHRR (Advanced Very High Resolution Radiometer)  -  , MODIS (Moderate Resolution Imaging Spectroradiometer)  -  , MERIS (Medium Resolution Imaging Spectrometer)  -  , SPOT (Satellite Pour l’Observation de la Terre)  -  and Landsat TM/ETM+  -  , for monitoring vegetation, detecting land cover change, identifying crops, mapping seasonal patterns, crop rotations and crop yield prediction. In Australia, Lee-Lovick and Kirchner  first examined the potential benefits of remote sensing for sugarcane and subsequently only a few studies have reported the use of time series observation of satellite imagery for sugarcane yield prediction       . Although, in recent times, regional yield forecasts from single within ?season remote sensing observations (SPOT 5) has produced high accuracies  , the timing of image capture and the influence of weather condition and other stresses (e.g., disease or nutrient deficiencies) has resulted in some incidences of high inaccuracy.
Sugarcane production depends on several environmental factors that vary seasonally and annually and the complex interactions among them. The advantage of time series satellite observations is that it provides a more accurate indication of crop response over time, whether it be the result of environmental changes occurring from rainfall distribution, drought, nutrient deficiency and other related factors or different human management practices such as fertilizer application, pests, weeds and diseases control etc.    . Moreover, time series observation can clearly define the time of year when the regional crop reaches its maximum growth vigour, and as such indicate when an image should capture for future single image based yield forecasting. Therefore, time series observation of satellite imagery has been suggested as a more appropriate method for achieving accurate yield forecasting in sugarcane industry      .
In addition to a pure remote sensing approach, several researchers have integrated crop models when predicting sugarcane yield. Duveiller et al.,  estimated sugarcane yield at the regional scale in Brazil using fraction of absorbed photosynthetically active radiation (fAPAR) derived from the VEGETATION sensor of SPOT 4 satellite images. Mulianga et al.,  reported a time integral and spatial aggregation of normalized difference vegetation index (NDVI) data to forecast sugarcane yield in western Kenya. On Reunion and Guadeloupe islands, Bégué et al.,  compared SPOT time series acquired maximum NDVI values and integrated NDVI values to estimate yield and found comparable RMSE of 13.2 and 15 t/ha respectively. Morel et al.,  compared integrated NDVI method, Kumar-Monteith efficiency model and forced coupling MOSICAS model to estimate sugarcane yield and revealed that the remote sensing based integrated NDVI method can provide the most accurate estimation of sugarcane yield on Reunion Island. In another study, Morel et al.,  successfully applied the force coupling technique using the fraction of intercepted photosynthetically active radiation (fIPAR) derived from time series SPOT satellite to MOSICAS model to improve the sugarcane yield forecasting technique. Although remote sensing data coupled with crop models produced reasonably satisfactory results under experimental conditions in different studies, the main limitation of these models under on-firm conditions is the ability to accurately incorporate the influences of abiotic and biotic influences on crop production  .
The normalized difference vegetation index (NDVI) from satellite imagery is one of the most widely used vegetation indices to determine crop phenology, biomass, and productivity in spatial resolution. Although it has been shown to be effective in reducing the influences of atmospheric attenuation and shading  , it is susceptible to saturation at higher leaf area index (LAI), an occurrence that regularly occurs in large biomass crops such as sugarcane. To overcome these issues green normalized difference vegetation index (GNDVI) has been employed in several studies to ascertain the temporal differences and to predict sugarcane yield   . Sequential observations of GNDVI can provide seasonal crop profiles that show the progression of crop canopy from emergence to senescence. These profiles reflect the crop performance based on environmental factors and are related to the final crop yields.
The objective of this study was to develop a novel model from time series Landsat TM/ETM+ data to predict sugarcane yield in Bundaberg region. The seasonal growth profiles over a period of 15 growing seasons were developed using time series GNDVI values derived from imagery captured between mid-November to July each year. The annual timing of peak growth was identified, therefore indicating the optimal timing for future yield prediction from a single image capture. The model derived highest GNDVI values were regressed against historical sugarcane yield data to estimate sugarcane yield at the regional level. In a recently cited article of Robson, et al.  average GNDVI from February to April produced a significant correlation (R2 = 0.91) to predict sugarcane yield in Bundaberg region for the six consecutive years from 2010 to 2015. Therefore, this study also compared the average GNDVI technique and model derived maximum GNDVI technique to predict sugarcane yield using the data from 2001 to 2015.
2. Materials and Methods
2.1. Study Area
The study was carried out at Bundaberg cane growing region (Figure 1) located in the South East of Queensland in Australia. The area is located between longitudes 151.91˚E and 152.49˚E, and latitudes 24.51˚S and 25.13˚S, covering an area of 20,700 ha. The climate of Bundaberg is subtropical with long hot summers and mild winters. The soil type in this area varied enormously due to climate, substance of parent material and topography. The mean annual rainfall of the area was recorded to 964 mm in 2015  .
2.2. Remote Sensing Data
In this study, 98 Landsat TM/ETM+ L1T images from 2001 to 2015 with cloud cover less than 40% over Bun-
Figure 1. Location of sugarcane grown area on Bundaberg in Queensland.
daberg region (Path 90 Row 77; Figure) were used. Images acquired between mid-November and July each year were selected as these were identified to encompass the sugarcane growth period in that region. Satellite images were downloaded from the US Geological Survey National Center for Earth Resources Observation and Science via the GLOVIS data portal (http://glovis.usgs.gov/).
Image pre-processing is one of the important steps to make all the images similar or nearly similar so that they can be considered to be captured in the same environmental conditions and by the same sensors  . Digital image processing software ENVI 5.1 and ArcGIS 10.2 were used for the image processing, analysis, and spatial data integration. In this study, all the Landsat images were geometrically referenced to UTM projection system “WGS 1984 UTM zone 56N” to match with the supplied land cover boundary images. Radiometric normalization was used for the acquired images to reduce reflectance variations between image dates due to atmospheric conditions and surface anisotropy. For the bulk corrections of atmospheric effects simple dark-object subtraction (DOS) method was applied  . Notably minor cloud, haze, shadow or bad data pixels were not considered while processing the images
Sugarcane boundary vector layers for each year were sourced from Bundaberg Sugar Limited. To ensure that the extracted spectral information was specific to sugarcane only, a 10 m internal buffer was applied to each field boundary. All the Landsat images were masked using these boundary layers and the green normalized difference vegetation index (GNDVI) derived, using the following equation 
where, RNIR and RGREEN are the reflectance measured in the near infrared and green spectral bands.
All available average GNDVI data from the 15 year period were plotted against the date of image acquisition and a quadratic model was fitted to visualize the progression of sugarcane crop growth. This model identified when maximum vigour of the annual crop was achieved, via the peak of the quadratic curve and, therefore, indicates when images should be captured to get maximum vigour of crop and ultimately predict yield in the season. The vertex form of the quadratic model as shown in Equation (2) was used to shift the vertical axis of the curve according to the acquired GNDVI value in maximum vigour period of a specific year.
Here, “a” is a value in the curve that indicates the curvature of the line, the sign (±) on “a” tells whether the quadratic opens up or opens down, the sign (+) indicates that the quadratic opens up and the sign (−) indicates the quadratic opens down; h is the horizontal shift of the curve from x = 0, for any standard quadratic equation
, and k is the vertical shift of the curve from x = 0, for any standard quadratic equation,.
2.4. Statistical Analysis
Linear regression analysis was performed to evaluate the relationship between the model derived maximum GNDVI value and sugarcane yield. The fifteen yield forecasts from 2001 to 2015 and yield observations at harvest were evaluated and compared using root means square error (RMSE).
Here, actual means actual yield data (t/ha), predicted means predicted yield data (t/ha) from measured GNDVI values and n is the number of observation. All analysis were performed using Microsoft Excel (version 10).
3. Results and Discussion
3.1. Time Profile of Sugarcane GNDVI
The average GNDVI time profile of the Bundaberg sugarcane growing region as a function of image acquisition date of year from 2001 to 2015 is illustrated in Figure 2. Although less number of cloud free images were obtained from November till March in most of the growing season, the evolution of GNDVI profile shows a distinct growth phases such as tillering, vegetative development and the senescence of sugarcane crop. Some infrequent variations in GNDVI values are related to the phenology of the crop and climatic conditions which was reported by Bégué et al.  as well.
3.2. Quadratic Model
All available GNDVI data over the 15 year period (2001 to 2015) were plotted against image acquisition date of year (Figure 3). A quadratic model was identified to best represent the growth profile of sugarcane crop with R2
Figure 2. Time series of Landsat GNDVI data over 15 year’s period (2001 to 2015) illustrating the sugarcane crop cycles in Bundaberg region. The images were captured only during the growing period from mid-November to July each year for sugarcane crop.
Figure 3. The average GNDVI data from 2001 to 2015 in the growing period of sugarcane (mid-November to July) against image acquisition date. Image acquired in the same date of different years were averaged to apt the trend line in the data.
= 0.72. The vertex form of the model shown in Figure 3 indicates that the GNDVI value reaches to its peak after 145 days of starting date (15th November) and about three months before harvest. This result is consistent with the previous findings   , whom reported that the best acquisition period of satellite images is about two months prior to the beginning of harvest for the prediction of sugarcane yield. The highest average GNDVI value from the model is 0.58.
3.3. Determination of Highest GNDVI from the Model
The model derived GNDVI values were plotted over the calculated GNDVI values from Landsat images (Figure 4). The model was shifted vertically in each year to pass through the calculated GNDVI value acquired near or at the maximum vigour period of sugarcane. The highest GNDVI value from the model was regressed against the final average crop yield measured in that year. In Figure 3, the GNDVI graph for Year 2001, 2005, 2009 and 2014 are shown as examples. From the model, the maximum GNDVI values for those years were 0.56, 0.58, 0.55 and 0.58 respectively. Although the GNDVI values at the senescent stage in 2009 and 2014 show higher values, the calculated GNDVI values near maximum vigour were considered to match with the model derived GNDVI values.
3.4. Prediction of Sugarcane Yield
A scatter plot of model derived maximum GNDVI against the annual harvested yield in terms of tonnes of cane per hectare (t/ha) from 2001 to 2014 is shown in Figure 5. The data point of 2013 is excluded, due to an extensive flood during 2013 that prevented the harvest of around 40% of final yield. A linear relation with good agreement (R2 = 0.69 and RMSE = 4.2 t/ha) was established between the model derived maximum GNDVI and the annual harvested yield (t/ha). The annual harvested yield of 2015 (87.3 t/ha) is used as a model validation data, which showed an overestimation of predicted yield by only 3.54 t/ha in that year.
Previous research by Robson et al.  identified that the time series method can produce a more accurate prediction than a single capture SPOT 5, when GNDVI values derived between February and April were averaged. Therefore, the average GNDVI data from the imagery acquired from February to April for all years were also regressed against the annual harvested yield according to Robson et al.  (Figure 6). Here also the data point of 2013 was not considered due to extensive flooding. The linear correlation was found with R2 = 0.48 and RMSE = 5.46 t/ha. The annual harvested yield of 2015 is over estimated by 5.30 t/ha using the average February to April GNDVI data.
Figure 4. The model derived GNDVI values (shown as lines) plotted over the calculated GNDVI values from Landsat images (shown as dots) for year 2001, 2005, 2009 and 2014 as examples. The calculated GNDVI values near maximum vigour of sugarcane were matched with the model derived GNDVI to get the maximum GNDVI.
Figure 5. The scatter plot of model derived GNDVI Vs annual harvested yield (t/ha) from 2001 to 2014. Data from 2013 was excluded, as the crops were not totally harvested in that year due to extensive flood.
The maximum GNDVI extracted from all Landsat images acquired around the critical April period for each year was also regressed against the annual harvested yield data (Figure 7), producing a correlation of R2 = 0.58
Figure 6. The scatter plot of average GNDVI (February to April) Vs annual harvested yield (t/ha) from 2001 to 2014. Data from 2013 was excluded, as the crops were not totally harvested in that year due to extensive flood.
Figure 7. The scatter plot of maximum GNDVI from all available captured images Vs annual harvested yield (t/ha) from 2001 to 2014. Data from 2013 was excluded, as the crops were not totally harvested in that year due to extensive flood. The image capture dates for each of the different years are shown as shadow text.
and RMSE = 4.89 t/ha. It is likely that the regression coefficient was slightly less than that achieved from the quadratic model as the timing of image capture for each year ranged between March 11 to May 9 and therefore not within the peak, early April growth period. Using the algorithm from this linear relationship, the predicted annual harvested yield for 2015 is over estimated by 5.09 t/ha. Note that in 2015, the maximum GNDVI was extracted from the image of 24th of April, which is very near to the optimum period of image capture.
This study identified how time series Landsat imagery could be effectively used for plotting the historic temporal pattern of sugarcane crop production in the Bundaberg growing region. In the simplest terms, these historic trends provide a benchmark to following year’s on how annual crop production should progress, in terms of canopy GNDVI value. Therefore, any deviation from this trend can indicate the onset of a widespread abiotic or biotic constraint.
In terms of yield forecasting, the maximum crop vigouror GNDVI value was historically achieved at 145 days from planting i.e. early April. This period is suggested as the optimal growth stage for the acquisition of satellite imagery to be used for regional yield forecasting. The development of a quadratic model that accurately depicts the growth profile of sugarcane has provided the opportunity for the maximum GNDVI to be calculated from imagery captured between November and March. Not only does this create the opportunity for predictions to be made earlier in a current growing season, but in the event that consistent cloud cover prevents the capture of a new image during the critical April period, the maximum GNDVI value for that season can still be calculated.
Although results from this study are highly encouraging, additional research is required to model temporal sugarcane growth across other domestic and global growing regions, as the influence of environmental conditions and cropping practices will likely vary the quadratic relationship between GNDVI and yield.
The authors acknowledge the receipt of an Endeavour Research Fellowship (Rahman) to conduct this study and also Sugar Research Australia for supporting the Australian yield forecasting project in which this study complimented. The authors gratefully acknowledge the assistance of Gavin Lerch (Bundaberg Sugar Limited) for providing the historical boundary layer map and the yield data.
 Duveiller, G., López-Lozano, R. and Baruth, B. (2013) Enhanced Processing of 1-km Spatial Resolution fAPAR Time Series for Sugarcane Yield Forecasting and Monitoring. Remote Sensing, 5, 1091-1116.
 Jensen, T.A., Baillie, C., Bramley, R.G.V. and Panitz, J.H. (2012) An Assessment of Sugarcane Yield Monitoring Concepts and Techniques from Commercial Yield Monitoring Systems. Proceedings of the Australian Society of Sugar Cane Technologists, 34, 7.
 Morel, J., Todoroff, P., Bégué, A., Bury, A., Martiné, J.-F. and Petit, M. (2014) Toward a Satellite-Based System of Sugarcane Yield Estimation and Forecasting in Smallholder Farming Conditions: A Case Study on Reunion Island. Remote Sensing, 6, 6620-6635.
 Goncalves, R.R.V., Nascimento, C.R., Zullo Jr., J. and Romani, L.A.S. (2009) Relationship between the Spectral Response of Sugar Cane, Based on AVHRR/NOAA Satellite Images, and the Climate Condition, in the State of Sao Paulo (Brazil), from 2001 to 2008. The Fifth International Workshop on the Analysis of Multi-Temporal Remote Sensing Images, Groton, 315-322.
 Nascimento, C.R., Goncalves, R.R.V., Zullo Jr., J. and Romani, L.A.S. (2009) Estimation of Sugar Cane Productivity Using a Time Series of AVHRR/NOAA-17 Images and a Phenology-Spectral Model. The Fifth International Workshop on the Analysis of Multi-Temporal Remote Sensing Images, Groton, 365-372.
 Mulianga, B., Bégué, A., Simoes, M. and Todoroff, P. (2013) Forecasting Regional Sugarcane Yield Based on Time Integral and Spatial Aggregation of MODIS NDVI. Remote Sensing, 5, 2184-2199.
 Xavier, A.C., Rudorff, B.F.T., Shimabukuro, Y.E., Berka, L.M.S. and Moreira, M.A. (2006) Multi-Temporal Analysis of MODIS Data to Classify Sugarcane Crop. International Journal of Remote Sensing, 27, 755-768.
 Bastidas-Obando, E. and Carbonell-Gonzalez, J. (2007) Evaluating the Applicability of MODIS Data for Forecasting Sugarcane Yields in Colombia. International Society of Sugar Cane Technologists (ISSCT), Durban.
 Si, Y., Schlerf, M., Zurita-Milla, R., Skidmore, A. and Wang, T. (2012) Mapping Spatio-Temporal Variation of Grassland Quantity and Quality Using MERIS Data and the PROSAIL Model. Remote Sensing of Environment, 121, 415-425.
 Zurita-Milla, R., Gomez-Chova, L., Guanter, L., Clevers, J.G.P.W. and Camps-Valls, G. (2011) Multitemporal Unmixing of Medium-Spatial-Resolution Satellite Images: A Case Study Using MERIS Images for Land-Cover Mapping. IEEE Transactions on Geoscience and Remote Sensing, 49, 4308-4317.
 Meng, J., Wu, B., Li, Q., Du, X. and Jia, K. (2009) Monitoring Crop Phenology with MERIS Data—A Case Study of winter Wheat in North China Plain. Electromagnetics Research Symposium, Beijing, 1225-1228.
 El Hajj, M., Bégué, A., Guillaume, S. and Martiné, J.-F. (2009) Integrating SPOT-5 Time Series, Crop Growth Modeling and Expert Knowledge for Monitoring Agricultural Practices—The Case of Sugarcane Harvest on Reunion Island. Remote Sensing of Environment, 113, 2052-2061.
 Gers, C. and Schmidt, E. (2001) Using SPOT 4 Satellite Imagery to Monitor Area Harvested by Small Scale Sugarcane Farmers at Umfolozi. Proceedings of the 75th South Africa Sugar Technologists’ Association (SASTA), Durban, 31 July-3 August 2001, 28-33.
 Lehmann, E.A., Wallace, J.F., Caccetta, P.A., Furby, S.L. and Zdunic, K. (2013) Forest Cover Trends from Time Series Landsat Data for the Australian Continent. International Journal of Applied Earth Observation and Geoinformation, 21, 453-462.
 Lu, D., Hetrick, S., Moran, E. and Li, G. (2012) Application of Time Series Landsat Images to Examining Land-Use/ Land-Cover Dynamic Change. Photogrammetric Engineering and Remote Sensing, 78, 747-755.
 Schmidt, E.J., Narciso, G., Frost, P. and Gers, C. (2000) Application of Remote Sensing Technology in the SA Sugar Industry. Review of Recent Research Findings. South African Sugar Technologists Association, 74, 192-201.
 Lee-Lovick, G. and Kirchner, L. (1990) The Application of Remotely Sensed (Landsat TM) Data to Monitor the Growth and Predict Yields in Sugarcane. Australian Society of Sugar Cane Technology, 65-72.
 Morel, J., Bégué, A., Todoroff, P., Martiné, J.-F., Lebourgeois, V. and Petit, M. (2014) Coupling a Sugarcane Crop Model with the Remotely Sensed Time Series of fIPAR to Optimise the Yield Estimation. European Journal of Agronomy, 61, 60-68.
 Bégué, A., Lebourgeois, V., Bappel, E., Todoroff, P., Pellegrino, A., Baillarin, F. and Siegmund, B. (2010) Spatio-Temporal Variability of Sugarcane Fields and Recommendations for Yield Forecast Using NDVI. International Journal of Remote Sensing, 31, 5391-5407.
 Benvenuti, F. and Weill, M. (2010) Relationship between Multi-Spectral Data and Sugarcane Crop Yield. Proceedings of the 19th World Congress of Soil Science and Soil Solutions for a Changing World, Brisbane, 1-6 August 2010, 33-36.
 Robson, A., Rahman, M.M., Falzon, G., Verma, N.K., Ohansen, K.J., Robinson, N., Lakshmanan, P., Salter, B. and Skocaj, D. (2016) Evaluating Remote Sensing Technologies for Improved Yield Forecasting and for the Measurement of Foliar Nitrogen Concentration in Sugarcane. 38th Australian Society of Sugar Cane Technologists, Mackay, in press.
 Hall, F.G., Strebel, D.E., Nickeson, J.E. and Goetz, S.J. (1991) Radiometric Rectification: Toward a Common Radiometric Response among Multidate, Multisensor Images. Remote Sensing of Environment, 35, 11-27.
 Chavez, P.S. (1988) An Improved Dark-Object Subtraction Technique for Atmospheric Scattering Correction of Multispectral Data. Remote Sensing of Environment, 24, 459-479.
 Gitelson, A.A., Kaufman, Y.J. and Merzlyak, M.N. (1996) Use of a Green Channel in Remote Sensing of Global Vegetation from EOS-MODIS. Remote Sensing of Environment, 58, 289-298.
 Rudorff, B.F.T. and Batista, G.T. (1990) Yield Estimation of Sugarcane Based on Agrometeorological-Spectral Models. Remote Sensing of Environment, 33, 183-192.
 Ueno, M., Kawamitsu, Y., Sun, L., Taira, E. and Maeda, K. (2005) Combined Applications of NIR, RS and GIS for Sustainable Sugarcane Production. Proceedings of the International Society of Sugar Cane Technologists (ISSCT), Guatemala, 204-210.