Integrated use of Sentinel-1 and Sentinel-2 data and open-source machine learning algorithms for land cover mapping in a Mediterranean region

ABSTRACT This paper aims to develop a supervised classification integrating synthetic aperture radar (SAR) Sentinel-1 (S1) and optical Sentinel-2 (S2) data for land use/land cover (LULC) mapping in a heterogeneous Mediterranean forest area. The time-series of each SAR and optical bands, three optical indices (normalized difference vegetation index, NDVI; normalized burn ratio, NBR; normalized difference red-edge index, NDRE), and two SAR indices (radar vegetation index, RVI; radar forest degradation index, RFDI), constituted the dataset. The coherence information from SAR interferometry (InSAR) analysis and three optical biophysical variables (leaf area index, LAI; fraction of green vegetation cover, fCOVER; fraction of absorbed photosynthetically active radiation, fAPAR) of the single final month of the time-series were added to exploit their correlation with the canopy structure and improve the classification. The random forests (RF) algorithm was used to train and classify the final dataset, and an exhaustive grid search analysis was applied to set the optimal hyperparameters. The overall accuracy reached an F-scoreM of 90.33% and the integration of SAR improved it by 2.53% compared to that obtained using only optical data. The whole process was performed using freely available data and open-source software and libraries (SNAP, Google Earth Engine, Scikit-Learn) executed in Python-script language.


