Soybean crop area estimation by Modis / Evi data

The objective of this work was to develop a procedure to estimate soybean crop areas in Rio Grande do Sul state, Brazil. Estimations were made based on the temporal profiles of the enhanced vegetation index (Evi) calculated from moderate resolution imaging spectroradiometer (Modis) images. The methodology developed for soybean classification was named Modis crop detection algorithm (MCDA). The MCDA provides soybean area estimates in December (first forecast), using images from the sowing period, and March (second forecast), using images from the sowing and maximum crop development periods. The results obtained by the MCDA were compared with the official estimates on soybean area of the Instituto Brasileiro de Geografia e Estatística. The coefficients of determination ranged from 0.91 to 0.95, indicating good agreement between the estimates. For the 2000/2001 crop year, the MCDA soybean crop map was evaluated using a soybean crop map derived from Landsat images, and the overall map accuracy was approximately 82%, with similar commission and omission errors. The MCDA was able to estimate soybean crop areas in Rio Grande do Sul State and to generate an annual thematic map with the geographic position of the soybean fields. The soybean crop area estimates by the MCDA are in good agreement with the official agricultural statistics.


Introduction
The globalized market of agricultural commodities requires precise and timely information on crop production for decision makers.The agricultural production of Brazil has an important role in the world commodities market.Therefore, it is relevant to develop objective methods that can provide precise and timely estimates of crop production.The Companhia Nacional de Abastecimento (Conab) and the Instituto Brasileiro de Geografia e Estatística (IBGE) are responsible for the country's official agricultural statistics.However, these estimates are based on subjective methods.
In Brazil, remote sensing satellite data have already shown adequate performance for some routine mappings of perennial and semi-perennial crops, such as sugarcane (Rudorff et al., 2005(Rudorff et al., , 2010) ) and coffee (Moreira et al., 2010).However, the high incidence of cloud cover during key identification periods of annual crops on Landsat-type images has hindered Pesq.agropec.bras., Brasília, v.47, n.3, p.425-435, mar. 2012 its operational use in providing agricultural statistics for those crops (Sano et al., 2007;Sugawara et al., 2008).A potential solution to overcome cloud cover is to increase the temporal resolution of orbital sensors, but there may be a reduction in spatial resolution due to technical constraints.Currently, the moderate resolution imaging spectroradiometer (Modis) sensor, on board of the Terra and Aqua satellites, combines an almost-daily revisit with a moderate spatial resolution of 250 m, which are favorable characteristics to map crops, such as soybean, cultivated in large scale.In addition, the geometric quality of the images allows the composition of time series data that guarantee the correct pixel geolocation (Justice et al., 2002;Wolfe et al., 2002).
Many studies have explored the potential of Modis imagery for agricultural crop surveys and monitoring.Lobell & Asner (2004) evaluated the impact of land use heterogeneity and soil cover on the classification of agricultural crop areas with Modis data and observed that the mean square error (MSE) decreases as the field size increases, tending to stabilize for fields larger than 500 ha.These authors also concluded that Modis images have considerable advantages over Landsat ones in the characterization of extensive agricultural crops, mainly due to their higher temporal resolution.Sakamoto et al. (2005) used time series of the enhanced vegetation index (Evi), obtained from Modis images, to monitor the phenology of rice crop in Japan and observed that the phenological stages can be determined within a mean error of 9 to 12 days.Doraiswamy et al. (2005), while assessing the quality of Modis data to provide information on both crop yield and area, estimated biophysical parameters from Modis images that were integrated into crop growth models, providing significant improvement in grain yield estimates.These authors also found that the overall mapping accuracy was 96.7%, when compared to a Landsat map.In the USA, Wardlow et al. (2007) investigated the applicability of Modis/Evi time series data to map agricultural lands and concluded that the 16-day composites of Modis images gave sufficient spatial, spectral, and temporal information to adequately separate crop fields from other land uses, and to express the phenology and climate characteristics of the region.
In Brazil, specifically in Rio Grande do Sul State, Rudorff et al. (2007) evaluated the potential of Modis images to identify and map soybean crops using the temporal-spectral response surface method developed by Vieira (2000).These authors observed that the mapping accuracy was dependent on the mean field size, in agreement with the results obtained by Lobell & Asner (2004).Also in Brazil, Epiphanio et al. (2010) used Modis in the same temporal-spectral response surface approach (Vieira, 2000) for mapping soybean in Mato Grosso State, in a region of plain relief and large fields, and found overall accuracy of 80%.D' Arco et al. (2007) used Modis vegetation indices to classify rice crop in the municipality of Santa Vitória do Palmar, RS, Brazil, and obtained best results with images from three specific crop growth periods.However, none of these studies have yet provided timely soybean area estimates.Furthermore, most of the researches are limited to only one crop year or region, indicating the potential of Modis data for crop forecast, but not its usefulness within a routine and systematic crop forecast system.
In this sense, the relative long-term Modis/Evi time series data now available, with its valuable spatial information on spectral-temporal crop growth profile, could be used to estimate cultivated soybean area for local and regional-scale operational crop forecasts.The Modis sensor on board of the Terra satellite completed, on 12/18/2009, a ten-year history of systematic global data acquisition with high-temporal, moderate-spatial, and improved spectral resolution (Justice et al., 2002), as well as correct geolocation of image pixels (Wolfe et al., 2002).
The objective of this work was to develop a procedure to estimate soybean crop areas in Rio Grande do Sul State, Brazil, based on temporal profiles of the Evi calculated from Modis images.

