Cropland area estimates using Modis NDVI time series in the state of Mato Grosso , Brazil

The objective of this work was to evaluate a simple, semi-automated methodology for mapping cropland areas in the state of Mato Grosso, Brazil. A Fourier transform was applied over a time series of vegetation index products from the moderate resolution imaging spectroradiometer (Modis) sensor. This procedure allows for the evaluation of the amplitude of the periodic changes in vegetation response through time and the identification of areas with strong seasonal variation related to crop production. Annual cropland masks from 2006 to 2009 were generated and municipal cropland areas were estimated through remote sensing. We observed good agreement with official statistics on planted area, especially for municipalities with more than 10% of cropland cover (R2 = 0.89), but poor agreement in municipalities with less than 5% crop cover (R2 = 0.41). The assessed methodology can be used for annual cropland mapping over large production areas in Brazil.


Introduction
Given the large extent of Brazilian agricultural activity and its economic importance in both national and international markets, the development and implementation of systematic crop area monitoring and mapping tools is of fundamental importance for the country.This system could be used to supply precise and up-to-date information for many agricultural related topics, such as commodity price definitions, crop forecasts, and the establishment of public policies.Moreover, systematic crop area monitoring could improve the acceptance of Brazilian products in the international market, since it would supply information on the environmental impacts of agricultural production.
Currently, annual Brazilian cropland statistics are based primarily on municipal-level data of the Instituto Brasileiro de Geografia e Estatística (IBGE), which is in contact with the main producers and agricultural Pesq.agropec.bras., Brasília, v.47, n.9, p.1270-1278, set.2012 agents of each municipality.This process requires an intricate interaction among several institutions, from local and municipal to state and national levels, in order to obtain valuable statistics, such as planted area, crop production, and monetary value for different crops at the municipal level (Produção agrícola municipal, 2010).Cropland mapping through remote sensing could complement official statistics, generating annual cropland maps along with planted area estimates.
The use of remote sensing for mapping land cover is well established in the literature, be it to evaluate natural vegetation cover, forest degradation and deforestation, cropland expansion, or land cover change (Pickup et al., 1993;Lu et al., 2003;Souza et al., 2003;Asner et al., 2004;Eva et al., 2004;Kastens et al., 2005;Rizzi & Rudorff, 2005;Lunetta et al., 2006;Panta et al., 2008).Many applications are based on moderately high spatial resolution images, such as the Landsat TM with 30 m pixels, and are restricted to relatively small areas or demand a large amount of work.The moderate resolution imaging spectroradiometer (Modis) produces near-daily images with 250 m spatial resolution, suitable for identifying large crop fields in regions with widespread use of mechanized agriculture (Wardlow et al., 2006;Ren et al., 2008).Composite images of worldwide Modis time-series are freely distributed by the U.S. Geological Survey (USGS) through the Land Processes Distributed Active Archive Center (LP-DAAC) with all necessary radiometric, atmospheric, and geometric corrections applied.Modis products are, therefore, well suited for monitoring land cover over large areas on an annual basis.
Much research has been carried out using time series of remotely sensed vegetation index (VI) data to detect land-use and land-cover change.Pioneering works relied on advanced very high resolution radiometer (AVHRR) products, but recently Modis has been used more often.Different techniques are used to perform VI time series analysis, such as Fourier transforms (Jakubauskas et al., 2002), Savitzky-Golay filtering (Jönsson & Eklundh, 2004), wavelet transforms (Sakamoto et al., 2005), mathematical functions for filtering and smoothing (Zhang et al., 2003), and probability distribution-based analysis (White & Nemani, 2006).Other researchers have focused on the use of spectral-temporal response surfaces, analyzing not only the VI time series but also the temporal variation of the different spectral bands in order to identify cropland areas (Rudorff et al., 2007;Epiphanio et al., 2010), or on the use of spectral linear mixture models applied to time series for vegetation cover identification (Silva et al., 2010).
The development of a systematic cropland mapping program in Brazil is challenging because of the high diversity of agricultural production systems at different technological stages, ranging from traditional small farmers to highly-mechanized commodity crop production operations.However, for some specific agricultural production areas where large-scale mechanized row-crop agriculture dominates the landscape, the development of cropland mapping tools is feasible.
In Mato Grosso, previous efforts related to VI time series analysis have focused on identifying the fate of recent forest clearings (Morton et al., 2006), mapping soybean cropland (Epiphanio et al., 2010), intensification of agricultural practices (Brown et al., 2007;Galford et al., 2008), and discrimination of Cerrado vegetation cover from natural or managed pastures (Silva et al., 2010).
The objective of this work was to evaluate a simple, semi-automated methodology for mapping cropland areas in the state of Mato Grosso, Brazil.

