Mapping Coastal Dune Landscape through Spectral Rao’s Q Temporal Diversity

: Coastal dunes are found at the boundary between continents and seas representing unique transitional mosaics hosting highly dynamic habitats undergoing substantial seasonal changes. Here, we implemented a land cover classiﬁcation approach speciﬁcally designed for coastal landscapes accounting for the within-year temporal variability of the main components of the coastal mosaic: vegetation, bare surfaces and water surfaces. Based on monthly Sentinel-2 satellite images of the year 2019, we used hierarchical clustering and a Random Forest model to produce an unsupervised land cover map of coastal dunes in a representative site of the Adriatic coast (central Italy). As classiﬁcation variables, we used the within-year diversity computed through Rao’s Q index, along with three spectral indices describing the main components of the coastal mosaic (i.e., Modiﬁed Soil-adjusted Vegetation Index 2—MSAVI2, Normalized Di ﬀ erence Water Index 2—NDWI2 and Brightness Index 2—BI2). We identiﬁed seven land cover classes with high levels of accuracy, highlighting di ﬀ erent covariates as the most important in di ﬀ erentiating them. The proposed framework proved e ﬀ ective in mapping a highly seasonal and heterogeneous landscape such as that of coastal dunes, highlighting Rao’s Q index as a sound base for natural cover monitoring and mapping. The applicability of the proposed framework on updated satellite images emphasizes the procedure as a reliable and replicable tool for coastal ecosystems monitoring. proxies of seasonality in vegetation biomass, presence of water surfaces, and bare surfaces.


Introduction
The increasing impact of human activities and the derived environmental transformations (i.e., climate and land-cover change, invasive species and habitat loss) are promoting changes in global biodiversity at an unprecedented rate in human history [1][2][3]. The effects of such changes are particularly severe on coastal dune landscapes [4,5], despite the fact that they host a highly specialized biodiversity [6] and provide essential benefits to society [7,8].
Coastal dunes are found at the boundary between continents and seas, representing unique transitional mosaics hosting highly dynamic habitats undergoing frequent and substantial changes in physical extent and environmental conditions. Along with this typical transitional condition,  [40,42]. (b) An example of coastal zonation. Reference system WGS84 UTM32 (epsg: 32632).

Methodology
We followed an unsupervised approach and classified coastal dune natural and semi-natural land cover types through a hierarchical cluster analysis. The classes identified in the clustering phase were then used as the response variable in a Random Forest (RF) model (an accurate learning method for discriminating differences among classes [43,44]), in order to quantify their accuracy and to predict their occurrence in the study area. In particular, the procedure was organized in the following  [40,42]. (b) An example of coastal zonation. Reference system WGS84 UTM32 (epsg: 32632).

Methodology
We followed an unsupervised approach and classified coastal dune natural and semi-natural land cover types through a hierarchical cluster analysis. The classes identified in the clustering phase were then used as the response variable in a Random Forest (RF) model (an accurate learning method for discriminating differences among classes [43,44]), in order to quantify their accuracy and to predict their occurrence in the study area. In particular, the procedure was organized in the following steps:

Sentinel-2 Imagery Selection
We used Sentinel-2 mission images as a sound support for land monitoring [45], with good spatial (10, 20 or 60 m), temporal (revisit time of 2-3 days at mid-latitudes) and spectral (13 bands ranging from 400 nm to visible to 2400 nm) resolutions [25]. Multispectral images of Sentinel-2 satellites were downloaded from Copernicus Open Access Hub (https://scihub.copernicus.eu/) at the level of bottom of atmosphere (Level-2A), and used to build a monthly temporal dataset for the year 2019. As the study area falls in two tiles (T33TVG, T33TWG), we selected 24 images (12 for tile) of

Sentinel-2 Imagery Selection
We used Sentinel-2 mission images as a sound support for land monitoring [45], with good spatial (10, 20 or 60 m), temporal (revisit time of 2-3 days at mid-latitudes) and spectral (13 bands ranging from 400 nm to visible to 2400 nm) resolutions [25]. Multispectral images of Sentinel-2 satellites were downloaded from Copernicus Open Access Hub (https://scihub.copernicus.eu/) at the level of bottom of atmosphere (Level-2A), and used to build a monthly temporal dataset for the year 2019. As the study area falls in two tiles (T33TVG, T33TWG), we selected 24 images (12 for tile) of each month with low cloud coverage (<5%, Supplementary Material). All images were atmospherically corrected by ESA through the "sen2cor" processor algorithm ( Figure 2, box 1) [46]. We specifically used Green band, with wavelength 560 ± 36 nm; Red band, with wavelength 665 ± 31 nm, and NIR band, with wavelength 833 ± 106 nm, with 10 m of resolution.

Spectral Indices Calculation and Temporal Variability Computation
For each image, we calculated three spectral indices: Modified Soil-Adjusted Vegetation Index 2 (MSAVI2), Normalized Difference Water Index 2 (NDWI2) and Brightness Index 2 (BI2), as they are good indicators of vegetation biomass, the presence of water surfaces, and bare surfaces, respectively (Figure 2 box 2 a, Table 1).
MSAVI2 is a vegetation index that quantifies the photosynthetic biomass based on Red and NIR bands (Table 1). MSAVI2 has been developed for mapping landscapes characterized by high percentages of bare surfaces [47,48] and its spatial behavior is a good support for land classification [49]. MSAVI2 ranges from −1 (absence of biomass vegetation) to 1 (maximum of biomass vegetation) with higher values indicating higher percentages of photosynthetic biomass [48].
NDWI2 is a water index that is useful to identify water surfaces, exploiting the Green and NIR bands (Table 1). NDWI2 is a remotely sensed index that is particularly efficient for identifying water surfaces and for mapping water-land transitions [50,51]. NIR band and Green band present opposite reflectance values behaviors and NDWI2 values range from −1 to 1, where values greater than 0 indicate water surfaces [52].
The BI2 index is sensitive to soil brightness, hence it quantifies the bare surfaces through the square root of brightness of each pixel (Table 1) [53]. BI2 accurately discriminates the bare surfaces from vegetation in heterogeneous environments [54,55]. The minimum value of BI2 is 0 and indicates the absence of bare surfaces, as growing positive values correspond to increasing percentages of bare surfaces.
We used ESA's Sentinel-2 toolbox-ESA Sentinel Application Platform 7.0 (SNAP) for index calculation. presence of bare surfaces [52] For each spectral index, we built an annual stack containing the 12 month values (i.e., MSAVI2 2019 , NDWI2 2019 , BI2 2019 ), standardized between 0 and 1 to make them comparable [56]. To summarize the within-year heterogeneity of ecological conditions (e.g., biomass, water surfaces and bare surfaces yearly variation), we calculated the temporal Rao's Q index for each annual stack ( Figure 2 box 2 b). The Rao's Q index has recently been borrowed from functional ecology and successfully applied in remote sensing contexts as an innovative measure of spectral heterogeneity [36]. Such a proposed version of temporal Rao's Q index accounts for both the relative abundances of the values assumed by a given pixel n throughout the temporal stack (e.g., 12 months), and the Euclidean distances among the pixel's numerical values.
The Rao's Q diversity index was applied on each annual stack (i.e., MSAVI2 2019 , NDWI2 2019 , and BI2 2019 ) according to the following formula (Equation (1)) [57,58]: where: Remote Sens. 2020, 12, 2315 6 of 18 Q index_n = Rao's Q quantifying the within-year variability of a spectral index (MSAVI2 2019 , NDWI2 2019 or BI2 2019 ) for the pixel n p = relative abundance of the index value assumed by the pixel n within a temporal stack i = month i j = month j d_n ij = distance between the month i-th and j-th index value of pixel n (d ij = d ji and d ii = 0). As similarly done in spectral diversity applications, Rao's Q adapted to detect temporal diversity quantifies the expected dissimilarity between two combinations of pixel values randomly selected within a pixel temporal stack ( Figure 3). Pixels representing highly seasonal cover types (e.g., deciduous or annual formations) are characterized by a pronounced within-year variability of biomass values and bare surfaces cover. Therefore, such pixels should assume high temporal Rao's Q values. On the contrary, pixels representing temporally stable coastal areas (e.g., open sand, pine wood) portraying weak or absent seasonal variations, should score low Rao's Q values. The temporal Rao's Q for each spectral index (Q MSAVI2 , Q NDWI2 , Q BI2 ) was computed through R package 'spacetimerao' 0.1 [59].
BI22019) according to the following formula (Equation (1)) [57,58]: where: Qindex_n = Rao's Q quantifying the within-year variability of a spectral index (MSAVI22019, NDWI22019 or BI22019) for the pixel n p = relative abundance of the index value assumed by the pixel n within a temporal stack i = month i j = month j d_nij = distance between the month i-th and j-th index value of pixel n (dij = dji and dii = 0). As similarly done in spectral diversity applications, Rao's Q adapted to detect temporal diversity quantifies the expected dissimilarity between two combinations of pixel values randomly selected within a pixel temporal stack ( Figure 3). Pixels representing highly seasonal cover types (e.g., deciduous or annual formations) are characterized by a pronounced within-year variability of biomass values and bare surfaces cover. Therefore, such pixels should assume high temporal Rao's Q values. On the contrary, pixels representing temporally stable coastal areas (e.g., open sand, pine wood) portraying weak or absent seasonal variations, should score low Rao's Q values. The temporal Rao's Q for each spectral index (QMSAVI2, QNDWI2, QBI2) was computed through R package 'spacetimerao' 0.1 [59]. In order to summarize the range of annual variation for each index, we also calculated their mean values (MBI2, MMSAVI2, MNDWI2), along with their 10 th and 90 th percentiles (10 th BI2, 10 th MSAVI2, 10 th NDWI2, 90 th BI2, 90 th MSAVI2, 90 th NDWI2, Figure 2 box 2 c).

Variables selection and Classification
The classification phase was implemented through consecutive cycles, each structured into three steps (Figure 2 box 3). Each step consists of a sequence of procedures: a) variables selection, b) pixel sampling and representativeness assessment, c) clustering and classification. As covariates for In order to summarize the range of annual variation for each index, we also calculated their mean values (M BI2 , M MSAVI2 , M NDWI2 ), along with their 10 th and 90 th percentiles (10 th BI2 , 10 th MSAVI2 , 10 th NDWI2 , 90 th BI2 , 90 th MSAVI2 , 90 th NDWI2 , Figure 2 box 2 c).