Materials and Methods
Soybean production was evaluated in Rio Grande do Sul state, Brazil, using a validation area of 223 municipalities in the north of the state (Figure 1).The region's climate is subtropical with four well-defined seasons.The average annual rainfall is 1,500 mm, relatively well distributed throughout the year, but with some dry periods.The sowing calendar for soybean goes from early October to late December, based on agricultural zoning for different soils, regions, and cultivars (Cunha et al., 2001).Depending on the sowing date, maximum plant growth is observed from late January to early March.
Pesq. agropec.bras., Brasília, v.47, n.3, p.425-435, mar. 2012 Three test sites located in intensive soybean production areas were used to define the sowing and maximum growth periods (Figure 1), which are key information to extract the correct parameters from the Modis/Evi values for the soybean area estimation procedure proposed in the present study.
Different types of data were used: yearly rainfall data (2000 to 2009), determined by the ten-day accumulated rainfall data of 16 meteorological stations, from September to October, from the Fundação Estadual de Pesquisa Agropecuária, which were used to refine the period of initial sowing; shuttle radar topography mission (SRTM) data, to generate a slope map with a 90-m spatial resolution, according to Rabus et al. (2003), in order to exclude areas improper for mechanization (slope >15%); a soil map (Instituto Brasileiro de Geografia e Estatística, 1986) in a 1:250,000 scale, to exclude the unfavorable soils for soybean cultivation (all Planosol and Plinthosol; BT1, BT2, BT3, BT4, BT5, BT6, BT7, and BT9 sub-classes of Chernosol; Ge, GHS2, GS1, GS2, GHe1, GHe2, and GHe3 sub-classes of Gleisol; Re1, Re2, Re6, Re7, Re9, Ae1, and Ae2 sub-classes of Neosol; and V1 and V3 sub-classes of Vertisol); annual soybean agricultural statistics, in state and municipality level (Instituto Brasileiro de Geografia e Estatística, 2009) for the entire study area, which were used to compare and evaluate the results obtained from the proposed soybean area estimation procedure; a soybean reference thematic map, available for the 2000/2001 crop year (Rizzi & Rudorff, 2005), obtained from multitemporal TM/Landsat images analysis, in a 30-m spatial resolution, to evaluate the Modis/Evi soybean thematic map developed in this procedure; Evi data from 2000 to 2009, derived from the Modis sensor on board of the Terra satellite, product MOD13Q1-collection 5; and TM images, obtained from the Instituto Nacional de Pesquisas Espaciais (2010).
The Evi data were chosen due to their potential ability to reduce atmospheric and soil background effects (Huete et al., 2002;Justice et al., 2002).The Evi is part of the MOD13Q1-V005 product, which comprises the best radiometric and geometric pixels selected within a 16-day period.The Evi formulation is 2.5(Nir -Red)/ (Nir + 6Red -7.5Blue + 1), in which: Nir, Red, and Blue are atmospherically or partially-atmospherically corrected surface reflectance of near infrared, red, and blue bands, respectively (Huete et al., 2002).The Modis images and its products were preprocessed by the National Aeronautics and Space Administration (Nasa) and are available at no charge at the National Aeronautics and Space Administration (2009).The SRTM data were used to generate the land slope variable in order to exclude all areas with slope greater than 15% from the analysis, since soybean is a highly mechanized crop and requires relative smooth land to allow the traffic of farm implements.Areas from hydromorphic soils were also excluded from the analysis because they are mainly recommended for flooded rice cultivation (D'Arco et al., 2007).
In the present study, the procedure developed for soybean classification was named Modis crop detection algorithm (MCDA), for which a diagrammatic flowchart was created (Figure 2).According to this procedure, a pixel is classified as soybean if it adheres to both A and B conditions.To establish these two conditions, the soybean Evi temporal profile is expected to have low values during the sowing period and high values at maximum crop development.
The set of initial parameters was obtained from the Modis/Evi images of three test sites of 100x100 pixels (Figure 1), during two specific periods: sowing (day of year, DOY, 273 to 337) and maximum crop development (DOY 1 to 65).The sowing period for each crop year was also defined based on rainfall data, in agreement with the soybean zoning provided by Brasil (2010).In the study area, the sowing period normally starts in October.
Two Modis/Evi images from the sowing period were averaged to obtain the minimum mean Evi image (Figure 3), which was used to define the lower (Lmin) and upper (Umin) thresholds of the Evi values for each crop year.Pixels below Lmin are typically associated to cloud shadows or water bodies.Umin was set at the convergence between minimum and maximum mean Evi images (Figure 3).Pixels above Umin will start to  1. Pixels with Evi values between Lmin and Umin were tagged as soybean, according to the A condition (Figure 2).
Following the sowing period, a rapid increase of the Modis/Evi values was observed due to intense plant growth, reaching maximum values in a relative short period (Wardlow et al., 2007).Three consecutive Evi images from the plant growth period (DOY 1 to 65; Figure 4) were used to obtain the maximum mean Evi image.The difference between the maximum and minimum mean Evi images was used to obtain the Evi amplitude image (Amp).The challenge is to obtain the best Amp value that includes not only pure soybean pixels with high values in the maximum Evi image and low values in the minimum Evi image, but also mixed pixels located at the border of soybean fields.The Amp value for each crop year was chosen at the convergence region between minimum and maximum Evi values in the scatterplots (Figure 3), i.e., the minimum difference between maximum and minimum mean Evi values to which a mixed soybean pixel can be tagged as soybean.The initial nine selected Amp values from 2000/2001 to 2008/2009 and the average value used in the MCDA procedure are shown in Table 1.Pixels with amplitude values greater than Amp were tagged as soybean, according to the B condition of the MCDA procedure (Figure 2).
After the first run of the MCDA procedure, Landsat satellite images (Table 2), in which soybean fields are easily identified, were used to interactively refine the most appropriate values of Lmin, Umin, and Amp in order to minimize omission and commission errors.
Since the soybean area was overestimated for the first run of the MCDA, when the Landsat images and the classification maps were overlapped, a second run was performed, adjusting the values of Umin and Amp.For each combination of Amp and Umin, new soybean classifications were generated, which were visually compared with the correspondent available Landsat images (Table 2).After several interactions, for several crop years, and comparing the results of each new MCDA classification faced with Landsat images, the combination with best performance to define the final values of the MCDA were 0.05, 0.47, and 0.21 for Umin, Lmin, and Amp, respectively.It should be emphasized that the parameters defined in the MCDA are the same for all of the analyzed crop years.
Soybean area estimation can be provided right after the maximum mean Evi image is available, which normally occurs in early February.A delay of about 20 days is expected in order to acquire the MOD13Q1 product, and, therefore, soybean estimation should be released no later than early March (second forecast of the MCDA).However, as an alternative to forecast the soybean area, a first estimate can be provided in mid-December of each crop year, based on the minimum mean Evi image of the current crop year and the maximum mean Evi image of a previous crop year (first forecast of the MCDA).If no usable images are found or a water deficit is observed (30 days without rainfall above 10 mm), images from previous normal crop years can be used to generate the maximum mean Evi image.For the 2000/2001 crop year, there was no first forecast of the MCDA due to unavailable Modis data.In order to generate the maximum mean Evi image in each crop year, 16-day Modis/Evi composites were selected (Table 3).The mean temporal Evi was determined for profiles of soybean, corn, and rice acquired over crop fields, which coexist in the evaluated area (Figure 4).
Soybean area was also estimated by municipality from 2000/2001 to 2008/2009 and compared to the official estimates provided by the Instituto Brasileiro de Geografia e Estatística (2009), using regression analysis.The performance of the MCDA was evaluated using a thematic soybean map elaborated from multitemporal Landsat images for the 2000/2001 crop year (Rizzi & Rudorff, 2005).The confusion matrix and the overall map accuracy were provided.

Results and Discussion
The MCDA soybean area estimates for the first and second forecasts agreed well with the IBGE estimates, with a maximum difference of 5.5%, from 2001/2002 to 2008/2009 (Figure 5).The MCDA soybean area estimates provided at a municipal level, due to the spatial distribution of the classified soybean area, were compared to the IBGE municipal statistics, which are only published about a year after the end of the soybean season.No early estimate could be provided for the 2000/2001 crop season, as this was the first crop year in which the Modis data became available.Therefore, compared to the current official methods for soybean area estimation in Brazil, the MCDA procedure represents a considerable improvement in time and in spatial information.
The linear least squares regression analysis was done for the municipal soybean estimates obtained from the MCDA (dependent variable) and from IBGE (independent variable) for the 2000/2001 to the 2006/2007 crop years, since the municipal information for 2007/2008 and 2008/2009 had not yet been released by IBGE (Table 4).The coefficients of determination (R 2 ) ranged from 0.91 to 0.94, indicating good agreement between the estimates.The test of b 0 = 0 indicated that b 0 was significantly different from zero (α = 0.05) for all tested crop years.Generally, this value was around -700 ha, indicating that municipalities with small soybean area are underestimated by the MCDA in relation to the IBGE statistics.The test of b 1 = 1 indicated that b 1 was significantly different from 1 (α = 0.05) for all tested crop years, with most values around 1.10, which indicates that the MCDA overestimated the soybean area in relation to IBGE for municipalities with large soybean areas.The root mean square error (RMSE) for both forecasts of the MCDA, in all crop years, was around 4,000 ha, indicating that the MCDA estimates are consistent.Chang et al. (2007)   Pesq.agropec.bras., Brasília, v.47, n.3, p.425-435, mar. 2012 in greater detail by Rizzi & Rudorff (2005), due to the much better spatial resolution of Landsat in comparison to Modis images (Figure 6 D).In fact, fields with smaller areas were more subject to errors than those with larger areas, in agreement with Lobell & Asner (2004) and Rudorff et al. (2007).
The comparison between the mapping provided by Rizzi & Rudorff (2005) and the MCDA procedure generated with Modis/Evi data, in the second forecast, resulted in a confusion matrix (Table 5).Despite the moderate spatial resolution of the Modis and the soybean field size heterogeneity, the overall map accuracy was approximately 82%, which indicates good agreement between those maps.The user's (commission) and producer's (omission) accuracy for the soybean class are similar, which indicates that these errors tend to cancel each other and that the estimates of total areas by municipality are near the actual value.Similar results were obtained by Rudorff et al. (2007) in the same study area.
It is important to point out that the parameters defined in the MCDA for the detection of crop areas are constant, as a fixed criteria, independently of the soybean crop area increases and of the dynamics of the crop years in Rio Grande do Sul State during the evaluated period, from the 2000/2001 to the 2008/2009 crop years.Pesq.agropec.bras., Brasília, v.47, n.3, p.425-435, mar. 2012 The MCDA can also be applied to corn and rice crops, but further adaptations to their respective agricultural calendars, parameters, and variables are necessary.In addition, soybean field size caused considerable differences in the classification results.Therefore, the MCDA should generate better results over the soybean producers of the Brazilian states where crop fields are much larger than the ones in the evaluated area (Lobell & Asner, 2004;Rudorff et al., 2007).

Conclusions
1.The Modis crop detection algorithm (MCDA) procedure is based on an adequate and objective methodology for estimating soybean crop areas, using moderate resolution imaging spectroradiometer (Modis) images and enhanced vegetation index (Evi) data.
2. The MCDA approach can assist the official soybean crop area estimates, providing reliable spatialized information, prior to the crop harvest period.

Figure 1 .
Figure 1.Rio Grande do Sul state, Brazil, and its 497 municipalities, indicating the validation area (223 municipalities) and the three test sites.

Figure 2 .
Figure 2. Flowchart of the MCDA classification based on the Modis/Evi images.

Figure 3 .
Figure 3. Scatterplot of 2000/2001 with minimum and maximum Evi values from the three test sites (30,000 pixels), showing symmetry and intersection lines to estimate the initial values of upper (Umin) and lower (Lmin) thresholds, and amplitude image (Amp), sorted downhill by difference.
compared soybean area estimates generated by Modis and by the National Agricultural Statistics Service of the United States Department of Agriculture and observed R 2 values ranging from 0.44 to 0.94 and RMSE varying from 41,465 to 120,955 ha for the entire USA.When comparing the MCDA estimates with the results obtained by Rizzi & Rudorff (2005), which used Landsat satellite images to estimate the soybean area in Rio Grande do Sul, Brazil, for the 2000/2001 crop year, it was noted that the two estimates are in good agreement with an overestimation of 5% by the MCDA.The soybean map obtained by Rizzi & Rudorff (2005) (Figure 6 A) and by the MCDA for the second forecast, for the 2000/2001 crop year (Figure 6 B), were overlapped (Figure 6 C and D), indicating the areas of agreement (green), omission (black) and commission errors (red).The major errors are associated with the field boundaries that were mapped

Figure 5 .
Figure 5.Comparison between the MCDA (first and second forecasts) and the IBGE soybean area estimates for Rio Grande do Sul State, Brazil.

Figure 6 .
Figure 6.Soybean map of the 2000/2001 crop year obtained by Rizzi & Rudorff (2005) (A) and by the MCDA procedure for the second forecast (B); validation between A and B (C); and enlargement of a validation sector (D).

Table 1 .
Annual values of lower (Lmin) and upper (Umin) thresholds, and amplitude image (Amp) extracted from the minimum and maximum Evi values from the three test sites (30,000 pixels), and average values initially used in the MCDA procedure.

Table 2 .
Dates of selected Landsat images for visual inspection during the interactive fitting phase.

Table 3 .
Selected Modis/Evi for each crop year.

Table 4 .
Regression analysis between the MCDA procedure forecasts and the IBGE data for the 2000/2001 to the 2006/2007 crop years.

Table 5 .
Confusion matrix resulting from the comparison between the MCDA and the Rizzi & Rudorff (2005) maps.