MAPPING OF SUMMER AGRICULTURAL CROPS IN THE ALTO PARAGUAI BASIN THROUGH EVI/MODIS TIME SERIES

The main land use and land cover (LULC) changes that a given area passes over the time can be evaluated by using spatial-temporal analysis of satellites images. Then, it is possible to identify the LULC changes, as well as the main causes of environmental impacts. The objective of this paper was to analyze the LULC changes of the main agricultural lands cultivated in the Alto Paraguai Basin (BAP). This paper focused on the summer crops (soybean and corn) and the analysis of agricultural expansion. The results, considering a16-year comparison, showed an increase of 40.60% in the expansion of agricultural areas. The evaluation of the accuracy showed the efficiency of the methodology of agricultural mapping, presenting a Kappa Index of 0.85 for the 2000/2001 and 0.86 for the 2015/2016 crop seasons.


INTRODUCTION
The states of Mato Grosso and Mato Grosso do Sul encompass the Alto Paraguai Basin (BAP), where the Pantanal is located, one of the most important and most conserved biomes in Brazil. The Pantanal has low aptitude for agricultural use, but its preservation depends on the type of land use and land cover (LULC) in the river basin (LIMA et al., 2009).
The states of Mato Grosso and Mato Grosso do Sul show great importance in the national scenario due to their agricultural and livestock aptitude. Due to Brazil's territorial extension, Brazil has agricultural areas for a greater agricultural expansion, being one of them favorable for the production of biofuels (GARDINI, 2007).
The report produced by the Food and Agriculture Organization of the United Nations (FAO), exemplifies the vision of land use in the Cerrado region of the BAP. In the scenario addressed, Brazil has increasingly benefited from low food production costs for livestock farming and is likely to remain a major producer of raw materials. The combination of land abundance and recent advances in infrastructure have been converting formerly remote areas, such as the state of Mato Grosso and the Cerrado, in the central zone of the country. These two regions have the lowest soybean and corn production costs (FAO, 2009).
However, the demand for new agricultural areas should lead to increased pressure on natural resources. The search for new arable land and even changes in the land use must be made within a sustainable concept of development in order not to aggravate even more the environmental imbalances of anthropic origin (LIMA et al., 2009).
The Pantanal and its surroundings are characterized by a rich biodiversity, which faces serious threats with the expansion of the agricultural frontier (IRIGARAY et al., 2011), being classified as the most threatened savannah in the planet (MITTERMEIER et al., 2004).
The change in land use is associated with the expansion of extensive livestock and high grain yield. In the BAP, the slaughterhouse industries intensified their performance, bringing the concept of agribusiness to the new economic contribution of the state of Mato Grosso (IRIGARAY et al., 2011). In Mato Grosso do Sul, mainly in the municipality of São Gabriel do Oeste, between 1984 and 2013, there was a relative increase of 25 to 30% in grain production (soybean and corn), 24% in pasture and a drop of 29% considering the native forest (BULLER et al., 2015(BULLER et al., , 2016.
Over the last few years, the Brazilian Pantanal region has undergone transformations that, from orbital remote sensing data, can be quantified due to the dimensions of the areas covered by high temporal and low spatial resolution sensors, which images are distributed free of charge.
One of the image processing techniques used in land use change analysis and agricultural monitoring is the calculation of vegetation indexes, which can be derived from the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor. Its high periodicity allows to analyze spatially and temporally the evolution of agricultural areas, as well as to carry out mappings and estimates of areas used for agriculture (JOHANN et al., 2012;COUTINHO et al., 2012;ARVOR et al., 2012;BROWN et al., 2013;ALVES et al., 2014;SOUZA et al., 2015;GRZEGOZEWSKI et al., 2016).
In order to verify the land use changes in agricultural areas, it was adopted the methodology of summer crops mapping proposed by Johann et al. (2012) in the state of Paraná. The proposed method uses the MODIS data, taking into account a period of maximum and minimum vegetative vigor in time-series of Enhanced Vegetation Index (EVI) images, from the analysis of sowing dates, vegetative peak dates and harvest of crop areas.
In this context, the aim of this study was to analyze the temporal dynamics of the main cultivated annual crops (soybean and corn) in the BAP by mapping and estimates of cultivated area from EVI/MODIS images between the 2000/2001 and 2015/2016 crop seasons.

MATERIALS AND METHODS
The study area comprises the states of Mato Grosso and Mato Grosso do Sul in the Center-West region of Brazil, with 65% of its territory located in the state of Mato Grosso do Sul and 35% in the state of Mato Grosso (SILVA; ABDON, 1998) ( Figure 1).

Figure 1 -Location of the Alto Paraguai Basin -BAP
The time-series of vegetation indexes was derived from the MODIS sensor, onboard the TERRA satellite. The images were obtained from the MODIS Products Database, a repository available by the Brazilian Agricultural Research Corporation (Embrapa) (ESQUERDO et al., 2010). Gamon et al. (1995) emphasize that the NDVI (Normalized Difference Vegetation Index), proposed by Rouse et al. (1973), may show saturation in the period of maximum vegetative development of agricultural crops. This makes it insensitive to the increase of the vegetal biomass, even when there is increase of the density of the crop canopy. In this way, we used the EVI proposed by Huete et al. (1997) because it is less susceptible to saturation and more sensitive to structure variation, canopy architecture and plant aspect (HUETE et al., 2002).
The 37 images 16-day compositions EVI of each crop season (2000/2001 and 2015/2016) were mosaics and clipped to the BAP. Then, these mosaics of images from different dates were stacked, resulting in a 3D cube of EVI images, where the "x" and "y" axes show the geographic coordinates and the "z" axis the temporal data (ABADE et al., 2015). This allowed the generation of the spectral-temporal profiles of EVI throughout each crop season.
After that, the Flat Bottom Smoother filter (WARDLOW et al., 2006) was applied to this cube of EVI images (WARDLOW et al., 2006) in order to filter cloud noise, as shown in the EVI temporal profile (Figure 2) obtained from a field with agricultural culture, when the beginning of crop development was found on November 1st, 2015 and the end on March 5th, 2016. The average temporal profile in the watershed indicates that the sowing of the summer agricultural crops starts in August and can be extended until November, with a vegetative peak having occurred between November/2015 and March/2016.

Figure 2 -EVI temporal profile of agricultural crops in the 2015/2016 crop season enhanced by Flat Bottom Smoother filtration
To define the range of images used to generate the minimum and maximum EVI compositions in the period of study, based on Johann et al. (2012), the Vegetation Temporal Analysis System (SATVeg), developed by Embrapa, was used. It consists in a Web tool that allows to visualize and to obtain spectral-temporal data of vegetation indexes from MODIS images (EMBRAPA INFORMÁTICA AGROPECUÁRIA, 2015).
Agricultural crop plots distributed throughout the BAP ( Figure 3) were selected to define the minimum EVI values, between 13/08/year1 and 17/11/year1, and the maximum EVI values, between 01/11/year1 and 22/03/year2 ( Figure 2). Johann et al. (2012) proposed this procedure of agricultural crops mapping, aiming to contemplate the variability of sowing season in the different regions of the state of Paraná. Thus, the minimum EVI period (lower EVI values among the images used) contemplates the pre-planting phase and the initial development of agricultural crops within the BAP. The maximum EVI period (highest EVI values among the images used) considers the maximum vegetative development of the agricultural crops or vegetative peak (final development, flowering and production formation) within the BAP area. This stage was operationalized through an image data extraction system developed in IDL (Interactive Data Language) (ESQUERDO et al., 2011).
The ENVI program was used to generate an RGB composition, following the proposition made by Johann et al. (2012). Thus, the maximum EVI image of each crop season was set in the R (red) channel and the minimum EVI image in the G (green) and B (blue) channels, generating a colored composition, where red areas (Figure 4) represent the agricultural crops in the BAP, and the other colors refer to other targets ( Figure 4).
The RGB composition with the original EVI values was then converted from 16 into 8 bits, with scales gray levels (GL) values ranging from 0 to 255, in order to normalize them for each crop season (JOHANN et al., 2012). In order to separate and to extract the red areas (Figure 4), a system developed in IDL language was used to obtain data of this RGB composition in gray levels (GL). The cut-off thresholds (0 to 255 GL) for each of the RGB channels and the studied harvest years were determined by simulation. Thus, the pixels defined according to these thresholds were extracted, resulting in the mapping of agricultural crops in BAP.

Figure 3 -Acquisition of MODIS pixels to evaluate the sowing and harvest period from the temporal profile of agricultural crops using SATVeg
As a post-classification procedure, the removal of isolated pixels was carried out. Preliminarily, the map images were converted from raster into vector format, using the Albers projection. Then, all the polygons were selected and dismembered independently. Thus, for each polygon, the area was determined and later, all isolated polygons with values smaller than 6.25 ha (1 pixel) were excluded.
The determination of sample points at a 95% confidence level was done according to Barbetta (2007) and 400 points were defined as the minimum sample set to evaluate the accuracy of the generated maps. As reference data for checking the results of the classifications, Landsat/TM/OLI images and the EVI temporal profiles of sample were used through the SATVeg. (1) where: E i = estimated summer crop area in the municipality; O i = summer crop area observed in the municipalities; = average summer crop area observed in the municipalities; n = number of municipalities.   Figure 5B), a total of 1,894,456 ha, which represented an agricultural expansion of approximately 40.6% (Figure 6). Similar results were obtained by Coutinho et al. (2016) in the mapping of annual croplands in the BAP through NDVI/MODIS time-series classification, once a growth of approximately 40% in grain agriculture with diffusion to other regions of the plateau in the period between 2001 and 2013 was obtained. This expansion is related to the growth of Brazilian grain exports as well as the increase of protein production in animal confinement systems (cattle and swine).

RESULTS AND DISCUSSION
Taking into account the study of WWF-Brazil and SOS Pantanal (2015) on land use and land cover mapping at BAP, which encompasses annual, perennial and semiperennial crops, using Landsat/TM, Resource-Sat/LISS III and Landsat/OLI, there was an increase in the participation of agricultural areas in BAP plateau from 20,164 km 2 to 21,590 km 2 between 2002 and 2008 corroborating with the growth observed in the mapping of this study with the MODIS sensor.
The methodology used in this study does not consider perennial and semiperennial crops (such as sugarcane), since it is designed to map annual crops, mainly soybean, corn and cotton.
The great advantage of the methodology used here is the low cost, speed and objectivity, allowing the spatialization of croplands and the knowledge of the temporal evolution of crops in the BAP.
The maps accuracy analysis showed a Kappa Index (Table 1) and Koch (1977), indicates an excellent thematic quality (KI ≥ 0.81). Regarding the overall accuracy, according to FOODY (2002), it is desirable that a classification achieves success rates higher than 85%, which was obtained in the presented maps produced by the RGB method in the BAP, with OA of 92. 5% (2000/2001) and 93.0% (2015/2016).
Areas that were not mapped as agricultural crops (Table 1) include forest remnants, scrublands, sugarcane and pasturelands, in addition to other cultivated crops that are present on a smaller scale in BAP.

Figure 6 -Expansion of agricultural areas in the BAP between the 2000/2001 and 2015/2016 crop seasons
The inclusion errors (IE) of 1.7% found for the "agricultural crops" class, considering the 2000/2001 crop season, occurred mainly due to confusions with pasturelands, followed by areas of forest remnants and scrublands. This suggests that the methodology misclassifies these areas as croplands. The same classification error was found for the 2015/2016 crop season with IE of 1.1% (Table 1). Considering the samples drawn on the "other targets" class, the IE were approximately 12% for the two studied crop seasons (Table 1). These errors are represented by croplands, which are not mapped by this methodology; however, most of these classification errors include sugarcane croplands.
The omission errors (OE) found for the "agricultural crops" class were 13.5% and 13.0%, respectively, for the 2000/2001 and 2015/2016 crop seasons (Table 1). This corroborates with Coutinho et al. (2012) in the mapping of the annual agricultural activity carried out in the state of Mato Grosso do Sul. These authors verified that the pixel contamination, caused by different land use and land cover classes, promotes a mixture of the spectral values, which consequently mischaracterizes the pattern of NDVI temporal profile in annual croplands.
Another factor that can contribute to increase the OE refers to the occurrence of small croplands, with areas bellow the MODIS pixel size of 6.25 ha (250 x 250 m), which, according to Wardlow and Egbert (2008), makes difficult the selection of pixels in the croplands boundaries during the classification process. Similar results were also found in the mapping of summer crop areas in the state of Paraná, compared to official statistical data (JOHANN et al., 2012;SOUZA et al., 2015;GRZEGOZEWSKI et al., 2016).
These are the main problems found concerning the mapping methods based on MODIS time-series data, especially in regions where the land structure and the spatial distribution pattern of agricultural activities define a mosaic occupation very heterogeneous (COUTINHO et al., 2012).
The spectral confusion between agricultural crops, pastures and forest remnants can be attributed to the similarity of the temporal profile, since the BAP vegetation is composed by a great variety of plants that compose landscapes of cerrado and forest, formations of high complexity of endemic species (GOODLAND;FERRI, 1979).
Based on the accuracy evaluation, it was verified that in 2000/2001, some areas that were occupied by pastures in 2015/2016 began to be occupied by agricultural crops in greater proportion (soybean and corn), as observed from the analysis done with Landsat images as auxiliary data. This finding corroborates with previous results of land use and land cover changes in BAP between 2002 and 2008 which increased in the plateau, with the conversion of 1.019 km 2 of pasturelands to croplands (WWF-Brazil et al., 2009).
The comparison of these results with the official IBGE data showed that there was an overestimation of crop area by this method in the two studied crop seasons. As the official data considers the total crop area per municipality, the municipalities partially included in the BAP (less than 90%) were not considered in this analysis. Then, from the 92 municipalities of BAP, the statistical analysis was performed in 47 of them. Thus, the greatest difference was found for the 2015/2016 crop season (41,139 ha) and the smallest one for the 2000/2001 crop season (38,199 ha), which represents, respectively, 7.11% and 9.89% of the BAP area with summer crops (Table 2).
Since the municipal areas (ha) of the official IBGE data and the estimates obtained by the mappings did not show normality, the Spearman correlation coefficient (r s ) was used to analyze the association between such data. Correlation values were obtained   (Figure 7). Thus, there was a trend of overestimation of summer crops cultivated area in the maps when compared to the official IBGE data. Grzegozewski et al. (2016) observed the same for soybean and corn crops. This fact may be related to soybean cultivation in larger areas, while corn cultivation may occur in areas smaller than the MODIS pixel size. This can also happen for rice and cotton crops, which are present in a lesser extent in BAP (IBGE, 2016).

CONCLUSIONS
The use of multitemporal EVI/MODIS spectral data allowed the mapping and estimation of cultivated area of the main summer crops in BAP.
The use of the time-series MODIS vegetation index data in the generation of minimum and maximum EVI compositions reduced the variability in the planting dates, allowing quickness in the mapping of agricultural areas by the RGB methodology, allowing at the same time the reproduction of the same period of images and cuts for other crop seasons.
Considering the two studied crop seasons, the OE were higher than the IE, which are represented by forest and pasture areas due to the spectral confusion of these targets with agricultural crops.
The analysis of the multitemporal evolution in the evaluated period indicated that there was a significant expansion of 40.6% in the area of summer agricultural crops in the BAP between the first and last season studied.
The use of the adopted methodology has the greatest advantage of being fast and efficient in the estimation of cultivated areas of summer crops, which allows in advance its determination of the harvesting season, with its proven efficiency by statistical methods, although the technique can be further improved in the future.