Evaluation of Sentinel-1 and 2 Time Series for Land Cover Classification of Forest – Agriculture Mosaics in Temperate and Tropical Landscapes

Monitoring forest–agriculture mosaics is crucial for understanding landscape heterogeneity and managing biodiversity. Mapping these mosaics from remotely sensed imagery remains challenging, since ecological gradients from forested to agricultural areas make characterizing vegetation more difficult. The recent synthetic aperture radar (SAR) Sentinel-1 (S-1) and optical Sentinel-2 (S-2) time series provide a great opportunity to monitor forest–agriculture mosaics due to their high spatial and temporal resolutions. However, while a few studies have used the temporal resolution of S-2 time series alone to map land cover and land use in cropland and/or forested areas, S-1 time series have not yet been investigated alone for this purpose. The combined use of S-1 & S-2 time series has been assessed for only one or a few land cover classes. In this study, we assessed the potential of S-1 data alone, S-2 data alone, and their combined use for mapping forest–agriculture mosaics over two study areas: a temperate mountainous landscape in the Cantabrian Range (Spain) and a tropical forested landscape in Paragominas (Brazil). Satellite images were classified using an incremental procedure based on an importance rank of the input features. The classifications obtained with S-2 data alone (mean kappa index = 0.59–0.83) were more accurate than those obtained with S-1 data alone (mean kappa index = 0.28–0.72). Accuracy increased when combining S-1 and 2 data (mean kappa index = 0.55–0.85). The method enables defining the number and type of features that discriminate land cover classes in an optimal manner according to the type of landscape considered. The best configuration for the Spanish and Brazilian study areas included 5 and 10 features, respectively, for S-2 data alone and 10 and 20 features, respectively, for S-1 data alone. Short-wave infrared and VV and VH polarizations were key features of S-2 and S-1 data, respectively. In addition, the method enables defining key periods that discriminate land cover classes according Remote Sens. 2019, 11, 979; doi:10.3390/rs11080979 www.mdpi.com/journal/remotesensing Remote Sens. 2019, 11, 979 2 of 20 to the type of images used. For example, in the Cantabrian Range, winter and summer were key for S-2 time series, while spring and winter were key for S-1 time series.


