Crop Yield Estimation in 2014 for Vojvodina Using Methods of Remote Sensing

Monitoring phenology of crops and yield estimate based on vegetation indices as well as other parameters such as temperature or amount of rainfall were largely reported in literature. In this research, MODIS Normalized Difference Vegetation Index (NDVI) was used as an indicator of specific crop condition; the other parameter was Land Surface Temperature (LST) which can indicate the amount of crop moisture. Trial years were 2011, 2012, and 2013. For those years sowing structure was acquired from agricultural organizations Nova Budućnost from Žarkovac and Sava Kovačević from Vrbas, both in Serbia. Also, satellite images with high and medium resolution for these areas and years were available. Multiple linear regression was used for crop yield estimate for Vojvodina Province, Serbia where the NDVI and LST were independent variables and the average yield for specific crop was the dependent variable. The results of crop yield estimate two months before harvest are presented (excluding wheat).


Introduction
Monitoring the condition of crops is important for the economic development of every country.Monitoring of plants using remote sensing has many advantages such as detection of very large areas (whole Vojvodina), or very easy derivation of vegetation and LST for near-real time monitoring of crops, which is the case with MODIS (Moderate resolution imaging spectroradiometer).Doraiswamy et al. (2007Doraiswamy et al. ( , 2003)), Prasad et al. (2005), Curnel & Oger (2007) reported that vegetation indices have direct relation with the condition of the plants, so they can be used for the purpose of yield estimation and creation of crop maps.The crop yields are of essential benefit for the economy of each country, especially if they were previously estimated.Satellite remote sensing data in combination with existing public data from agricultural companies and D. Jovanović* • F. Sabo • M. Govedarica University of Novi Sad, Faculty of Technical Sciences, 6 Trg Dositeja Obradovića, 21000 Novi Sad, Serbia e-mail: dusanbuk@uns.ac.rsB. Marinković University of Novi Sad, Faculty of Agriculture, 8 Trg Dositeja Obradovića, 21000 Novi Sad, Serbia Ratar. Povrt. 51:3 (2014) [145][146][147][148][149][150][151][152][153] government statistics agencies, represents a viable solution for early assessment of yield.Doraiswamy et al. (2007) used MODIS data to predict yield for soybean and corn and they compared the results with official results.The predictions were within a 20% standard deviation of the official estimates.Prasad et al. (2005) used AVHRR (Advanced very high resolution radiometer) NDVI with 8x8 km spatial resolution together with rainfall data (RF), land surface temperature (LST) and soil moisture (SM) in order to estimate yield for corn and soybean.Their method involved nonlinear regression (Quasi-Newton method).They mention that variations of NDVI, SM, LST and RF data do not follow linear combination with respect to crop yield.That study was carried out based on the 19 years datasets, while this study observed only three years, addressed local area, and plots were not as large as in Prasad et al. (2005) study.
This study tested the usefulness of the MODIS NDVI and LST for several cultivated plants.NDVI dynamics is directly related to biomass of certain plants, while LST indicates the amount of plant humidity.Observed cultivated plants in this study were: wheat, mercantile and seed corn, sugar beet and soybeans.
The main study objectives were: a) to identify the crop of interest using accurate data and highresolution images on a wider area; b) to determine the beginning and the end of the season for each cultivated plant by tracking phenophases of cultivated plants using NDVI; c) to obtain yield estimate for 2014 for each cultivated plant by using multiple regression with monitored parameters (NDVI and temperature).

