1. 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.
Remote sensing time-series are widely used to map agricultural land use and dynamics [
18,
19,
20,
21]. Vegetation index time-series show consistent results in the classification of pasture and crop areas [
22,
23,
24,
25,
26,
27]. It is also used to study land-use changes, detecting different land-use patterns in the same year (e.g., summer-crop and winter-crop) and over several years [
28,
29,
30,
31,
32,
33].
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.
2. Materials and Methods
2.1. 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.
2.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.
2.3. 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.
2.4. 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].
2.5. 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].
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 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.
2.5.2. 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. 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.
2.5.3. 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.
2.6. 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.
2.6.1. 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.
2.6.2. 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.
4. 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 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.
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-2A, and Sentinel-2B are not sufficiently long-lasting to cover an ICL period of analysis. In the next few years, there is the potential for integrated systems to use these products individually or with fusion techniques, as several studies have shown good results for mapping and monitoring agriculture [
20,
37,
51,
52]. Despite the different characteristics of the products, a united source of Landsat-8 and Sentinel-2 data can provide 10–30 m information and produce more frequent observations for monitoring [
20]. The TWDTW method had certain applications with Landsat-8 OLI and Sentinel-2 products [
37,
49,
53]. The exploration of the TWDTW method with other classification techniques offers an opportunity for future studies. Belgiu and Csillik [
49] showed good results when they combined object-oriented classification and the TWDTW method.