Introduction
Mapping forest-agriculture mosaics is essential for understanding effects of land cover changes on economic (e.g., agricultural and silvicultural production, timber extraction, ecotourism) and non-market ecosystem services (e.g., carbon sequestration, water quality, stream flow, species conservation) [1].Cropland and pastures are the largest use of Earth's ice-free land [2] and have particularly high biodiversity potential [3].However, this capacity refers to landscapes with heterogeneous semi-intensive agricultural systems and semi-natural elements (hedgerows, wetlands, forests) [4][5][6].In many parts of the world, land cover changes affect landscapes, causing habitat loss and fragmentation [7].Habitat fragmentation decreases the size of habitats and populations and increases distances between patches of habitat, which may result in species loss [8].
Each species uses different kinds of habitats and requires different amounts of them [7].Many species perceive landscapes in ways more complex than as a "habitat-matrix" dichotomy (i.e., dividing the landscape into "habitat", which contains all necessary resources, and "matrix", which is hostile) and use resources from different cover types [5].Thus, landscape management can be considered a complex process that requires characterizing the landscape composition and structure suitable for the target species.Habitat selection and resistance models are widely used as a decision support tool, since they can apply to multiple groups of species (e.g., insects, birds, herpetofauna, mammals, plants).Landscape maps are essential for developing these models, and most environmental features used in the models are derived from remote sensing imagery [9].Remotely sensed imagery is a common tool for straightforward land cover classification.The main difficulty is that the land cover classes in forest-agriculture mosaics are distributed along a landscape gradient, which results in misclassification.Images with high and very high spatial resolutions provide greater spatial detail and precise information [10][11][12][13].However, the low temporal resolution of such data does not enable monitoring vegetation dynamics to better discriminate land cover classes that are spectrally similar.In contrast, satellite time series provide great opportunities for this purpose.
MODIS and Landsat optical time series have been used extensively for land cover classification since the 1970s.Despite a 1-2-day revisit time, however, land cover derived from MODIS has low local accuracy due to the maximum spatial resolution of 250 m.The Landsat constellation has an 8-day revisit time but has a higher spatial resolution (30 m) than MODIS.The S-2 satellite has even higher capacities, with a five-day revisit time and a 10 m spatial resolution [14].Despite these high capacities, optical images sometimes cannot be used (e.g., lack of cloud-free periods) or are at the limit of their utility (e.g., optical reflectance provides information only about the top layer of vegetation).In this context, synthetic aperture radar (SAR) images provide a reliable alternative to optical images because they are not greatly influenced by atmospheric conditions and can be acquired during the day or night.Several scattering mechanisms, including single bounce, double bounce, and volume scattering can contribute to SAR backscatter over distributed targets such as agricultural fields [15], depending on crop properties (e.g., biomass, architecture, height) [16][17][18] and soil conditions (e.g., roughness, moisture) [19][20][21].Unlike optical data, however, SAR data are rarely used to map land cover.Some studies have focused on the complementarity of SAR and optical data, concluding that using them together provided better results than using them separately [22].The SAR signal is sensitive to the geometry (e.g., roughness, texture, internal structure) and wetness of the observed targets, while optical reflectance is influenced by their physiology.Thus, these two spectral domains provide complementary information.
The recent SAR Sentinel-1 (S-1) and optical Sentinel-2 (S-2) time series provide a great opportunity to monitor forest-agriculture mosaics due to their high spatial and temporal resolutions, with a five-day revisit time and decametric resolution.They are freely available under an open license, unlike most SAR data (e.g., TerraSAR-X, ALOS2, RADARSAT-2, COSMO-Sky Med) and optical data (e.g., SPOT, Quickbird, WorldView, Geo-Eye, Ikonos).Recent studies have demonstrated the potential of S-2 data for mapping 6-12 land cover classes from a single-date image [23][24][25][26][27]. Studies that use S-2 time series focus on one or a few specific land cover classes, usually on homogeneous agricultural or forested landscapes located in the same climatic region [28][29][30][31][32][33][34].Use of S-1 data has been limited to combining them with S-2 or Landsat data, which increases classification accuracy [35][36][37][38][39]. Hence, the use of S-1 time series alone for mapping land cover classes has not yet been evaluated.
The objective of this study was to assess the potential of S-1 data alone, S-2 data alone, and combined S-1 and 2 data for mapping land cover in forest-agriculture mosaics.We focused on two landscapes with contrasting vegetation gradients: a temperate mountainous landscape in the Cantabrian Range (Spain) and a humid tropical forested landscape in Paragominas (Brazil).Although these forest-agriculture mosaics differ greatly, they perform similar ecological functions, such as conserving high levels of biodiversity and storing carbon.
Satellite images were classified using an incremental procedure based on the rank of importance of input features derived from S-1 and S-2 time series.The method automatically selects the relevant features (spectral bands and/or vegetation indices) and time periods to use to best classify land cover types in the forest-agriculture mosaics.

Study Area
The Spanish study area is located in northwestern Spain, extending from eastern Lugo to western Cantabria.Covering the entire Cantabrian Range, it has an area of 35,700 km 2 (Figure 1a).The region has an elevation of 0-2468 m, an Atlantic climate, and remnant forest fragments that are distinct from a non-forest matrix composed mainly of pastures and heathlands in abandoned meadows [40,41].Croplands constitute 22% of the study area and are composed of vines and cereals crops including wheat, barley, rye, sunflower and oats [42].These crops are located in the south of the study area (Castile-Leon), most of them are irrigated [43].In northern slopes of the Cantabrian range, oaks (Quercus robur and Q. petraea), beeches (Fagus sylvatica) and chestnuts (Castanea sativa) are the main species found in forested areas, while southern slopes are dominated by semi-deciduous oaks (Q.pyrenaica and Q. faginea) and evergreen oaks (Q.ilex), Pinus sp. and Eucalyptus globulus being relatively abundant [42].The major decrease in the forest area dates back to 3000 years BP, with the transformation of large patches of natural forests into pasturelands (traditional cattle grazing and selective timber extraction).More recently, additional forests were lost due to road construction, surface mining, increased frequency of anthropogenic fires, and creation of forest plantations.In the past 20 years, reserves were established that impose legal restrictions on new land uses [40].Many studies focused on habitat quality for brown bear and wood grouse, two native species of the Cantabrian Range, indicating that the forest cover is essential to their presence [41,44].Brown bear home ranges include forests, shrublands, rocky areas, and grasslands [41].
The second study area is the municipality of Paragominas, Pará State, in the eastern Brazilian Amazonia, 217 km south of Belem (Figure 1b).The municipality extends for 19,342 km 2 , has an elevation up to 190 m, and has a tropical climate.It experienced a continuous period of deforestation from the 1960s to 2010, due mainly to cattle ranching, soybean and corn cultivation, and creation of forest plantations.Overlogging and fires degraded the remaining forests greatly [45,46].The municipality was the main timber-producing region in Brazil in 1990.A biodiversity survey performed in the municipality showed that both landscape and within-forest disturbances contributed to biodiversity loss, with the greatest negative effects on species of high conservation and functional value [47].In 2005, in the context of the new federal policy to fight deforestation, Paragominas was asked to establish measures against deforestation.The municipality established a new local governance model, called "Municipio Verde" [48], which aimed to fight deforestation and strengthen the ability of local institutions to develop specific environmental policies based on the Brazilian forest code.Since then, annual deforestation rates have decreased by up to 80%.Paragominas represents a case of remarkable deforestation that became a national reference for municipal-level anti-deforestation policies.