Materials and Methods
The state of Mato Grosso is located in the Central-West region of Brazil, with an area of 903 thousand km 2 (18 o to 8 o S and 62 o to 50 o W).The state is known for its highly mechanized agriculture and is responsible for 27% of soybean and 46% of cotton planted areas in Brazil (Companhia Nacional de Abastecimento, 2011).In 2009, soybean accounted for 67.8% of all planted area in the state, followed by corn (19.3%), and cotton (4.1%) (Produção agrícola municipal, 2010).Cropland masks were estimated for the 141 municipalities of the state.
Time series of 16-day composite Modis NDVI (normalized difference vegetation index) images with 250 m spatial resolution (Huete et al., 2002), downloaded from the USGS LP-DAAC, were used for the classification.Overall, the procedure consisted of obtaining and organizing the images, applying a simple local minimum filter, performing a Fourier transform to obtain the amplitude images of the Pesq.agropec.bras., Brasília, v.47, n.9, p.1270-1278, set.2012 harmonic components, performing an unsupervised classification, and extracting crop areas.
Images from the Modis product MOD13Q1 version 5, which provides 16-day composite VI data, were obtained for tiles h12v10, h13v10, h11v09, h12v09, and h13v09, from 2005-2009.The NDVI bands for each date and tile were extracted, mosaicked, and clipped to the extent of the state.The mosaics were then stacked according to the crop calendar, beginning July 12 of the previous year (day of the year, DOY,193)  Though the 16-day composite Modis NDVI data are largely free from cloud contamination, some problems of this nature cannot be completely eliminated and appear as strong negative spikes in the NDVI time series.To mitigate these effects, a simple filtering method, similar to that of Wardlow et al. (2006), was applied, substituting all local minima in the time series with the closest value that occurred immediately before or after that observation.
All pixel-level NDVI profiles from the multitemporal image stacks were analyzed using the Fourier transform, which converts a signal measured in the time domain to a signal in the frequency domain (Sakamoto et al., 2005).The Fourier transform decomposes a signal into an additive constant term (which is the mean value of the signal) and into a series of harmonic terms (cosine waves), each having a unique amplitude and phase value (Jakubauskas et al., 2002).The first harmonic of a Fourier transform is a sinusoid with one complete cycle over the length of the time series, whereas the second harmonic is a sinusoid with two complete cycles, and so on.
A Fourier transform was applied to each multitemporal image stack, resulting in amplitude images for the first two harmonic components, along with mean annual NDVI (component 0).The first harmonic of the NDVI time series represents the annual variation in the signal and is captured by component 1.The second harmonic (component 2) reflects the six-month (biannual) variation present in the NDVI signal, whereas the third harmonic reflects the variation on the four month scale, and so on.The images were then classified using the Isoseg unsupervised classification algorithm (Richards & Jia, 1999), which was chosen due to its simplicity and availability in most remote sensing analysis software.
The resulting classes were analyzed based on their spatial patterns in order to determine which classes corresponded to cropland areas.
Cropland area estimates obtained using Modis NDVI were compared to official government estimates of crop planted area from IBGE, for the 141 municipalities of Mato Grosso, from 2006 to 2009.The IBGE estimates are part of the Municipal Agricultural Production ("Produção Agrícola Municipal", PAM) survey, which is based on interviews with local technicians, producers, farmers' associations, and others.Municipal-level data are officially released in November of each year and can be obtained from the IBGE web site (Instituto Brasileiro de Geografia e Estatística, 2012).Total cropland area reported by IBGE takes into consideration the fact that many crops, especially corn, are planted as double crops.In fact, in 2010, 93% of the corn in the state was planted in the second growing season, which accounts for 20% of all cultivated cropland in Mato Grosso, and 89% of the beans were planted in the second or third growing season.Therefore, for this comparison, cropland area for each municipality was estimated as the sum of soybean, first season corn, first season bean, and cotton planted areas.