Study area
Area of interest included wider territory of Žarkovac (municipality of Ruma) and wider territory of Vrbas (Figure 1) in the Province of Vojvodina, Serbia.For these study areas high-resolution images and sowing structure of cultivated plants for the years 2011, 2012, 2013, and 2014 were available.Ratar. Povrt. 51:3 (2014)  On the basis of accurate data from the agricultural companies, plants were identified in the wider area, not only in the parcels of agricultural company for which data were available, but also for other, smaller parcels.This is achieved by using the spectral characteristics of Landsat 5, 7 and 8 as well as RapidEye (RE) satellite sensors.In this way, monitoring sampling area was expanded, so we have achieved representative NDVI and temperature parameters.The area of the municipality of Vrbas (Fig. 1a) was representative for Bačka region and Žarkovac for Srem region (Fig. 1a), and yield estimate was calculated for the entire territory of Vojvodina.

Data
MODIS's 16-day composite vegetation indices (MOD13Q1 label) and MODIS's 8-day LST (MOD11A2) were used for sampling and monitoring of cultivated crops.Landsat mission (5, 7, and 8) and RapidEye data together with data from agricultural company (farms) were used in order to identify the agricultural crops of interest.RapidEye has five bands, while Landsat missions have a slightly higher number, but mainly we have used the following bands: blue, green, red, infrared, mid-infrared and RapidEye's red edge band.The spatial resolution of the MODIS NDVI is 250 m (231.66 m), MODIS LST spatial resolution is 1 km, Landsat 30 m, 15 m pan sharpened images (Landsat 7 and 8), whereas RapidEye has the highest spatial resolution of 5 m.Projection for all images was UTM, zone 34, WGS-84 ellipsoid.MODIS's default projection is sinusoidal, so it was necessary to reproject it in conformal transverse cylindrical projection (UTM), for MODIS's data this is done using MRT-web online program.All images were imported into Erdas Imagine software in standard imagine format.
In this study we have used a total of 4 RapidEye images, 15 Landsat missions images, 75 composite MODIS NDVI images, and 36 LST images all for the period from January 2011 until June 2014.
The NDVI is one of the most popular and surely one of the most used indices in remote sensing.
In theory the NDVI values range from -1 to +1, in practice this range is narrower.Barren rock, sand areas have very low NDVI values, around 0.1.Areas with scarce vegetation have greater NDVI values, from 0.2 to 0.5.High NDVI values are reserved for dense vegetation areas, forests and cultivated plants in peak of their season.NDVI was first used by Rouse et. al (1973).Formula for NDVI is as follows: (1) NIR represents Near-infrared reflection, MODIS band 2. RED is red reflection i.e. reflection of MODIS band 1.
It is very important to understand the difference between air temperature and LST.LST is the temperature that would be felt if the Earth surface was to be touched in a particular location.This surface can include rooftops, bare ground, top of vegetation cover, water, etc. Remote sensing sensors record long-wave infrared (LWIR) energy in their thermal bands and this data together with emissivity values and water vapour content can be used to derive LST.MODIS generates the LST by using split -window algorithm (Wan 1999) and bands 31 and 32 (LWIR).