Reference Data
A set of 828 polygons was manually discriminated on the Spanish study area using aerial photographs and the Forest Map of Spain [49], which was developed from 1997-2006 and is the official forest inventory of the country.Its minimum mapping unit is 2.25 ha.The map is composed of polygons containing land cover and land use attributes derived from interpretation of aerial photographs and field inventory data.The 8 land cover classes discriminated in the sample data were permanent bare soils, artificial surfaces, water bodies, forested areas, shrublands, permanent herbaceous vegetation (herbaceous vegetation containing chlorophyll throughout the year), summer herbaceous vegetation, and winter herbaceous vegetation.The total area of samples per class ranged from 260-290 ha.
A set of 328 GPS points was recorded during a 3-week field mission in Paragominas in September 2017.The samples were produced by calculating polygons around the GPS points and using aerial photographs for photo-interpretation.The seven land cover classes discriminated were bare soils, artificial surfaces, water bodies, forested areas, croplands, pastures, and "young secondary forests" (fallow land with dense but low vegetation representing early regeneration stages of forests after abandonment of agriculture or pasture ).The total area of samples per class ranged from 14 ha (bare soil) to 530 ha (forested area).

Sentinel-1 Time Series
S-1 Ground Range Detected (GRD) products were acquired in Interferometric Wide Swath mode and delivered with VV and VH polarizations and an incidence angle ranging from 30.6-46.0 0 (Table 1).The GRD products provide multilook intensity (5 and 1 looks according to slant range and azimuthal direction, respectively) and are projected to a ground range based on an Earth ellipsoid model [50].The range and azimuth spatial resolutions were 20 and 22 m, respectively, and the pixel spacing was 10 ˆ10 m (Table 1).Sets of 66 and 42 S-1A images were downloaded for the Spanish and Brazilian study areas, respectively.S-1A images were acquired on two orbits (five-day and 12-day revisit time for Spain and Brazil, respectively), and 3 images were required to cover each study area.In total, 22 mosaics were produced from December 2016 to September 2017 for the Spanish study area, while 14 mosaics were produced from November 2016 to December 2017 for the Brazilian study area.

Sentinel-2 Time Series
The S-2 sensor acquired optical images during the same periods, with spatial resolutions of 10 and 20 m, and a spectral resolution of 10 bands (Table 2).For the Spanish study area, a series of 6 S-2 mosaics was acquired from December 2016 to August 2017.The tiles were downloaded in level 2A, which provides Top Of Canopy (TOC) reflectances and a cloud and shadow mask [51].Eight S-2 tiles were necessary to cover the entire study area (7-43 days between two tiles).Tiles with minimum cloud cover were selected.The mosaics were projected onto the ETRS89/ETRS-TM30 system (EPSG 3042).
For the Brazilian study area, only one mosaic was produced due to the frequent heavy cloud cover (July 2017).The tiles were downloaded in level 1C, which provides Top Of Atmosphere reflectances and orthorectified images [51].Nine S-2 tiles were necessary to cover the entire study area (no more than 12 days between 2 tiles).The mosaic was projected onto the WGS84/UTM 22S system.Table 2. Main characteristics of Sentinel-2 L1C/ L2A images [51].