Introduction
Mapping the vegetation composition, besides providing information on the quantitative and qualitative status of the area under study, is a necessary early step in the analysis and monitoring protocols of the state of the vegetation and ecosystems' responses affected by various environmental disturbances (Choudhury et al., 2021;Grabska et al., 2020;Monroe et al., 2020;Pollino & Modica, 2013;Rodman et al., 2021;Semeraro et al., 2019), including wildfires (Gitas et al., 2012), storms (Giannetti et al., 2021;Hamdi et al., 2019), deforestation (Nicolau et al., 2021), forest degradation (Modica et al., 2015), desertification (Hill et al., 2008) and climate change effect (Yang et al., 2013). Therefore, mapping the composition of forest vegetation is fundamental for the concrete implementation of sustainable land management policies at any scale, regional to global (e.g., the REDD + activities; Gulinck et al., 2018;Nicolau et al., 2021).
In the context of vegetation mapping and monitoring, several remote sensing techniques based on different types of multispectral sensors have been developed and successfully used over the years (Grabska et al., 2019;De Luca et al., 2019;Modica et al., 2016;Morin et al., 2019;Praticò et al., 2021;Solano et al., 2019). The use of spectral signatures, temporally differentiated following the phenological cycles of the various seasons, allows a better spectral separability of the investigated vegetation types and, therefore, their recognition and characterization (Aragones et al., 2019;Aubard et al., 2019;Grabska et al., 2020Grabska et al., , 2019Morin et al., 2019;Praticò et al., 2021). Grabska et al. (2019) used a Sentinel-2 (S2) time-series to map forest composition showing the effectiveness of seasonal phenology variations in improving spectral discrimination between species, achieving better accuracy results than using single images. Moreover, the spectral vegetation indices (VIs) enhance the sensibility of single-bands spectral signals to the variability of the bio-physical state of plant tissues, the photosynthetic activity, and leaf productivity (Aragones et al., 2019;Marzialetti et al., 2019;Praticò et al., 2021;Semeraro et al., 2019). Strong correlations were found between specific regions of the electromagnetic spectrum and species-specific physiological characteristics useful in estimating forest cover, especially using VIs based on infrared wavelengths: the normalized difference vegetation index (NDVI; Marzialetti et al., 2019;Spadoni et al., 2020), the normalized difference red-edge index (NDRE; Evangelides & Nobajas, 2020), the normalized burn index (NBR; Praticò et al., 2021;Shaun et al., 2020). The red-edge, near infra-red (NIR) and short-wave infra-red (SWIR) regions have, respectively, a long-established correlation to the leaf pigments content, vegetation net primary productivity, and leaf water content, being very effective in vegetation monitoring (Arevalo-Ramirez et al., 2020;Delegido et al., 2011;Eitel et al., 2011;Knipling, 1970).
The free availability of the higher temporal and spatial resolution Copernicus S2 mission multispectral data (ESA Sentinel Homepage, 2021), provided by the European Space Agency (ESA), improved the accuracy of forest cover classification maps and allowed for the launch of several successful monitoring studies at a higher scale of detail (Grabska et al., 2020;Immitzer et al., 2016;Inglada et al., 2017;Praticò et al., 2021;Solano et al., 2019).
In addition to the use of multispectral data, several authors studied the applicability of active synthetic aperture radar (SAR) systems for mapping land cover (Lapini et al., 2020;Nicolau et al., 2021;Perko et al., 2011;Waske & Braun, 2009). Besides the all-day and all-weather operational capability, these sensors provide different and complementary physical information helpful for improving the spectral data when combined with optical imagery (Nicolau et al., 2021;Spracklen & Spracklen, 2021;Stroppiana et al., 2015). The total signal backscattered from forest vegetation results from the combination and interaction of the canopy and ground backscatters (Lapini et al., 2020;Saatchi, 2019;Yu & Saatchi, 2016). This backscatter response is affected by implicit sensor variables, such as wavelength and polarization, and by some vegetation features as cover shape, structure, and orientation, moisture content, geometric and dielectric property of the surface (Lapini et al., 2020;De Luca et al., 2021a).
The Copernicus mission provides two polar-orbiting SAR satellite platforms belonging to the Sentinel-1 (S1) constellation (S1-A and S1-B) carrying a C-band sensor (wavelength of 5.6 cm) with both cross-polarized (VH) and co-polarized (VV) polarization (ESA Sentinel Homepage, 2021). At these wavelengths, the backscatter is mainly due to the leaves, needles, and small branches of the upper canopy and presents lower penetration power than longer wavelengths (Lapini et al., 2020). Potentially, the information from the upper canopy could allow the discrimination between forest and non-forest areas.
Referring to forest applications, recently, Nicolau et al. (2021) assessed the potential of S1 time-series for land use/land cover (LULC) purposes in tropical forests, while Numbisi et al. (2019) utilized S1 time series to discriminate agroforests environments in a heterogeneous savannah-forest transition zone. On the other hand, Lapini et al. (2020) assessed the multi-frequency approach for Mediterranean forest classification, discriminating forest from non-forest areas and broadleaved from coniferous forests, using data from different SAR sensors (X-, C-and L-band). These authors concluded that the L-band is better for the first purpose, but C-band and X-band performed better for distinguishing coniferous and broadleaves.
The utility of SAR signal in forest vegetation discrimination can be also explicated by its particular sensitivity to the forest stand height (Deutscher et al., 2013;Li et al., 2020;Perko et al., 2011;Siqueira, 2019). The simple SAR backscatter is indirectly and empirically related to the forest stand height since its value increases with a high presence of canopy scattering elements, proportional to forest height (vertical distribution) and canopy density, as a function of wavelength and polarization. Moreover, there is a geometric relationship between the SAR signal and the heights of the objects on the Earth's surface, estimable through SAR interferometry (InSAR) models (Siqueira, 2019). The InSAR technique exploits the phase information of the radar signal to obtain information about the topography and height of the Earth's surfaces (Ferretti et al., 2007;Ghosh et al., 2020). The S1 constellation observes the same scene at two different times, applying the repeat-pass InSAR. The amount of temporal phase decorrelation occurring between two passes is one of the models used to estimate the forest stand height. The temporal decorrelation is assumed to be higher the greater the height of the canopy due to a more significant presence of small scatter elements (Siqueira, 2019). The interferometric coherence can represent the temporal phase decorrelation: the higher is the time phase decorrelation, the lower is the resulting coherence. Several authors (Deutscher et al., 2013;Ghosh et al., 2020;Perko et al., 2011;Siqueira, 2019) applied empirical models to estimate the forest stand height from the interferometric coherence measure, with levels of accuracy that can vary greatly depending on various factors. For this reason, in this study, the coherence relationship with the forest stand height was exploited to discriminate the presence of standing forest concerning the other surrounding LULC classes.
Considering the research experiences mentioned above, the combined use of both optical and SAR data would further improve the identification of forest cover, as confirmed by several authors (Ienco et al., 2019;Morin et al., 2019;Polychronaki et al., 2014;Spracklen & Spracklen, 2021;Zhang et al., 2019). Spracklen and Spracklen (2021) used S1 and S2 timeseries to distinguish natural and plantation forests in a tropical monsoon climate zone, concluding that the different sensitivity of these two sensors makes them complementary in analyzing the investigated vegetation surface. In particular, while the SAR backscatter depends on the vegetation's physical properties, the optical signal is correlated to the biochemical state of vegetation. Several studies use combined data S1 and S2 data for vegetation cover purposes. Ienco et al.
(2019) proposed a Convolutional Neural Network (CNN) architecture, combining S1 and S2 time-series for LULC mapping in tropical regions, obtaining satisfactory results. Zhang et al. (2019) used the differences in temporal signatures between vegetation cover types, combining three temporal information of S1 and S2, for distinguishing woody canopy from the herbaceous canopy in savanna ecosystems using the Support Vector Machines (SVM) classifier. Morin et al. (2019) combined the use of S1, S2 and ALOS-PALSAR data to estimate forest structure parameters and the aboveground biomass in maritime pine plantations. In Mediterranean environment, besides the literature concerning the single use of optical sensors (e.g., Aragones et al., 2019;Aubard et al., 2019;Modica et al., 2016;Praticò et al., 2021), there seems to be a lack of studies exploiting the SAR (Lapini et al., 2020) and its integration with optical imagery for LULC classification. Polychronaki et al. (2014), integrated the Système Pour l'Observation de la Terre (SPOT) optical data with the European Remote Sensing (ERS) C-band VV for LULC object-based classification affected by a fire in a Mediterranean landscape. They observed that the use of SAR backscatter improved the accuracy, reducing the commission errors related to forest land cover class and the misclassification between forest and shrub classes in shadowed areas. Chust et al. (2004) assessed the performances of combining ERS and SPOT images for Mediterranean LULC discrimination, resulting in a slight improvement of the obtained accuracy. Chatziantoniou et al. (2017) evaluated the combined use of Sentinel-1 and Sentinel-2 data for a regional-scale Mediterranean wetlands classification, and they concluded that SAR data did not significantly improve classification accuracy. Some other authors focused on Mediterranean crop detection (Campos-Taberner et al., 2017;Lobo et al., 1996;Villa et al., 2015). Considering this, more studies should be carried out to explore the potential of integrating multispectral and SAR sensors on mapping heterogeneous Mediterranean ecosystems and to investigate how the single information contributes to the obtained accuracy.This work aimed to develop a supervised classification procedure by integrating S1 and S2 data for forest cover mapping in a Mediterranean area of southern Portugal. The obtained results are essential to fulfilling the main research framework in which this work was conducted, based on applying remote sensing methodologies to analyze and monitor wildfires' effects in Mediterranean forest ecosystems. The first analysis on this study area was carried out to develop an unsupervised classification of the burned areas and based only on SAR S1 data and reported in De Luca et al. (2021a), while the combined use of SAR S1 and optical S2 data allowed to map the spatial distribution of burn severity (De Luca et al., 2021b).
The initial goal was to create a binary map to distinguish the forest cover from other LULC classes (pastures/shrubs, urban, agricultural, etc.). Afterward, we decided to implement this workflow further by subdividing the forest ecosystems into three forests LULC classes better representing the territory: Eucalyptus, Pine, and native broadleaf forest (Quercus suber, Q. ilex, etc.). In the same way, the main surrounding noforest LULC classes were classified individually (pastures/shrubs, bare soil, urban and agricultural). For this purpose, we provided an original and open workflow (i. e., implemented using diverse open-source software and freely available upon a reasonable request to the authors) based on an advanced coupling of SAR (S1) and multispectral (S2) time-series imagery. In particular, the time-series of SAR S1 backscatter (both VH and VV polarizations) and two derived indices, radar vegetation index (RVI) and radar forest degradation index (RFDI), were combined to the time-series of the optical S2 bands and three derived VIs: NDVI, NDRE, and NBR. In order to optimize the classification procedure, the coherence measure coming from InSAR analysis of different pairs of dates in July 2018 was also added as additional information, as well as the optical-based biophysical variables fraction of green vegetation cover (fCOVER), the fraction of absorbed photosynthetically active radiation (fAPAR), and the leaf area index (LAI) calculated for the same month. The well-known Random Forests (RF) machine learning algorithm (Breiman, 2001) was applied for classification through the use of open-source and Python-based libraries (The Python Language Reference, 2021). One of the main problems in applying machine learning classification algorithms is choosing the optimal values of the model's hyperparameters. In this direction, another original contribution of this study was to provide an open workflow in which a thorough grid search approach automatically sets the optimal hyperparameters. The feature importance was performed during the RF classification process to evaluate each input variable's contribution to the final mapping performance.