(a) Variables Selection
During each cycle, the 12 covariates were checked for their multicollinearity (Graham 2003) using the Variance Inflation Factor (VIF, Figure 2 box 3 a). VIF quantifies the multiple correlation of a variable with respect to all the other variables through linear regressions [60]. As VIF values greater than 5 indicate multicollinearity problems [61], we retained, for the analysis, the variables with VIF < 5 [62].

b) Pixel Sampling and Representativeness Assessment
To reduce the computational costs, we ran the classification process on a subset of pixels randomly extracted from the study area. In particular, classification was carried out on a 20%-pixel sample and the results were extrapolated throughout the entire study area using an RF algorithm. The degree of extrapolation on values of covariates lying outside the RF calibration range was assessed through the Multivariate Environmental Similarity Surface, (MESS, Figure 2 box 3 b) [63]. MESS quantifies the similarity of all the pixels in the study area with respect to the covariates of the extracted 20% sample used for classification. Low MESS values indicate a high similarity between covariate values of calibration and prediction pixels, suggesting a high representativeness of the former. In each classification cycle, we indicated the percentage of highly dissimilar pixels (MESS value < 0).

(c) Clustering and Classification
The random sample of pixels was classified through a hierarchical cluster analysis using the Ward's minimum variance method [64], which minimizes the cluster's internal variability using the sum-of-squares [65,66]. We used as distance the Bray-Curtis dissimilarity metric [67,68] measured as follows (Equation (2)) where BC pq is the Bray-Curtis dissimilarity between pixels (p, q), x pi are the values of the n selected remotely sensed variables (e.g., Q MSAVI2 , M NDWI2 , 90 th BI2 ) on pixel p and x qi the variables value on pixel q. BC pq ranges from 0 when the pixels are identical, to 1 when the pixels are completely different [69].
On each cycle, we identified the optimal number of classes (form from 2 to 10 classes; Supplementary Material) by calculating three indices (Silhouette index [70]; Calinski-Harabasz index [71] and Davies-Bouldin index [72]) on the just built hierarchical cluster. The selected indices summarize two cluster characteristics: the compactness of classes (e.g., how closely pixels are grouped inside a class), and the separation between classes (e.g., how the classes are different from each other).
The classes identified in the clustering phase were used as response variables in a RF model (R package 'caret' 6.0-85; Supplementary Material) [73][74][75] in order to evaluate their accuracy, to extrapolate the classification throughout all the pixels of the study area, and to evaluate the variables' contribution in defining such classes. To optimize RF parameters, we set a high number of uncorrelated decision trees (Ntree = 1000) [76][77][78], while testing different combinations of the number of variables randomly selected at each node (Mtry parameter) and split rules [79][80][81][82], then choosing the combination that yielded the highest Kappa statistic value. Specifically, we tested Mtry values ranging from 2 to the total number of variables in the cycle, and the Gini index and Extra-Trees algorithm as possible split rules (Supplementary Material).
Once the optimal RF model was identified, it was used to predict the membership of all the study area pixels to one of the classes identified in the clustering phase. Moreover, we estimated the relative variables' importance (Supplementary Material) [83].
At the end of each cycle, we repeated the three steps (variables selection, pixel sampling and representativeness assessment, clustering and classification) inside the classified land cover classes.
The final classes predicted by the classification phase were then interpreted in terms of coastal cover types through an expert-based approach based on field detection and visual interpretation of high resolution aerial images (~1 m).

