Crop Yield Forecasting on the Canadian Prairies by Satellite Data and Machine Learning Methods by Michael David Johnson B.Sc., The University of Winnipeg, 2011 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Atmospheric Science) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2013 ? Michael David Johnson, 2013 ii Abstract The production of grain crops plays an important role in the economy of the Canadian Prairies and early reliable crop yield forecasts over large areas would help policy makers and grain marketing agencies in planning for exports. Forecast models developed from satellite data have the potential to provide quantitative and timely information on agricultural crops over large areas. The use of nonlinear modeling techniques from the field of machine learning could improve crop forecasting from the linear models most commonly used today. The Canadian Prairies consist of the provinces of Alberta, Saskatchewan and Manitoba and three of the major crops in this region are barley, canola and spring wheat. The agricultural land on the Canadian Prairies has been divided into Census Agricultural Regions (CAR) by Statistics Canada. A clustering model was applied to the crop yield data to group the CARs for the development of forecast models. The normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) derived from the Moderate Resolution Imaging Spectro-radiometer (MODIS), NDVI derived from the Advanced Very High Resolution Radiometer (AVHRR) and several climate indices were considered as predictors for crop yields. A correlation analysis between crop yield and the time series of each potential predictor was performed to determine which variables showed the most forecasting potential and at what time during the growing season their values were most correlated to crop yield. Various combinations of MODIS-NDVI, MODIS-EVI and NOAA-NDVI were used to forecast the yield of barley, canola and spring wheat. Multiple linear regression as well as nonlinear Bayesian neural networks and model-based recursive partitioning forecast models were developed using the various sets of predictors. The models were trained using a cross-validation method and the forecast results of each model were evaluated by calculating the skill score from the mean absolute error, with 95% confidence intervals for the skill scores calculated using a bootstrap method. The results were compared in an effort to determine the optimal set of predictors and type of forecast model for each crop. iii Preface Currently no part of this thesis has been published, however a paper based on this thesis will be submitted for publication to peer-reviewed journals. This thesis contains research conducted by the candidate, Michael David Johnson, under the supervision of Dr. William Hsieh, Dr. Alex Cannon and Dr. Phil Austin. The satellite and agricultural data sets used in this study were provided by Dr. Andrew Davidson (Agriculture and Agri-Food Canada) and Dr. Fr?d?ric B?dard (Statistics Canada Agriculture Division). Monthly climate indices for the El Ni?o/Southern Oscillation, North Atlantic Oscillation, Arctic Oscillation, Pacific/North American Pattern and Antarctic Oscillation were downloaded from the website of the Climate Prediction Center, NOAA (http://www.cpc.ncep.noaa.gov/). Monthly values of the Pacific Decadal Oscillation were obtained from the Joint Institute for the Study of the Atmosphere and Ocean, University of Washington (http://jisao.washington.edu/pdo/PDO.latest). The monthly values for the North Pacific Gyre Oscillation were obtained from Di Lorenzo et al., (2008) (http://www.o3d.org/npgo/npgo.php). Figures 2.3 and 3.1 were reproduced from Mkhabela et al. (2011) with permission from the publisher. The development of the machine learning/statistical models and the analysis of the results were primarily the work of the candidate. The supervisory committee provided the original research topic, direction and support in the development of the models and critical feedback on the research methods and presentation of results. iv Table of Contents Abstract .......................................................................................................................................... ii Preface ........................................................................................................................................... iii Table of Contents ......................................................................................................................... iv List of Tables ................................................................................................................................ vi List of Figures .............................................................................................................................. vii Acknowledgements ........................................................................................................................x Chapter 1: Introduction ................................................................................................................1 1.1 Background ..................................................................................................................... 1 1.2 Research Objectives ........................................................................................................ 3 1.3 Organization of Thesis .................................................................................................... 3 Chapter 2: Literature Review .......................................................................................................4 2.1 Introduction ..................................................................................................................... 4 2.2 Remote Sensing .............................................................................................................. 4 2.3 Prairie Climate and Climate Indices ............................................................................... 6 2.4 Machine Learning Techniques ........................................................................................ 8 2.5 Crop Yield Forecasting ................................................................................................. 13 2.6 Summary ....................................................................................................................... 16 Chapter 3: Data ............................................................................................................................17 3.1 Study Area .................................................................................................................... 17 3.2 Data Sets ....................................................................................................................... 17 Chapter 4: Methods .....................................................................................................................20 v 4.1 Predictor Selection ........................................................................................................ 20 4.2 Clustering Analysis ....................................................................................................... 24 4.3 Multiple Linear Regression........................................................................................... 25 4.4 Bayesian Neural Networks ........................................................................................... 26 4.5 Model-Based Recursive Partitioning ............................................................................ 30 4.6 Model Evaluation .......................................................................................................... 31 4.7 Summary ....................................................................................................................... 32 Chapter 5: Results and Discussion .............................................................................................33 5.1 Barley ............................................................................................................................ 33 5.2 Canola ........................................................................................................................... 39 5.3 Spring Wheat ................................................................................................................ 45 5.4 Summary ....................................................................................................................... 49 Chapter 6: Conclusion .................................................................................................................54 6.1 Research Objectives ...................................................................................................... 54 6.2 Future Research ............................................................................................................ 55 Bibliography .................................................................................................................................57 vi List of Tables Table 4.1: Predictor combinations for MLR. ................................................................................ 26 Table 4.2: Predictor combinations for BNN. ................................................................................ 29 Table 4.3: Regression and partitioning variable combinations for MOB. .................................... 30 Table 5.1: Skill score of MLR forecast models for barley. .......................................................... 35 Table 5.2: Skill score of BNN forecast models for barley. ........................................................... 35 Table 5.3: Skill score of MOB forecast models for barley. .......................................................... 35 Table 5.4: Skill score of MLR forecast models for canola. .......................................................... 41 Table 5.5: Skill score of BNN forecast models for canola. .......................................................... 41 Table 5.6: Skill score of MOB forecast models for canola. .......................................................... 41 Table 5.7: Skill score of MLR forecast models for spring wheat. ................................................ 46 Table 5.8: Skill score of BNN forecast models for spring wheat. ................................................ 46 Table 5.9: Skill score of MOB forecast models for spring wheat. ............................................... 46 vii List of Figures Figure 2.1: The general structure of a neural network (Hsieh and Tang, 1998). ............................ 9 Figure 2.2: A diagram illustrating the problem of over-fitting. The dashed curve shows a good fit to noisy data (indicated by the squares), while the solid curve illustrates over-fitting, where the fit is perfect on the training data (squares), but is poor on the test data (circles) (Hsieh and Tang, 1998). ............................................................................................................................................ 10 Figure 2.3: Agro-Climatic Zones of the Canadian Prairies (Mkhabela et al., 2011. Reproduced with permission from publisher). .................................................................................................. 15 Figure 3.1: Canadian Prairie census agriculture regions (CAR) (Mkhabela et al., 2011. Reproduced with permission from publisher)............................................................................... 17 Figure 4.1: The evolution of correlation between vegetation indices and barley, canola and spring wheat. d) Correlation between MODIS-EVI and MODIS-NDVI throughout the growing season........................................................................................................................................................ 21 Figure 4.2: Correlation of barley, canola and spring wheat with climate indices. ....................... 22 Figure 4.3: Correlation of barley, canola and spring wheat with AAO for the years 1987-2011. 23 Figure 4.4: Dendrograms from hierarchical clustering of CARs. ................................................. 25 Figure 4.5: Illustration of the double cross-validation approach used to train and test the BNN models. The outer loop (CV1) is used to train and then test the forecast skill of the model. In CV1 the test data is shown in red and the training data in white. In this loop the test data segment is moved repeatedly in 1-year increments from the start to the end of the data record, so the forecast skill is tested on the entire data record. The CV2 loop is used to determine the optimal number of hidden neurons to use for the BNN model trained in CV1. The training data segment from CV1 is further segregated into validation data (in green) and a smaller training data segment (in white). In CV2 BNN models are trained using various numbers of hidden neurons and then evaluated on the validation data. The number of hidden neurons from the BNN model with the highest forecast skill on the validation data is selected to be used in the BNN in CV1. In the CV2 loop the validation data segment is moved repeatedly in 1-year increments from the start to the end of the training data record, so the forecast skill is evaluated on the entire training data record. .................................................................................................................................... 29 Figure 5.1: Grouping of CARs from hierarchical clustering model based on barley and by agro-climatic zones................................................................................................................................ 34 Figure 5.2: Skill score of MLR forecast models for barley with 95% confidence intervals. ....... 36 Figure 5.3: Skill score of BNN forecast models for barley with 95% confidence intervals. ........ 36 viii Figure 5.4: Skill score of MOB forecast models for barley with 95% confidence intervals. ....... 37 Figure 5.5: Skill scores of the best performing models for barley................................................ 38 Figure 5.6: Grouping of CARs from hierarchical clustering model based on canola and by agro-climatic zones................................................................................................................................ 40 Figure 5.7: Skill score of MLR forecast models for canola with 95% confidence intervals. ....... 42 Figure 5.8: Skill score of BNN forecast models for canola with 95% confidence intervals. ....... 42 Figure 5.9: Skill score of MOB forecast models for canola with 95% confidence intervals. ...... 43 Figure 5.10: Skill scores of the best performing models for canola. ............................................ 44 Figure 5.11: Grouping of CARs from hierarchical clustering model based on spring wheat and by agro-climatic zones. ...................................................................................................................... 45 Figure 5.12: Skill score of MLR forecast models for spring wheat with 95% confidence intervals........................................................................................................................................................ 47 Figure 5.13: Skill score of BNN forecast models for spring wheat with 95% confidence intervals........................................................................................................................................................ 47 Figure 5.14: Skill score of MOB forecast models for spring wheat with 95% confidence intervals. ........................................................................................................................................ 48 Figure 5.15: Skill scores of the best performing models for spring wheat. .................................. 49 Figure 5.16: Scatter plots illustrating the relationship between MODIS-NDVI and barley yield for the entire Canadian Prairie region. The data points from the year 2002 are marked by a black circle. The CAR means for both yield and NDVI have been removed from the data. The climatology predictions are shown as well as the forecasts for a) MLR, b) BNN and c) MOB models using only MODIS-NDVI as a predictor. ........................................................................ 50 Figure 5.17: Scatter plots illustrating the relationship between MODIS-NDVI and canola yield for the entire Canadian Prairie region. The data points from the year 2002 are marked by a black circle. The CAR means for both yield and NDVI have been removed from the data. The climatology predictions are shown as well as the forecasts for a) MLR, b) BNN and c) MOB models using only MODIS-NDVI as a predictor. ........................................................................ 51 Figure 5.18: Scatter plots illustrating the relationship between MODIS-NDVI and spring wheat yield for the entire Canadian Prairie region. The data points from the year 2002 are marked by a black circle. The CAR means for both yield and NDVI have been removed from the data. The climatology predictions are shown as well as the forecasts for a) MLR, b) BNN and c) MOB models using only MODIS-NDVI as a predictor. ........................................................................ 52 ix List of Abbreviations Abbreviation Meaning AAO Antarctic Oscillation AB Alberta ACZ Agro-climatic zones ANN Artificial neural networks AO Arctic Oscillation AVHRR Advanced Very High Resolution Radiometer BNN Bayesian neural network CAR Census agricultural region CART Classification and regression tree DRI Drought research initiative ENSO El Ni?o/Southern Oscillation EVI Enhanced vegetation index MAE Mean absolute error MB Manitoba MLR Multiple linear regression MOB Model-based recursive partitioning MODIS Moderate Resolution Imaging Spectro-radiometer NAO North Atlantic Oscillation NASA National Aeronautics and Space Administration NDVI Normalized difference vegetation index NIR Near infrared NPGO North Pacific Gyre Oscillation PAR Photosynthetically active radiation PDO Pacific Decadal Oscillation PNA Pacific/North American Pattern SK Saskatchewan SRPG Soil rating for plant growth SS Skill score SST Sea surface temperature x Acknowledgements I would like to thank Prof. William Hsieh for his never-ending guidance and support. The time that was put into this research could not have been done without his coordination and exceptional knowledge of Atmospheric science and machine learning methods. He is who made it possible for me to not only do my schooling and research at UBC but to also work at Earthrisk Technologies as an intern where I was able to put my skills to productive use and acquire added knowledge. William Hsieh?s constant advisement and considerate attitude is what allowed the completion of this research without any external pressure or tension. I would also like to thank Dr. Alex Cannon for his commitment to learning and advisement for this research. I would also like to acknowledge and thank the members of the UBC climate prediction group; Carlos Gaitan-Ospina, Aranildo Rodrigues and Andrew Snauffer. Not only were they my fellow students who helped and supported, but they were also my friends. I would like to thank Prof. Phil Austin and Prof. Susan Allen, who taught me during my research, and this knowledge was able to contribute to my research. Additionally I would like to thank Dr. Andrew Davidson and Dr. Fr?d?ric B?dard for providing data necessary for this research. Lastly I would like to thank Prof. Jacqueline Binyamin who was my professor at the University of Winnipeg and is who encouraged me to go to grad school and obtain my Master?s degree. Without her constant support and encouragement, I would have never continued on to UBC to get my masters. 1 Chapter 1: Introduction 1.1 Background The major grain crops from the Canadian Prairies include barley, canola and spring wheat. The production of grain crops plays an important role in the economy of the Canadian Prairies and crop yield is a key element for rural development and an indicator for national food security (Li et al., 2007). Early and reliable crop yield forecasts over large areas would help policy makers and grain marketing agencies in planning for exports. Having a better estimate of the crop yield would help the grain industries to be more organized during harvest time. They will be able to set prices based on the expected yield and organize transportation for the crops much earlier than before. The development of forecast models from remotely sensed data is ideal because remote sensing techniques have the potential to provide quantitative and timely information on agricultural crops over large areas and many different methods have been developed to estimate crop conditions (Li et al., 2007). Crop yield can be influenced by many factors such as precipitation and temperature during growing season (Guo, 2012), soil conditions (Alvarez, 2009), disease, and anthropogenic factors like irrigation or fertilizers (Prasad et al., 2006). Climatic data are often readily available, but some of the other factors can be difficult or impossible to quantify. Remotely sensed vegetation indices can directly measure crop growth and account for many of these factors. Normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) derived from the Moderate Resolution Imaging Spectro-radiometer (MODIS) on the Terra and Aqua satellites have been used to monitor crop conditions in many locations including Canada. Mkhabela at al. (2011) showed that MODIS-NDVI could be used to effectively predict crop yields on the Canadian Prairies one to two months before harvest. In the past crop forecasts relied on vegetation indices developed from the Advanced Very High Resolution Radiometer (AVHRR) which has a resolution of only 1 km. The MODIS sensor has a higher spatial resolution (up to 250m) than AVHRR and could lead to more accurate forecasts. Prasad et al. (2006) developed a linear model that used AVHRR-NDVI along with soil moisture, surface temperature and rainfall to predict corn and soybean yield for Iowa and it is suggested that higher resolution data could 2 further improve the model. Precipitation and temperature data sets can be difficult to work with because they are often only measured at certain sites. Interpolation or spatial averaging techniques would be required to apply such data sets to crop yield prediction problems (Vincente-Serrano et al., 2003). It is possible that climate indices could be included as crop yield predictors to help account for the large-scale meteorological conditions. Hsieh et al. (1999) found a nonlinear relationship between the Pacific sea surface temperature and Canadian Prairie wheat yield and it is possible that more teleconnections such as this exist. Forecast models with climate indices included as predictors could have higher skill than models using vegetation indices alone. The Canadian Prairies cover a very large geographical area that has been divided up into Census Agricultural Regions (CARs) by Statistics Canada. The conditions over this large area can vary greatly and developing a single forecast model for the entire region may not be appropriate. Mkhabela at al. (2011) categorized the CARs based on soil type into three distinct agro-climatic zones including sub-humid, semi-arid and arid then a forecast model was developed for each agro-climatic zone. A potential problem with this method is that some CARs could have multiple soil types and this would make it difficult to classify such a CAR into the appropriate agro-climatic zone. It is possible that a clustering algorithm could be applied to the yield data to better group the CARs for the purpose of crop yield forecasting. Most current crop yield forecasting uses linear statistical models, which cannot capture nonlinear relations in the data. The cropland ecosystem is complex and many of the processes involved are non-linear. This causes problems with traditional statistical models, such as multiple linear regression (Jiang et al., 2004). The results of Mkhabela at al. (2011) show a possible nonlinear relationship between MODIS-NDVI and the yield of barley, canola, field peas and spring wheat. Nonlinear neural networks have been successfully used to predict crop yield using remotely sensed vegetation indices. Li et al. (2007) and Kaul et al. (2005) both used artificial neural networks to forecast corn and soybean yields in the USA. This study attempts to use fully nonlinear models from the field of machine learning, such as model-based recursive partitioning (MOB) and Bayesian neural networks (BNN) to forecast crop yields. These machine learning 3 methods have been used successfully in the past in speech and handwriting recognition, computer vision, robotics, etc. 1.2 Research Objectives This study is focused on improving crop yield forecasting on the Canadian Prairies and therefore the research goals are to: i) Compare the effectiveness of several vegetation indices derived from satellite data (MODIS-NDVI, MODIS-EVI and AVHRR-NDVI) in forecasting crop yields on the Canadian Prairies. ii) Compare linear and nonlinear statistical/machine learning models to determine which is most appropriate for forecasting crop yields on the Canadian Prairies. iii) Determine if the introduction of climate indices as predictors will improve the skill of crop yield forecasting models for the Canadian Prairies. The ultimate goal is that in completing the research objectives an effective crop yield forecast model for the Canadian Prairies will be developed. 1.3 Organization of Thesis Chapter 2 provides a literature review covering topics related to crop yield forecasting, remote sensing and machine learning techniques. Background theory on various machine learning and remote sensing topics will be covered. The results from the reviewed crop forecasting studies as well as the modeling techniques used will be discussed with emphasis on how they can be applied to this research. Chapter 3 describes in detail the study area and data sets obtained for this study. Chapter 4 outlines the methods used to conduct this research and describes the developed forecast models and evaluation methods. In Chapter 5 the results from all developed forecast models for each crop are discussed in detail. The thesis concludes with Chapter 6 where the original research objectives are addressed and recommendations for future research are made. 4 Chapter 2: Literature Review 2.1 Introduction Agriculture is a major industry and it is vital to the economy of many countries. Due to the importance of agriculture there has been a wide variety of studies conducted in this field. The purpose of this literature review is to justify the research objectives of this study in light of previous work by investigating past crop yield prediction studies and determining where further research is needed. Literatures related to remote sensing, prairie climate and the various types of machine learning methods used in this study are reviewed. Machine learning theory and past applications are examined to show that these methods will likely perform well in a crop forecasting application. The influence of this material on the field of crop yield forecasting will be discussed and how it can be applied to this specific project is explained. 2.2 Remote Sensing Remotely sensed vegetation indices from the MODIS sensor are an essential aspect of this study. In the past vegetation indices were a product of the AVHRR, which is a sensor aboard the National Oceanic and Atmospheric Administration (NOAA) weather satellites. Now vegetation indices are usually a product of the MODIS sensor which is aboard the National Aeronautics and Space Administration (NASA) satellites, Terra and Aqua. The MODIS sensor has been gathering data for about 12 years and captures more spectral bands at a higher resolution than the AVHRR. Huete et al. (2002) discusses two types of MODIS vegetation indices and explains the theory behind them. The first is the normalized difference vegetation index (NDVI). The theory of NDVI is based on the fact that live plants absorb solar radiation in the photosynthetically active radiation (PAR) spectral region, which they use as a source of energy in the process of photosynthesis. Plants have also evolved to scatter solar radiation in the near-infrared spectral region and therefore live plants appear relatively dark in the PAR and relatively bright in the near-infrared spectral region. NDVI is the normalized ratio of the near infrared (NIR) and red MODIS bands: , (2.1) 5 where ?NIR and ?red are the surface bidirectional reflectance factors for their respective MODIS bands (Huete et al., 2002). The second vegetation index is the enhanced vegetation index (EVI) which was developed to optimize the vegetation signal with improved sensitivity in high biomass regions and improved vegetation monitoring through a de-coupling of the canopy background signal and a reduction in atmospheric influences (Huete et al., 2002): , (2.2) where ? are atmospherically corrected or partially atmospherically corrected (Rayleigh and ozone absorption) surface reflectances, L is the canopy background adjustment that addresses non-linear, differential NIR and red radiant transfer through a canopy. C1 and C2 are the coefficients of the aerosol resistance term, which uses the blue band to correct for aerosol influences in the red band (Huete et al., 2002). Both EVI and NDVI are unit-less values ranging from -1 to 1, with healthy green vegetation having high positive values. Rocks and bare soil have values near zero and surfaces such as water, snow or clouds usually have negative values. The performance and validity of the MODIS vegetation indices have been demonstrated by comparison of the vegetation indices to measured test sites that represented semi-arid grass/shrub, savannah, and tropical biomes (Huete et al., 2002). Chen et al. (2006) assessed the suitability of using MODIS vegetation indices for agricultural monitoring and found that the results closely resembled actual crop conditions. It has also been shown that the vegetation indices derived from MODIS data are more effective at monitoring vegetation than the vegetation indices from AVHRR (Huete et al., 2002). Since MODIS NDVI and EVI have been shown to be effective tools for assessing vegetation condition it is valid to use them to predict crop yield. However it is not clear whether NDVI or EVI will be a better predictor of crop yield in the case of this study. Li et al. (2010) used MODIS NDVI and EVI to monitor several different vegetation cover types in the northern Hebei Province of China and compared the results. In general NDVI was found to outperform EVI in this region, but the results could be different for the Canadian Prairies. 6 2.3 Prairie Climate and Climate Indices Climate indices are values that represent the phases of various teleconnections or oscillations in the climate system. The phase of these teleconnections can affect weather in Canada which will in turn affect crop yield. Many crop yield models from previous studies use meteorological variables and not just vegetation indices (Prasad et al. 2006; Kaul et al., 2005) and therefore using climate indices as additional predictors could improve the skill of the forecast models developed in this study. Historical climate index data is available from the National Oceanic and Atmospheric Administration (NOAA) and for this study the indices for the El Ni?o/Southern Oscillation, North Atlantic Oscillation, Arctic Oscillation, Pacific/North American Pattern, Pacific Decadal Oscillation, North Pacific Gyre Oscillation and the Antarctic Oscillation were investigated. The El Ni?o/Southern Oscillation (ENSO) is a large-scale ocean-atmosphere oscillation in the tropical Pacific that influences climatic conditions around the globe including Canada. ENSO consists of a warm and a cold phase called El Ni?o and La Ni?a, respectively, and they occur on average every two to seven years (Bonsal and Shabbar, 2011). The North Atlantic Oscillation (NAO) is the dominant mode of atmospheric variability over the North Atlantic and is most pronounced during the cold season. The NAO index is a measure of the difference in the strength of the Icelandic Low and the Azores High and has been shown to affect climate over Europe and eastern Canada (Wettstein and Mearns, 2002). The NAO is highly correlated to the Arctic Oscillation (AO) which represents atmospheric circulation variability over the extra-tropical Northern Hemisphere where sea-level pressures over the polar regions varies in opposition with that over middle latitudes (about 45?N) (Thompson and Wallace, 1998). The Pacific/North American (PNA) pattern describes the variation of atmospheric circulation patterns over the Pacific Ocean and North America (Bonsal and Shabbar, 2011). The Pacific Decadal Oscillation (PDO) is a measure of the sea surface temperatures in the North Pacific and is spatially somewhat similar to ENSO (but relatively stronger in the extratropics) with cycles of 20 to 30 years (Bonsal and Shabbar, 2011). The North Pacific Gyre Oscillation (NPGO) is the second dominant mode of sea level pressure variability in the North Pacific (Di Lorenzo et al., 2010). 7 Finally, the Antarctic Oscillation (AAO) is a low frequency mode of atmospheric variability in southern hemisphere (Song et al., 2009). All of the climate indices investigated in this study were considered as potential predictors of crop yields because they can all have an effect on the Canadian Prairie climate. The strongest relationship that occurs between large scale teleconnections and the Canadian climate occurs with cold season temperature. The positive phases of ENSO, PNA and PDO are associated with above average winter temperatures across western Canada, while the negative phase are associated with colder than normal winters (Bonsal and Shabbar, 2011). Although the relationship is not as strong as with ENSO, colder than normal winters are associated with positive NAO and AO values and vice versa (Bonsal and Shabbar, 2011). Large scale teleconnections are also related to cold season precipitation in Canada. For example, El Ni?o events are associated with a pattern of below normal precipitation stretching from southern British Columbia, through the Prairies, and into the Great Lakes region and La Ni?a events lead to an opposite response (Bonsal and Shabbar, 2011). Although not as strong as in the winter, relationships do exist between teleconnections and the summer climate in Canada. Shabbar and Skinner (2004) found that El Ni?o events led to a summer moisture deficit in western Canada, while La Ni?a events produced an abundance of moisture, mainly in extreme western Canada. Although it may seem unlikely, even the Antarctic Oscillation (AAO) could have an effect on crop yields on the Canadian Prairies. It has been shown that the AAO accounts for a large fraction of extra-tropical variability in the Southern Hemisphere. It has also been found that the AAO can have an influence on the Southern Hemisphere subtropical/tropical latitudes and can even extend into the subtropics of the Northern Hemisphere as well. It has been well-proven that the Northern Hemisphere winter is very sensitive to tropical forcing. This means that the AAO could impact the Northern Hemisphere about 25-40 days later through the tropics in the boreal winter (Song et al., 2009). This does not imply that the AAO index will have any predictive ability for crop yields, but it shows that further investigation is justified. Hsieh et al. (1999) showed the potential for using teleconnections to predict wheat yield on the Canadian Prairies. El Ni?o?Southern Oscillation (ENSO) was shown to have a nonlinear relationship to wheat yields. Further research into the potential for using climate indices to predict crop yields is justified 8 because it is possible that such a relationship exists between ENSO or another climate index and other crops on the Canadian Prairies. Droughts are an important part of the Prairie ecosystem and they can have an extreme effect on crop yield and therefore it is vital for crop yield forecast models to be able to accurately capture these effects. Better prediction of droughts is extremely important because agriculture is a major part of the Canadian economy and a drought can have devastating financial effects. The Drought Research Initiative (DRI) was established to study the recent drought on the Canadian Prairies that began in 1999 and ended in 2004/05 (Stewart et al., 2011). The DRI?s main goals were to better predict droughts in Canada, and to assess whether there will be a drying of the continental interior in the future. The DRI was able to relate global winter sea surface temperature (SST) patterns to droughts and extreme wet seasons on the Canadian Prairies. Some of these patterns closely resemble ENSO and the Pacific Decadal Oscillation (PDO) and since there is as much as a six-month lag time between the winter and growing season, this means that these SST patterns could be used to predict droughts on the Prairies (Shabbar et al., 2011). Since a drought would directly affect vegetation it is likely that the MODIS vegetation indices will be able to account for drought conditions, but based on the finding of Shabbar et al. (2011) there is reason to believe that ENSO, PDO and possibly other climate indices may help improve forecasts in drought years, especially at longer lead times. 2.4 Machine Learning Techniques Most current crop yield forecasting on the Canadian Prairies uses linear or other simple regression statistical models (Mkhabela et al., 2011). Artificial neural networks are nonlinear models capable of capturing nonlinear relationships in data. Neural networks were originally developed from investigations into human brain function and it is an adaptive system that changes as it learns (Hsieh and Tang, 1998). There are two main types of machine learning: supervised learning and unsupervised learning. In supervised learning the model has a data set with the predictors and corresponding dependent variables given and the model is trained on this data to develop a statistical relationship between the variables. Examples of this are regression (predicts a continuous value output) and classification (discrete value output e.g. 1 or 0). Crop 9 yield prediction problems are another example of supervised learning. A forecast model would be given training data containing both the predictors and the crop yields from previous years. If for example a neural network is used, a nonlinear relationship will be determined from this data. Many types of neural network models exist, but in general the basic structure of a neural network is as shown in Figure 2.1. Figure 2.1: The general structure of a neural network (Hsieh and Tang, 1998). The input xi are mapped to the hidden layer of neurons hj by: ? , (2.3) and then onto the output yk by: ? ? ? , (2.4) where f and g are activation functions, and ? are weight parameter matrices and and ? are offset parameters and their optimal values are learned by the model during training (Hsieh 10 and Tang, 1998). In this study f is the hyperbolic tangent function (Eq. 2.5) and g is the identity function (Eq. 2.6), i.e. , (2.5) , (2.6) The presence of a nonlinear activation function at the hidden layer is what allows the neural network to have nonlinear modeling capabilities (Hsieh and Tang, 1998). One of the main challenges in developing a neural network model is how to address the problem of over-fitting. An over-fitted neural network could fit the data very well during training, but produce poor forecast results during testing (Hsieh and Tang, 1998). Over-fitting occurs when a model fits the noise in the data and it will not generalize well to new data sets as shown in Figure 2.2. Figure 2.2: A diagram illustrating the problem of over-fitting. The dashed curve shows a good fit to noisy data (indicated by the squares), while the solid curve illustrates over-fitting, where the fit is perfect on the training data (squares), but is poor on the test data (circles) (Hsieh and Tang, 1998). 11 Typically regularization (i.e. the use of weight penalty to avoid over-fitting from the use of large weights) is used to prevent over-fitting. This requires some of the training data to be used as validation data to determine the optimal regularization parameter to prevent over-fitting. Machine learning is a quickly evolving field and many types of neural networks and different techniques for optimizing their performance have been developed. For example, Bayesian neural networks (BNN) have become quite popular. A BNN can address the problem of over-fitting and does not require validation data to determine the regularization parameter. Using Bayesian probability theory the probability distribution of the network weights given the training data is calculated. For a data set, rather than selecting the single optimal weight vector, the probability distribution of the weights from the given data set is determined. Bayesian neural networks have been compared to standard neural networks and multiple linear regression techniques and have been found to outperform linear models, and in most cases, outperform standard neural network models (Bhat and Prosper, 2005; Rasmussen, 1996). Bayesian neural network packages are available in many popular code languages including MATLAB and due to the fact that they do not require validation data Bayesian neural networks are especially useful in situations where there is a limited amount of data available for training. Neural networks are not the only machine learning method that use supervised learning. An alternative and often more interpretable model is a classification and regression tree (CART). Classification trees are used when the predictand takes a discrete number of values whereas regression trees are used when the predictand is continuous. Regression trees can represent non-linear functions in a piecewise, discretized manner (Gardner and Dorling, 2000). The output from a CART model consists of a binary tree structure and associated with each split in the tree is a rule involving one or more predictor variables (Gardner and Dorling, 2000). The model-based recursive partitioning (MOB) package available in R is a nonlinear tree based method. The MOB model partitions the data like a standard CART model, but then a linear model is developed for the data in each partition. The basic steps are: (1) fit a parametric model to a data set, (2) test for parameter instability over a set of partitioning variables, (3) if there is some overall parameter instability, split the model with respect to the variable associated with the highest instability, (4) repeat the procedure in each of the child nodes (Zeileis et al., 2008). 12 Zeileis et al. (2008) successfully used a MOB model in several nonlinear applications in the field of healthcare and real-estate. Clustering is a common machine learning technique that can have applications to crop yield forecasting because it can be useful for grouping geographical areas or data sets before a forecast model is applied. Clustering is an example of unsupervised learning where only the input data are provided and the model must find patterns or structure in the data on its own. A common clustering method is k-means clustering which starts by guessing initial random locations for the cluster centers and then the following two steps are repeated until convergence: i) For each data point the closest cluster center is determined based on Euclidean distance and the data point is assigned to belong to that cluster. ii) For each cluster, the cluster center is changed to be the mean position of all the data points belonging to that cluster. Often the optimal number of clusters is not known and it can be difficult to determine. There is a method of clustering analysis called hierarchical clustering (Wilks, 2006) which builds a hierarchy of clusters. This method begins with each data point in its own cluster and then merges points as the number of clusters decreases until there is just a single cluster. The results of hierarchical clustering are usually presented in a dendrogram which is basically a tree diagram that shows how the data is merged into various numbers of clusters with the vertical axis representing the distance between merged clusters. At the start of the process when similar clusters are merged the distances are small and usually increase slightly as the merging continues. Later in the process there will be a small number of clusters separated by large distances. If a point in the process can be discerned where the distances between merged clusters changes drastically, the process can be stopped before these distances become too large (Wilks, 2006). 13 2.5 Crop Yield Forecasting There have been many studies done about crop yield prediction in various parts of the world and in this section relevant crop forecasting studies will be reviewed. The relationship between the reviewed studies and this research will be discussed with particular emphasis on the types of models developed and the predictor variables used. Many crop prediction studies use only linear techniques, e.g. Prasad et al. (2006) which focused on predicting corn and soybean yield for Iowa using AVHRR-NDVI, soil moisture, surface temperature and rainfall data. Piece-wise linear regression models were developed and found to fit both the corn and soybean data well. However, the models were not tested on independent data so no actual forecasts were made. Even still, the fact that a piece-wise regression model fit the yield data well is suggestive that a MOB model may also perform well in crop forecasting applications as the output is also the results of several linear fits. It is suggested that other factors such as pest, diseases and human activities may cause variations in crop yield that may limit the forecast ability of this model. However, the use of NDVI in the model may account for the loss of crops due to these factors if they directly affect the vegetation and it is also suggested that higher resolution data could further improve this model, which gives motivation for using MODIS-NDVI as a predictor in this study. Li et al. (2007) is similar to Prasad et al. (2006) except that nonlinear techniques are included and higher resolution NDVI was used. Li et al. predicts corn and soybean yield using both AVHRR-NDVI and MODIS-NDVI. Neural network and MLR models at the county level were developed for the ?corn belt? of the USA which includes the states of North Dakota, South Dakota, Nebraska, Kansas, Minnesota, Iowa, Missouri, Michigan, Wisconsin, Illinois, Indiana and Ohio. It was concluded that the neural network models using the higher resolution MODIS-NDVI values as predictors were the most effective at predicting corn and soybean yields. The only potential problem with the MODIS-NDVI is that there is not a long historical data record. This means that there will be less training data available and MODIS-NDVI might not be able to be incorporated into models using meteorological variables which would have a much longer historical record. 14 Kaul et al. (2005) is also focused on predicting the yield of corn and soybeans. In this case the study takes place in Maryland and both multiple linear regression and artificial neural network models are developed. Historical yield, rainfall data and soil rating for plant growth (SRPG) values were used as predictors and although no remotely sensed variables were used as predictors in Kaul et al. (2005) there are still aspects that relate to this current research. The predictive abilities of each type of model were compared at various spatial scales, the ANN models were found to outperform the MLR models in every case and the yield predictions improved as the area being modeled decreased. It was stated that larger areas will have larger variability of cropping conditions and the models are only applicable to the conditions for which they were developed (Kaul et al., 2005). Kaul et al. (2005) is a good example of the advantage of using nonlinear modeling techniques for agricultural applications and it also shows that developing a single crop forecast model for a large area is not the best approach. In most cases dividing the study area into smaller regions, such as CARs in the case as the Canadian Prairies, should lead to better crop yield prediction. Jiang et al. (2004) predicted winter wheat in China and the study is similar to Kaul et al. (2005) in that MLR and ANN models are compared and the study area is divided up into smaller regions to improve crop yield prediction. The difference is that Jiang et al. (2004) focused on using remotely sensed data to predict the yield of winter wheat in China. In addition to vegetation indices water stress and surface temperature indices were also included as predictors. The study found that it was possible to predict winter wheat yields using the ANN models developed and that it performed better than the traditional MLR approach. The success of the model was attributed to the non-linear capabilities of the ANN, the use of remotely sensed data, and the precise division of the study area. Mkhablea et al. (2011) developed various types of regression models to predict the yields of barley, canola, field peas and spring wheat on the Canadian Prairies using MODIS-NDVI as a predictor and is therefore closely related to this research. The goal of Mkhablea et al. (2011) was to determine if it is possible to use MODIS-NDVI to forecast crop yield on the Canadian Prairies and to determine the best time during the growing season to make the forecast. The study area 15 consisted of Manitoba, Saskatchewan and Alberta and it was divided up into census agricultural regions (CAR). Most of the Canadian Prairies consists of brown, dark brown and black and gray chernozemic soils. The soils are generally associated with the agro-climatic zones, i.e., the black and gray soils are dominant in the sub-humid zone, the dark brown soils in the semi-arid zone and the brown soils in the arid zone (Mkhabela et al., 2011). To account for the different soil types the CARs were then divided up into the 3 agro-climatic zones based on their annual precipitation and separate models were developed for each agro-climatic zone. Figure 2.3: Agro-Climatic Zones of the Canadian Prairies (Mkhabela et al., 2011. Reproduced with permission from publisher). To determine the best time during the growing season to make the forecast a correlation analysis between crop yields and weekly composite MODIS-NDVI values was performed. The times of the year that had the highest correlation coefficients were determined to be the best times to make a crop yield prediction. The best time to make the prediction varied slightly for each region and crop, but it was determined that accurate forecasts could be made as early as late June to late July. The MODIS-NDVI values from these times of the year were used as predictor variables in a regression analysis to predict crop yield and it was concluded that MODIS-NDVI is effective at forecasting the yields of barley, canola, field peas and spring wheat on the Canadian Prairies. A 16 power function regression was found to fit the yield data best which suggests that a nonlinear approach like a neural network could improve crop forecasting in this region. 2.6 Summary Studies from the fields of remote sensing, climate and teleconnections, machine learning and crop yield forecasting were reviewed and it is clear that much effort has been put into agricultural research and crop forecasting. The material covered in this literature review strongly justifies the approach taken in this study and shows the motivation behind the research objectives. It is clear that MODIS vegetation indices have the potential of being effective predictors of barley, canola and spring wheat on the Canadian Prairies. The importance of dividing large study areas up into smaller regions is clear because in most cases crop forecast models were found to have higher skill when developed for smaller areas. The advantage of nonlinear models was shown and the applications of BNN and MOB models in this research should provide some advantage over linear methods based on the results of the previous studies covered in this review. 17 Chapter 3: Data 3.1 Study Area The Canadian Prairies consists of the provinces of Alberta (AB), Saskatchewan (SK) and Manitoba (MB) and collectively these provinces have about 30 million ha of crop land. The Prairies extend northward from approximately 49?N to 54?N latitudes and westward from approximately 96?W to 114?W longitudes (Mkhabela at al., 2011). Statistics Canada has divided the agricultural land in Canada into Census Agricultural Regions (CARs). There are 40 CARs in total with 8 in AB, 20 in SK and 12 in MB. Figure 3.1: Canadian Prairie census agriculture regions (CAR) (Mkhabela et al., 2011. Reproduced with permission from publisher). 3.2 Data Sets This study requires several data sets including crop yield data, remotely sensed vegetation indices and climate indices. The Agricultural Division of Statistics Canada collects detailed crop 18 yield data every year for all crops across Canada. To account for the fact that not all agriculture land is seeded each year, crop yield is defined as the amount of crop harvested per seeded area (kg/ha). The yield of barley, canola and spring wheat by CAR (kg/ha) was obtained for the period of 2000-2011. Some of the yield data was classified as confidential or missing, but for both barley and canola there was less than 0.9% of the data missing and only 4.8% missing for spring wheat so the missing data points were simply ignored. The mean yield for each CAR was removed from all yield data belonging to that CAR to account for possible difference between the CARs that could affect crop yield such as elevation, climate, etc. The vegetation index data by CAR for the period 2000-2011 was obtained from Statistics Canada. The data set included EVI and NDVI derived from the MODIS sensor as well as NDVI derived from the AVHRR. The NDVI data was made up of a time series of 7-day composite values that span the length of the growing season (approximately April to September) for each CAR. The EVI data was similar, except it contained 16-day composite values. An agricultural land mask was applied to all the vegetation index data sets to ensure that the values were not affected by pixels not containing agricultural land, such as residential areas or forests. As with the yield data the mean value of NDVI and EVI for each CAR was removed to account for possible differences between CARs. The NDVI and EVI time series were also smoothed to reduce noise in the data. A four point moving average was applied to the NDVI data and a two point moving average was applied to the EVI data making them both essentially monthly average values. There were several climate indices considered as potential predictors of crop yields in this study and the historical data records for each of them were available for download online. Monthly climate indices for the El Ni?o/Southern Oscillation, North Atlantic Oscillation, Arctic Oscillation, Pacific/North American Pattern and Antarctic Oscillation were downloaded from the website of the Climate Prediction Center, NOAA (http://www.cpc.ncep.noaa.gov/). Monthly values of the Pacific Decadal Oscillation were obtained from the Joint Institute for the Study of the Atmosphere and Ocean, University of Washington (http://jisao.washington.edu/pdo/PDO.latest). The monthly values for the North Pacific Gyre 19 Oscillation were obtained from Di Lorenzo et al., (2008) (http://www.o3d.org/npgo/npgo.php). The climate index data was smoothed by a 3 month moving average. 20 Chapter 4: Methods 4.1 Predictor Selection As previously mentioned, the vegetation index data sets are in the form of a time series of 7-day composite values for NDVI and 16-day composite values for EVI that span the growing season and as in Mkhabela at al. (2011) it was necessary to determine the best time to make an accurate yield prediction. However, the lead time must still be large enough for the forecast to be useful and therefore to ensure a significant lead time only data from before August (before weeks 31-35) were considered for predictors. Figure 4.1 shows the preliminary correlation analysis between the yield of each crop and the smoothed MODIS-EVI, MODIS-NDVI and NOAA-NDVI time series for the growing season. The product term MODIS-EVI x MODIS NDVI was considered as input to the multiple linear regression models and the correlation between this cross term and each crop is also shown. The highest correlation between EVI or NDVI and the yield of barley, canola or spring wheat was found to be during the month of July (weeks 27 - 30) which is comparable to the results of Mkhabela at al. (2011). Hence, the predictors used in the forecast models were the values from July. The correlation between MODIS-EVI and MODIS-NDVI is shown in Figure 4.1d, where the correlation drops around weeks 27 - 30. Thus, at this point in the growing season EVI and NDVI may be able to provide different information to the model. 21 Figure 4.1: The evolution of correlation between vegetation indices and barley, canola and spring wheat. d) Correlation between MODIS-EVI and MODIS-NDVI throughout the growing season. Monthly averages of the climate indices for ENSO, NAO, AO, PNA, PDO, NPGO and AAO were obtained from the NOAA website. These data sets were time series for each climate index that spanned from December of the previous year to the end of the year. Due to the large oscillation periods and scales of these climate indices the data was smoothed to a 3 month moving average and a similar correlation analysis was performed for the years 2000-2011 between each climate index and the yields of barley, canola and spring wheat (Fig. 4.2). 22 Figure 4.2: Correlation of barley, canola and spring wheat with climate indices. 23 Based on the correlation analysis it was clear that the climate indices did not show much potential for predicting crop yields with the exception of AAO. Although it is possible that there is a relationship between AAO and crop yields on the Canadian Prairies it is unlikely. To further investigate this relationship the crop yield and AAO data records were extended and the correlation analysis was repeated for the longer period of 1987-2011 (Fig. 4.3). Figure 4.3: Correlation of barley, canola and spring wheat with AAO for the years 1987-2011. The previous correlation seen between crop yields and AAO was clearly a spurious relationship because when the data record is extended the high correlation does not exist. It was still possible that there could be a nonlinear relationship between some of the climate indices and crop yield and to confirm that none of the climate indices would be useful as predictor variables some preliminary nonlinear models were run using the various climate indices as predictors and the models were found to have little to no skill. J F M A M J J A-0.0500.050.10.150.20.25AAO Correlation With Crop Yieldmonthcorrelation BarleyCanolaSpring Wheat24 4.2 Clustering Analysis The Canadian Prairies cover a large area and the geography and climate can change from one location to the next. It does not make sense to develop one crop yield forecast model for all the Canadian Prairies and developing one model for each CAR will likely be unreliable due to the small number of data points used to train the model. The best approach is to group CARs which behave similarly with respect to yearly crop yield. Mkhabela et al. (2011) grouped the CARs based on soil type into three agro-climatic zones: sub-humid, semi-arid and arid because it was thought that crops in each region would behave similarly and this would improve the accuracy of the forecasts. In this study a clustering model for each crop was developed in an attempt to determine if there is a way to group the CARs that will lead to higher forecast skills. A clustering model groups the data into a number of subsets called clusters, so that the data in each cluster are more closely related to each other than data from other clusters. The clustering model grouped the CARs based on their crop yield time series over the period 2000-2011 and a forecast model was then developed for each cluster. A common problem with clustering analysis that was encountered in this study is that the optimal number of clusters is not known. To address this issue a method of clustering analysis called hierarchical clustering was used, which starts with each data point in its own cluster and then the data is merged as the number of clusters decreases. The Ward algorithm for computing the distance between clusters was selected (Murtagh and Legendre, 2011). 25 Figure 4.4: Dendrograms from hierarchical clustering of CARs. After reviewing the dendrograms (Figure 4.4) it was estimated that the optimal number of clusters was in the range of 1 to 6 because for each crop the distance between cluster centers appears to noticeably increase at around 6 clusters. Each type of model was developed with data divided up into 1 to 6 clusters and the results were compared to determine the optimal number of clusters. 4.3 Multiple Linear Regression Multiple linear regression forecast models using various combinations of predictors were developed in MATLAB. During training a multiple linear regression model finds a linear relationship between the target data and multiple predictor variables xi (i 1,?,m): , (4.1) 26 where is the model output, ai (i 0, ?, m) are the regression parameters determined during training. The model was trained using a leave one year out cross-validation scheme. There are 12 years of data spanning from 2000-2011, the model was trained on 11 years of data with 1 year left unseen by the model to be used as independent test data. This process was repeated until every year has been used as test data. The MLR model was developed with several different combinations of predictors which are shown in Table 4.1, where for MLR model 5, the product of MODIS-EVI and MODIS-NDVI served as the third predictor. Cross terms such as this are used in regression equations to capture the interaction between two variables and the addition of this term may allow the model to gain additional information from the vegetation indices. Table 4.1: Predictor combinations for MLR. MLR model Predictor 1 Predictor 2 Predictor 3 1 EVI 2 MODIS-NDVI 3 AVHRR-NDVI 4 EVI MODIS-NDVI 5 EVI MODIS-NDVI EVI x MODIS-NDVI 4.4 Bayesian Neural Networks Bayesian neural network forecast models were developed in MATLAB using the Neural Network Toolbox. A neural network with a single output is trained using a data set , , where x are the predictors and the predictand, by adjusting network parameters or weights (including offset or bias parameters) w so as to minimize an objective function (also called cost or loss or error function), such as, , (4.2) ? ? , ? , (4.3) 27 ? ? , (4.4) where is the number of training data points, is the network output, is the target data, is an error function measuring how close the network output is to the target, is the weight penalty term to prevent over-fitting (i.e. fitting to the noise in the data) and ? and ? are hyperparameters as they control the distributions of other parameters, i.e. the weights w of the network (van Hinsbergen, 2009). A BNN treats the network weight parameters as random variables with an assumed prior distribution and will give an estimation of the optimal weight penalty parameter. This allows the problem of over-fitting to be addressed without the need for validation data. Using the training data the prior distribution is updated to a posterior distribution by Bayes? theorem. The posterior probability distribution of the weights (w) given the training data , is denoted by | , | | , (4.5) where is the normalization factor, | represents a noise model on the target data and corresponds to the likelihood function, and is the prior probability of the weights (van Hinsbergen, 2009). The simplest and most commonly used forms of the prior and likelihood function are Gaussian and are as follows: , (4.6) | , (4.7) where ? and ? are normalizing constants. The expression for the posterior is then given by: 28 | , , (4.8) where , ? is a normalizing constant. The maximum of the posterior distribution is the most probable value of the weight vector and can be found by minimizing the negative logarithm of Eq. (4.8), which is equivalent to minimizing Eq. (4.2) (van Hinsbergen, 2009). In summary, the BNN procedure starts by constructing different networks with different initial weight values, then for a model, initial weight values for the hyperparameters are drawn from their priors. Then the networks are trained by the scaled conjugate gradient algorithm and for every step of the algorithm, values for ? and ? are re-estimated (van Hinsbergen, 2009). A common problem in the development of a neural network is determining the optimal number of hidden neurons. Based on the number of data points and preliminary model runs it was estimated that the optimal number of hidden neurons for this application was between 1 and 3. A more complicated double cross-validation scheme was using to train the BNN so that the optimal number of hidden neurons could be determined with more accuracy. 29 Figure 4.5: Illustration of the double cross-validation approach used to train and test the BNN models. The outer loop (CV1) is used to train and then test the forecast skill of the model. In CV1 the test data is shown in red and the training data in white. In this loop the test data segment is moved repeatedly in 1-year increments from the start to the end of the data record, so the forecast skill is tested on the entire data record. The CV2 loop is used to determine the optimal number of hidden neurons to use for the BNN model trained in CV1. The training data segment from CV1 is further segregated into validation data (in green) and a smaller training data segment (in white). In CV2 BNN models are trained using various numbers of hidden neurons and then evaluated on the validation data. The number of hidden neurons from the BNN model with the highest forecast skill on the validation data is selected to be used in the BNN in CV1. In the CV2 loop the validation data segment is moved repeatedly in 1-year increments from the start to the end of the training data record, so the forecast skill is evaluated on the entire training data record. Figure 4.5 shows the rigorous double-cross validation scheme. The inner loop (CV2 in Fig. 4.5) determines the optimal number of hidden neurons, which is then used in the BNN model trained in the outer loop (CV1 in Fig. 4.5). The final output for the BNN model was an ensemble forecast of 30 model runs and as with the MLR model various different combinations of predictors were used and are shown in Table 4.2. Table 4.2: Predictor combinations for BNN. BNN model Predictor 1 Predictor 2 1 EVI 2 MODIS-NDVI 3 AVHRR-NDVI 4 EVI MODIS-NDVI 30 4.5 Model-Based Recursive Partitioning The MOB forecast models were developed using the party package in R. The MOB model partitions the data like a standard regression tree, but then a linear model is developed for the data in each partition region. The MOB function uses the following algorithm: i) Fit a linear model ( ? ) for the observations in the current node. ii) Assess the stability of the model parameters with respect to each of the partitioning variables ( , , ). If there is some overall instability, select the variable associated with the highest instability, otherwise stop. iii) Compute the split points(s) for the selected variable that optimize the objective function. iv) Repeat the procedure in each of the child nodes (Zeileis et al., 2008). The MOB requires a hyperparamter called minsplit that gives the minimum number of data points that can be in any terminal node. For this application the values considered for minsplit were 10%, 20%, 30%, 40% and 50% of the training data. The MOB model was trained using the same double cross-validation method as used with the BNN model (Fig. 4.5), except that the inner loop (CV2 in Fig. 4.5) was used to determine the value for minsplit instead of the number of hidden neurons. The MOB model was developed using various different combinations of regression and partitioning variables as shown in Table 4.3. Table 4.3: Regression and partitioning variable combinations for MOB. MOB model Predictor 1 Predictor 2 Partitioning Variable 1 Partitioning Variable 2 1 EVI EVI 2 MODIS-NDVI MODIS-NDVI 3 AVHRR-NDVI AVHRR-NDVI 4 EVI MODIS-NDVI EVI MODIS-NDVI 5 EVI MODIS-NDVI MODIS-NDVI 6 MODIS-NDVI EVI MODIS-NDVI 31 4.6 Model Evaluation As previously mentioned, the forecast models were evaluated using a cross-validation method so the model forecast skill was evaluated over the entire data record. The model was evaluated on the test data segment by calculating the mean absolute error (MAE) and the skill score (SS). MAE is calculated as follows: ? | | , (4.10) where n is the number of data points, is the predicted value and is the measured value. The skill score is used to determine the skill of a forecast model by comparing it to a base or reference model such as climatology or persistence. In the case of this study the skill score is calculated as follows: ( ) , (4.11) where MAE from using the MLR, BNN or MOB model, MAE of the climatology model (i.e. average crop yield for a given CAR) and 0 MAE of the perfect model. A model with SS > 0 means that it has greater skill than the base model, whereas SS < 0 means less skill than the base model. The skill score makes it easy to compare different models and to determine if a model has any forecast skill. In addition, 95% confidence intervals for the skill scores were calculated using the bootstrap confidence interval method of Gilleland (2010). Bootstrapping is a technique in which data points are resampled with replacement from the sample data to create a large number of bootstrap samples for the purpose of approximating the sampling distribution of a statistic (Singh and Xie, 2008). In this case the skill score is computed for each of the bootstrap samples to determine the distribution of the skill score which can then be used to determine the confidence intervals. 32 4.7 Summary In summary, a correlation analysis was performed to select the optimal time during the growing season to make a crop yield prediction and determine the predictor variables. A clustering model was applied and the CARs were grouped into 1 to 6 clusters based on their historical crop yield. Multiple linear regression, Bayesian neural networks and model-based recursive partitioning forecast models were developed using various sets of predictors and each type of model was trained on data from individual clusters determined by hierarchical clustering and by the Agro-Climatic Zones scheme of Mkhabela et al. (2011). Finally, the MAE and skill score was calculated to evaluate the models. 33 Chapter 5: Results and Discussion 5.1 Barley Multiple linear regression, Bayesian neural networks and model-based recursive partitioning models to predict the yield of barley were developed for the entire Canadian Prairie region, for clusters of CARs determined by a hierarchical clustering model and for each agro-climate zone used in Mkhabela et al. (2011). The hierarchical clustering was based on the historical barley yield of each CAR and the grouping of the CARs can be seen in Figure 5.1 which illustrates how the CAR grouping changes as the number of clusters varies. In general the clustering model grouped CARs that are close together which makes sense because CARs with similar geography will likely behave similarly with respect to barley yield. It should also be noted that when using 3 clusters (Fig. 5.1b) the grouping from the clustering model is significantly different than the grouping based on the three agro-climatic zones (Fig. 5.1f) from Mkhabela et al. (2011). Using each possible scheme for grouping CARs a forecast model was developed for each cluster using various combinations of predictors. In order to accurately compare the results of models with different numbers of clusters the skill score of each model was calculated over the entire data set. 34 Figure 5.1: Grouping of CARs from hierarchical clustering model based on barley and by agro-climatic zones. 35 Tables 4.1, 4.2 and 4.3 show the combinations of predictors used for each type of model. Tables 5.1, 5.2 and 5.3 show the skill scores of the MLR, BNN and MOB models for various predictor sets and methods used to group the CARs, with the highest skill score in bold. Table 5.1: Skill score of MLR forecast models for barley. # of Clusters: 1 2 3 4 5 6 ACZ Mean MLR 1 SS: 0.090 0.086 0.088 0.088 0.082 0.087 0.083 0.086 MLR 2 SS: 0.222 0.220 0.206 0.202 0.199 0.196 0.221 0.209 MLR 3 SS: 0.052 0.052 0.044 0.041 0.032 0.040 0.047 0.044 MLR 4 SS: 0.223 0.220 0.238 0.236 0.244 0.240 0.219 0.231 MLR 5 SS: 0.226 0.235 0.242 0.234 0.243 0.236 0.227 0.235 Mean: 0.163 0.163 0.164 0.160 0.160 0.160 0.159 Table 5.2: Skill score of BNN forecast models for barley. # of Clusters: 1 2 3 4 5 6 ACZ Mean BNN 1 SS: 0.088 0.077 0.072 0.078 0.076 0.076 0.073 0.077 BNN 2 SS: 0.216 0.223 0.197 0.197 0.185 0.180 0.230 0.204 BNN 3 SS: 0.045 0.039 0.038 0.032 0.000 0.041 0.042 0.034 BNN 4 SS: 0.237 0.248 0.215 0.213 0.203 0.200 0.230 0.221 Mean: 0.147 0.147 0.131 0.130 0.116 0.124 0.144 Table 5.3: Skill score of MOB forecast models for barley. # of Clusters: 1 2 3 4 5 6 ACZ Mean MOB 1 SS: 0.088 0.084 0.083 0.088 0.093 0.097 0.083 0.088 MOB 2 SS: 0.225 0.202 0.225 0.211 0.208 0.195 0.230 0.214 MOB 3 SS: 0.047 0.047 0.039 0.035 0.030 0.043 0.047 0.041 MOB 4 SS: 0.231 0.238 0.248 0.211 0.203 0.190 0.232 0.222 MOB 5 SS: 0.235 0.236 0.248 0.206 0.196 0.189 0.234 0.221 MOB 6 SS: 0.215 0.220 0.211 0.204 0.213 0.199 0.240 0.215 Mean: 0.174 0.171 0.176 0.159 0.157 0.152 0.178 36 Figures 5.2, 5.3 and 5.4 compare the skill scores of models with different predictor sets and show the 95% confidence intervals developed using the bootstrap confidence interval method of Gilleland (2010). Figure 5.2: Skill score of MLR forecast models for barley with 95% confidence intervals. Figure 5.3: Skill score of BNN forecast models for barley with 95% confidence intervals. MLR 1 MLR 2 MLR 3 MLR 4 MLR 5-0.0500.050.10.150.20.250.30.350.4MLR Skill Score BarleySS 123456ACZBNN 1 BNN 2 BNN 3 BNN 4-0.1-0.0500.050.10.150.20.250.30.350.4BNN Skill Score BarleySS 123456ACZ37 Figure 5.4: Skill score of MOB forecast models for barley with 95% confidence intervals. It is clear that models developed using only MODIS-EVI (MLR 1, BNN 1 and MOB 1) or only AVHRR-NDVI (MLR 3, BNN 3 and MOB 3) as predictors have little to no forecast skills because for MLR, BNN and MOB none of the models based on these predictors had a skill score above 0.1 and the lower range of the 95% confidence interval is not above 0. However, all models which include MODIS-NDVI as a predictor have significant forecast skill and although using MODIS-EVI alone as a predictor was not successful, adding it to MODIS-NDVI as an additional predictor appears to improve the forecast skill in most cases. The best results for the MLR models was a skill score of 0.244 for MLR 4 (MODIS-EVI and MODIS_NDVI as predictors) developed with 5 clusters. In general MLR 4 and MLR 5 (EVIxNDVI predictor added) showed the most forecast skills for the MLR models, which means that addition of MODIS-EVI and possibly the EVIxNDVI term as predictors seemed to improve forecast skill over using just MODIS-NDVI (MLR 2). For the Bayesian neural network models, BNN 4 (MODIS-NDVI and MODIS-EVI as predictors) developed with 2 clusters had the highest skill score of 0.248. The highest skill obtained by a MOB model was 0.248 for MOB 4 and 5 developed with 3 clusters. MOB 4 uses MODIS-NDVI and MODIS-EVI as both MOB 1 MOB 2 MOB 3 MOB 4 MOB 5 MOB 6-0.0500.050.10.150.20.250.30.350.4MOB Skill Score BarleySS 123456ACZ38 regression and partitioning variables and MOB 5 uses MODIS-NDVI and MODIS-EVI as regression variables and just MODIS-NDVI as a partitioning variable. Figure 5.5: Skill scores of the best performing models for barley. Figure 5.5 compares the best performing models and shows that the forecasting of barley does not appear to be very sensitive to the type of model used which implies a mainly linear relationship between barley yield and the vegetation indices. The best results were obtained by BNN and MOB models, but there was not a large difference in skill from the best linear models. It is not clear what the optimal number of clusters to group the CARs into for barley prediction is, but it is clear that performing a clustering analysis before making a forecast is useful because in this case for almost every model at least one of clustering schemes lead to better results than developing the model from CARs grouped by the agro-climatic zones. MLR 4 MLR 5 BNN 4 MOB 4 MOB 500.050.10.150.20.250.30.350.4Skill Scores of Best Models for BarleySS 123456ACZ39 5.2 Canola A hierarchical clustering model was applied to the canola yield data and Figure 5.6 shows how the CAR groupings evolved as the number of clusters varied. The CAR grouping is similar to the results for barley (Fig. 5.1) with much of Alberta in one cluster and most of Manitoba in another, the main difference being the grouping in Saskatchewan, especially as the number of clusters increases. 40 Figure 5.6: Grouping of CARs from hierarchical clustering model based on canola and by agro-climatic zones. Tables 5.4, 5.5 and 5.6 show the skill scores of the MLR, BNN and MOB models for various predictor sets and methods used to group the CARs, with the highest skill score in bold. 41 Table 5.4: Skill score of MLR forecast models for canola. # of Clusters: 1 2 3 4 5 6 ACZ Mean MLR 1 SS: 0.037 0.035 0.033 0.04 0.056 0.059 0.035 0.042 MLR 2 SS: 0.165 0.163 0.170 0.170 0.169 0.167 0.157 0.166 MLR 3 SS: 0.059 0.065 0.059 0.085 0.082 0.083 0.053 0.069 MLR 4 SS: 0.159 0.158 0.159 0.157 0.192 0.189 0.144 0.165 MLR 5 SS: 0.152 0.153 0.156 0.151 0.185 0.191 0.139 0.161 Mean: 0.114 0.115 0.115 0.121 0.137 0.138 0.106 Table 5.5: Skill score of BNN forecast models for canola. # of Clusters: 1 2 3 4 5 6 ACZ Mean BNN 1 SS: 0.03 0.012 0.002 0.006 0.049 0.042 0.020 0.023 BNN 2 SS: 0.147 0.159 0.154 0.151 0.147 0.147 0.159 0.152 BNN 3 SS: 0.056 0.053 0.033 0.063 0.077 0.075 0.049 0.058 BNN 4 SS: 0.146 0.123 0.127 0.128 0.167 0.155 0.136 0.140 Mean: 0.095 0.087 0.079 0.087 0.110 0.105 0.091 Table 5.6: Skill score of MOB forecast models for canola. # of Clusters: 1 2 3 4 5 6 ACZ Mean MOB 1 SS: 0.033 0.035 0.033 0.04 0.056 0.068 0.032 0.042 MOB 2 SS: 0.136 0.152 0.155 0.162 0.172 0.172 0.143 0.156 MOB 3 SS: 0.052 0.053 0.055 0.075 0.064 0.059 0.044 0.057 MOB 4 SS: 0.120 0.121 0.133 0.148 0.182 0.178 0.128 0.144 MOB 5 SS: 0.123 0.139 0.149 0.152 0.190 0.181 0.129 0.152 MOB 6 SS: 0.139 0.135 0.166 0.149 0.149 0.161 0.131 0.147 Mean: 0.101 0.106 0.115 0.121 0.136 0.137 0.101 Figures 5.7, 5.8 and 5.9 illustrate the effectiveness of different predictor sets for each type of model by comparing the skill score and 95% confidence intervals for each set of predictors. 42 Figure 5.7: Skill score of MLR forecast models for canola with 95% confidence intervals. Figure 5.8: Skill score of BNN forecast models for canola with 95% confidence intervals. MLR 1 MLR 2 MLR 3 MLR 4 MLR 5-0.1-0.0500.050.10.150.20.250.30.35MLR Skill Score CanolaSS 123456ACZBNN 1 BNN 2 BNN 3 BNN 4-0.15-0.1-0.0500.050.10.150.20.250.3BNN Skill Score CanolaSS 123456ACZ43 Figure 5.9: Skill score of MOB forecast models for canola with 95% confidence intervals. The results for MLR canola forecast models again show that only models containing MODIS-NDVI as a predictor (MLR 2, MLR 4 and MLR 5) have significant forecast skills with MLR 4 (MODIS-EVI and MODIS-NDVI as predictors) developed with 5 clusters having the highest skill score of 0.192. MLR 5 tended not to perform as well as MLR 2 or 4, which means the EVIxNDVI cross term was not as useful in this case. Both BNN 2 (MODIS-NDVI as predictor) and BNN 4 (MODIS-EVI and MODIS-NDVI as predictors) had forecast skill for canola, but the BNN models did not perform as well as MLR for canola. The best BNN model was BNN 4 developed with 5 clusters, which had a skill score of only 0.167. The MOB model outperformed the BNN for canola and was more comparable to MLR. MOB 2, 4, 5 and 6 all had forecast skill with the highest skill score being 0.19 for MOB 5 (MODIS-EVI and MODIS-NDVI for regression and MODIS-NDVI for partitioning) developed with 5 clusters. MOB 1 MOB 2 MOB 3 MOB 4 MOB 5 MOB 6-0.1-0.0500.050.10.150.20.250.3MOB Skill Score CanolaSS 123456ACZ44 Figure 5.10: Skill scores of the best performing models for canola. In general it appears that MLR and MOB models are more effective for canola prediction than BNN models with the highest overall skill of 0.192 from MLR 4. As with barley, the application of a clustering algorithm will likely improve canola forecasting and based on Figure 5.10 it can be estimated that the optimal number of clusters was 5 or 6. MLR 2 MLR 4 MLR 5 BNN 2 BNN 4 MOB 4 MOB 5-0.0500.050.10.150.20.250.30.35Skill Scores of Best Models for CanolaSS 123456ACZ45 5.3 Spring Wheat The results from the hierarchical clustering model for spring wheat are shown in Figure 5.11. The grouping is not exactly the same as for barley and canola, but once again the CARs appear to be grouped largely based on geographical location. Figure 5.11: Grouping of CARs from hierarchical clustering model based on spring wheat and by agro-climatic zones. 46 Tables 5.7, 5.8 and 5.9 show the skill scores of the MLR, BNN and MOB models for various predictor sets and methods used to group the CARs, with the highest skill score in bold. Table 5.7: Skill score of MLR forecast models for spring wheat. # of Clusters: 1 2 3 4 5 6 ACZ Mean MLR 1 SS: 0.074 0.069 0.065 0.047 0.045 0.045 0.070 0.059 MLR 2 SS: 0.224 0.230 0.232 0.243 0.236 0.232 0.223 0.231 MLR 3 SS: 0.067 0.066 0.061 0.057 0.053 0.056 0.063 0.060 MLR 4 SS: 0.228 0.227 0.229 0.267 0.261 0.257 0.224 0.242 MLR 5 SS: 0.239 0.244 0.230 0.269 0.267 0.263 0.247 0.251 Mean: 0.166 0.167 0.163 0.177 0.172 0.171 0.165 Table 5.8: Skill score of BNN forecast models for spring wheat. # of Clusters: 1 2 3 4 5 6 ACZ Mean BNN 1 SS: 0.073 0.075 0.062 0.042 0.038 0.040 0.056 0.055 BNN 2 SS: 0.216 0.222 0.222 0.217 0.214 0.212 0.235 0.220 BNN 3 SS: 0.063 0.057 0.046 0.033 0.018 0.036 0.047 0.043 BNN 4 SS: 0.229 0.236 0.228 0.219 0.223 0.221 0.245 0.229 Mean: 0.145 0.148 0.140 0.128 0.123 0.127 0.146 Table 5.9: Skill score of MOB forecast models for spring wheat. # of Clusters: 1 2 3 4 5 6 ACZ Mean MOB 1 SS: 0.074 0.072 0.062 0.031 0.028 0.028 0.050 0.049 MOB 2 SS: 0.209 0.185 0.191 0.227 0.220 0.216 0.217 0.209 MOB 3 SS: 0.067 0.056 0.051 0.029 0.029 0.041 0.058 0.047 MOB 4 SS: 0.221 0.196 0.219 0.222 0.219 0.209 0.234 0.217 MOB 5 SS: 0.221 0.198 0.212 0.225 0.221 0.211 0.245 0.219 MOB 6 SS: 0.195 0.193 0.171 0.244 0.237 0.225 0.222 0.212 Mean: 0.165 0.150 0.151 0.163 0.159 0.155 0.171 Figures 5.12, 5.13 and 5.14 illustrates the effectiveness of different predictor sets for each type of model by comparing the skill score and 95% confidence intervals for each set of predictors. 47 Figure 5.12: Skill score of MLR forecast models for spring wheat with 95% confidence intervals. Figure 5.13: Skill score of BNN forecast models for spring wheat with 95% confidence intervals. MLR 1 MLR 2 MLR 3 MLR 4 MLR 5-0.100.10.20.30.40.50.6MLR Skill Score Spring WheatSS 123456ACZBNN 1 BNN 2 BNN 3 BNN 4-0.2BNN Skill Score Spring WheatSS 123456ACZ48 Figure 5.14: Skill score of MOB forecast models for spring wheat with 95% confidence intervals. As with barley and canola, for spring wheat the models using only MODIS-EVI (MLR 1, BNN 1 and MOB 1) or only AVHRR-NDVI (MLR 2, BNN 2 and MOB 2) as a predictor showed poor forecast skill. The MLR models for spring wheat using MODIS-NDVI as a predictor (MLR 2, MLR 4 and MLR 5) all show significant forecast skills with the highest skill score being 0.269 for MLR 5 (MODIS-EVI, MODIS-NDVI and EVIxNDVI as predictors) developed with 4 clusters. BNN 2 and BNN 4 showed forecast skill, with the highest skill score of 0.245 being for BNN 4 (MODIS-EVI, MODIS-NDVI as predictors) developed using the agro-climatic zone groupings for the CARs. Similarly all MOB models using MODIS-NDVI (MOB 2, MOB 4, MOB 5 and MOB 6) had significant forecast skill with the highest skill of 0.245 obtained for MOB 5 (MODIS-EVI and MODIS-NDVI for regression and MODIS-NDVI for partitioning) developed using the agro-climatic zone groupings for the CARs. MOB 1 MOB 2 MOB 3 MOB 4 MOB 5 MOB 6-0.3-0.2-0.100.10.20.30.40.5MOB Skill Score Spring WheatSS 123456ACZ49 Figure 5.15: Skill scores of the best performing models for spring wheat. The MLR model performed slightly better than the nonlinear models, but in general the forecasting of spring wheat does not appear to be very sensitive to the type of model used, implying a mainly linear relationship between spring wheat yield and vegetation indices. The optimal CAR grouping method varied, but for the nonlinear models the highest skill was usually obtained with the agro-climatic zone method. 5.4 Summary For all crops many of the forecast models developed showed forecast skill. For barley the highest skill score was 0.248 which was obtained by BNN 2 developed with 2 clusters and both MOB 4 and MOB 5 developed with 3 clusters. The highest skill score for canola was 0.192 for MLR 4 developed with 5 clusters and for spring wheat the highest skill score was 0.269 for MLR 5 developed with 4 clusters. In general skill scores for barley and spring wheat were higher than for canola and the type of model used did not appear to have a large effect on the skill score. Barley was the only case where the nonlinear models showed a slight forecasting advantage over multiple linear regression. The reason for this is likely because the models were mainly based on MODIS-NDVI and the relationship between each crop yield and MODIS-NDVI is essentially linear. MLR 4 MLR 5 BNN 4 MOB 4 MOB 500.050.10.150.20.250.30.350.40.45Skill Scores of Best Models for Spring WheatSS 123456ACZ50 Figure 5.16: Scatter plots illustrating the relationship between MODIS-NDVI and barley yield for the entire Canadian Prairie region. The data points from the year 2002 are marked by a black circle. The CAR means for both yield and NDVI have been removed from the data. The climatology predictions are shown as well as the forecasts for a) MLR, b) BNN and c) MOB models using only MODIS-NDVI as a predictor. 51 Figure 5.17: Scatter plots illustrating the relationship between MODIS-NDVI and canola yield for the entire Canadian Prairie region. The data points from the year 2002 are marked by a black circle. The CAR means for both yield and NDVI have been removed from the data. The climatology predictions are shown as well as the forecasts for a) MLR, b) BNN and c) MOB models using only MODIS-NDVI as a predictor. 52 Figure 5.18: Scatter plots illustrating the relationship between MODIS-NDVI and spring wheat yield for the entire Canadian Prairie region. The data points from the year 2002 are marked by a black circle. The CAR means for both yield and NDVI have been removed from the data. The climatology predictions are shown as well as the forecasts for a) MLR, b) BNN and c) MOB models using only MODIS-NDVI as a predictor. Figures 5.16 - 5.18 show actual crop yield, the climatology predictions and the model predictions for the entire data set and illustrate the relationship between crop yield and MODIS-NDVI which is fairly similar for all three crops. Most of the data points are closely grouped together and the relationship appears to be close to linear for those points, but when the yield data points in the extreme low and high ends are considered an overall slightly nonlinear relationship can be seen. In this study most of the low yield data came from one year (2002) during which parts of the Canadian Prairies suffered a severe drought. This lack of low yield data points outside of year 2002 could explain why the BNN and MOB models failed to capture this apparent nonlinear 53 relationship. In cross-validation when this extreme low yield year is used as the test data, the models have not seen such low values during training and must extrapolate. In most cases the BNN models tapered off towards the extremes and were not able to accurately predict high or low yield data. The MOB models performed better than the BNN at the extremes, but the MLR was the best model for predicting the extreme values of 2002. The BNN model behaves differently from the MLR model during extrapolation due to the bounded hyperbolic tangent activation function used by BNN. The BNN model output will remain bounded even if the predictor variables become extremely large in magnitude, while the MLR output will not (Hsieh, 2009, p.304). Therefore when the low yield data points from the year 2002 were used as test data the MLR model was able to better extrapolate to these extreme values not seen by the model during training than the BNN. If the data record was expanded and more years with high or low yield data were added, the BNN and MOB functions would likely be able to pick up the slight nonlinear relationship and better forecast the extreme values. 54 Chapter 6: Conclusion 6.1 Research Objectives Machine learning methods were used to develop crop yield forecast models for the Canadian Prairies. This study focused on the prediction of barley, canola and spring wheat using remotely sensed vegetation indices and compared the effectiveness of MODIS-NDVI, MODIS-EVI and AVHRR-NDVI. The results of Mkhablea et al. (2011) were confirmed as MODIS-NDVI was shown to be an effective predictor of crop yield. However, the lower resolution AVHRR-NDVI was found to have little to no forecast skills which agrees with the results of Li et al. (2007) where MODIS-NDVI was shown to outperform AVHRR-NDVI for corn and soybean yield prediction. Li et al. (2010) compared MODIS-NDVI and MODIS-EVI and found that in general the NDVI outperformed EVI in monitoring several different vegetation cover types. Similarly in this study MODIS-NDVI was found to be more effective than MODIS-EVI as forecast models developed using only MODIS-EVI showed poor crop forecasting ability. From this study it is clear that MODIS-NDVI was the most effective vegetation index for forecasting crop yields on the Canadian Prairies. However, the use of MODIS-EVI as an additional predictor tended to increase the skill score of models using just MODIS-NDVI, so EVI appeared to supply complementary information useful for crop yield prediction. As shown in Figure 4.1d the correlation between MODIS-EVI and MODIS-NDVI drops around the time the crop forecast is made, so the EVI is likely measuring different properties of the crops than NDVI. Many current crop forecasting methods use only linear techniques (Prasad et al. 2006) which could miss nonlinear relationships in the data. In many cases such as Jiang et al. (2004), Kaul et al. (2005) and Li et al. (2007) nonlinear neural networks were found to outperform linear techniques for crop yield prediction. In this study multiple linear regression, neural networks and model-based recursive partitioning forecast models were developed and evaluated on the same barley, canola and spring wheat data sets. In general the results from all three types of models were very similar. For barley the nonlinear models slightly outperformed the MLR models, for both canola and spring wheat the MLR models showed the highest skill, but the MLR skill scores were only slightly higher than those for BNN or MOB models. The results of this study 55 were not able to demonstrate the advantages of nonlinear techniques and this was likely due to the mainly linear relationship between crop yields and vegetation indices and/or the data record being too short to provide sufficient training data for the BNN and MOB models to learn the nonlinear relations. With only one extreme drought year in the data record (2002), when the yield data from 2002 were the test data during cross-validation, since the models were trained without extreme drought data, they had to extrapolate. Unfortunately, BNN models are structured to extrapolate more gently than linear models and that lead to BNN under-predicting the severity of the extreme low yields. Based on the literature review in Chapter 2 there was sufficient evidence to suggest that climate indices could be used to predict crop yields on the Canadian Prairies. For example, Hsieh et al. (1999) found a nonlinear relationship between the Pacific sea surface temperature and Canadian Prairie wheat yield. It was thought that climate indices could help account for the large-scale meteorological conditions affecting crop yields and improve forecasting skills, however based on the correlation analysis (Chapter 4.1) the climate indices showed little to no forecasting potential. Models with climate indices as extra predictors were tested (not shown), but the results were not useful and therefore were not included as final models for this study. It seems whatever climate effects have on the crop yield, the satellite vegetation indices may have already incorporated the climate variability information, so adding climate indices as extra predictors was not fruitful. 6.2 Future Research This study clearly demonstrates the potential for using MODIS-NDVI to forecast crop yield on the Canadian Prairies and the addition of more predictor variables to these models could improve upon these results. Many studies including Jiang et al. (2004), Kaul et al. (2005), Prasad et al. (2006) used more than just vegetation indices to predict crop yields. Meteorological variables such as precipitation and surface temperature could be added to the models developed in this study to further improve crop yield forecasting on the Canadian Prairies. The addition of these new variables will make the problem more complex and possibly better demonstrate the advantages of nonlinear BNN and MOB models. 56 The main limitation for this research was the limited data record for MODIS vegetation indices which spanned from 2000-2011. Meteorological variables like precipitation and temperature have much longer data records and if these variables could be added to the models in this study it might be possible to extend the data record to better train the model. If more predictors were added the model might be less sensitive to the MODIS-NDVI variable and possibly AVHRR-NDVI could be used to fill in for MODIS-NDVI in earlier years of the data record. The best way to group the CARs for the development of crop forecasting models was not clear because the optimal number of clusters was not easily determined. However, in most cases several of the CAR groupings from the clustering algorithm led to models with higher skills than those developed from the agro-climatic zone CAR grouping. If a longer data record was available it might be possible to use some of the data for validation to better estimate the optimal number of clusters. 57 Bibliography Alvarez, R., 2009: Predicting average regional yield and production of wheat in the argentine pampas by an artificial neural network approach. European Journal of Agronomy 30, 70-77. Bhat, P. C., Prosper, H. B., 2005: Bayesian neural networks. Proceedings of statistical problems in particle physics. Astrophysics and Cosmology, Oxford, UK, 12?15, 151-154. Bonsal, B., Shabbar, A., 2011: Large-scale climate oscillations influencing Canada, 1900-2008. Canadian Biodiversity: Ecosystem Status and Trends 2010, Technical Thematic Report No. 4. Canadian Councils of Resource Ministers. Ottawa, ON. iii + 15 pp. http://www.biodivcanada.ca/default.asp?lang=En&n=137E1147-0 Chen, P., Fedosejevs, G., Tiscareno-Lopez, M., Arnold J. G., 2006: Assessment of MODIS-EVI, MODIS-NDVI and vegetation-NDVI composite data using agricultural measurements: an example at corn fields in western Mexico. Environmental Monitoring and Assessment 199, 69-82. Di Lorenzo, E., Cobb, K.M., Furtado, J.C., Schneider, N., Anderson B.T., Bracco, A., Alexander, M.A., Vimont, D.J., 2010: Central Pacific El Ni?o and decadal climate change in the north Pacific Ocean. Nature Geoscience 3, 762-765. Di Lorenzo, E., Schneider, N., Cobb, K. M., Franks, P. J. S., Chhak, K., Miller, A.J., McWilliams, J. C., Bograd, S.J., Arango, H., Curchitser, E., Powell, T. M., Riviere, P., 2008: North Pacific Gyre Oscillation links ocean climate and ecosystem change. Geophysical Research Letters, Vol. 35, L08607, doi: 10.1029/2007GL032838. Gardner, M. W., Dorling, S. R., 2000: Statistical surface ozone models: an improved methodology to account for non-linear behaviour. Atmospheric Environment 34, 21-34. 58 Gilleland, E., 2010: Confidence Intervals for Forecast Verification. http://citeseerx.ist.psu.edu/ viewdoc/download?doi=10.1.1.174.5002&rep=rep1&type=pdf Guo, W. W., Xue, H., 2012: An incorporative statistic and neural approach for crop yield modeling and forecasting. Neural Computing and Applications 21, 109-117. Hsieh, W. W., Tang, B., 1998: Applying neural network models to prediction and data analysis in meteorology and oceanography. Bulletin of the American Meteorological Society 79, 1855-1870. Hsieh, W. W., Tang, B., Garnett E. R., 1999: Teleconnections between Pacific sea surface temperatures and Canadian Prairie wheat yield. Agriculture and Forest Meteorology 96, 209-217. Hsieh, W. W., 2009: Machine Learning Methods in Environmental Sciences. Cambridge University Press, Cambridge, 349 pp. Huete, A., Didan, K., Miura, T., Rodriguez, E.P., Gao, G., Ferreira, L.G., 2002: Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment 83, 195-213. Jiang, D., Yang, X., Clinton, N., Wang, N., 2004: An artificial neural network model for estimating crop yields using remotely sensed information. International Journal of Remote Sensing 25, No. 9, 1723?1732. Kaul, M., Hill, R. L., Walthall, C., 2005: Artificial neural networks for corn and soybean yield prediction. Agricultural Systems 85, 1-18. 59 Li, A., Liang, S., Wang, A., Qin, J., 2007: Estimating crop yield from multi-temporal satellite data using multivariate regression and neural network techniques. Photogrammetric Engineering and Remote Sensing 73, No. 10, 1149-1157. Li, Z., Li, X, Wei, D., Xu, X., Wang, H., 2010: An assessment of correlation on MODIS-NDVI and EVI with natural vegetation coverage in northern Hebei Province, China. Procedia Environmental Sciences 2, 964-969. Mkhabela, M.S., Bullock, P., Raj, S., Wang, S., Yang, Y., 2011: Crop yield forecasting on the Canadian Prairies using MODIS NDVI data. Agricultural and Forest Meteorology 151, 385-393. Murtagh, F., Legendre, P., 2011: Ward?s hierarchical clustering method: clustering criterion and agglomerative algorithm. arXiv preprint arXiv:1111.6285, http://arxiv.org/abs/1111.6285. Prasad, A. K., Chai, L., Singh, R. P., Kafatos, M., 2006: Crop yield estimation model for Iowa using remote sensing and surface parameters. International Journal of Applied Earth Observation and Geoinformation 8, 26?33. Rasmussen, C. E., 1996: A practical Monte Carlo implementation of Bayesian learning. Advances in Neural Information Processing Systems 8, 598-604. Shabbar, A., Bonsal, B. R., Szeto, K., 2011: Atmospheric and oceanic variability associated with growing season droughts and pluvials on the Canadian Prairies. Atmosphere- Ocean, 49, 339-355. Shabbar, A., Skinner, W., 2004: Summer drought patterns in Canada and the relationship to global sea surface temperatures. Journal of Climate 17, 2866-2880. 60 Singh, K., Xie, M., 2008: Bootstrap: A Statistical Method. http://www.stat.rutgers.edu/home/mxie/rcpapers/bootstrap.pdf Song, J., Zhou, W., Li, C., Qi, L., 2009: Signature of the Antarctic oscillation in the northern hemisphere. Meteorology and Atmospheric Physics 105, 55-67. Stewart, R., Pomeroy, J., Lawford, R., 2011: The drought research initiative: a comprehensive examination of drought over the Canadian Prairies. Atmosphere-Ocean, 49, 298-302. Thompson, D.W. J., Wallace, J.M., 1998: The Arctic oscillation signature in the wintertime geopotential height and temperature fields. Geophysical Research Letters 25, 1297- 1300. van Hinsbergen, C., van Lint, J., van Zuylen, H., 2009: Bayesian committee of neural networks to predict travel times with confidence intervals. Transportation Research Part C: Emerging Technologies 17, 498 ? 509. Vincente-Serrano, S. M., Saz-Sanchez, M. A., Cuadrat, J. M., 2003: Comparative analysis of interpolation methods in the middle Ebro Valley (Spain): application to annual precipitation and temperature. Climate Research 24, 161-180. Wettstein, J.J., Mearns, L.O., 2002: The influence of the North Atlantic?Arctic Oscillation on mean, variance, and extremes of temperature in the northeastern United States and Canada. Journal of Climate 15, 3586?3600. Wilks, D. S., 2006: Statistical Methods in Atmospheric Sciences. 2nd edition. Elsevier, Amsterdam. 627 pp. 61 Zeileis, A., Hothorn, T., Hornik, K., 2008: Model-based recursive partitioning. Journal of Computational and Graphical Statistics 17, 492-514.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Crop yield forecasting on the Canadian Prairies by...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Crop yield forecasting on the Canadian Prairies by satellite data and machine learning methods Johnson, Michael David 2013
pdf
Page Metadata
Item Metadata
Title | Crop yield forecasting on the Canadian Prairies by satellite data and machine learning methods |
Creator |
Johnson, Michael David |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | The production of grain crops plays an important role in the economy of the Canadian Prairies and early reliable crop yield forecasts over large areas would help policy makers and grain marketing agencies in planning for exports. Forecast models developed from satellite data have the potential to provide quantitative and timely information on agricultural crops over large areas. The use of nonlinear modeling techniques from the field of machine learning could improve crop forecasting from the linear models most commonly used today. The Canadian Prairies consist of the provinces of Alberta, Saskatchewan and Manitoba and three of the major crops in this region are barley, canola and spring wheat. The agricultural land on the Canadian Prairies has been divided into Census Agricultural Regions (CAR) by Statistics Canada. A clustering model was applied to the crop yield data to group the CARs for the development of forecast models. The normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) derived from the Moderate Resolution Imaging Spectro-radiometer (MODIS), NDVI derived from the Advanced Very High Resolution Radiometer (AVHRR) and several climate indices were considered as predictors for crop yields. A correlation analysis between crop yield and the time series of each potential predictor was performed to determine which variables showed the most forecasting potential and at what time during the growing season their values were most correlated to crop yield. Various combinations of MODIS-NDVI, MODIS-EVI and NOAA-NDVI were used to forecast the yield of barley, canola and spring wheat. Multiple linear regression as well as nonlinear Bayesian neural networks and model-based recursive partitioning forecast models were developed using the various sets of predictors. The models were trained using a cross-validation method and the forecast results of each model were evaluated by calculating the skill score from the mean absolute error, with 95% confidence intervals for the skill scores calculated using a bootstrap method. The results were compared in an effort to determine the optimal set of predictors and type of forecast model for each crop. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-10-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0103360 |
URI | http://hdl.handle.net/2429/45281 |
Degree |
Master of Science - MSc |
Program |
Atmospheric Science |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2013_fall_johnson_michael.pdf [ 1.64MB ]
- Metadata
- JSON: 24-1.0103360.json
- JSON-LD: 24-1.0103360-ld.json
- RDF/XML (Pretty): 24-1.0103360-rdf.xml
- RDF/JSON: 24-1.0103360-rdf.json
- Turtle: 24-1.0103360-turtle.txt
- N-Triples: 24-1.0103360-rdf-ntriples.txt
- Original Record: 24-1.0103360-source.json
- Full Text
- 24-1.0103360-fulltext.txt
- Citation
- 24-1.0103360.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0103360/manifest