Study area
The study area (Figure 1) extends over the Serra de Monchique mountain range located in the southern region of Portugal, Algarve (37° 18ʹN; 08° 30ʹW). Part of the study area is a Special Area of Conservation (SAC) falling within the European Natura 2000 network (Natura 2000 Site Code: PTCON0037). The territory is characterized by the typical heterogeneous and fragmented Mediterranean mountain landscape. The forest cover was mainly composed of Eucalyptus plantations (Eucalyptus globulus, Labill. 1800), mixed Mediterranean indigenous deciduous forests (Quercus suber L., Quercus ilex L., and other secondary Mediterranean native species), and coniferous plantations, composed by Pinus pinea L. and Pinus pinaster Aiton. A part of the autochthonous oak forest cover can be associated with the typical semi-natural landscape of the Iberian peninsula (dehesa and montado): woodlands and agro-forestry systems used for cork harvesting and grazing. A large part of the territory was covered by non-forest LULC classes represented by heathlands, sclerophyllous shrublands, pastures, bare soil, general uncultivated lands (e.g., derived from harvested forest plantations, agricultural and urban lands) (Sistema Nacional de Informação Geogrãfica, SNIG, 2021).

Materials and methods
The implemented procedure ( Figure 2) was carried out using different free and open-source software solutions, starting from download to the classification output and testing and exploiting their interoperability. Most of the Sentinel images were unavailable on the official Copernicus Open Access Hub platform due to the Long-Term Access policy adopted by the ESA (Copernicus Long Term Archive Access, 2021). Therefore, we adopted other alternative ways to speed up the data download phase. The S1 images were downloaded using the Alaska Satellite Facility (ASF) interface (ASF, 2021), while the Google Earth Engine (GEE) Python API (Google Earth Engine Guides, 2021) was employed to pre-process and then download the S2 dataset. The data pre-processing for S1 was carried out using the Sentinel-1 Toolboxes implemented in the SNAP v.8.0.3 open-source software (ESA SNAP Homepage, 2021) provided by ESA. The classification algorithms were processed using the modules integrated into the Scikit-learn Python library (Pedregosa et al., 2011).

Sentinel-1 dataset and pre-processing
The SAR dataset was composed of a time-series of S1-A/B ground range detected (GRD), acquired in interferometric wide (IW) mode, for each of the two available polarizations: co-polarized VV and cross-polarized VH. The time series comprised the period from April 2017 to July 2018, immediately prior to the fire event of August 2018, and it was composed of 82 images deriving from both ascending (42 images) and descending (40 images) flight paths. In order to perform the InSAR analysis for coherence extraction, a total of six Single Look Complex (SLC) format images (three for ascending and three for descending flight path), presenting respective same orbit-path and covering the month of July 2018, were downloaded. This month was chosen as it was the closest to the fire event, and, presumably, the vegetation conditions were more consistent with those at the time of the event. The images were in IW mode, which contains both amplitude and phase information of the backscattered SAR signal.

Sentinel-1 ground range detected (GRD) preprocessing
Concerning the GRD imagery, the application of the auto-downloaded orbit file and the thermal noise removal started the S1 pre-processing workflow. Subsequently, the images were first radiometric calibrated (β0) and then applied the radiometric terrain correction (RTC). The reduction of geometric and topographic errors was carried out by applying the radiometric terrain flattening and the terrain correction processes using the shuttle radar topography mission (SRTM) digital elevation model (DEM), characterized by a spatial sampling of 1 arc-second.
The correction of topographical errors, in particular the ground flattening process, is an essential step to map the land cover because it reduces the error due to the difference in the structure and shape of the ground surface and, therefore, in the backscatter values (Mendes et al., 2019;Small, 2011).
The bilinear interpolation resampling method was used for both DEM and output image resampling. All the S1 images were stacked using the product geolocation as the initial offset method. Since the speckle noise is an impactful error for land cover mapping purposes (Lapini et al., 2020)  5 × 5 pixel window size was applied, followed by a backscatter time (monthly) averaging (BTA). This last step further improves the image signal concerning speckle noise and minimizes the effects of environmental and seasonal variables (Lapini et al., 2020;De Luca et al., 2021a).