Results and Discussion
Cropland areas in Mato Grosso are dominated by large continuous fields, managed using highly mechanized agricultural practices, resulting in a land cover mosaic of large, uniform objects that are frequently discernible in Modis NDVI images.These cropland areas frequently follow a double cropping pattern, with soybean as a first crop, followed by a second crop, such as millet for land cover, or corn or cotton.This results in a distinct VI time series pattern relative to other covers, such as pasture, forest, or Cerrado (Figure 1).
Amplitude images for harmonic components 0, 1, and 2 were extracted (Figure 2).The mean annual NDVI captured by the zero component is high (light gray) in areas with dense vegetation cover throughout most of the year (e.g., gallery forests and Cerrado with dense vegetation).Pasture and Cerrado have lower values Pesq.agropec.bras., Brasília, v.47, n.9, p.1270-1278, set.2012 than those of areas with dense, perennial vegetation cover, but not as low as croplands (dark gray to black), which are covered with sparse or no vegetation during multiple months outside of the cropping season.
Annual variation, which is captured by the first harmonic, is low (dark gray to black) in some areas, such as forests, forested Cerrado, and also in open waters and urban areas (Figure 2).Vegetation greenness in Cerrado and pasture follows more closely the seasonal cycle of precipitation and, therefore, has higher (light gray) annual variation than forest, showing intermediate values of amplitude in the first harmonic.As for croplands, the short phenological cycles of commercial and cover crops are fairly well represented by the sum of harmonic components 1 and 2. This is because the soybean cropping patterns (either double cropping or single crop followed by cover) in Mato Grosso have both strong annual and semestral variation.The stronger amplitude of component 2 in some croplands could be an indication of areas where the second-season crop has a more vigorous green cover, which could potentially be used to distinguish cover crops (millet) from the commercial second crop (corn or cotton).A false color composite, using average annual NDVI and the first two harmonic components, is also shown to illustrate how cropland (bright magenta) is readily separable from other vegetation covers.
Comparison of different classification results indicated that using ten classes was adequate, capturing the different land covers of the area without fragmenting the cropland area.When comparing the classified maps to higher resolution remote sensing images, it was possible to identify which class represented cropland, allowing for the creation of annual cropland masks.A visual inspection of the 2009 cropland mask showed good agreement with crop areas identified in higher resolution images (Figure 3), though a more thorough spatial evaluation is desirable.However, total cropland area estimated through remote sensing showed good agreement with the reported planted area.
Percentage of cropland areas for municipalities in Mato Grosso for each of the 2006-2009 crop years, estimated using Modis NDVI time series, showed good agreement with percentage of planted area calculated using statistics obtained by IBGE {(soybean + first season corn + first season bean + cotton) / municipal area}, and a significant linear regression was observed (significant F<0.01%, R 2 = 0.96).However, a closer inspection revealed that the cropland area estimates obtained using Modis had the best agreement with IBGE estimates for municipalities with large planted areas (R 2 = 0.89).For comparison, determination coefficients for municipalities with crop covers ranging from 5 to 10% and below 5% are 0.12 and 0.41, respectively (Figure 4).This indicates that, despite the significant linear regression found with all municipalities, the method is more appropriate for regions with greater densities of cropland.
For regions with low cropland fractions, in particular those where agricultural fields are smaller, the moderate spatial resolution of the pixel (250 m) can result in mapping errors due to the occurrence of mixed vegetation covers.In addition, errors in the time series, such as sensor problems or cloud contamination not screened out by the filtering procedure, could lead to the incorrect identification of croplands (false positives).The local minima filtering procedure applied to the images proved to be a simple method for reducing cloud contamination in the time series.Several other filtering methods are available, some of which were reviewed by Hird & McDermid (2009).Those authors found that methods based on the double logistic and asymmetric Gaussian function-fitting techniques are superior to other methods tested, though caution should be taken when using noise reduction procedures, since they may influence the information being extracted from the time series.The local minima method used  Other studies identifying cropland areas using Modis, which mainly focused on soybeans, have shown that the moderate spatial resolution of the sensor is adequate for the identification of large crop fields; however, the classification of smaller plantations is less efficient (Rudorff et al., 2007;Epiphanio et al., 2010).Working specifically over the state of Mato Grosso, Epiphanio et al. (2010) used the spectral-temporal response surface method to successfully identify soybean plantations.However, the procedure proved to be highly computationally intensive and its application over an entire Modis tile was considered impracticable (Epiphanio et al., 2010).The method evaluated in the present study, based on Fourier transforms, has the advantage of being simple and fast, not demanding great computational resources.As indicated by IBGE estimates, presented in Table 1, there are close to 30 municipalities in Mato Grosso with more than 10% of their area covered by croplands, which accounts for approximately 80% of total planted area in the state.By applying the described method, it was possible to estimate cropland areas in these municipalities, with results showing good agreement with official cropland statistics (significant F <0.01%, R 2 = 0.89).Furthermore, the method is able to map cropland areas and is an important tool to evaluate the spatial dynamics of croplands over large regions.Pesq.agropec.bras., Brasília, v.47, n.9, p.1270-1278, set.2012 Conclusions 1.The Fourier transform of the normalized difference vegetation index (NDVI) time series is able to extract the strong seasonal variations related to crop phenology and management.
2. Modis-NDVI time series can be used to identify croplands over the state of Mato Grosso, Brazil, and to estimate total cultivated area, with results comparable to those reported by the official Instituto Brasileiro de Geografia e Estatística (IBGE).
3. Modis based estimates of croplands show better agreement with IBGE data for municipalities with large cultivated areas, but poor agreement in municipalities with low cropland densities.