Accuracy Assessment
The RF predictive accuracy was assessed by an internal 10-fold cross validation and an independent validation based on 300 random checkpoints. We selected 300 points as to assure a minimum standard number of 20 control points for each land cover class [84]. In the independent accuracy assessment, Remote Sens. 2020, 12, 2315 8 of 18 we compared the assignment of the checkpoints according to the land cover classes obtained by satellite classification with their description obtained from field observations and the visual inspection of Google Earth images. Because of the differences in spatial resolution between Sentinel-2 (10 × 10 m) and Google Earth images (~1 m), we focused our visual inspection on circular buffer areas of 5-m radius (10 m' diameter) around each check point [85][86][87][88]. In both accuracy assessments, we calculated the same performance metrics. We built a confusion matrix and calculated the percentages of overall accuracy, producer's accuracy, user's accuracy and Kappa statistic [84]. Moreover, we calculated the Matthews Correlation Coefficient, both overall and for each land cover class [89]. This coefficient is widely used to evaluate the accuracy of classifications [86,90], as it proved reliable to evaluate the accuracy of classification when the classes had different sizes; its range of −1 to 1 and values close to 1 represent a perfect accuracy [87,91].

Unsupervised Classification
We obtained seven land cover classes after three classification cycles, which are organized in two hierarchical levels, with marked differences in temporal diversity values and in the range of the spectral indices (Figures 4 and 5). During the first classification cycle, we identified the first hierarchical level including three classes (see Table 2, Figure S1): Water (W), Sand (S), and Vegetation (V, Figure 4, Table 2). The second and third classification cycles established the second hierarchical level in which sand and vegetation clusters were split into two and four classes, respectively (see Table 2). Particularly, the second RF cycle applied to Sand class (S) identified two categories (Table 2, Figure S2): Water Edge (WE) and Open Sand (OS; Figures 4 and 5). The WE class represents the sea-land transition including the inter-tidal area, while OS includes dry sand dunes partially covered by sparse annual vegetation. In the third cycle, RF split the class Vegetation (V) into four categories (Table 2, Figure S3 Figures 4 and 5). Since we focused on terrestrial cover types, we did not explore the presence of sub-classes inside the Water class (Table 3). For each cycle, RF showed extremely low MESS values, suggesting no extrapolation effect on models' predictions throughout the study area. These outcomes evidenced that the selected pixels for each classification cycle are representative of the study area.

Accuracy Assessment
In all the three cycles, the RF results achieved very high accuracy values according to both crossvalidation and field assessments ( Table 3).
For the first cycle, the overall accuracy, Kappa statistic and MCC indicate an almost perfect agreement under both assessments (Table 3). Similar results were also obtained by performance metrics by land cover classes. In the field accuracy assessment, both user's and producer's accuracy values in all classes are high, with the Sand class accuracy resulting lower than 90% (Table S3). The second classification cycle reported cross-validation performances of 99% for overall accuracies, 96% for Kappa statistic, and 0.98 for Matthews Correlation Coefficient, while the performance assessed through the field assessment is slightly lower (overall accuracy: 87%, Kappa statistic: 72%, Matthews Correlation Coefficient: 0.66, Table 3). All performance metrics by land cover classes calculated through cross-validation indicate an almost perfect agreement, whereas the comparison of the natural and semi-natural land cover map between the field and visual inspection shows a decrement

Accuracy Assessment
In all the three cycles, the RF results achieved very high accuracy values according to both cross-validation and field assessments ( Table 3).
For the first cycle, the overall accuracy, Kappa statistic and MCC indicate an almost perfect agreement under both assessments (Table 3). Similar results were also obtained by performance metrics by land cover classes. In the field accuracy assessment, both user's and producer's accuracy values in all classes are high, with the Sand class accuracy resulting lower than 90% (Table S3). The second classification cycle reported cross-validation performances of 99% for overall accuracies, 96% for Kappa statistic, and 0.98 for Matthews Correlation Coefficient, while the performance assessed through the field assessment is slightly lower (overall accuracy: 87%, Kappa statistic: 72%, Matthews Correlation Coefficient: 0.66, Table 3). All performance metrics by land cover classes calculated through cross-validation indicate an almost perfect agreement, whereas the comparison of the natural and semi-natural land cover map between the field and visual inspection shows a decrement of accuracy (Table S4) Table S4). The Open Sand class displays a higher agreement than the previous class, and, only for the Matthews Correlation Coefficient, the agreement drops to under 0.70 (Table S4). Finally, the performance of the third classification cycle has very high results under cross-validation assessment, with all metrics above 95% (overall accuracy: 97%, Kappa statistic: 96%, Matthews Correlation Coefficient: 0.96), while in the field accuracy assessment, the overall metrics show a substantial agreement (Table 3). Similar to the previous cycles, all land cover classes were accurately predicted, as showed by cross-validation results (Table S5). In the field accuracy assessment, the user's accuracy of all classes shows substantially good performances, with the metrics values resulting higher than 75%, especially in the Mobile Dune Herbaceous Vegetation class (Table S5). Equally, in the producer's accuracy, all classes show a substantial agreement, and the values are higher than 70%; indeed, the values range from the minimum of substantial agreement (73%) in Deciduous and Humid Herbaceous Vegetation to the maximum of almost perfect agreement (93%) in Fixed Dune Herbaceous Vegetation with Shrubs (Table S5).