Sentinel-1 SLC InSAR pre-processing
For each of the two S1 platforms (A and B), all the pair combinations of the three SLC images were subjected to TOPS InSAR processing. Before applying the coregistration (Back-Geocoding) and Network Enhanced Spectral Diversity (NESD) optimization for each coregistered S1pair, we split the single images (TOPS-Split) and added the auto-downloaded orbit file. Then, we extracted the coherence information, characterized by a range value from 0 (minimum coherence) to 1 (maximum coherence). A final terrain correction was applied, using the 1 arc-second digital elevation model (DEM) derived from the Shuttle Radar Topography Mission (SRTM; Farr et al., 2007), where the pixel spacing was resampled to 10 m x10 m. The three coherence maps resulting for each of the two flight paths were averaged to reduce the speckle noise, obtaining one coherence map for ascending and one for descending paths.

Sentinel-2 dataset and pre-processing
The optical S2 Level-2A (Bottom-Of-Atmosphere, BOA) multispectral time-series was composed of 64 images, considering the same time period of the S1 dataset (April 2017 -July 2018). Subsequently, a monthly average for each band was performed, while a pixel size resampling of 10 m x10m was performed semi-automatically by the system during the download step by setting the scale parameter (pixel resampled size = 10) and the projection system (EPSG: 32,629), using the default nearest neighbor resampling algorithm. Each single-date image was masked by clouds before this last step, which generated an image for each month and each band from April 2017 to August 2018 (as the final SAR dataset). In this case, we used the S2-Cloud Probability product available in the GEE data catalog as a mask image upsampled to 10 m spatial resolution. The cloud probability mask product is created with the Sentinel2cloud-detector package (s2cloudless; "Sentinel Hub's cloud detector repository," 2021) for automated pixelbased cloud detection and it is based on the LightGBM machine learning library (LightGBM documentation, 2021), developed by Sentinel Hub's research team (Sentinel Hub Homepage, 2021). Each pixel contains a value between 0 and 100, representing the probability that the pixel is cloudy. Higher values are more likely to represent dense clouds or highly reflective surfaces but may omit less dense clouds. Lower values, although able to detect all clouds, could increase the risk that medium-high reflective surfaces could be mistaken for clouds (false positives). For our purpose, the value of 10 was considered the pixel's threshold value, greater than which a pixel is considered as a cloud. No-data values replaced the masked pixels of S2, filled by applying a temporal-linear interpolation between consecutive images. Water surfaces were also masked since three static water bodies characterized the study area (Barragem de Odelouca, Barragem do Funcho, and Barragem do Arade), and these were not targeted LCLUs. Effectively, the stable spectral features of the water surfaces, characterized by a significant absorption of most of the NIR wavelength radiation (Donchyts et al., 2016;Gao, 1996;Schwatke et al., 2019), made it simple to detect them by masking pixels of the S2 B8 band image (Jul2018) using the threshold <0.09.
Where t represents one of the months constituting the time-series (from April 2017 to July 2018). The final SAR dataset used in the classification process was composed of the BTA VV, BTA VH , RVI, RFDI timeseries, and the two coherence maps (ascending and descending) for July 2018.
It is expected that the SAR vegetation indices potentially improve the discrimination of land cover because of the combination of VH polarization, related to the scattering elements of the canopy, and the VV polarization, sensitive to the topographic and morphological caratheristics of the ground (Nicolau et al., 2021).
Moreover, three additional biophysical variables were automatically computed for July 2018 using the SNAP toolbox biophysical variables processor: LAI, fCOVER, and fAPAR. The criterion of this month's choice for biophysical variables estimation was the same as that already explained for the InSAR analysis. The quantitative-qualitative conditions of the vegetations should be more similar to those at the time of the event (occurred in August 2018). These biophysical variables have been found helpful in integrating the information of the interferometric coherence considering their direct correlation with the structural characteristics of the vegetation canopy (Ghosh et al., 2020). The final S2 dataset resulted was composed of all the 10 m resampled bands, the relative four vegetation indices for each month, and the three biophysical variables referred to only July 2018. The three biophysical variables calculated for July 2018 were clipped on the same area and stacked to the S1 dataset (Section 3.3.1) using the latter dataset as the master extent (S1+ S2).

Image classification
Supervised pixel-based image classification of forest cover was carried out using the RF algorithm (Breiman, 2001;Cutler et al., 2007). This, based on a set of decision trees, is a machine learning model widely used in land cover mapping, forest classification and the estimation of other forest structural parameters (Ghosh et al., 2020;Li et al., 2020;De Luca et al., 2019;Morin et al., 2019;Numbisi et al., 2019;Praticò et al., 2021). The RF algorithm was performed using the RFClassifier function from the Scikit-learn Python library package. In this study, the optimal RF hyperparameters were set using an exhaustive grid search approach implemented in Scikit-learn (GridSearchCV), based on a cross-validation analysis between all the possible combinations of a given set of hyperparameters values (Table 1) and for a given training input set.
The image classification has been implemented according to the following seven LULC classes: Eucalyptus, Euc; Pinus, Pin; Autochthonous Forest, AuFor; Soil; Pasture and/or Shrubs, Past/Shr; Urban, Urbe; Agriculture, Agri. The training pixels were selected on 950 regions of interest (ROIs) with a square size of 4 × 4 pixels. The ROIs, scattered over the entire study area and with a balanced distribution among the seven LULC classes (Euc,298;Pin,78;AuFor,131;Soil,95;Past/Shr,240;Urbe,52;Agri,56), were manually drawn by visual interpretation supported by the use of Google Earth very high-resolution satellite images (Google Earth Homepage, 2021).
Numerous scholars (Kattenborn et al., 2019;Nicolau et al., 2021;Zhang et al., 2019) stated that training data for land use classification could be based on high-resolution imagery instead of field observations, allowing to achieve comparable results a saving time and costs, and adequate representativeness of the full range of environmental characteristics present in those large portions of territory, which are more challenging to catch by punctual field observations.