Identifying culture and phenology of crops
Monitoring phenology with remote sensing is not a novelty.This method of involving remote sensing into something new was first introduced by Tucker et al. (1979).Since then, many studies have used NDVI as the phenology indicator.Some case studies were based on Enhanced Vegetation Index (EVI) from MODIS (Son et al. 2013, Sakamoto et al. 2005).Apart from just monitoring phenology, several studies dealt with methods of determining exact phenophases of crops (Curnel & Oger 2007, White et al. 2002).On the other hand, the main goal of our study was to estimate yields based on ground truth data, auxiliary data and based on farming practices.Farming practice experience helped us in more precise phenology detection.Based on all this data we were able to detect phenophases for each crop specifically and for each year.Because of the spatial resolution of MODIS's indices (250 m), in order to detect one pixel of cultivated plant on the image it was necessary to identify the parcels that are larger than 6 ha.Because of the fact that one pixel (6 ha) will be mixed, heterogeneous pixel (consists of different cultivated plants and other land cover, e.g.corn, soybean and road or smaller lake), it is desirable to monitor parcels larger than 10 ha.This study has shown that parcels larger than 10 ha can be discriminated with MODIS spatial resolution, whilst parcels less than 10 ha result in mixed pixels and must be excluded from further analysis.
Agricultural companies "Sava Kovačević" in Vrbas and "Nova Budućnost" in Žarkovci own large parcels which can be monitored with this resolution (50 ha average size).White et al. (2002) mentioned that the average size of their areas of research was 32 ha or larger, we had larger areas, so no mixed pixels were included.In order to avoid mixed pixels, it was necessary to select training polygons in that way that the entire pixel belongs to the observed crops (Figure 2).Poor selection of polygons will include pixels that do not belong to the observed crop, the so-called boundary pixels.
In addition to the good selection of pixels, it was necessary to extend the sampling area using Landsat and RapidEye's spectral properties.One spectral profile was taken as a reference (exact planting structure), and based on that reference profile the area of monitored crops was extended for each crop and for each year.Figure 3 (territory of Žarkovci) is an example of spectral profile which has been used for identification of soybean in nearby area.This procedure was applied to all monitored cultivated crops, for all available images (Landsat 5, 7, 8, and RapidEye).In Figure 4b it can be seen that due the plant growth, NDVI values increase gradually because of the higher NIR reflectance.With the increase of NIR reflectance values over time, the NDVI values reaches its peak value, this period is defined as maximum number of "green" days, growing days.After that, NDVI values begin to gradually decline, this represents the aging days of plants or "brown" days.Figure 4a and 4c shows two different stages of corn, June 19 and August 9, 2013.In June, corn is slowly starting to grow, while in July, and August, reaches its maximum vegetation period, which primarily depends on the weather conditions.Figure 4d and 4f shows wheat for two periods, May 21 and June 19 from the remote sensing aspect.Wheat reaches a maximum maturity of vegetation in May and in June slowly begins to change colour (Figure 4e).At the end of June wheat is the most characteristic part of the satellite image, it becomes "scorched" (yellow-brown), and is readily distinguished from other cultivated plants.Figure 5 depicts how remote sensing (NDVI data) notices the growth stages of wheat in 2012.Start of "green" period for wheat is on 65th day of the year (March 5) and the maximum growing period is reached on 129th day of the year (May 8).From that day wheat gradually begins to age and on 177th day of year (June 25th) reaches its minimum vegetation period when the harvest is performed.It is possible that the harvest is performed before or after that date, but certainly just close to that day.6 is approximately 0.6, and after that in mid-July 2014 NDVI has its peak value.Figure 6 shows that at the beginning of June, corn is still small and not dense enough, but at mid-July, when the density is high and corn reaches height of approximately 190 cm (measured on field), NDVI increases to a value between 0.8 and 0.9.Such monitoring of crops can discriminate almost all the plants; the curve in Figure 6 is different for all plants.In this way we can identify the beginning and the end of season for each plant individually.Ratar. Povrt. 51:3 (2014)  Figure 6 is actually a spectral temporal profile of mercantile corn, where x-axis is the day of the year, and y-axis is the mean NDVI value of days per year.This profile is obtained by calculating the mean value of NDVI for each 16 day composite period for the entire year (which gives a total of 22 composite images).The same principle is applied to all other plants that are monitored in this paper.The gradual increase of NDVI value up to the peak is seen as a growth of observed plants (period when the plants are greenest).Gradual reduction of NDVI value from the peak value to the date when the harvest performed is characterized as aging period of plants ("brown-yellow" days of certain plants).Table 1 shows the total area of monitored crops for each year separately.This monitored area is based on ground truth area, which was increased as possible.Plantation structure of seed corn for 2011 year was not available so we excluded it from analysis.
Statistical analysis of the data From Figure 6 it can be concluded that the value of NDVI curve from mid-July is predicted in order to estimate yields for the current year.Growing season for mentioned plants ends in late September or early October, depending on the plant and weather conditions.Results of this paper were created on 17 July 2014.We have predicted NDVI values for corn, sugar beets, seed corn and soybeans and those results were used to estimate the yield of 2014.
It was noted that the first half of the season 2014 (from 7 April to 26 June) for the observed plants are in a high correlation with the first half of the season 2013.Thus the values of NDVI for the first half of the season of 2013 were used to obtain a single linear regression and to predict values for the rest of the season for 2014.Value of NDVI from the first half of the season 2014 was dependent variable Y, independent variable X was NDVI value for first half of the season 2013.
Based on the equations from Table 2 and on the basis of the value of NDVI from the second half of the season in 2013 for the above observed plants, the behaviour of plants in 2014 was predicted.
Figure 7.A strong linear relationship between the two seasons for mercantile corn Such prediction was made for all plants using plants from 2013 and linear regression (Figure 7).For the purpose of yield estimate, we used multiple linear regression.NDVI and LST values were used as independent variables, and the average yield as the dependent variable.
(2) where, Y is the average yield, parameter α is constant value, β1 and β2 are regression coefficients.α, β1 and β2 are derived from multiple linear regression by using the Least Squares Method.
NDVI parameter is summarized for the period from the beginning of greening to the end of the aging for every plant.The parameter LST is summarized for the period from April 6 until to July 3, since we could not predict these values, the correlation was low for any indicator.NDVI and LST parameters for each crop for the period from 2011 to 2013 (trial years) are then used in equation ( 1) with the average yield Y, after which coefficients α, β1, β2 for each crop are calculated.Based on these coefficients the equation for predicting yield for 2014 was formed.
Using seasonal values for NDVI and LST for 2014 as in the previous years and derived coefficients of regression, crop yield estimation was performed.
Average yields for the previous years were supplied by the Republic Statistical Office and these values were used for the Y (Table 3).Table 3 shows the low average yields for the drought year of 2012.The NDVI values for that year were also very low, NDVI curve reached lower values than in other years, and summarized LST was higher by about 15°C.