Variables Importance
Among the most important variables in the first RF cycle, the 10th percentile of NDWI2 (10 th NDWI2 ) showed a 63% importance, followed by the 10th percentile of BI2 (10 th BI2 ) with 22%, the temporal Rao of NDWI2 (Q NDWI2 , 11%), and the temporal Rao of BI2 (Q BI2 ,4%). Water class (W) is characterized by high values of 10 th NDWI2 and Q NDWI2 (Figure 4) and includes the sea, rivers, channels and wetland. In this class, the 10 th BI2 shows overall low values and variability, while Q BI2 and Q MSAVI2 are close to zero. These features are also evident when inspecting the annual trend of each of the three spectral indices (MSAVI2, NDWI2 and BI2; Figure 4). For the Water class, NDWI2 values are above 0.5 in all months. The BI2 and the NDWI2 values are always lower than 0.25 and 0, respectively.
The most important variables in the second RF cycle are temporal Rao of MSAVI2 (63%, Q MSAVI2 ) followed by temporal Rao of BI2 (26%, Q BI2 ) and the 90th percentile of BI2 (11%, 90 th BI2 , Table 2). Water Edge and Open Sand showed similar trends in annual spectral indices (Figures 4 and 5), only differing in the upper values of brightness (90 th BI2 ), which is higher in Open Sand. Open Sand shows higher BI2 values in summer (Figure 4), and lower variability in the upper values of biomass (90 th MSAVI2 , Figure 4). The most important variables in the third RF cycle are the lower values of bare surfaces (10 th BI2 , 33%) and the temporal variability of biomass (Q MSAVI2 , 32%) and of bare surfaces (Q BI2 , 26%, Table 2 Figure 4). The most important variables in the third RF cycle are the lower values of bare surfaces (10 th BI2, 33%) and the temporal variability of biomass (QMSAVI2, 32%) and of bare surfaces (QBI2, 26%, Table 2

Discussion
In this study, we used information from annual fluctuations of vegetation biomass, water surface, and bare soil combined into the temporal Rao's Q index, in an effort to identify and map seven land cover classes that clearly depicted the sequence of coastal dune natural and semi-natural ecosystems occurring along the sea-inland gradient.
Our modelling approach showed high levels of classification accuracy, substantially confirming the usefulness of time-series analysis for semi-natural and natural land cover mapping [16,92,93]. Moreover, we pointed out the high potential of integrating such spectral indices as MSAVI2, NDWI2 and BI2 to describe land cover seasonality in such heterogeneous environments as those of coastal dunes.
The classification framework that made it possible to identify and map seven coastal semi-natural cover types was organized in two hierarchical levels. For each set of cover types, different variables emerged as the most important for classification. In the first RF cycle implemented on the overall coastal landscape, water surface was split from the other cover classes because of the typical behavior of water-related spectral indices. Specifically, the minimum of NDWI2 (10 th NDWI2 ) and its temporal diversity (Q NDWI2 ) assumed particularly high values and emerged as the most important variables in the first RF cycle. As previously observed, the spectral responses of water and dryland are very different [94,95]. Such differences clearly emerged at coarser classification levels on the transition area we analyzed here. However, the percentage of bare surfaces also appeared as an important variable in the first cycle, revealing the presence of open sand and sparse vegetation in the coastal dune mosaic. In the second and third cycles performed on non-water areas, the temporal diversity of biomass (Q MSAVI2 ) showed an important role in land cover classification, confirming the importance of phenology in coastal landscape mapping [16]. In particular, the second cycle, which focused on open sand and sparse vegetation areas, reported the phenology (Q MSAVI2 ) along with the annual fluctuations on bare soil (Q BI2 ) as important variables discriminating bare areas from complex mosaics of pioneer vegetation and sand [96,97]. The third RF cycle focused mainly on vegetated areas hosting herbaceous and woody land cover classes, reporting the temporal variability of vegetation biomass (Q MSAVI2 ) and the low percentages of bare soil (10th BI2 ) as the most important predictors for classification. In addition, low temporal variation on bare soil cover (low Q BI2 ) helped to discriminate seasonal deciduous formations from the contiguous evergreen vegetation [38,96].
The classification of Rao's Q temporal diversity measured on spectral indices made it possible to identify and map the mosaic characterizing the different sectors of coastal dune zonation with high accuracy levels. However, such transition classes as Water Edge and Fixed Dune Herbaceous Vegetation with Shrubs evidenced the lowest values, perhaps due to the spatial dimension of Sentinel-2 images (10 m), which is too coarse to describe this intricate transition mosaic.
The seven identified classes clearly depicted the sequence of coastal dune natural and semi natural cover classes occurring from the sea to the inland, each characterized by a specific seasonal pattern. Indeed, classes ranged from the Water class including the sea close to the seashore, to Deciduous formations and Humid Herbaceous Vegetation growing in the inner dune slacks [16,41,92]. As for the Water Edge class, it represented the inter-tidal area, with higher presence of water during the winter months possibly caused by storms. The Open Sand class included the dry sand beach with the sparse annual vegetation on the drift line, while the embryonic shifting dune formations, such as the Dense Herbaceous Vegetation type, represented the mobile dunes with Ammophila arenaria (see also [16]). The Fixed Dune Herbaceous Vegetation with Shrubs class included a fine mosaic of herbaceous and shrub vegetation occurring on transition areas, as well as fuzzy edges between herbaceous and woody vegetation on dune sectors towards the sea and shrubs, and ruderal formations on inland disturbed areas [16,96]. As for the Evergreen Woody Vegetation class, we reported low temporal diversity values, both in vegetation biomass (Q MSAVI2 ) and bare surfaces (Q BI2 ), confirming the already known absence of seasonality of this cover class [16,98]. The Deciduous and Humid Herbaceous Vegetation class enclosed riparian deciduous woody vegetation and vegetation of humid environments. This class was characterized by a similar seasonality of vegetation biomass and presented two peaks in spring and summer that alternate with lower values in autumn and winter [99,100].
The results we obtained in this study clearly identified remotely sensed seasonality and Rao's Q temporal diversity as an effective source of information for landscape mapping in coastal areas. Using Sentinel-2 images with accurate spatial, temporal and spectral resolutions [25], we derived sound temporal heterogeneity variables for land cover mapping summarizing spectral indices previously implemented separately (e.g., using MSAVI2 [48], or NDWI2 [52], or BI2 [53]).
Land cover classification based on multi-temporal remotely sensed data analysis has improved here for landscape mapping on Mediterranean coasts by relying on a set spectral indices rarely analyzed together, as well as by the introduction of Rao's Q for summarizing their temporal variability. Furthermore, the implemented RF classification organized in sequential cycles ensured accurate and ecologically coherent results. For instance, water presence and seasonality evidenced by the yearly behavior of NDWI2 were important variables in the first cycle of classification implemented on overall coastal area extent, which confirmed its potential for discriminating water from emerged areas [94,95] and suggested its role when mapping transition systems between terrestrial and marine realms. Similarly, our results confirmed that biomass and bare soil indexes and their temporal variability are important variables for mapping terrestrial cover classes [16,28,92,101]. Based on our results, it possible to affirm that for remotely sensed monitoring and mapping, the inclusion of variables depicting bare soil, water and biomass seasonality are highly advisable.

Conclusions
The proposed procedure for classifying and mapping natural and semi-natural land cover effectively summarized the dynamic nature of coastal systems, capturing all of the main components of coastal dune mosaics together with their seasonal variation (i.e., vegetation, bare surfaces and water surfaces) by combining spectral indices that are rarely used together (MSAVI2, NDWI2, or BI2).
The information provided by Rao's Q offered a sound base for natural and semi-natural land cover monitoring and mapping. In particular, the here proposed implementation of Rao's Q temporal heterogeneity allowed the accurate depiction of the seasonality of the different cover types conforming the coastal dune mosaic along the sea-inland gradient (e.g., biomass phenology, water seasonality, yearly variation on bare surfaces). Furthermore, the applicability of the proposed framework on the available updated sentinel images emphasized the procedure as a promising tool for cover monitoring and reporting. Even more interestingly, the integration of other remotely sensed data with higher spatial resolutions derived by satellite, such as Planet images, UAV, or LiDAR data, may further improve the classification of coastal zonation.
From an applied perspective, the natural and semi-natural land cover map provided in this study yields relevant knowledge for coastal monitoring and management; therefore, we hope new studies exploring increasingly larger areas will be analyzed to further test the proposed classification and, at the same time, to provide homogeneous information for coasts in the Mediterranean.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/14/2315/s1, Table S1: Sentinel-2 dataset indicating for each image the month, day, platform (Sentinel-2A or Sentinel-2B) and the hour of acquisition, the cloud percentage, and the relative tile (T33TVG for Molise north, T33TWF for Molise south). Table S2: The optimal number of classes for each cycle was identified by Silhouette, Calinski-Harabasz, and Davies-Bouldin indices calculated on a specific hierarchical cluster produced through Ward's minimum variance. These three indices identify the optimal number of classes by comparing the intra-class compactness (degree of aggregation between pixels inside each class) and the separation between classes (degree of differences between classes). We fixed as the optimal number of classes, the most frequent output produced by the three indexes. Figure S1: Graphical representation of the Silhouette, Calinski-Harabasz, and Davies-Bouldin indices used for identifying the optimal number of classes in the 1 • cycle of Random Forest. All three of the indices identified three as the optimal number of classes. Figure S2: Graphical representation of Silhouette, Calinski-Harabasz, and Davies-Bouldin indices used for identifying the optimal number of classes in the 2 • cycle of Random Forest. Two (Silhouette and Calinski-Harabasz) of the three indices have established two as the optimal number of classes. Figure S3: Graphical representation of Silhouette, Calinski-Harabasz, and Davies-Bouldin indices used for identifying the optimal number of classes in the 3 • cycle of Random Forest. Two (Silhouette and Calinski-Harabasz) of the three indices suggested four as the optimal number of classes. Table S3: The percentages of accuracy assessment values by land cover classes for first classification cycle: Water (W), Sand (S), Vegetation (V). In particular, User's accuracy (U ACC), Producer's accuracy (P ACC), and Matthews Correlation Coefficient (MCC). In the Random Forest accuracy assessment are indicated the mean values and standard deviation for the performance metrics. Table S4: The percentages of accuracy assessment values by land cover classes for second classification cycle: Water Edge (WE), Open Sand (OS). In particular, User's accuracy (U ACC), Producer's accuracy (P ACC), and Matthews Correlation Coefficient (MCC). In the Random Forest accuracy assessment are indicated the mean values and standard deviation for the performance metrics. Table S5