Accuracy assessment
The accuracy assessment of the classified map was performed using a validation dataset formed by 658 ROIs of different random sizes (from 1 × 1 to 20 × 15 pixels) drawn with the same criteria used for training ROIs and balancedly distributed among the seven classes. The confusion matrix was implemented and, following the calculation of the producer's and user's accuracies (Congalton & Green, 2019), the single-class F-score (Fscore i ; eq. 6) and the multi-class F-score (F-score M ; eq. 7; Modica et al., 2021;Sokolova & Lapalme, 2009) were computed. The F-score i was calculated from the producer's and user's accuracies, representing the harmonic mean (De Luca et al., 2021b;Modica et al., 2021). .
where i represents the single class. The number of trees in the RF model max_depth 10,20,30,40,50,60,70,80,90,100,110,300,500,800,1000 The maximum depth of the tree min_samples_split 2, 5, 10 The minimum number of samples required to split an internal node min_samples_leaf 1, 2, 4 The minimum number of samples required to be at a leaf node max_features "auto", "None", "log2" The number of features to consider when looking for the best split The F-score M was calculated using the multi-class producer's M (eq. 8) and user's M (eq. 9) metrics, derived averaging the respective single-class producer's and user's accuracies for all the classes and expressed as: where n represents the total number of classes.

Dataset combinations results comparison
In addition to the integrated S1+ S2 dataset (main target of this study), the entire classification process and subsequent accuracy assessment was applied using each of the two single optical and SAR datasets, in order to compare the results and evaluate how effectively the integration of the two types of information can improve the outcomes. The optimal RF parameter values and the feature importance were also computed for each S1 and S2 datasets. Table 2 shows the hyperparameters set for the RF algorithm using the exhaustive grid search optimization approach. The parameter optimization was carried out for each of the three dataset combination: the integrated SAR and optical (S1+ S2), the SAR (S1) and the optical (S2). Among the tested ones, the values resulted as optimal were: 1200 (S1+ S2, S2) for the number of trees; 110 (S1+ S2, S2) and 800 (S1) the maximum depth of a tree; the default values 2, 1, and "auto" (S1+ S2, S1, S2) resulted for the min_samples_split, min_samples_leaf and the max_features, respectively. The value "auto" implies that the maximum number of features considered equals the square root of the total number of features. The initial out-of-bag (OOB) error, expressing a predicted accuracy performance estimated by the RF model during the training step, resulted in 98.13% (S1+ S2), 88.98% (S1) and 97.40% (S2) using the respective optimized hyperparameters set.

Results
Figure 3 (top) shows the land cover map resulted from the RF classification process applied to the S1 + S2 dataset and covering the entire study area. They are also magnified (bottom) three sample areas showing the details of the obtained classification. The total surface occupied by forest classes resulting from the classification equals 673.59 km 2 (excluding water surface that accounts for 9.81 km 2 ).
The resulted distribution among the forest classes showed in Figure 3 (pie chart) is 139.13 km 2 for Eucalyptus, 11.13 km 2 for Pinus, 29.71 km 2 for Autochthonous Forest. The classes corresponding to the unforested surfaces present an area equal to 344.77 for Pasture/Shrubs, 79.26 km 2 for Soil, 64.80 for Agriculture km 2 , 6.79 km 2 for Urban.
The accuracy level of the land cover map was assessed for each class deriving the producer's and user's metrics and their harmonic mean (F-score) from the confusion matrix showed in Figure 5.