Results and Discussion
Comparing yield values from Table 3 and Table  4, good average yields for 2014 can be noticed.NDVI values were higher than in previous years with the sampled areas.By inserting the summarized values for NDVI and LST in the equations from Table 4, we can obtain crop yields for 2014.The yield of corn seed was not calculated due to lack of crop planting structure for one year.It is very difficult to distinguish seed and mercantile corn based on only spectral characteristics.Based on spectral-temporal characteristics of NDVI this is possible (Figure 9).It should be noted that these are early estimates of yield and that by increasing the sampling areas, i.e. supplying ground trough data from other conglomerates in Vojvodina results can be improved.Also the increasing number of monitored years can improve results.Using the average yields for 2014 and total harvested area of crops, total yields for the current year can be calculated.
In order to assess the yields, this paper discussed a method that predicts the behaviour of the plants for the end of 2014.We can conclude that this is a good support of remote sensing for yield assessment, due to the fact that images cover very large territory and because of direct relationship with the condition of plants throughout the years.Crop yield results are based on the sampled areas "Nova Budućnost" and "Sava Kovacevic", which certainly have better yields due to better cultivation of crops than conventional farmers.Also, we used images of high and medium resolution combined with accurate plantings for the period from 2011 to 2014.To minimize the impact of "good crops", sampling area was extended as much as possible, but the parcels had to be large in order to detect them in the images with 250 meters resolution.Also, in order to improve the determination, sampling pixels belonged to a particular plant, i.e. mixed pixels were avoided.After that, phenology for each plant was identified and NDVI and LST value were summarized for the seasonal period so that we could be able to estimate yields using multiple linear regression.The contribution of this study is to support early crop yield estimate in Vojvodina by methods of remote sensing together with ground truth data.This method can be used for different geographical areas, not only Serbia.By monitoring phenology it was proven that it is possible to discriminate all cultivated crops throughout all of the monitored years.Furthermore, it was noted that extreme seasons can be noticed easily with NDVI phenology, season 2012.Also, the behaviour of crops was simulated because of early crop yield predictions.Regarding estimated crop yields we conclude that this method demands more sampling areas (more pixels in analysis).The weakness regarding directly to the MODIS LST is in very low spatial resolution (1 km).Sampling areas for LST must be picked carefully in order to avoid water areas or mixed pixels (roads, urban areas).This weakness can be bypassed if LST from Landsat platforms is used.The quality of MODIS images is critical to this method; the algorithm method uses Maximum Value Composite filtering to obtain a composite NDVI.For a period of 16 days, this algorithm takes into account only pixels with the highest value in order to eliminate the interference of clouds.Based on the fact that MODIS produces the best possible products which are geometrically, radiometrically and atmospherically corrected (Huete et al. 2002, Vermote et al. 2002, Wolfe et al. 2002), we can conclude that MODIS can provide exact results if a proper method is used and if proper number of sampling areas is included.Results presented in this paper for average yield crops for 2014 are not compared with the results of the Statistical Institute, because there is no assessment of the yield before harvest.The results will be compared when the official results of the Statistical Office are published.
Figure 1.Vojvodina's location map with study areas of Vrbas (a) and Ruma (b).Study areas are shown with RapidEye images, true colour display Figure 2. Example of good and bad selection of training polygon.
Left part of Figure 2 shows the RapidEye image overlaid with MODIS's grid (250 m) with good (dotted purple) and bad (solid red) polygon selection.Right part of Figure 2 represents the correspondent MODIS pixels.

