Framework for Mapping Integrated Crop-Livestock Systems in Mato Grosso , Brazil

Integrated crop-livestock (ICL) systems combine livestock and crop production in the same area, increasing the efficiency of land use and machinery, while mitigating greenhouse gas emissions, and reducing production risks, plant diseases and pests. ICL systems are primarily divided into annual (ICLa) and multi-annual (ICLm) systems. Projects such as the “Integrated crop-livestock-forest Network” and the “Livestock Rally” have estimated the ICL areas for Brazil on a state or regional basis. However, it remains necessary to create methods for spatial identification of ICL areas. Thus, we developed a framework for mapping ICL areas in Mato Grosso, Brazil using the Enhanced Vegetation Index time-series of Moderate Resolution Imaging Spectroradiometer and a Time-Weighted Dynamic Time Warping (TWDTW) classification method. The classification of ICL areas occurred in three phases. Phase 1 corresponded to the classification of land use from 2008 to 2016. In Phase 2, the ICLa areas were identified. Finally, Phase 3 corresponded to the ICLm identification. The framework showed overall accuracies of 86% and 92% for ICL areas. ICLm accounted for 87% of the ICL areas. Considering only agricultural areas or only pasture areas, ICL systems represented 5% and 15%, respectively.


Introduction
The need to increase food production and the restriction of agricultural expansion by environmental policies led to the adoption of production systems with higher productivity and lower environmental impacts, such as integrated crop-livestock (ICL) systems [1][2][3][4][5].These systems employ strategies that aim to increase the land production such as, the integrated use of land as pastures for livestock and crop activities in succession, rotation, or intercropping [6][7][8].In Brazil, ICL systems are associated with an improvement in pasture conditions and an intensification in land use [3,9,10].The Brazilian Agriculture Ministry is aiming for four million new ICL areas by 2020 as part of the Brazilian Low Carbon Agriculture Plan-ABC Plan [11].
The interaction between crop and livestock makes the nutrient cycle more efficient, improves soil quality, and increases overall productivity [12][13][14].Thus, integrated systems increase land and machinery efficiency, mitigate greenhouse gas emissions, and reduce production risks, plant diseases, and pests [15,16].The ICL systems can be divided into two groups, annual and multi-annual [6].In the annual group, land is used for grain cultivation in summer and pasture cultivation for livestock in winter.In the multi-annual ICL systems, usually for pasture reform, crop and pasture are cultivated in different years in rotation or succession.
The "Integrated crop-livestock-forest network" project ("Rede Integração-Lavoura-Pecuária") conducted interviews to estimate the extent of the ICL system areas in Brazil on a state scale [17].With the support of the "ABC Plan", the "Livestock Rally" project performs fieldwork campaigns to understand and quantify livestock production and ICL systems in select regions in Brazil.Gil et al. [10] quantified farms with different ICL system combinations on a municipal scale in the Brazilian state of Mato Grosso.Despite these efforts to quantify ICL areas throughout Brazil, it is necessary to develop a framework to identify the spatial distribution of ICL areas within the municipalities.Remote sensing is an alternative approach for the systematic identification of ICL areas.
Many time-series classification methods have been developed to identify land use through remote sensing time-series [20,[34][35][36].Among these methods, Maus et al. [29] developed a Time-Weighted Dynamic Time Warping (TWDTW) method that improves the land-use classification in different years.
To improve the land-use classification within the same year, they used a Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation index time-series.The TWDTW method involves comparing two time-series, finding the optimal time-alignment, and then generating a measure of the dissimilarity (TWDTW distance), in which lower values indicate a greater similarity in the time-series temporal behavior.The use of a vegetation index times-series and a TWDTW classification method showed effectiveness compared with other classification methods, such as random forests [37].
Therefore, this study aims to develop a framework for mapping ICL system areas in a region of Brazil using vegetation index time-series from MODIS and a TWDTW classification method.

Study Area
The study area is located in the Sinop, Cláudia, and Santa Carmém municipalities in the Mato Grosso state of Central-Western Brazil, shown in Figure 1.Mato Grosso is among the states that have the largest ICL areas in Brazil [10,17].Since 2001, intensification in land use has occurred in this region, moving from low-production pasture for cattle ranching to crop production or higher livestock production [38].
In order to understand the distribution and arrangements of ICL systems in the study area, three fieldwork campaigns were conducted.We visited properties that had adopted ICL systems and attended technical meetings at research institutes in order to understand the integration of grain and livestock production in the region.We gathered information about land management, historical land use, and ICL systems.
In addition, to develop and understand the production systems, ground truth information of ICL systems, crops (summer and winter crops), pasture, forest, and other land-use types were collected for the classification framework during the fieldwork campaigns.Information gathered in the first (April 2016) and second (December 2016) field campaigns was used for training the classifier and validation, while information collected from the third (June 2017) campaign was used only for validation.

MODIS Data
We chose to use the Enhanced Vegetation Index (EVI) product of the MODIS/Terra sensor as it had images recorded throughout the study period and adequate time resolution for seasonal agricultural identification [39][40][41].Among the MODIS products available, we selected EVI from MOD13Q1, with a 16-day composition and 250 m spatial resolution [42].Despite the absence of a consensus regarding the best vegetation index for classification [43][44][45], an EVI time-series is more sensitive for high biomass than the NDVI and had been used in previous studies for the temporal characterization and identification of pasture and crop [24,46].
In the ICL classification framework, we used a total of 210 EVI images, from 31 June 2007 to 31 September 2016.The Savitzky-Golay filter was applied to reduce the time-series noise across all periods [47], as shown in Figure 2.

MODIS Data
We chose to use the Enhanced Vegetation Index (EVI) product of the MODIS/Terra sensor as it had images recorded throughout the study period and adequate time resolution for seasonal agricultural identification [39][40][41].Among the MODIS products available, we selected EVI from MOD13Q1, with a 16-day composition and 250 m spatial resolution [42].Despite the absence of a consensus regarding the best vegetation index for classification [43][44][45], an EVI time-series is more sensitive for high biomass than the NDVI and had been used in previous studies for the temporal characterization and identification of pasture and crop [24,46].
In the ICL classification framework, we used a total of 210 EVI images, from 31 June 2007 to 31 September 2016.The Savitzky-Golay filter was applied to reduce the time-series noise across all periods [47], as shown in Figure 2.

MODIS Data
We chose to use the Enhanced Vegetation Index (EVI) product of the MODIS/Terra sensor as it had images recorded throughout the study period and adequate time resolution for seasonal agricultural identification [39][40][41].Among the MODIS products available, we selected EVI from MOD13Q1, with a 16-day composition and 250 m spatial resolution [42].Despite the absence of a consensus regarding the best vegetation index for classification [43][44][45], an EVI time-series is more sensitive for high biomass than the NDVI and had been used in previous studies for the temporal characterization and identification of pasture and crop [24,46].
In the ICL classification framework, we used a total of 210 EVI images, from 31 June 2007 to 31 September 2016.The Savitzky-Golay filter was applied to reduce the time-series noise across all periods [47], as shown in Figure 2.

Field Data
Two types of in situ data were collected across the study area.The first type of ground truth was the history of agricultural practices for individual fields.The focus was to record the location and years where ICLa and ICLm systems occurred.A data acquisition protocol in which farmers/farm managers participated in an interview was established.The participants indicated the property fields by drawing polygons over a Landsat 8 OLI image and described the previous year's land use in each area.If at least one pixel from MODIS-250 m fitted inside the polygon, a point in the region was created with an unmixed MODIS-250 m pixel.Otherwise, the polygon and the field information was discarded.Thus, during the three fieldwork campaigns, 185 fields in 22 properties were recorded and validated.In the first fieldwork campaign, the field locations were indicated by local researchers and, in the second and third campaigns, the randomly selected points in the classification validation step guided the field visits.
The second type of ground truth was collected in the first fieldwork campaign (April 2016).It was conducted by local field observation for double-crop, single-crop, pasture, and forest.With a Geographic Information System software, Global Positioning System information, Landsat 8 OLI image, and MODIS-250 m pixel grid, fields were visited, and the points were recorded in the center of a homogeneous field which size fitted in one unmixed MODIS 250 m pixel.As the campaign was in April, it was possible to differentiate between single-and double-crop by the observation of winter crops or fallow in crop areas.In total, 280 points were recorded, see Table 1.From the field points, ten points were selected for each class, creating the training data set.The selection of training data set points for ICLa, ICLm, Double-Crop (DC), and Pasture, were from the interview data.The rest of the points were allocated in the test data set.

TWDTW Method
We selected a Time-Weighted Dynamic Time Warping (TWDTW) method for time-series identification [29].This method is an adaptation of the Dynamic Time Warping (DTW) method, and the R software in the "dtwSat" package were used.Bagnall et al. [48] studied time-series classification and reported that DTW-based methods were among the best performing time-series classification methods.The DTW method involves the comparison of two time-series and subsequent generation of a measure of the dissimilarity, in which lower values indicate a greater similarity in the time-series temporal behavior.However, the DTW method is sensitive to temporal seasonality.In the TWDTW method, Maus et al. [29], adapted the DTW method to find the optimal time-alignment and then generate a measure of dissimilarity (TWDTW distance).Time variations occur in phenological cycles of natural or planted vegetation.
In land-use classification, the identification of time-series was accomplished by defining temporal patterns of land-use classes and generating TWDTW distances through the comparison of temporal patterns with the time-series for each pixel.Each pixel was classified using the temporal pattern with the lowest TWDTW distance, k-Nearest Neighbor where k = 1.The classification method was sufficiently robust to classify single-crops, double-crops, forest, pasture, and certain winter crops using an EVI/MODIS time-series over several years [29,49].
The temporal patterns definition for TWDTW distance generation was created from training samples obtained in the field campaigns.For each class created, the EVI time-series from the training samples were grouped and the Generalized Additive Models (GAM) time-series filter was applied, generating one temporal pattern for each class.The GAM time-series filter was chosen as it was able to generate one time-series from a data set and could be successfully applied to remote sensing data [29,50].

Classification
The classification framework was designed to detect the two ICL systems, ICLa and ICLm, in the study area of 2016.The ICLa systems had detectable characteristics in one year (crop and pasture in the same agricultural year), requiring classification only for the 2016 agricultural year (September 2015 up to August 2016).The ICLm systems present multi-year characteristics, and for this reason, classification was necessary for each year between 2008 and 2016, as shown in Figure 2.
The following thematic classes were established for all classification processes: Forest, Pasture, SC, DC, ICLa, and ICLm.During the phases of classification, described below, the thematic classes may have been grouped or divided into subclasses to achieve better ICL differentiation, see Table 2.The classification of ICL areas occurred in three phases, as shown in Figure 3.The first corresponded to the classification from 2008 to 2016 of Forest, Pasture, SC, and DC classes, using a temporal period of one agricultural year (September to August).In the second phase, ICLa was identified in a new period of analysis (June to August).In the third phase, an ICLm-class was identified through the land-use change analysis, based on Phase 1 classification.Based on the "Terrclass 2014 project" urban areas, bodies of water, and mining areas were grouped in the class, Others [19].
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 18 pattern with the lowest TWDTW distance, k-Nearest Neighbor where k = 1.The classification method was sufficiently robust to classify single-crops, double-crops, forest, pasture, and certain winter crops using an EVI/MODIS time-series over several years [29,49].The temporal patterns definition for TWDTW distance generation was created from training samples obtained in the field campaigns.For each class created, the EVI time-series from the training samples were grouped and the Generalized Additive Models (GAM) time-series filter was applied, generating one temporal pattern for each class.The GAM time-series filter was chosen as it was able to generate one time-series from a data set and could be successfully applied to remote sensing data [29,50].

Classification
The classification framework was designed to detect the two ICL systems, ICLa and ICLm, in the study area of 2016.The ICLa systems had detectable characteristics in one year (crop and pasture in the same agricultural year), requiring classification only for the 2016 agricultural year (September 2015 up to August 2016).The ICLm systems present multi-year characteristics, and for this reason, classification was necessary for each year between 2008 and 2016, as shown in Figure 2.
The following thematic classes were established for all classification processes: Forest, Pasture, SC, DC, ICLa, and ICLm.During the phases of classification, described below, the thematic classes may have been grouped or divided into subclasses to achieve better ICL differentiation, see Table 2.The classification of ICL areas occurred in three phases, as shown in Figure 3.The first corresponded to the classification from 2008 to 2016 of Forest, Pasture, SC, and DC classes, using a temporal period of one agricultural year (September to August).In the second phase, ICLa was identified in a new period of analysis (June to August).In the third phase, an ICLm-class was identified through the land-use change analysis, based on Phase 1 classification.Based on the "Terrclass 2014 project" urban areas, bodies of water, and mining areas were grouped in the class, Others [19].In the previous analysis, ICLa and DC classification was tested in one step using the annual pattern, which resulted in a great misclassification among these two classes.Therefore, we decided to group all the crops for the two seasons (summer and winter crops) in one class in the first classification phase, and then classify DC and ICLa in the second classification phase.This approach was taken in other studies, where the two-seasons classes were classified separately [24,37,40,46].The annual temporal pattern for the classes, as shown in Figure 4a, was determined from the training data collected in the field campaigns.The TWDTW distance was calculated for each class from 2008 until 2016, as shown in Figure 4b.By comparison at the pixel level, the class with the lowest value of TWDTW distance was assigned to the pixel.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 18 In the previous analysis, ICLa and DC classification was tested in one step using the annual pattern, which resulted in a great misclassification among these two classes.Therefore, we decided to group all the crops for the two seasons (summer and winter crops) in one class in the first classification phase, and then classify DC and ICLa in the second classification phase.This approach was taken in other studies, where the two-seasons classes were classified separately [24,37,40,46].The annual temporal pattern for the classes, as shown in Figure 4a, was determined from the training data collected in the field campaigns.The TWDTW distance was calculated for each class from 2008 until 2016, as shown in Figure 4b.By comparison at the pixel level, the class with the lowest value of TWDTW distance was assigned to the pixel.

Phase 2
In Phase 2, the areas classified as DC + ICLa in 2016 were separated into DC and ICLa.The winter crop season (June to August) was defined as the period of analysis.After the first tests with TWDTW distance for DC and ICLa separation, the misclassification of cotton cultivation areas as ICLa was observed.Therefore, another class, Cotton, was created at this stage and later reincorporated into DC.Based on field training samples, the temporal patterns for the three classes in this stage (ICLa, Cotton, and DC) were established, and the TWDTW distances for each pattern were generated, see Figure 5.
In the areas classified as DC + ICLa in Phase 1 for the year 2016, the TWDTW distances for the three thematic classes were compared.An area was considered to belong to the class ICLa when the TWDTW distance of ICLa was lower than 2.2, and was the smallest among the other two classes (DC and Cotton), see Figure 6

Phase 3
In order to detect areas where ICLm occurred as either rotational systems or pasture reform with crop integration, we analyzed changes in the land-use classification between crops and pasture, during 2008-2016.
To be classified as ICLm, pixels were required to have had at least twice the land-use changes-Pasture-Crop (DC or SC), Pasture or Crop (DC or SC), or Pasture-Crop (DC or SC).At least two consecutive years of pasture permanence were also required to avoid including areas that were mistakenly classified as pasture areas in Phase 1. Further, to exclude areas of recent conversion to crop in ICLm, areas classified as crop (SC or DC) since 2014 were not regarded as ICLm, as shown in Figure 6.The classification for the final year, 2016, was evaluated by merging the results generated in Phases 2 and 3.If the same pixel was classified as ICLa and ICLm, the pixel was classified as ICLa.

Classification Validation
We validated the ICL classification and Phase 1 classification in two different manners; on randomly selected control points and on field data.This was necessary because our field data does not cover the entirety of the study area and is limited to the 2016 classification.The ICL ground truth is difficult to record; therefore an interview protocol is necessary.The randomly selected control points approach made the validation of all study areas over all time periods in Phase 1 possible.

Validation of the MODIS Classification on Randomly Selected Control Points
Phase 1 validation was carried out through analysis of 400 randomly distributed pixels for each year from 2008 to 2016, generating the confusion matrix.The reference for each validation pixel was obtained based on the EVI time-series interpretation.For 2016 classification, 52 points from the test data set were used for validation.

Phase 3
In order to detect areas where ICLm occurred as either rotational systems or pasture reform with crop integration, we analyzed changes in the land-use classification between crops and pasture, during 2008-2016.
To be classified as ICLm, pixels were required to have had at least twice the land-use changes-Pasture-Crop (DC or SC), Pasture or Crop (DC or SC), or Pasture-Crop (DC or SC).At least two consecutive years of pasture permanence were also required to avoid including areas that were mistakenly classified as pasture areas in Phase 1. Further, to exclude areas of recent conversion to crop in ICLm, areas classified as crop (SC or DC) since 2014 were not regarded as ICLm, as shown in Figure 6.The classification for the final year, 2016, was evaluated by merging the results generated in Phases 2 and 3.If the same pixel was classified as ICLa and ICLm, the pixel was classified as ICLa.

Classification Validation
We validated the ICL classification and Phase 1 classification in two different manners; on randomly selected control points and on field data.This was necessary because our field data does not cover the entirety of the study area and is limited to the 2016 classification.The ICL ground truth is difficult to record; therefore an interview protocol is necessary.The randomly selected control points approach made the validation of all study areas over all time periods in Phase 1 possible.

Validation of the MODIS Classification on Randomly Selected Control Points
Phase 1 validation was carried out through analysis of 400 randomly distributed pixels for each year from 2008 to 2016, generating the confusion matrix.The reference for each validation pixel was obtained based on the EVI time-series interpretation.For 2016 classification, 52 points from the test data set were used for validation.
For ICL classification validation, 114 stratified pixels were generated randomly within the area classified as ICL, and 86 stratified pixels were generated randomly outside this area.Pixels in areas classified as ICLa were validated with a visit to the area to collect land-use information.This was necessary, as areas of ICLa, Cotton, and DC are not clearly distinguishable using only the visual interpretation of EVI time-series.For pixels in ICLm areas and outside the ICL areas, verification was based on an EVI time-series visual interpretation and the conditions used in the Phase 3 classification.Visits were made for collecting information on the land-use history of 20 ICLm and five double-crop areas.

Validation of the MODIS Classification on Field Data
The validation of field data was performed for ICL and Phase 1 classification for the year 2016, using all 405 points of the test data set.For Phase 1 2016 classification, a test data set was reordered for the four classes, DC + ICLa, SC, Pasture, and Forest.For ICL classification, the ICL classes merge into one class with non-ICL classes in another.

Validation of the MODIS Classification on Randomly Selected Control Points
In Phase 1, the classifications for DC + ICLa, SC, Pasture, and Forest over the time range of 2008-2016 were generated, as shown in Figure 7.By analyzing the confusion matrix of the Phase 1 classification on randomly selected control points, shown in Table 3, it was found that the overall classification accuracy varied between 86-91%.The crop classes, DC + ICLa and SC, presented a user´s accuracy (UA) between 0.60-0.88,and a producer´s accuracy (PA) between 0.61-0.87.The higher amount of errors observed for SC were in the DC + ICLa and Pasture, and for DC + ICLa, were in the SC and Pasture classes.Apart from the year 2016, misclassification within DC + ICLa and SC did not influence the ICLm classification.This is because, in the analysis on conversion between crop and pasture in Phase 3, crop classes were grouped.The Pasture class had a UA of 0.65-0.81and a PA of 0.68-0.84,and errors mainly occurred in classes, SC and Forest.Forest had the largest separation from the other classes with a UA and PA ranging from 0.94-0.99.The classification of DC + ICLa for the year 2016 was made based on the ICLa identification and because the highest UA and PA values, 0.87 and 0.88, respectively, were recorded in this year.
Using the proposed ICL classification framework and validation on randomly selected control points, an overall accuracy of 86% was achieved, as shown in Table 4.The UA and PA were 0.80 and 0.94, respectively.Evaluation of the ICL-class classifications, ICLa and ICLm, showed that the crop classes (SC and DC) generated the greatest errors in Phases 1 and 2 of the classification, see Table 5.For the ICLa class, the greatest errors occurred in the differentiation of the targets in Phase 2, with difficulty in the separation of ICLa and Cotton.As shown in Figure 5, these two classes had a similar temporal pattern.The highest misclassification errors for ICLm arose when differentiating between Pasture and SC.This led to errors in Phase 3 of the classification as the detection of a land-use change for ICLm generation was compromised.Errors in differentiating between SC and DC occurred in the year 2016, due to a long period of drought during the winter harvest season.However, such errors did not affect the classification of ICL areas significantly.

Validation of the MODIS Classification on Field Data
The confusion matrix generated from the classification validation on field data for the Phase 1 2016 classification had an accuracy of 0.93, see Table 6.The classes, DC + ICLa, Pasture, and SC presented a UA between 0.79-0.97and PA between 0.84-0.94.A larger quantity of errors was found in the misclassification of DC + ICLa and Pasture, and DC + ICLa and SC.The ICL classification for

Validation of the MODIS Classification on Field Data
The confusion matrix generated from the classification validation on field data for the Phase 1 2016 classification had an accuracy of 0.93, see Table 6.The classes, DC + ICLa, Pasture, and SC presented a UA between 0.79-0.97and PA between 0.84-0.94.A larger quantity of errors was found in the misclassification of DC + ICLa and Pasture, and DC + ICLa and SC.The ICL classification for the year 2016, achieved an overall accuracy of 92%, as shown in Table 7.The UA and PA were 0.80 and 0.94, respectively.

Integrated Crop-Livestock Systems Area
The ICL area in 2016 was 219.18 km 2 , see Table 8, with ICLa and ICLm areas accounting for 28.28 km 2 and 190.90 km 2 , respectively.Multi-year systems (ICLm) accounted for 87% of the ICL areas.The area covered by forests is the largest in this region, followed by agriculture and pasture.The ICL systems represented 5% of the total agricultural area found within the study area.If only the pasture is considered, ICL areas represent 15% of the total area.Indicating that a significant number of livestock producers in the region intensified production by using integrated systems.

Discussion and Conclusions
The framework presented herein provided a good accuracy overall in land-use classification and ICL classification, shown in Figure 8.The largest errors for the classes ICLm and ICLa, were distinguishing cotton cultivation and the misclassification between pasture and crop, respectively.The classification techniques proved to be efficient for the characteristics of the land uses in the region.The ICL systems made up 5% of the total agricultural area and 15% of the total pasture area in the study area.The proposed framework could detect ICL areas within the Sinop, Mato Grosso region.The detection of ICL areas in other regions of Mato Grosso or Brazil will require methodological adjustments due to the possible presence of other land uses, such as sugarcane, natural pastures, and double-crop combinations.However, in regions with the exclusive presence of forest, grains, pasture, and similar ICL systems, only a few adjustments related to the locality may be necessary, such as different periods for establishing the EVI temporal patterns and the creation of new classes in classification Phases 1 and 2.
The overall ICL classification accuracies, in the two validations, were similar to those reported by previous studies that classified summer and winter crops in Mato Grosso using a MODIS vegetation index.Using MODIS/NDVI, Chen et al. [40] classified six crop arrangements; the soyapasture (ICLm) class had a UA of 0.81 and PA of 0.74.Maus et al. [29] demonstrated the functionality of the TWDTW method for property-scale classification with an overall accuracy of 90%.The classification by Arvor et al. [46] was in two phases; first, the land uses of Cerrado, forest, agriculture, and pasture in the state of Mato Grosso, obtaining an accuracy rate of 86%; second, five crop arrangements obtaining an accuracy rate of 74%.
The most difficult step in this study was the validation step.The reference data obtained from the collected field information was fundamental for validating the accuracy of the classification.However, in ICL areas, the process for obtaining reference information was more complex.This is because, in most cases, it was not possible to determine whether ICL systems were adopted with a single visit to the study area; therefore, either more than one visit was required during the study period or it was necessary to determine the land-use history by interviewing the landowners.Despite the differences in the two validation methods, the results for the Phase 1 2016 classification did not have large differences and obtained similar values of UA and PA.For ICL classification the validation on field data had a higher accuracy and consequently, more accurate UA and PA values.
The choice of MODIS/Terra in this study was due to the temporal resolution and images available from 2007 to 2016.However, satellite-derived products, such as Landsat-8 OLI, Sentinel- The proposed framework could detect ICL areas within the Sinop, Mato Grosso region.The detection of ICL areas in other regions of Mato Grosso or Brazil will require methodological adjustments due to the possible presence of other land uses, such as sugarcane, natural pastures, and double-crop combinations.However, in regions with the exclusive presence of forest, grains, pasture, and similar ICL systems, only a few adjustments related to the locality may be necessary, such as different periods for establishing the EVI temporal patterns and the creation of new classes in classification Phases 1 and 2.
The overall ICL classification accuracies, in the two validations, were similar to those reported by previous studies that classified summer and winter crops in Mato Grosso using a MODIS vegetation index.Using MODIS/NDVI, Chen et al. [40] classified six crop arrangements; the soya-pasture (ICLm) class had a UA of 0.81 and PA of 0.74.Maus et al. [29] demonstrated the functionality of the TWDTW method for property-scale classification with an overall accuracy of 90%.The classification by Arvor et al. [46] was in two phases; first, the land uses of Cerrado, forest, agriculture, and pasture in the state of Mato Grosso, obtaining an accuracy rate of 86%; second, five crop arrangements obtaining an accuracy rate of 74%.
The most difficult step in this study was the validation step.The reference data obtained from the collected field information was fundamental for validating the accuracy of the classification.However, in ICL areas, the process for obtaining reference information was more complex.This is because, in most cases, it was not possible to determine whether ICL systems were adopted with a single visit to the study area; therefore, either more than one visit was required during the study period or it was necessary to determine the land-use history by interviewing the landowners.Despite the differences in the two validation methods, the results for the Phase 1 2016 classification did not have large differences and obtained similar values of UA and PA.For ICL classification the validation on field data had a higher accuracy and consequently, more accurate UA and PA values.

Figure 2 .
Figure 2. Examples of annual integrated crop-livestock (ICLa) and multi-annual integrated crop-livestock (ICLm) Enhanced Vegetation Index (EVI) temporal behavior, non-filtered, and filtered by Savitzky-Golay method.The brackets highlight land-use changes between crop and pasture across multiple agricultural years, and within an agricultural year, represented by ICLm and ICLa, respectively.

Figure 3 .
Figure 3.The flowchart of the proposed method for mapping integrated crop-livestock systems.Note: ICL, integrated crop-livestock.

Figure 3 .
Figure 3.The flowchart of the proposed method for mapping integrated crop-livestock systems.Note: ICL, integrated crop-livestock.
2.5.1.Phase 1 In this step, the classes, Forest, Pasture, SC, and DC + ICLa were classified from 2008 to 2016.The class, DC + ICLa was created by grouping classes, DC and ICLa as they have similar EVI temporal patterns, as shown in Figure 2. In both cases, cultivation of summer crops (usually soya) and winter cultivation (grains for DC and pasture for ICLa) occur, generating a similar annual pattern of vegetation index.
2.5.1.Phase 1 In this step, the classes, Forest, Pasture, SC, and DC + ICLa were classified from 2008 to 2016.The class, DC + ICLa was created by grouping classes, DC and ICLa as they have similar EVI temporal patterns, as shown in Figure 2. In both cases, cultivation of summer crops (usually soya) and winter cultivation (grains for DC and pasture for ICLa) occur, generating a similar annual pattern of vegetation index.
. The definition of values lower than 2.2 occurred based on classification tests.Areas classified as DC + ICLa in Phase 1 and not classified as ICLa in Phase 2 were considered to belong to DC.

Figure 5 .
Figure 5. (a) Temporal patterns and (b) TWDTW distance subset for classes in Phase 2.Figure 5. (a) Temporal patterns and (b) TWDTW distance subset for classes in Phase 2.

Figure 5 .
Figure 5. (a) Temporal patterns and (b) TWDTW distance subset for classes in Phase 2.Figure 5. (a) Temporal patterns and (b) TWDTW distance subset for classes in Phase 2.

Figure 6 .
Figure 6.The flowchart of the mapping framework for integrated crop-livestock systems, Phase 1-3.

Figure 6 .
Figure 6.The flowchart of the mapping framework for integrated crop-livestock systems, Phase 1-3.

Table 1 .
Number of field data/points collected by interviews and field observation, for 2016.

Table 2 .
Description of land cover classes.
DescriptionICLaIntegrated crop-livestock in succession in the same area and year ICLm Integrated crop-livestock in rotation or succession in the same area and in different

Table 2 .
Description of land cover classes.

Table 3 .
Confusion matrices and accuracy statistics for Phase 1 classification, validation on randomly selected control points.Note: P, pasture; F, forest; O, other uses; UA, user´s accuracy; PA, producer´s accuracy; OA, overall accuracy.

Table 4 .
Confusion matrices and accuracy statistics for integrated crop-livestock classification, validation on randomly selected control points.

Table 6 .
Confusion matrices and accuracy statistics for Phase 1 classification, validation on field data.

Table 7 .
Confusion matrices and accuracy statistics for integrated crop-livestock classification, validation on field data.

Table 8 .
Classification area, percentage of total agricultural areas and percentage of pasture area.