Figure 1 .
Figure 1.Temporal pattern of the normalized difference vegetation index (NDVI) for different land covers.NDVI values for each land cover were extracted from pixels where land cover was identified in the field.
will only create new values at down spikes, preserving most of the original information in the time series.

Figure 2 .
Figure 2. Mean annual normalized difference vegetation index (NDVI) for the state of Mato Grosso, Brazil, in 2009 (A); annual amplitude obtained from the first harmonic component (B); semestral amplitude from the second harmonic component (C); false color composite using the three images with component 2 in red, component 0 in green, and component 1 in blue (D).Lighter colored areas in A through C indicate higher amplitude values.Cropland areas appear as bright magenta in D.

Figure 3 .
Figure 3. Cropland mask for crop years 2006-2009, along with close-up views of the 2009 cropland mask, compared to higher resolution images from ArcGIS World Imagery map service (15 m spatial resolution, acquisition date not specified).

Figure 4 .
Figure 4. Moderate resolution imaging spectroradiometer (Modis) estimates of municipal crop cover (percentage of area) against estimates from the Municipal Agricultural Production (PAM) survey carried out by the Brazilian Institute of Geography and Statistics for all municipalities, municipalities with more than 10% crop cover, between 5 and 10%, and with less than 5% cropland in the 2006 crop year.
to July 11 of the nominal crop year (DOY 192); for example, crop year 2006 extends from 7/12/2005 to 7/11/2006.This resulted in four multi-temporal image stacks with 23 bands each, covering the 2006-2009 crop years.

Table 1 .
Soybean crop statistics for the state of Mato Grosso, Brazil(1).