Figure 3 .
Figure 3. Example of spectral match of two same crops, soybean, near Žarkovac.Landsat 8 OLI in 6 5 4 false colour display

Figure 4 .
Figure 4. Phenophases example of wheat and maize (b, e) and RapidEye and OLI images for May, June and August (a, c, d, f)

Figure 5
Figure 5 illustrates period of growth, the aging period, the beginning and the end of green vegetation, as well as the entire length of the season, which is equal to the sum of the growth period and the aging period.The following values are obtained:  the beginning of green vegetation= 65th day of the year (on March 5)  the end of green vegetation = 177th day of the year (on June 25)  the length of the season = 177-65 = 112 days  growth period = 64 days  aging period = 48 days  maximum NDVI value is on 129th day of the year

Figure 8 .
Figure 8. Spectral-temporal profile for mercantile corn in 2014, current and predicted values of NDVI

Figure 9 .
Figure 9. Monitored phenology of crops in 2014 Conclusion This paper presents method for yield estimating by using parameters obtained from the MODIS sensor, NDVI and LST (T) data for the territory of Vojvodina for 2014.Crop yield results are based on the sampled areas "Nova Budućnost" and "Sava Kovacevic", which certainly have better yields due to better cultivation of crops than conventional farmers.Also, we used images of high and medium resolution combined with accurate plantings for the period from 2011 to 2014.To minimize the impact of "good crops", sampling area was extended as much as possible, but the parcels had to be large in order to detect them in the images with 250 meters resolution.Also, in order to improve the determination, sampling pixels belonged to a particular plant, i.e. mixed pixels were avoided.After that, phenology for each plant was identified and NDVI and LST value were summarized for the seasonal period so that we could be able to estimate yields using multiple linear regression.The contribution of this study is to support early crop yield estimate in Vojvodina by methods of remote sensing together with ground truth data.This method can be used for different geographical areas, not only Serbia.By monitoring phenology it was proven that it is possible to discriminate all cultivated crops throughout all of the monitored years.Furthermore, it was noted that extreme seasons can be noticed easily with NDVI

Table 1 .
Monitored years for each crop, together with total number of MODIS pixels and area

Table 2 .
Equation for predicting the 2014 and the coefficient of determination R 2

Table 3 .
Average yields for Vojvodina (source-Republic Institute for Statistics)