Discussion
The main objective of this study was to obtain a map of the forest cover of an area around the municipality of Monchique, in southern Portugal, where a high severe fire occurred in August 2018 (De Luca et al., 2021a, 2021b. The vegetation cover mapping and, especially, the distinction between forest and non-forest vegetation before a disturbance event is a decisive purpose for improving the analysis of its effects on the ecosystem and their monitoring from the short to the long term (Chu & Guo, 2013). The integrated use of S1 and S2 time-series for forest tree species classification was evaluated to deal with this aim. Several SAR vegetation indices (RVI, RFDI) and optical vegetation indices (NDVI, NFRE, NBR) were computed for each month of the time series and included. The coherence information of the last month of the time series (July 2018) was calculated from the InSAR process and implemented as additional information for image classification. The classification approach was carried out by training the RF machine learning algorithm.

Accuracy and uncertainty
The good overall accuracy achieved by the integrated optical and SAR datasets in this approach, represented by an F-score M equal to 90.33%, is in line with the outcomes obtained from other studies where the combination of optical and SAR data was used (Spracklen & Spracklen, 2021;Zhang et al., 2019). This value should have been higher considering that most of the single classes exceeded the 90% threshold of F-score i , with the class relating to the pasture cover (Past/Shr) reaching the value of 95.67%. The lowest F-score i value was obtained by the AuFor (79.58%). It was caused by the wrong classification of pixels belonging to the AuFor class as Eucalyptus cover. These errors of commission of the AuFor class are represented by an user's accuracy equal to 70.32%, in contrast with the producer's accuracy value, equal to 91.65% in the same class. However, it is noticeable how the integration of the two sensors has improved the accuracy of this LULC class, considering that the Fscore i value for the optical single dataset was 71.95%. Probably, the origin of these errors relies on the combination of two main factors. First, the AuFor class comprises a mix of different Mediterranean broadleaved species, therefore not constituted by a specific and univocal spectral signature. Second, several isolated nuclei of Eucalyptus are scattered within the mixed Authocthonous forest. At the same time, within some gaps in the areas covered by Eucalyptus, there may be Table 2. The adopted Random forests (RF) parameters values for each dataset combination tested (integrated SAR and optical, S1+ S2; only SAR, S1; only optical, S2), set using the exhaustive grid search approach.

Parameter name
Values set S1+ S2 S1 S2  n_estimators  1200  1750  1200  max_depth  110  800  110  min_samples_split  2  2  2  min_samples_leaf  1  1  1  max_features "auto" (Square root of the total number of features) "auto" (Square root of the total number of features) "auto" (Square root of the total number of features) some small nucleus of species belonging to native vegetation, challenging to detect based on satellite imagery with 10 m of spatial resolution. It has to be considered that this spatial resolution of the Sentinel sensors could involve errors and uncertainties in the classification output to identify covers/crowns smaller than this dimension. Zhang et al. (2019) supposed that this could simultaneously conduct two opposite scenarios for each pixel: overestimation or underestimation/missing of the actual cover of smaller covers/trees. Furthermore, Zhang et al. (2019) pointed out some other elements that could cause uncertainty in its study, based on the classification of wooded areas using S1 and S2 data: the method of selection of training data, which may not be optimally representative of the actual conditions of the study area depending on their number and spatial distribution; a temporal miss-match between the interpreted Google Earth VHR data and Sentinel images; the presence of flooding, residual clouds, and shadows; other disturbances (e.g., fire, deforestation, etc.) that could affect the spectral signatures. Another source of errors could derive from edge pixels between different land classes (Lapini et al., 2020). However, this does not spoil the excellent effectiveness of these data, considering the small extent of these errors on the final accuracy values and, above all, the free availability of the images and software to process them. In fact, the precision obtained in the classification of areas not covered by forest is to be considered optimal, despite the complexity of the cover background that composes it (crops, orchards, towns, pastures, shrubs, grasses, soil, rock, etc.). On the other hand, this result was expected since the use of interferometric coherence is strongly correlated with the structural characteristics of forest arboreal vegetation (Ghosh et al., 2020;Siqueira, 2019). Aspect demonstrated by the relevance resulting from the feature importance analysis of this data.

Coherence
As the tree canopy height is a dendrometric measure describing the vegetation structure and it is a significant indicator of the aboveground biomass of forest areas (Ghosh et al., 2020), in this study, we empirically exploited the geometric relationship between the InSAR coherence and the forest stand height as additional information to discriminate the forest cover. Through the InSAR analysis, it is possible to use the phase information of the SAR signal to characterize the topography and height of the Earth's surface and Figure 4. The figure shows the feature importance (Gini importance) expressed by the first fifteen image layers with the highest importance (left column) and by the last fifteen image layers with the lowest importance (right column), calculated for each dataset combination tested (integrated SAR and optical, S1+ S2; only SAR, S1; only optical, S2). the objects present on it. Phase decorrelation, with consequent loss of coherence values, occurs where the surfaces' conditions change during the two moments of acquisition of the pair of interferometric images. The relationship between this decorrelation and the height of the forest surface is based on the empirical assumption that as the height of the trees increases, the volume and quantity of canopy increases, causing greater decorrelation. Therefore, a decrease in coherence could indicate an increase in canopy height (Ghosh et al., 2020;Perko et al., 2011;Siqueira, 2019). In this study, although the time gap was not large between the two acquisition dates of the SAR SLC images the average of the three coherence layers for each of the two flight paths was used to reduce the adverse effects of decorrelation. When observations are made at different times, the targets within a SAR resolution cell (e.g., small branches and leaves) may have moved, causing an error in measuring the trigonometric gaze angle and, therefore, a reduction in the interferometric coherence (Siqueira, 2019). The use of shorter wavelength SARs (X-band, C-band) involves a higher temporal decorrelation due to interaction with smaller objects, even for a time interval of just one day (Ghosh et al., 2020). Concerning the coherence information in this study, the two polarizations resulted equally influent in the dataset classification. About this aspect conflicting testimonies were found in literature: Ghosh et al. (2020) used only the co-polarized (VV) coherence information for forest stand height estimation due the higher noise produced by the volumetric effect of canopy scatter elements on cross-polarized SAR signal; on the other hand, Siqueira (2019) affirms that the cross-polarized coherence signature is in general more appropriate for characterizing forest structure as it is more correlated to the multiple volume scattering of vegetation canopy.

Optical and SAR layer integration; feature importance
In this study, temporal spectral signatures for each optical band were exploited as a valuable expression of the seasonal phenological and photosynthetic activity variations. In particular, the use of optical vegetation indices is essential in the phenological discrimination and observation of the different types of forest vegetation, such as Mediterranean evergreen conifers and deciduous forests (Aragones et al., 2019). However, the temporal spectral signature may also have limitations. Although their efficiency is proven, in Grabska et al. (2019), it is reported that the species-specific characteristics do not only and uniquely influence the spectral reflectance, but often there are intra-specific variations in reflectance caused by the age of the trees, by states of stress or disease, or by local site conditions such as level of coverage, the effects of the surface texture of the roof and the resulting shadow effects, the type of soil present which in conditions of sparse foliage contributes to the reflected signal. These factors can cause as much significant deviations of the typical spectral signature values as to create overlaps between species and make it difficult for classifiers to discriminate. Vegetation indices improve these aspects, better characterizing the species-specific temporal dynamics of the various forest cover types. Aragones et al. (2019) used the NDVI time-series to characterize the phenological changes of five Mediterranean Pinus species for species discrimination. They observed how the index is subject to a significant decrease during the summer drought. Zhang et al. (2019), focusing on a savanna ecosystem, encountered such high differences in NDVI time signatures between woody and nonwoody as to allow their distinction.
However, the feature importance analysis shows that not all the spectral bands/indices contributed equally to the forest cover classification. Whether the integrated dataset or the single dataset is used, the feature importance shows the SWIR-based index (NBR), the red-edge-based (NDRE) and the B8A, B8, B7 and B12 bands as the most significant optical layers for RF prediction. This confirmed the findings of other studies concerning the high efficiency of the SWIR and red-edge bands for vegetation mapping (Grabska et al., 2019;Immitzer et al., 2016). The red-edge bands are more sensitive to the photosynthetic pigments (chlorophyll a and b) levels and their variations from the biophysical and biochemical points of view. In contrast, the SWIR wavelengths are sensitive to the water content of the surface, since their optical absorbance increases with the increase of water content. Spracklen and Spracklen (2021) observed how the SWIR-based indices were the most important in distinguish plantation from natural forest. This is probably due to the different water content in plant tissues, with a more significant presence in plantations with a high growth rate. Observing the feature importance of respect to the imagery's seasonality, the autumn-winter images seem to have most influenced the algorithm performance, followed by the summer months. In these periods, the phenological difference between species is higher, contributing to the discrimination of broadleaved species (Grabska et al., 2019). The fCOVER and fAPAR S2 biophysical variables also appeared among the most influential image layers. It should be considered that Ghosh et al. (2020), combining the use of S1 coherence and S2 biophysical variables (LAI and fraction of vegetation cover, FVC) and modeling using the RF regressor, found an excellent correlation with the canopy height. They concluded that since coherence information could potentially be affected by any object on the Earth's surface greater than the SAR wavelength, the biophysical variables support this gap by indicating the presence and status of the vegetation cover, as well as by a proven direct correlation with the height of the canopy.
Cloud cover, which is the main obstacle in dealing with time-series of optical data, was effectively addressed in this study by applying the S2-Cloud Probability mask product available in GEE, based on an automated machine learning pixel-based cloud detection, followed by linear interpolation for filling the missing pixel values. Other errors could derive from discontinuities of time series due to cloud covers or other artifacts. In general, the gap-filling methods may impact the quality of the images at various levels. However, the consistency of these errors is relatively not relevant, and, generally, they do not weigh the statistical quality metrics, although they may be visually distinguishable (Inglada et al., 2017). If, on the one hand, the optical data require few and consolidated pre-processing steps, the SAR data are more complicated to manage during the processing and interpretation of the information due to various factors intrinsic to the characteristics of the signal/ sensor and their interaction with the characteristics of the affected surface (De Luca et al., 2021a;Tanase et al., 2020). In order to reduce speckle noise, both a monthly BTA and a speckle filter were applied. In fact, as Lapini et al. (2020) stated, although the multitemporal averaging inherently introduces a reduction in speckle-noise, a performance improvement was observed when a filter was applied, mainly if the time series consists of a relatively low number of images. Moreover, these authors pointed out how the execution of the backscatter temporal average of the whole SAR time series reduces seasonality effects (soil moistures, presence/absence of leaves, trees water content, etc.) for forest classification purposes, preserving radiometric discrimination between classes.
Recently several other authors (Lasko, 2019;Lehmann et al., 2015;De Luca et al., 2021b;Morin et al., 2019;Spracklen & Spracklen, 2021;Stroppiana et al., 2015;Zhang et al., 2019) have demonstrated that the combined use of optical and SAR data optimize both usage potential, filling the gaps of each other's and increasing the accuracy of the returned products. However, excluding the case of InSAR coherence already discussed, no SARderived images are present among the most significant, and a co-polarized BTA image achieved the lowest value of feature importance. This indicates that, in the present case, the information given by the time-series of the SAR data was not actually decisive in classifying the forest cover. This can be confirmed by comparing the accuracy values, where the integration of the two sensors resulted in an improved F-score M of only 2.53%, in line with what was stated by Chust et al. (2004) and Chatziantoniou et al. (2017). However, the accuracy improvements of more than 7% (F-score i ) found in some classes (e.g., Urbe), probable positive effect of the use of InSAR coherence on surfaces with lower decorrelation, should not be underestimated.
Concerning the SAR information, there are also much more pronounced limitations. Nicolau et al. (2021) demonstrated that classes with similar ground cover and similar backscatter have a lower separability potential than classes with distinct scattering mechanisms. Finally, he concluded that using dual-polarimetric SAR indices, relating the two different types of polarization (VV and VH), improved the spectral separability between classes. In the present study, we calculated the RVI and RFDI for each month, using the respective monthly BTA. The use of both polarizations, combined through dual-polarimetric indices, integrates their respective information (De Luca et al., 2021a;Nicolau et al., 2021). The crosspolarized backscatter is more sensitive to the distribution of volume scatters of canopy elements (leaves, branches, etc.) than the co-polarized signal. The latter is more associated with the underlayer soil backscatter (Meyer, 2019). Indeed, in the Mediterranean context, Lapini et al. (2020) better individuated the forest land cover when cross-polarization was used, even in a simple visual RGB SAR composite. Nevertheless, other authors such as Nicolau et al. (2021), using S1 time series for land cover classification, observed that the single polarization (VV or VH) did not achieve good results for multiple land cover classes classification. For this reason, besides the simple backscatter values (VV and VH), they used SAR indices such as the VV/VH ratio and the modified RFDI as additional spectral information to train the classifier.
As mentioned above concerning coherence information, the sensitivity of the SAR signal to the structure of forest vegetation also depends on the wavelength used. Shorter wavelengths, such as C-band (5 cm) and the Xband (3 cm), are more sensitive to the canopy surface characteristics because comparable with the dimensions of needles and leaves, and therefore more helpful in differentiating between coniferous and broad lives (Lapini et al., 2020). However, due to confusion with non-forest or secondary vegetation cover, misclassification is encountered (Lapini et al., 2020;Nicolau et al., 2021). Another factor to consider is the ease of saturation that the C-band has towards dense forest cover due to lower penetration capacities than longer wavelengths (Lapini et al., 2020;Meyer, 2019). Longer wavelengths (L-band and P-band) are proven to be better to estimate forest parameters (Morin et al., 2019), and therefore to better discriminate the forest vegetation (Lapini et al., 2020). Lapini et al. (2020), using only SAR data, proved that combining more wavelengths (multi-frequency approach), coming from different sources, leads to the highest accuracy of forest cover map.
It is expected that the imminent (in 2022) launch of the BIOMASS mission (Le Toan et al., 2011), consisting of a P-band polarimetric SAR satellite, will further improve forests biomass estimations. This could also concern other aspects such as forest type classification and other forest parameters estimation. The future objective will be to evaluate these new data as integrative information of the already consolidated optical spectral response.
Meanwhile, the approach proposed in this study demonstrated the effective potential of the combined use of S1 and S2 imagery in classifying forest cover in a fragmented and heterogeneous Mediterranean landscape. The applicability of remote sensing in these conditions has always been a complex task due to the anthropogenic influence on the landscape and the natural variability of plants structure and response to the climaterelated typical events (e.g., summer aridity, higher fire severity, etc.; Lapini et al., 2020). The situation is worsened by the topographic aspect of the study area. The roughness and irregular Earth's surfaces are more complicated to classify correctly than flat and homogeneous areas, mainly if SAR data are used, due to the effect of the slope exposure on both the sensor signal and on the characteristics of the vegetation cover (Inglada et al., 2017).
With this research, we demonstrated the reliability of our open workflow based on diverse free and opensource software to accurately map the forest cover, distinguishing the non-forest vegetation among the others LULC classes. Moreover, our workflow has been able to classify the vegetation cover more in-depth, discriminating between Eucalyptus, pine, and other forest types. This information could also be helpful in the framework of programmes aimed at long-term mapping and monitor of natural (or semi-natural) and planted forest cover. Future analysis dealing with sustainable forest management, habitat and biodiversity monitoring, carbon cycle estimation, and forest inventory could further exploit the reliability of our proposed research.

Conclusions
The forest and plant composition mapping is essential to analyze ecosystems' quantitative and qualitative characteristics to facilitate various monitoring applications of their state and condition. In this context, remote sensing techniques and tools demonstrated to be very efficient, especially with the advent of easily accessible and open-source solutions (data and software). This study explored the potential of the combined use of SAR Sentinel-1 SAR and optical Sentinel-2 band and indices time-series, integrated the InSAR coherence measure and the optical biophysics variables to classify forest cover and discriminate it from the surrounding non-forest land covers. A new cloud detection tool, provided by the Sentinel Hub team, was implemented to optimize the cloud masking, followed by a temporal linear interpolation for gaps filling. This allowed us to effectively manage one of the main problems encountered in the analysis of optical images. Among the layers that have had greater importance in obtaining good results are the NBR and NDRE, mainly from the autumn, and interferometric coherence. This study aimed to optimize a supervised classification procedure for the quantitative and qualitative analysis of the forest vegetation cover in the Mediterranean region in the time period immediately prior to a wildfire occurred in August 2018. Obtaining a map of the vegetation cover, with a good level of accuracy achieved (> 90%), is a significant advantage to improve future monitoring and analysis of the study area and to be able to carry out an effective and more targeted management on specific ecosystems. This aspect is even more interesting if we consider that the approach presented by this study has been implemented with the combination of free imagery and open-source software, proving the efficiency of the interoperability of the various web platforms and open-source libraries, from download to final process (GEE, ASF, ESA SNAP, Scikit-learn, etc.). The spatial resolution offered by these satellites, although allowing classification at a good scale of detail, still does not allow the detection of small landscape elements, such as narrow roads and small patches of land cover, and still creates confusion where the land cover is not pure (e.g., commission errors were found where small nuclei of Eucalyptus settled in the autochthonous forest vegetation). In the distinction and classification of vegetation cover, the contribution of Cband backscatter was found not to be as decisive as the InSAR coherence itself or some optical bands, as already observed in the literature. It is expected that the imminent launch of new SAR satellite platforms, such as the SAR L-band mission BIOMASS, will improve the contribution of this type of information. Future studies should test and validate the proposed approach on different Mediterranean study areas, equally contributing to the poor state of the art on the use of multisensor data LULC classification in this biome.The results obtained in this work will be fundamental to set up a more focused investigation on the integrated use of SAR and optical data for monitoring over time the study area so as to assess the ecosystems' response to wildfires.

Acknowledgments
Giandomenico De Luca was supported by the European Commission through the European Social Fund (ESF) and the Regione Calabria (POR Calabria FESR FSE 2014-2020). João Silva was funded by the Forest Research Centre, a research unit funded by Fundação para a Ciência e a Tecnologia I.P. (FCT), Portugal (UIDB/00239/2020). The authors would like to thank the two reviewers for their insightful comments on the earlier version of the manuscript, contributing to improving it significantly

Data availability statement
The data that support the findings of this study are available from the corresponding author, upon reasonable request. These data were derived from the following resources available in the public domain, ESA Sentinel Homepage, https:// sentinel.esa.int/web/sentinel/home.

Disclosure statement
No potential conflict of interest was reported by the author(s).