Methodology
The method was developed to assess the potential of S-1 data alone, S-2 data alone, and combined S-1 and 2 data for mapping land cover of forest-agriculture mosaics (Figure 2).We focused on automatically selecting the relevant features (spectral bands and/or vegetation indices) and time periods to classify land cover types.The method used to select and classify features into land cover types was based on the incremental classification developed by Inglada et al. [35].

Samples Selection
A set of 30 pairs of training and validation polygons were randomly selected from the reference polygons with a 50/50 ratio.The 30 random pairs were performed in order to obtain 30 trainings and 30 corresponding validations, allowing to compute average performances with confidence intervals [35].The 50/50 ratio allows a more reliable comparison between training and validation samples than a ratio with a lower proportion of validation samples.For each set of selected training polygons, 300 and 1000 pixels per class were randomly selected for the Brazilian and Spanish study areas, respectively.The same random selection of pixels was applied to each set of associated validation polygons.This selection was performed to avoid over-representing large polygons at the expense of small polygons in the training and validation samples and to have the same number of samples for each class.

Preprocessing
The S-1 images were preprocessed using the Sentinel-1 Toolbox (Sentinel-1 toolbox, ESA, http: //step.esa.int/main/toolboxes/sentinel-1-toolbox/).The images were first radiometrically calibrated to extract the backscattering coefficients (σ 0 VV, σ 0 V H) (Figure 2).A speckle filter was then applied (Lee refined, with a 7 ˆ7 window, [52]).The images were geo-coded using Shuttle Radar Topography Mission data to correct topographic deformations (geometric correction accuracy < 1 pixel) and were projected onto the ETRS89-TM3/ETRS0 system (EPSG 3042) and the WGS84/UTM 22S system (EPSG 32722) for the Spanish and the Brazilian study areas, respectively.A ratio channel was produced (VV/VH) from the backscattering coefficient images.All images were then converted from linear to decibel (dB).
Spectral bands and indices were derived from S-2 images.For the Brazilian study area, the Sen2Cor application (Sen2Cor, ESA, http://step.esa.int/main/third-party-plugins-2/sen2cor/) was used to transform L1C tiles to level L2A (TOC).Four vegetation indices commonly used in remote sensing were calculated for each S-2 mosaic for both study areas: EVI (Enhanced Vegetation Index), NDVI (Normalized Difference Vegetation Index), NDWI (Normalized Difference Water Index), and SAVI (Soil-adjusted Vegetation Index) (Table 3).Each mosaic was composed of 10 spectral bands and 4 vegetation indices.

Feature Selection and Classification
Incremental classification was performed using a random forest (RF) classifier.Two output analyses-mean rank of importance of each input feature and mean kappa index as a function of the number of features-were used to select the relevant features to use in the final classification.Incremental classification analyzes mapping quality as the type and number of features used are added, thus, determining at what combination and number of features the classification reaches an acceptable quality [35].For each training and validation set (30 pairs), a RF algorithm was first applied to all of the features (dates and bands) to rank them in order of importance based on the mean decrease in the Gini index.Calle and Urea [57] demonstrated more robust results for exploring ranking stability using the mean decrease in the Gini index instead of the mean decrease in accuracy.The mean rank of importance of each feature was derived from the 30 ranks obtained from the 30 pairs of training and validation samples.We then ran as many RF algorithms as the number of features, starting with the two most important features and then adding the less important features until all of the features were processed.The kappa index was calculated for each classification to estimate the accuracy of the land cover classifications [58].Rosenfield and Fitzpatrick-Lins [59] advised using the kappa index to measure classification accuracy.
We chose the RF algorithm [60] due to its advantages of having low sensitivity to feature selection, simple parameterization, and short calculation time [61,62].We used it to select and classify features.In both cases, the number of trees was set to 100 to reduce calculation time sincePelletier et al. [63] demonstrated that it can be set to 100 without a major loss in accuracy.The experiments were performed using the randomForest package of R software (randomForest package version 4.6-14, Andy Liaw, https://cran.r-project.org/web/packages/randomForest/index.html).We used the OrfeoToolbox to calculate the final RF prediction (OrfeoToolbox version 6.6.1.,CNES, https://www.orfeo-toolbox.org/).

Percentage of Pixels Confused
For both study areas, the features selected from analyses of S-1 data alone and S-2 data alone were used to calculate 30 RF classifiers.For each pair of classes, we calculated the percentage of confused pixels in relation to the total number of misclassified pixels.The following equation: where i and j are the classes in the pair, was applied to the 30 confusion matrices derived from the 30 RF classifications using the 30 pairs of training and validation samples.The mean of the percentages was calculated for each pair of classes.We then calculated the mean percentage of confused pixels per pair of classes in relation to the total number of misclassified pixels using the most important S-1 features and most important S-2 features.

Comparison of Classifications
We used McNemar's X 2 test to analyze the significance of the difference or resemblance between the classifications derived from S-1 data alone vs. S-2 data alone, S-1 data alone vs. the combined S-1 and 2 data, and S-2 data alone vs. the combined S-1 and 2 data.This test, based on a binary 2 ˆ2 contingency matrix, shows the proportion of pixels correctly and incorrectly classified in two classification runs and allows the use of non-independent samples.A X 2 value lower than 3.14 means that the two classifications compared are non significantly different [64].

Contribution of Sentinel 1 & 2 Time Series to Map Land Cover
To estimate the contribution of the optical and SAR Sentinel time series, we analyzed the mean kappa index as a function of the number of input features.The combined use of S-1 and 2 data provided the best results for both study areas (Figure 3).However, for the standard deviation of the kappa index, results of the combined use of S-1 and 2 data were similar to those of S-2 data alone.The improvement in classification due to combining S-1 and 2 data rather than using S-2 data alone was perceptible up to 30 input features for the Spanish study area.For the Spanish study area, using the top ten S-2 features out of the 84 total S-2 features increased the mean kappa index by ca.0.25 (i.e., from 0.58 to 0.83 ); the mean kappa index continued to increase slightly up to 40 input features and then stabilized.For the Brazilian study area, using the top 7 S-2 features out of the 14 total S-2 features increased the mean kappa index by ca.0.15 (i.e., from 0.70 to 0.85); the mean kappa index continued to increase slightly up to 9 input features and then stabilized.Using S-1 data alone yielded the lowest accuracy for both study sites.The mean kappa index was higher for the Spanish area than for the Brazilian area, as was the maximum kappa index (ca.0.73 and 0.60, respectively).Compared to the mean kappa indices of S-2 data alone, those of S-1 data alone were lower by 0.08-0.20 for the Spanish study area and by 0.28-0.40for the Brazilian study area.

Importance of Input Features
When using S-2 data alone, bands 11 and 12 (the short-wave infrared (SWIR) domain) were the most important features for both study areas (mean rank of importance of 1.77 (B11) and 4.07 (B12), and of 1.43 (B12) and 1.63 (B11) for the Spanish and Brazilian study areas respectively) (Figure 4).For the Spanish area, the vegetation indices were relevant, with 11 of them among the top 20: 4 SAVI and 4 NDVI (1 each in January, March, June, and July), 2 NDWI (in January and March) and, although less relevant for classifying forest-agriculture mosaics, EVI.For the Brazilian study area, NDWI was the most relevant index (ranked 3rd), while other spectral bands were more relevant than EVI, NDVI, or SAVI.
When assessing the relevance of time periods using S-2 data alone (Figure 4a), January was the most important month for discriminating the 8 land cover classes, since it was the month of the top 6 features (mean rank of importance = 1.77-11).In addition, March and July appeared among the top ten features.The least important months were August and April.When using S-1 data alone, the VV/VH ratio was not useful for mapping land cover in these two forest-agriculture mosaics.For both study areas, VV polarization was the most important feature.For the Spanish and Brazilian study areas, the VV mean ranks of importance were of 2.03 and 2.13, respectively (Figure 5).The top 20 features, however, had a nearly equal number of VV and VH features (10 VV and 9 VH polarization dates for the Spanish study area, 11 VV and 9 VH polarization dates for the Brazilian study area).When assessing the relevance of time periods using S-1 data alone, the months from December to May were relevant for the Spanish study area.In contrast, months from both dry and wet seasons were among the top ten features (January, March, April, July, September, November, December) for the Brazilian study area.
When combining S-1 and 2 data, the top 20 important features were all S-2 features for the Spanish study area and mostly S-2 features (the top 14) for the Brazilian study area (Figure 6).As a trade-off between accuracy and processing time, we decided to maintain a minimum but sufficient number of features for classification (Table 4).In the following experiments, we choose to keep the 20 most important S-1 features and 10 S-2 features, and the 10 most important S-1 features and 7 S-2 features on the Spanish and Brazilian study areas, respectively (Table 4).Due to the large majority of S-2 features in the top 20 important features of combined S-1 and 2 data (Figure 6), we decided to combine the previously selected features of S-1 data alone and S-2 data alone to predict the land cover classes combining S-1 and 2 data (Table 4).

Confusion between Classes
Using the S-1 features for the Spanish study area, forested areas vs. shrublands, shrublands vs. permanent herbaceous vegetation, and bare soils vs. winter vegetation were the pairs of classes most often confused (Figure 7a).Using the S-2 features, bare soils vs. artificial surfaces, bare soils vs. winter vegetation, and forested areas vs. shrublands were the most confused pairs of classes (Figure 7b).Forested areas and shrublands were not well discriminated using S-1 or S2 features.Shrublands vs. permanent vegetation were better discriminated using S-2 data than S-1 data.Using the S-1 features for the Brazilian study area, most forested areas and young secondary forests pair of classes was confused, unlike the other pairs of classes, including forested areas vs. pastures (Figure 8a).Using the S-2 features, cropland vs. bare soils, pastures vs. young secondary forests, and forested areas vs. young secondary forests were the most confused pairs of classes (Figure 8b).Forested areas vs. young secondary forests were better discriminated by S-2 than S-1 features.However, pastures and young secondary forests were better discriminated by S-1 than S-2 features.On both study areas, bare soils vs. artificial surfaces were discriminated better by S-1 data than S-2 data features.

Prediction of Selected Features
Misclassification errors in the land cover map of forest-agriculture mosaics of the Spanish study area-due to cloud confusions using S-2 data alone and slopes using S-1 data alone-were partly corrected using combined S-1 (Figure 9).The combined S-1 and 2 classification of the Brazilian study area (Figure 10) was more accurate than classification by S-2 data alone.Misclassification of cropland, which was confused with artificial surfaces using S-2 data alone, was partly corrected with combined S-1 and 2 data.

McNemar χ 2 Test Results
According to McNemar's χ 2 test, all land cover class predictions differed significantly (χ 2 > 3.14), and S-2 data alone vs. combined S-1 and 2 data for the Brazilian study area is not significant (p-value = NS) (Table 5).Thus, we concluded that all classifications differed.

Relative Contributions of S-1 and S-2 Data to Map Land Cover of Forest-Agriculture Mosaics
The land cover classes that comprise forest-agriculture mosaics are distributed along a landscape gradient, with transition classes such as shrublands to forest or pasture to young secondary forests, which result in misclassifications.S-1 data alone were least accurate for mapping land cover of forest-agriculture mosaics, with a best mean kappa index of 0.73, vs. 0.87 and 0.89 for S-2 data alone and combined S-1 and 2 data, respectively (Figure 3).For example, S-2 data discriminated forested areas from young secondary forests better than S-1 data for the Brazilian study area (Figure 8).This may be because SAR sensors are sensitive to vegetation structure, while optical sensors are sensitive to chlorophyll content.Thus, when observed at a 10m resolution, broad vegetation classes were discriminated more easily by their physiology than their physical structure.Also, VV and VH polarizations were the most important S-1 features for discriminating land cover classes, while the VV/VH ratio was not useful (Figure 5).While it is known that VH polarization is more sensitive to vegetation than VV polarization [65] (Patel et al., 2006), results show that both polarizations had an equivalent low contribution to classify forest-agriculture mosaics, as well as the ratio VV/VH.We can conclude that in this case the C-band of Sentinel-1 was not relevant to classify vegetation classes.In general, the C-band SAR is less suitable than L-Band for forest change monitoring due to a lower penetration depth and a rapid saturation of the signal [66]).For example, Patel et al. [65]) demonstrated that L-band is more sensitive to plant density than C-band, and that C-band interacts mostly with the primary branches of the canopy while L-band within the vegetation canopy.We used S-1 GRD products that record backscattering coefficients (VV and VH) because they require less processing time than Single Look Complex products, which preserve the phase information.However, using texture, coherence, and polarimetric indices derived from full-polar radar data, such as RADARSAT-2 data, could improve classification based on S-1 data alone [35,36,[67][68][69][70]. Thus, although Bagdhadi et al. [71] concluded that using polarimetric indices (Shannon Entropy (SE) and span) derived from polSAR RADARSAT-2 images did not improve estimates of soil moisture and vegetation parameters, Betbeder et al. [72] demonstrated the greater utility of the SE index than backscattering coefficients for mapping wetland vegetation using dual-pol TerraSAR-X time series.The best mean kappa index achieved with S-1 data was higher for the Spanish study area than for the Brazilian study area (0.73 vs. 0.60, respectively) (Figure 3).This could be due to differences in the classification methods (different numbers of classes and land cover types) and to the number of S-1 images used to classify land cover types (22 dates for the Spanish study area vs. 14 dates for the Brazilian study area).For the Spanish study area, the classification results achieved using S-1 data alone show that misclassification errors are most often located on mountain slopes.These errors are due to the acquisition of SAR images in slant range geometry that causes layover and foreshortening effects [15].Moreover, it is known that soil moisture and roughness influence the radar backscatter depending on frequency, polarization and the incidence angle of the incoming microwave.Holah et al. [73] found that the sensitivity to surface roughness is more important using HH and HV polarizations than VV polarization.However, Baghdadi et al. [74] demonstrated that the soil moisture was not very dependent on polarization.The effects of soil moisture on forest-agriculture classification are rather limited using C-band due to a low penetration depth of the microwaves compared to L-band [75], however C-band is more sensitive to roughness [76].Lastly, the higher the incidence angle, the more important the sensitivity of the SAR to surface roughness [77].
However, S-1 data sometimes discriminated better than S-2 data-for bare soils and artificial surfaces-(Figures 7 and 8) because the SAR signal reacts differently to these two land cover classes, with double-bounce on buildings and single-bounce on flat bare soils [15], while their spectral signatures in the optical domain are similar (high reflectance values in green, blue and red bands).When using S-2 data alone, SWIR bands were the most important features for discriminating land cover classes (Figure 4).The importance of S-2 SWIR bands has been demonstrated for mapping vegetation and forests [23,78].Vegetation indices derived from S-2 data (e.g., SAVI, NDVI, NDWI) discriminated better than spectral bands for the Spanish study area, while NDWI was the most important vegetation index for the Brazilian study area.The sensitivity of the blue band to water contained in vegetation could explain why NDWI is well suited to tropical vegetation, while NDVI saturates at high biomass values [55,79].The EVI, which was developed for MODIS, was a less important feature than other vegetation indices for both study areas.It could be interesting to calculate other vegetation indices that are used for crop discrimination such as S2REP (Sentinel-2 Red-Edge Position), IRECI (Inverted Red-Edge Chlorophyll Index), MTCI (MERIS Terrestrial Chlorophyll Index) or SRI (Simple Ratio Index) [80,81].
The best mean kappa index using combined S-1 and 2 data (0.88) did not differ greatly from that using S-2 data alone (0.86) (Figure 3).Thus, the use of S-2 data alone was sufficient to discriminate the land cover classes with high accuracy.However, the backscattering intensity from S-1 data alone provided additional relevant information, since predictions of land cover classes using S-1 data alone, S-2 data alone, and combined S-1 and 2 data differed according to McNemar's χ 2 test (Table 5).For example, the misclassification between cropland and artificial surfaces for the Brazilian study area when using S-2 data alone was partly corrected using combined S-1 and 2 data (Figure 10).

Using S-1 and S-2 Data to Identify the Key Time Periods for Classifying Land Cover
Based on the classification results using S-1 data alone, S-2 data alone, and combined S-1 and 2 data, classification accuracy was strongly related to the number of dates.Classification accuracy increased when dates were added to the RF classifier (Figures 3, 4a, 5 and 6).The relevance of high temporal resolution underlines the importance of describing and taking into account the phenology of vegetation to map forest-agriculture mosaics.For the Spanish study area, January, March, June, and July features are among the top 20 important S-2 features (Figure 4a).This may have been because the S-2 sensor is sensitive to chlorophyll content in vegetation, and these months are key periods of chlorophyll activity.The features of the months from December to May were among the top 20 S-1 features for this study area, while the summer period was not relevant (Figure 5a).This is probably because SAR sensors, which are sensitive to the internal structure of elements, penetrate vegetation better during the leaf-off period.Only one S-2 date was used for the Brazilian study area, which precluded studying the importance of time period.In addition, no clear time periods emerged from the importance ranks of the S-1 time series (Figure 5b).Thus, no intra-annual time period seems more important than others for S-1 time series for discriminating land cover types in Paragominas.Unlike the Cantabrian Range, where the seasons differ, temperatures and precipitation are relatively constant in Pará, with a drier season from July to October [82].In addition, most of the dynamics of human practices in pastures, forested areas, and plantations are inter-annual rather than intra-annual [83].

The Robustness of the Method for Different Landscapes
The results show that the method can be applied to forest-agriculture mosaics in two landscapes with contrasting climates (i.e., tropical, temperate) and vegetation types (e.g., tropical forests, shrublands, coniferous and deciduous forests).S-1 and S-2 features were selected according to their respective importance in each study area, which highlighted relevant features and time periods for the two areas.Land cover maps of the two study areas were produced with a high degree of accuracy and a short processing time.Time series are used to map the land cover of forest-agriculture mosaics, but processing them requires high computational capacity.The method enabled focusing on specific time periods and features to reduce the time required for image processing.

Conclusions
An incremental procedure based on ranking the importance of the input features derived from S-1 and S-2 time series was used to discriminate land cover classes in forest-agriculture mosaics.The method automatically selected relevant features (spectral bands and/or vegetation indices) and time periods to classify land cover types in such landscapes.S-2 data alone were more relevant than S-1 data alone for mapping the land cover of forest-agriculture mosaics, and combining S-1 and 2 data slightly improved results compared to those of S-2 data alone.Using polarimetric indices, such as the SE index, which has already shown potential for characterizing vegetation, can improve predictions of S-1 data alone.Indeed, SAR data are useful in cloudy regions; the high cloud cover in the S-2 time series was the main source of misclassification errors.

Figure 1 .
Figure 1.Location of (a) the Spanish study area ( c EuroGeographics for the country boundaries) and (b) the Brazilian study area ( c 2018 GADM for the state boundaries).

Figure 3 .
Figure 3. Mean kappa index of the incremental classifications for Sentinel-1 data alone (green), Sentinel-2 data alone (red), and combined Sentinel-1 & 2 data (blue) as a function of the number of input features for (a) the Spanish study area and (b) the Brazilian study area.

Figure 4 .
Figure 4. Mean rank of importance of the most important features for Sentinel-2 data alone for (a) the Spanish study area and (b) the Brazilian study area.

Figure 5 .Figure 6 .
Figure 5. Mean rank of importance of the most important features for Sentinel-1 data alone for (a) the Spanish study area and (b) the Brazilian study area.

Figure 7 .
Figure 7. Mean percentage of pixels confused per pair of classes in relation to the total number of misclassified pixels using (a) the top 20 Sentinel-1 features and (b) the top ten Sentinel-2 features of the Spanish study.

Figure 8 .
Figure 8. Mean percentage of pixels confused per pair of classes in relation to the total number of misclassified pixels using (a) the top ten Sentinel-1 features area and (b) the top 7 Sentinel-2 features of the Brazilian study area.

Figure 9 .
Figure 9. Classification of the Spanish study area using combined Sentinel-1 & 2 data.

Figure 10 .
Figure 10.Classification of the Brazilian study area using combined Sentinel-1 and 2 data.

Funding:
This research was funded through the 2015-2016 BiodivERsA COFUND call for research proposals, with the national funders ANR, MINECO, and BELSPO.A. Mercier has a Ph.D. grant from MESRI (Ministry of Higher Education, Research and Innovation of France), PhD grant 2017.The study in Brazil was sponsored by the ODYSSEA Project "Observatory of the dynamics of interactions between societies and environment in the Amazon", financed by the European Commission, Horizon 2020, and supported by the French Space Agency CNES and TOSCA (Terre, Océan, Surfaces continentales, Atmosphère) program CASTAFIOR project (Caractérisation et dynamique des Agro-écoSystèmes Tropicaux Amazoniens par Fusion d'Image Optique et Radar).

Table 4 .
Number of top important features selected for predictions for the Spanish and Brazilian case studies.

Table 5 .
McNemar's χ 2 test results for the Spanish and Brazilian case studies.The χ 2 value less than 3.14 indicates that the two classifications were not significantly different.