Over the past two decades, the population in Pakistan has increased exponentially, and, this rapid growth has caused two major problems among many others. The two significant issues of importance here are i) increased demand for food; and ii) the decrease in arable land. Food production, both in terms of higher productivity and improved nutrient quality, has not kept pace with the needs of the expanding human population, leading to alarming levels of food insecurity. Extreme weather conditions, such as floods and droughts, associated with climate shifts across the globe are projected to cause even greater adverse effects on the food supply (Briscoe and Qamar, 2006). Accurate and timely information obtained from monitoring observations of crop growth and development and forecasting yield at large scale is instrumental in addressing these concerns.
Agriculture is central to Pakistan’s economy, and wheat is the main agricultural product which is grown on a large scale over a wide spectrum of climatic and soil conditions. Wheat production comprises about 9.9% of the total agricultural productivity and contributes 2.0% to GDP (Khan et al., 2009a). Nearly 80% of the farmers are directly involved with the wheat production and this crop covers 40% of the total cropped area in Pakistan (Faruqee et al., 1997). Wheat production is acknowledged as a prime source of livelihood in rural and agrarian sectors of Pakistan. In 1975-76, about 19.5% household income in the agricultural sector of Pakistan was generated through wheat (Cornelisse and Opdam, 1982). Food security is the primary requirement for the stability of a nation. Thus, it is essential to develop, increase, and sustain cultivable land and improve the nation’s agriculture system. Understanding and following worldwide farming practices are crucial to cope with both short-term and long-term threats to sustainable food production. The expected increase in deficient food supplies is attributed to variations in environmental and meteorological factors, water scarcity, and the rapidly growing population. Accurate estimation of potential yield prior to harvest can significantly aid the policy guidance and contributes to better and timely decisions (Ali, 2001).
The traditional approach of measuring the crop yield is based on ground visits or field surveys. This type of data is crucial since it provides thorough information, and because it embodies direct observations in a point of interest typically associated with high quality measurements. On the other hand, the modern approach is based on remotely sensed satellite-based vegetation index observations. Both approaches have their advantages and limitations, however, the integration of these techniques could lead to the accurate estimation of the crop yield. These two kinds of data complement each other in several ways in the yield estimation process. The evaluation of crop yield output often depends on the traditional data gathering techniques, which, while generally accurate at the point source, are time intensive, and susceptible to bias errors (Reynolds et al., 2000). In many countries, the limited or late access to the available data inhibits the assertive actions that can be taken to prevent potential food shortages. The traditional structure operated by the provincial government for recording and archiving agricultural output is reliant on manual and time-consuming ground data surveying in randomly sampled villages. The field survey observations from inaccessible areas are also not likely recorded. Most of the data are obtained from field surveys and enumeration methods, and data synthesis and assimilation are finalised a number of months after the harvest. Hence, the application of such data for prompt decision-making and satisfactory planning for wheat surplus or shortages is often limited (Ibrar, 2012). This is especially true if information and communication technology is not sufficiently advanced to disseminate timely information for decision makers.
In contrast to this, remotely sensed satellite data provide timely coverage to large areas even to those regions where the ground survey is not possible. The evolution in remote sensing and the availability of real-time satellite imagery have provided wide access to spatial variability of the Earth features (Schuler, 2002). Satellite data have the potential to identify not only the crop types but also the crop yield estimates (Mohd et al., 2009). Hence, remote sensing can be used to monitor and detect crop growth stages and to predict crop yield (Sawasawa, 2003). The phenological stages and the crop health can be identified through the spectral responses of the crop by using multi-spectral data. Normalised difference vegetation index (NDVI) is a measure of vegetation health. Healthy plants show high values of NDVI since they have high reflectance in infrared wavelength and somewhat low reflectance in red wavelength. Crop yield estimation at regional scale using MODIS data through data assimilation method is described in (Liang and Fang, 2004). In this study, low-resolution MODIS data product including the enhanced vegetation index (EVI) and leaf area index (LAI) has been used for yield estimates. However, this is also a fact that the low or coarse resolution of images results in the generalisation of the crop status and yield estimates. The coarse resolution data combine the crop and other non-vegetation observations that contribute to NDVI data, which could later be correlated to the final crop yield recorded (Sawasawa, 2003). Several studies conducted for vegetation mapping and crop growth monitoring have reported a high correlation between vegetation indices (VI) and crop yield (Groten, 1993; Rasmussen, 1997). However, most of these studies were done at field levels or in a research environment, which involves very small plots as an area of interest (AOI) with spectral data being collected either from low altitude aerial platforms or from some ground based platforms. Such an approach provides an efficient level of control over many extraneous factors and quality of data, which eventually result in an exceptional correlation between the remotely sensed and physically measured data (Staggenborg and Taylor, 2000). It has been proposed in Quarmby et al. (1993) that NDVI can be a key factor for the estimation of crop yield. According to (Sawasawa, 2003), for the accurate estimation of crop yield under variable agro-meteorological conditions, the potential of regression modelling has been efficiently established. In (Ren et al., 2007), the work on winter wheat yield prediction has been described at a regional scale, which uses crop biomass estimation by employing data from multiple sources (i.e. soil moisture data, meteorological data, and remote sensing data) for the Huanghuaihai plain in North China. It was concluded through this research that the application of data from various sources was effective and showed satisfactory results for crop yield estimation.
In (Doraiswamy et al., 2007a), the work on operational assessment of corn and soybean crop yield using MODIS data for Iowa and Illinois States has been described. It was concluded from this study that the differences between the National Agricultural Statistics Service (NASS) State level yield and the regression model estimates for both Iowa and Illinois were less than 4 and 2 bushels per acre (b/ac) in 2006 season for corn and soybeans respectively.
In (Bernardes et al., 2012), it is suggested that the correlation between the biennial coffee crop yield and VI is quite significant. Even in the presence of low spatial resolution of MODIS data and the fact that yield is a composite factor and is sensitive to several conditions, the VI were able to establish the relationship between leaf biomass and coffee crop and could be valuable for the estimation of biennial yield.
In (Kalubarme and Sharma, 2014), the work on wheat production forecasting in Himachal Pradesh using remote sensing data and GIS is presented. Since wheat is cultivated mainly as a rainfed crop in these districts, the accuracy of wheat acreage estimates also varies to some extent with normal and unusual rainfall conditions. However, the performance of agro-meteorological yield models in general, is comparable with Department of Agriculture estimates, except for very unusual drought conditions. The Relative Deviations (RD) of wheat acreage estimates of remote sensing based district level was within 0.2 to 17% in comparison with Department of Agriculture estimates for all districts.
The national space agency of Pakistan, Space and Upper Atmosphere Research Commission (SUPARCO), started developing crop area estimation procedures and crop yield models in 2005. The crop acreage was estimated by land stratification (in terms of agriculture and non-agriculture land) based on visual interpretations followed by extensive ground surveys and supervised classification using SPOT-5 satellite data. For supervised classification, Gaussian Maximum Likelihood method was applied. For yield prediction, SPOT vegetation data was used and crop phenology was mapped along with the related NDVI values (Ahmad et al., 2014). Although an accuracy of 85-95% was achieved in the classification stage but the use of supervised learning method and extensive ground surveys made the overall process slow and very much region specific and also practically infeasible to scale on the national level due to its expensive and time-consuming nature.
In Dempewolf et al. (2014) an approach has been proposed for the prediction of the wheat yield of Punjab Province of Pakistan, by using moderate (Landsat 4-5 and 7) and low-resolution MODIS data at different stages of wheat growing season. The land cover classification of was performed using a supervised learning approach i.e. bagged decision tree method. The training stage of classification required a former understanding of wheat growing areas and data for classification were selected by visual interpretation of the Landsat time series, based on the phenology of wheat. The regression analysis was performed for yield forecast using the time series based vegetation indices of MODIS and Landsat data and their performance comparison was done with statistical measurements provided by Crop Reporting Services (CRS) for Punjab province. The drawback of this approach is linked to the supervise learning where training data were selected through visual interpretation of optical images. The quality of training samples was dependent on optical images with no or low cloud cover. This has made the approach more region specific and reduced the feasibility of transferring the methodology to other, more cloud prone regions. It can be inferred from above that among different approaches used to predict crop yield from remotely sensed data; the most common approach is to develop a regression model based on the direct empirical relationship between NDVI measurements and crop yield. The research work presented in this paper is also inspired by this widely accepted approach to crop yield prediction; however, prior to yield forecast, the primary focus of this work is establishing a relationship between the meteorological parameters (temperature and rainfall) and wheat growing stages considering 14 years of data (2001-14). The impact of rainfall and temperature on wheat development stages has been investigated by examining the data recorded on the daily and monthly basis. In addition to this, the rainfall and temperature normals (1981-2010) have also been incorporated in the analysis. The regions of interest are classified in terms of rainfed and irrigated areas so that the effect of meteorological parameters can be analysed on wheat yield based on different cropping system. For land cover stratification, unsupervised classification has been adopted followed by wheat area mapping using crop calendar. In addition to this, parabolic least square fitting (LSF) has been applied considering the non-linear response of wheat crop in terms of NDVI. The application of parabolic LSF was instrumental in approximating an appropriate value for NDVI maxima, which resulted in obtaining improved crop yield estimates.
The scope of this research work is confined to rainfed and irrigated croplands i.e. Chakwal and Faisalabad districts of Punjab province, however, once the study is established, it can be scaled up to any crop land. The selection of the agro-meteorological attributes is limited to temperature and rainfall. Apart from precipitation and temperature, there are many other human-induced factors such as the seed quality, selection of soil type, crop type, inflow water, and fertilisers used, which can affect the yield of the crop. These factors have not been incorporated in this study mainly due to their inaccessibility and non-availability. Also, most of these factors are hard to apply over a large study area as their properties vary from region to region. Such factors are suitable to apply if the study is conducted in a small region with known land properties. This study will be helpful for our policy makers to formulate an appropriate food security policy for rainfed and irrigated areas. This would also allow the policy makers and agriculture planners to ascertain the amount of crop to export or import while conserving sufficient resources for nationwide food security.
Materials and methods
The methodology is divided into three major tasks. The first step is to generate the wheat area maps. For this purpose, the land cover mapping is performed using unsupervised classification of satellite data and then the classified images are refined for wheat area mapping using the crop calendar. The second step is to analyze the effect of meteorological parameters on wheat yield and identifying any relationship between them. This has been accomplished by first generating wheat profiles of the study area based on NDVI data and then the effect of variation in rainfall and temperature on wheat yield for the crop season has been analysed. The final step is to estimate the wheat yield based on a regression model. This has been achieved by regressing the NDVI data against the reported yield and then validated the estimated data by comparing it with the recorded yield. Figure 1 represents the framework of entire methodology.
Chakwal (rainfed cropland)
Chakwal (Figure 2) is located in Potohar plateau at 32:9303o N, 72:8556o E and an altitude of 498m in Punjab province, Pakistan, with a total area of about 6,524 km2. Agriculture in Potohar region is largely dependent on weather conditions. The reported average annual rainfall for Chakwal is about 880mm. In summer weather is hot while in winter it is cold and dry. The average temperature during winter is 8oC, which rises up to 42oC during summer. The major crops are Wheat, Ground Nut, Sorghum, Rape/Mustard Seed, Millet (Bajra) and Gram. Wheat is grown on most of the land, which is then used for both commercial and personal purposes.
Faisalabad (irrigated cropland)
Faisalabad (Figure 2) district is geographically located at 31.4187° N, 73.0791° E and an altitude of 184m in Punjab province, Pakistan, with an area of 5,856 Km2. It has a population of approximately 3.5 million and is the second most populous district in Punjab (Aser, 2012). The district has extreme meteorological conditions. In summers, the average maximum and minimum temperature reach 39 and 27oC respectively, whereas in winter, average maximum and minimum temperature approaches to 21 and 6oC, respectively. The winter season begins in November and stays till March while the summer season begins in April and lasts until October. The cropping system of Faisalabad district is dependent on canal-based irrigation or cultivated through tube wells or is rainfed. The district of Faisalabad is incomparable for its agricultural production. The main Kharif crops of the region are millet, maize, sugar cane, and rice and winter crops include wheat, gram and barley (itspakistan, 2012).
For crop type mapping, impact assessment of meteorological parameters on crop, and yield estimation, different types of data sets were acquired from multiple sources (Table 1).
For trend analysis of wheat yield in relation to meteorological variations, MOD13Q1 NDVI product of MODIS sensor of TERRA satellite was acquired for the year 2001-14 from (https://ladsweb.nascom.nasa.gov/) with 16 days of temporal resolution and 250 meters of spatial resolution. The MODIS VI products (MOD13) give regular, temporal and spatial contrasts of vegetation conditions on a global scale. The phenological stages and the crop health can be identified through the spectral responses of the crop by using multi-spectral data. In this context, NDVI is known as an excellent measure of vegetation health. The NDVI is the normalised values of the ratio (between -1 and +1) of NIR and Red bands and is commonly expressed as:
where ρNIR is the reflectance of the near-infrared wavelength band and ρRed is the reflectance of the red wavelength band (Hamlyn and Robin, 2010).
It can be further used for monitoring various photosynthetic activities that take place in plants and also their phenological conditions at different stages of their growth. Since NDVI is a ratio, it has the advantage of minimising the atmospheric noise among the correlated bands. The ratio is also helpful for reducing the calibration or instrument associated errors. The extent to which the ratio can minimise the noise is linked to the correlation of noises between the Red and NIR bands (Solano et al., 2010).
Crop calendar data was obtained from National Agriculture Research Council (NARC), Islamabad. The crop calendar states the key agronomic observations in cropping system per year which varies from area to areas such as sowing and harvesting period. Most of the original data or records were gathered by the United Nations Food and Agriculture Organization (FAO) (FAO, 2016b).
Meteorological data was obtained from Pakistan Meteorological Department (PMD), Islamabad. This data comprised of total monthly rainfall and mean monthly temperatures, and rainfall and temperature normals (1981-2010) for 2001-14.
Other ancillary data, such as wheat yield statistics for the years 2001-13 were provided by NARC and Pakistan Bureau of Statistics, Islamabad (Table 2).
Land cover classification and wheat area mapping
Major land covers were classified by using unsupervised classification technique. In contrast to supervised classification, unsupervised classification is considered to be simple as no prior signatures or samples are required. It allows determining the parameters, which the computer utilises as guidelines to reveal spectral patterns in the data. Iterative self-organising (ISO) data clustering was applied in this research among several unsupervised classification techniques. It provides an understanding of the particular spectral-temporal clusters in the multi-temporal NDVI data and their relative landscape positions. The resultant classes of unsupervised classification are spectral in nature. As the classes are solely based on the natural or implicit clustering in the image data so their identity will not be known in prior. ISO data clustering is a threestep process. In the beginning, it randomly assigns preliminary cluster vector. In the following step, it classifies every pixel based on its closest cluster. In the last step, the subsequent cluster mean vectors are computed on the basis of the entire pixels comprising one cluster. The last two steps are recurrent until the change or variation between the iterations is negligible or very small. The variation can be determined either by calculating the distances of the mean cluster vectors from iteration to iteration or alternatively by measuring the ratio of pixels that have altered across the iterations. The ISO data algorithm is further modified by dividing and combining the clusters. The clusters are fused if either the number of pixels in a cluster is below a certain limit or threshold or if the proximity of clusters centres exceeds a predefined value (Jensen, 1986). The clusters are divided into two different clusters if the cluster standard deviation breaches a pre-set value and the number of pixels is doubled the threshold for the minimum number of members. By using ISO data clustering in ERDAS Imagine 2013, stacked NDVI profiles of 14 years were classified to a range of pre-set classes (5 to 100). The maximum number of iterations was set to 100 and convergence threshold was defined as 0.999 respectively. As a result, in total, 20 classified images were generated using batch processing. Subsequently, the mean statistics of all stacked NDVI layers was computed and NDVI profiles were generated for the study areas and mainland covers were identified. At the end, the crop calendar was also used to refine NDVI profiles for wheat.
Investigating the effect of meteorological parameters on wheat yield
In order to analyze the effect of meteorological variations, the max NDVI values of the wheat profiles have been extracted from the preceding step. The NDVI profile of wheat represents a parabolic shape. This non-linear behaviour is owing to the green pigment in the crop. The NDVI curve changes with the growth of the crop as it is linearly related to total chlorophyll content. When the crop is in the preliminary stage, the green pigment is not in high density, but as the crop grows in size, the green colour becomes more prominent and influences the magnitude of NDVI (Reddy et al., 2001). The max NDVI is observed at the centre of the curve and corresponds to the development stage of the crop. When the crop reaches its maturation stage and turns golden, the NDVI curve starts to decline. The wheat NDVI profiles (2001-14) of the study areas are then correlated with the total monthly rainfall and mean monthly temperature data. The rainfall and temperature normals (1981-2010) are also used where thorough analysis was required. The normals (1981-2010) are latest three-decade averages of climatological variables, including temperature and precipitation. They are used as a reference or benchmark in climatological and agricultural studies. The analysis was categorised in terms of rainfed and irrigated croplands so that the effect of meteorological variation can be discerned under different cropping systems.
Wheat yield estimation using linear and non-linear Regression
The final step is to estimate the wheat yield based on a regression model. This has been achieved by regressing the NDVI data against the reported yield and then validated the estimated data by comparing it with the recorded yield. The non-linear regression has also been applied to wheat NDVI profiles in order to improve the yield estimates.
Results and discussion
In land cover mapping, the unsupervised classification has been performed to classify the land. In order to select the most suitable classified image, which can best represent the different classes, the divergence statistics were evaluated for each test image using the Separability feature of Erdas Imagine. Divergence measures the spectral distance between the generated cluster’s signatures. The minimum divergence represents strong similarity among test classes while average divergence demonstrates the existence of distinct classes.
The unsupervised classified image with 60 classes was selected by analysing the divergent curve made through divergent maxima values of all outputs of the classified images (Figure 3). The selected file of 60 classes was used to separate the pixels associated vegetation and non-vegetation lands based upon NDVI values.
After unsupervised classification, mean of each layer was computed by using statistics tool in Erdas Imagine to generate NDVI profiles for the crops. On the basis of these profiles, the major land covers were identified, i.e. i) vegetation and ii) nonvegetation lands. The vegetation cover consists of several crops with overlapping growth cycle. In order to perform wheat area mapping, phenology curves were plotted for each class in the selected image. From phenology curve dates corresponding to peaks and dips were noted and these dates were cross-checked with crop calendar to identify the crop type. The NDVI profiles, which matched the wheat crop calendar, were extracted and labeled as Wheat and rest of the profiles were labeled as vegetation and non-vegetation land. The seasonal profiles of wheat for Chakwal and Faisalabad districts for the years (2001-14) are shown in Figure 4.
From these profiles, the wheat crop phenology at different stages can be determined over time. The monitoring of this NDVI pattern is also a good indicator of the crop yield pattern. In Chakwal and Faisalabad, the wheat crop shows the maximum NDVI peak from mid of February to early March. At harvesting stage, the NDVI values of the crop start to decline and reach a minimum value approaching 0.2. In Figure 4, it is observed that the NDVI curves are very low for the year 2009-10 which was the drought year for Chakwal (Kazmi and Rasul, 2012). In contrast to this, the years 2004-05 (wet year), 2006-07 and 2012-13 profiles have shown the maximum trend as they were the high-yielding years (Table 2).
As Chakwal is a rainfed area, wheat yield was severely affected by the dry weather. This is evident in Figure 5 also.
Impact of meteorological parameters on wheat crop
Wheat is a Rabi crop, i.e. a crop that is sown in winter and harvested in the spring season. In Pakistan, sowing of wheat takes place from October to December and harvesting during the month of March to May. In Punjab, sowing starts in November and goes up to the mid of December whereas harvesting period is April and May (PMD, 2012). However, in Sindh province, wheat crop sowing period is November-December and harvesting period is March-April (PMD, 2015). Normally the wheat crop takes 100-120 days to get fully matured (PMD, 2012). In this section, the relation of crop yield with NDVI and the impact of meteorological parameters [in comparison with the normals (1981-2010)] on wheat crop growth and its yield have been discussed.
Relationship among wheat yield, rainfall, and normalised difference vegetation index
The rain is the key parameter impacting the crop growth. The water requirement of the vegetation cover varies in terms of composition and the climate of their origin (Morison and Baker, 2008). The effect of rainfall is more pronounced in rainfed areas as compared to irrigated croplands.
In Faisalabad (Irrigated Cropland), most of the area is irrigated with runoff water from rainfall in upland regions as an additional water source. Rain affects the wheat yield with the difference in precipitation during its growing stages. Wheat requires a rainfall of 50 cm to 100 cm during the growing season but for irrigated lands a rainfall of 40 cm to 50 cm is sufficient. The excessive rainfall beyond this threshold would damage the crop and insufficient rainfall may parch the grain (Jaydeep, 2002; Bose, 2013). Annual variation in rain and snowfall over plains and hills of the Pakistan also cause the proper supply of canal irrigation water.
The growth season of wheat crop consists of several phenological stages (Table 3). These stages need water for the whole development period, but there are some stages, which are more vulnerable to water shortage. Any water shortage during this period may result in serious yield losses. The highest amount of rainfall during the winter season is received in the month of February followed by March and January. Rain before sowing and shooting to grain formation stage triggers yield production, but with suitable intervals in required amount of the respective phenological phase as continuous rainy and moist atmosphere can adversely affect the normal growth of the wheat crop. Rain before full maturity phase also contributes to normal or above the normal yield of the crop. Rain just after sowing causes decrease in the number of germinated wheat seeds and also at the time of full maturity or harvesting negatively affects the crop growth and ultimately reduces the yield (PMD, 2012).
The wheat crop is irrigated three times during the entire season before full maturity. First irrigation is made to the wheat crop 21 days after sowing. During the tillering stage, second irrigation is given and the third irrigation is given during the heading stage. The shortage of irrigation water at tillering, shooting and early grain fill period results in significant yield losses. In order to understand the impact of precipitation on crop growth, the total rainfall has been plotted with yield and NDVI in Figure 6.
The recorded rainfall seems high in 2004-05 and 2007-08, however, the yield has increased in 2004-05 and dropped in the latter year. The difference in yield was owing to the amount of rainfall received during subsequent stages of wheat growing season in these years. The observed rainfall in 2007-08 was below normals (1981-2010) during the months of February and March (Figure 7). These are the months, which ought to receive the maximum amount of the seasonal rain. During this year only 18.2 mm rainfall occurred in early days of February and no rainfall occurred during heading to flowering stage i.e. from last week of February to 2nd week of March. It is very important that the crop receives maximum rain during this period. The 3rd week of April also received 24 mm of rain during full maturity which negatively affected the crop resulting in reduced yield. So the lack of rainfall in critical stages and excessive precipitation during the stages when rainfall was not required subsequently declined the yield. In 2004-05 considerable amount of rainfall occurred at appropriate times across the growing season, which has ultimately maximised the crop yield. This year was called as a wet year (National Coordinator Wheat Plant Sciences Division, Pakistan Agricultural Research Council, 2013).
The total rainfall of 2002-03 and 2006-07 is almost the same but there is a difference in their final wheat yield (Figure 6). The low yield in 2002-03 was attributed to more rainfall during the flowering stage. In the wake of third irrigation, a heavy spell of 54.7 mm in 1st week of March had damaged the crop. The observed rainfall was also higher than normal during shooting, heading and flowering stages (Figure 7). However, in December 2006, the rainfall of 46.2 mm was received which proved to be beneficial for the crop as it happened after the emergence phase. The month of February and March received sufficient rainfall and no rainfall was observed during full maturity, which was good for the crop. As the crop received water at all critical stages so the crop yield was higher relative to 2002-03.
The crop years 2003-04, 2005-06 and 2010-11 have the total rainfall almost on the same scale but they differ in their final wheat yield (Figure 6). In 2003-04, a rainfall of 25 mm after flowering and 34.3 mm during mature stage badly affected the crop and subsequently declined the yield (Figure 7). In 2005-06, the rainfall was below normal in most of the months except March. Insufficient water during heading and flowering stages affected the crop. No rainfall was received during first 13 days of March. As, the rainfall did not occur when it was highly required, subsequently resulted in the reduction of wheat yield.
In 2010-11, there is a rise in yield owing to the amount of rainfall (34 mm) received during shooting and heading stages. A rainfall of 17 mm was received before full maturity and no rainfall was observed during full maturity, which was good for the crop and ascended its yield.
The final comparison was made between 2008-09 and 2012-13. The total rainfall in these years was nearly the same but again the difference was in their final yield (Figure 6). In 2008-09, the rainfall was above the normal during December and significantly below the normal during March (Figure 7). During emergence, the third leaf, tillering and shooting stages, a considerable amount of rainfall (31.4 mm) was experienced which fulfilled the crop water requirement. During March, precipitation amounting to 12 mm [<normals (1981-2010)] was experienced, which affected the yield to some extent.
In 2012-13, the crop yield was relatively higher owing to the increased crop area but the actual yield was relatively less with reference to the area utilised for crop production. The low yield may partly be attributed to the unfavourable weather conditions and partly to the late sowing (3rd week of November) of the crop than the recommended period. Moreover, more rain in the month of February (59.5 mm) and April (25.4 mm) caused a disease in wheat crop. This unfavourable weather conditions may be responsible for below than normal yield of wheat crop (PMD, 2012).
The wheat yield and rainfall relation can also be discerned from NDVI spectral and temporal profiles that provide wide uses of remote sensing data for crop yield prediction. The total seasonal rainfall over the study period and the corresponding wheat NDVI trend are plotted in Figure 6. NDVI curves showed a positive correlation with the precipitation (where rainfall is received at all critical stages in the required amount) and negative correlation (where excessive rainfall has occurred at inappropriate periods). Hence, wheat yield and NDVI showed an identical pattern when analysed with total rainfall of the wheat season (Figure 6).
In the case of Chakwal, a rainfed area, the yield is directly proportional to the required rainfall, the higher the precipitation the larger the yield. Chakwal district is constrained by the monsoon range. In addition to the seasonal rainfall, there are two other rainy periods. The first season is triggered by the monsoon winds inset, which originates from the Bay of Bengal over the mid of July to mid-September. The second rainy period is caused by the Mediterranean winds and occurs in the last two weeks of December and the first two weeks of January. The average net rainfall is 558.8 to 635 mm. The sub-division, Choa Saidan Shah, has the highest rainfall in the region (Afzal et al., 2015). In Chakwal, there is a great variation in rainfall due to its geo-location. The wheat crop requires almost 450 to 650 mm (FAO, 2016a) of water per growing season. During 2001-14, mostly the amount of yearly rainfall was more than the average annual rainfall.
The relationship of total rainfall with crop yield can be seen in Figure 8. The year 2009-10 showed the minimum yield as it was the drought year and has low precipitation value of the decade which resulted in extremely low crop yield (Kazmi and Rasul, 2012). During this year, the environment became harsh and dry with a limited supply of rains and crop sowing was also delayed. The high-temperature extremities and water stress resulted in low crop yield (Ahmed et al., 2012).
The year 2011-12 also faced the drought conditions. During this year the rainfall was below the normals (1981-2010) (EPD, 2012), except for the month of October and April (Figure 9). The rainfall was not received at critical periods of growing season and instead, it occurred at full maturity stages, which has drastically scaled down the crop yield. The year 2005-06 was also considered as a low yielding year. During this year, the rainfall was below normals (1981-2010) in all months except March. The insufficient water during heading and flowering stages affected the crop. A low total rainfall of 135mm was received across the season which declined wheat yield (Figure 9).
Although the total rainfall of 2002-03 and 2007-08 was higher than 2003-04, a decline in yield was observed in these years when compared to the yields of 2003-04 (Figure 8). For 2007-08, this decline was owing to poor germination and frost damage (Subhani, 2011). The rainfall was below normal in most of the months except January and April and only 3 mm of rainfall received in March (Figure 9). It is pertinent to note here that an insufficient rainfall at critical stages of the crop can also affect the crop yield adversely (Jaydeep, 2002; Bose, 2013). The yields of 2004-05 and 2006-07 are very close but the rainfall received during 2004-05 was more than 2006-07 (Figure 8). The prime reason behind the affected yield in 2004-05 is the occurrence of 40mm rainfall at full maturity stage of the crop (Figure 9). From the above analysis, it can be inferred that the rainfall was the major parameter that affected the yield over the study period of 2001-14 in Chakwal. The inter-annual variability of the rainfall was also captured in the remotely sensed NDVI data as shown in Figure 8. The wheat yield and NDVI profiles showed strong positive correlation with the total precipitation across the study period as shown in Figure 10, where correlation coefficient R approached to 0.80 and 0.94 respectively. This observation is true for rainfed areas as the cropping system is totally reliant on the seasonal rainfall. In irrigated areas (like Faisalabad), the cropping system is supported by various water sources other than rainfall. Therefore, in such areas, it is difficult to establish such correlations.
Relationship among wheat yield, temperature, and normalised difference vegetation index
Temperature plays a critical role in the growth of the wheat crop. The photosynthesis process of the plant is sensitive to temperature and it has been found from a number of studies in South Asia and China that a 1oC increase in temperature during the wheat development season may result in a decline of 3% in crop yield (Hussain, 2010). During growth season, wheat requires low temperature and moisture while harvesting period is constrained by dry and warm weather. Wheat has absolute minimum temperatures ranging from 2 to 5oC, optimum minimum from 15 to 20oC, optimum maximum from 23 to 33oC, and absolute maximum from 27 to 38oC (Doraiswamy et al., 2007b; Hollinger and Angel, 2011).
The mean temperature required during sowing is 10oC to 15oC and during harvesting, it requires a higher temperature, but sudden rise of temperature is harmful and can adversely affect its growth (Jaydeep, 2002; Bose, 2013). The higher temperatures usually speed up the wheat development cycle so that it matures sooner, thus shortening the duration of grain filling period (Ashfaq et al., 2011). The increased heat stress often intensifies the pressure on water resources that are essential for crop growth.
As the meteorological data for this study was not publically available so its acquisition was a serious problem. The owners of meteorological data are government organisations that induce several administrative procedures and overheads prior to data access. Sometimes the data is also not available due to the nonexistence of the recording stations. Owing to the above problems, there was no access to the Chakwal district temperature profiles and only Faisalabad district temperature profiles have been analysed.
The crop seasonal mean temperature, yield and NDVI of Faisalabad district for the years 2001-13 have been plotted in Figure 11 and NDVI data and wheat yield were found to exhibit identical behaviour when observed with the mean temperature data.
The temperature profiles of the years 2003-04, 2005-07, 2008-11, and 2012-13 were favourable across the season. The temperature across the growing season was in compliance with the temperature normal (1981-2010) as shown in Figure 12. Lower temperatures were observed during early growing stages and higher at full maturity. The decline in yield in any of these years was owing to the erratic behaviour of rainfall. The rainfall was not either received in sufficient amounts when required or occurred in excess at odd periods.
On the other hand, the temperature profiles of the years 2001-03, 2007-08, and 2011-12 showed slightly high temperature than temperature normal (1981-2010) during early stages of the crop (Figure 12). In 2001-02, during tillering and shooting stages the temperature was high and heat stress hampered the crop growth. Similarly in 2002-03 high temperatures in early stages of crop affected the crop growth. In 2007-08 and 2011-12, the temperature was high during emergence stage and the crop did not receive rain when it was highly required thus resulting in low yield.
The year 2004-05 was the wet year and its temperature profile was considered optimal for the crop growth. The satisfactory rainfall at all critical stages and low temperature at early stages of crop produced high yield (Figure 12).
The crop seasonal mean temperature and yield have been regressed linearly to determine the relationship with them as shown in Figure 13. The correlation coefficient, R, is approximately 0.5 because the temperature is not the only factor affecting yield in irrigated areas. Apart from precipitation and temperature, there are many other human-induced factors such as the seed quality, selection of soil type, crop type, and fertilisers used, which can affect the yield of the crop.
Statistical analysis for yield prediction
The correlation between meteorological factors, wheat phenology (NDVI), and yield was analysed statistically. The behaviour or patterns were analysed for mean monthly meteorological telemetry to identify whether a substantial drift in the trend of temperature, precipitation, and wheat growth stages followed over the study period.
The reported yield and max NDVI were plotted together to discern the relationship between them (Figure 14). It is observed that max NDVI of the study areas followed the trend of the reported yield over time.
Considering the positive correlation of NDVI with the reported yield for both rainfed and irrigated croplands, the estimates of the wheat yield of the study areas were obtained by regressing the reported crop yield against the NDVI (See linear regression in Figure 15).
The regression equations obtained for both districts are as follows:
In equations, 2 to 4 x corresponds to the maximum NDVI. By using these equations, the yields of 2001-13 were estimated as shown in Figure 15.
These figures illustrate that the predicted yield has successfully captured the fluctuations in the wheat yield over time and is in compliance with the reported data. The correlation between the estimated and reported wheat yield of 2001-13 was also calculated with R2 to be 0.88 and 0.69 for Chakwal and Faisalabad respectively. The standard error for the yield estimation equations (2 to 4) was found to be 118.23 and 118.11 for Chakwal and Faisalabad districts respectively.
Improving yield estimates using parabolic least square fit
The wheat NDVI profiles have exhibited a parabolic behaviour across the Rabi season as shown in Figure 4. It has been discussed earlier in section 2.4 also that this non-linear behavior is owing to the green pigment in the crop. The max NDVI is observed at the centre of the parabolic curve and relates to the development stage of the crop. The non-linear pattern has manifested that wheat NDVI profile across the season has geometric properties of a parabola. For this purpose, the LSF of NDVI curves were computed. In the previous section, the relationship between NDVI and yield has been established using the max NDVI. But if the NDVI curves are analysed carefully then it is observed that the power of NDVI is distributed among the neighboring peaks also. So in order to fine-tune the NDVI curves and to uniformly distribute the max NDVI values, LSF has been applied. The other reason for applying LSF is owing to the temporal resolution of NDVI data. The NDVI data has a temporal resolution of 16-days so the curves represented in Figure 4 cannot represent true max NDVI and therefore an approximation was required to bridge the resolution gap and quantify the max NDVI value. As the R2 of Chakwal district seems quite suitable and is already closer to 1 so the parabolic LSF of NDVI curves were computed for Faisalabad district only (Figure 16) due to its relatively small R2.
Hence, the estimation equation was reformulated (as shown below) with Parabolic LSF NDVI as a parameter of yield prediction (Figure 17).
By using the equation 5, the crop yield for the years 2001-13 was re-estimated for Faisalabad. The reported yield and estimated yield (with and without parabolic LSF) are shown in Figure 17, and it is evident that in contrast to estimated yield without parabolic LSF, the estimated yield using parabolic LSF has shown high correlation and good agreement with the reported yield. Figure 17 shows the scattered plot of the reported and predicted wheat yield (using the parabolic LSF) from 2001-13 with R2 to be 0.73 and standard error was 110, indicating that the estimated results obtained after applying parabolic LSF have a high degree of reliability and adherence with the recorded telemetry.
Conclusions and Recommendations
This research has demonstrated that MODIS imagery captured at high temporal resolution provide an optimal way of monitoring crops and estimating yields. In this study, the effect of temperature and rainfall on wheat yield has been analysed and the relationship between NDVI and crop yield has been established. The meteorological variables affect wheat crop at different stages of production including germination, vegetation, and maturity and have a positive effect if they occur in the proper amount and at critical stages of growth.
The effect of rainfall, temperature and other associated factors on the wheat crop is mapped using NDVI profiles. It is a fact that in rainfed areas like the Potohar region, the crops are totally reliant on the frequency and spatial distribution of rainfall, therefore, the excess or deficiency of rain may adversely affect the yield, especially at critical growth stages of wheat. It is concluded from this study that the total rainfall is not the only factor that can trigger the increase or decrease in yield, but an appropriate amount of rainfall at various crop development stages is a requisite for maximising the yield. Different crop growth stages have variable sensitivity levels of development to water stress. The most critical time is when the wheat crop is approaching its final maturation stage. At this stage, the wheat crop requires warmer temperature for dryness and the excessive rainfall during this stage will decline its yield. In fact, persistent rain during the germination stage is also detrimental for crop yield. It does not provide enough gaps to the soil suitable for seed sowing. On the other hand, if rain is delayed then deficient soil moisture fails to sustain the seed germination process. The analytical results have demonstrated that in addition to crop monitoring, NDVI has the potential to predict crop yield as well. This has been realised from the positive correlation found between NDVI and the reported crop yield. The wheat yield has been calculated for the years 2001-14 by applying the regression analysis and the forecasted yield was validated by comparing with the reported yield statistics (Provided by NARC and Pakistan Bureau of Statistics, 2009). The R2 between the reported and predicted yield was found to be 0.88 and 0.69 for Chakwal and Faisalabad respectively. It has been observed that wheat NDVI profiles exhibit parabolic shape due to the variation in the green pigment of the crop across the growth period. A least square parabolic fit has been applied to these wheat profiles to fine tune the NDVI maxima. The least square parabolic fit proved to be instrumental in improving the wheat yield estimates and the R2 between the reported and predicted yield approached to 0.73 for Faisalabad respectively. This clearly indicates that the method is capable of providing competitive accuracies that are statistically reliable with the already identified pattern of vegetation in the region. It can then be concluded that the regression method has the potential to predict the crop yield with acceptable accuracy. However, the accuracy of the yield estimates is highly dependent on the temporal and spatial resolution of images, the texture or type of land cover and the size of historical data used for regression analysis. Based on the entire analysis and conclusion of this work, the following are the recommendations for further work.
The empirical results can be additionally improved by acquiring high-resolution spatial-temporal imagery like Landsat 30m data. The R2 value can be improved by increasing the size of the reported yield data to build yield estimation equation.
In this research, the crop yield is only analysed against temperature, precipitation, and area. Other associated parameters like fertilisers, soil type, pesticides, and irrigation applied, and seeds etc. are not addressed. It is recommended that for precise yield prediction, these factors should be incorporated in the yield estimation equation. Regression analysis for yield estimation can be further refined by applying nonlinear or polynomial based models.
It is recommended that the implemented technique should be applied on regions having distinctive climatic conditions so the potential of this research can be realised and accepted for diverse climatic conditions.
In addition to the above recommendations, this work can be also extended to explore the impact of the recent practice of multiple cropping on remote sensing analysis. The multiple cropping practices are beneficial for maximising the food productivity and make effective use of inputs such as soil, water, and fertilisers. The impact of changing cropping practices on remote sensing techniques is not trivial. For instance, we require high-resolution spatial data for crop acreage estimation and hyperspectral data for the classification of the crop in the same field. However, the remote sensing techniques would not be directly affected by the variation in crop growing cycle as crop cycle is dependent on environmental changes. If there is a change in crop growing cycle then the crop calendar needs to be updated and the requisite data would be obtained accordingly. Multiple cropping techniques have benefits, which are manifold but the lack of information, research, resource, and skills are some of the reasons for its low adoption. Keeping in view the economic benefits of multiple cropping, there is a need to promote it among the farming community.