*Article* **Mapping Coastal Dune Landscape through Spectral Rao's Q Temporal Diversity**

**Flavio Marzialetti <sup>1</sup> , Mirko Di Febbraro <sup>1</sup> , Marco Malavasi 2,\* , Silvia Giulio <sup>3</sup> , Alicia Teresa Rosario Acosta <sup>3</sup> and Maria Laura Carranza <sup>1</sup>**


Received: 27 June 2020; Accepted: 17 July 2020; Published: 18 July 2020

**Abstract:** 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 classification approach specifically 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 classification 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., Modified Soil-adjusted Vegetation Index 2—MSAVI2, Normalized Difference Water Index 2—NDWI2 and Brightness Index 2—BI2). We identified seven land cover classes with high levels of accuracy, highlighting different covariates as the most important in differentiating them. The proposed framework proved effective 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.

**Keywords:** coastal habitats; ecosystem monitoring; land cover mapping; random forest algorithm; Sentinel-2; modified soil-adjusted vegetation index 2–MSAVI2; normalized difference water index 2–NDWI2; brightness index 2–BI2

#### **1. 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–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,

seasonality also plays a pivotal role in such ever-changing nature, affecting both abiotic (e.g., magnitude and intensity of weather and sea conditions) and biotic (e.g., phenology) factors, whose interaction finally gives rise to an extremely dynamic and complex mosaic of psammophilous plant communities, bare surfaces and water surfaces [9,10]. Moreover, being highly endangered, coastal dunes are of conservation concern worldwide and, as such, they need updated monitoring and mapping protocols [11,12]. Such protocols, traditionally based on field biodiversity data and habitat photointerpretation of aerial imagery [13,14], have improved in the last decades with the support of remotely sensed images and biomass spectral information (e.g., [15]). Nonetheless, in order to better discriminate between different types of standing biomass, remote sensing approaches for habitat mapping and monitoring in coastal areas are conventionally performed through classification of images captured at the peak of the plant growing season (see [16]). That said, developing an efficient mapping protocol that accounts for the dynamic nature of coastal systems, thus capturing all of the main components of coastal dune mosaics together with their seasonal variation (i.e., vegetation, bare surfaces and water surfaces) [17], still represents a challenge for delivering more accurate and viable maps.

At present, the availability of remotely sensed data, at various spatial, spectral and temporal resolutions, offers a great potential to carry out accurate and cost-effective studies accounting for ecosystem phenology and seasonality [18–20], thus reinforcing this research field, which traditionally relies on ground-based observations [21,22]. The analysis of remotely sensed phenological trends may effectively support land cover classification and mapping [23,24]. In addition, this potential is steadily growing with the continuous improvement of the temporal resolution of orbiting satellites [25,26].

The remotely sensed analysis of seasonality, traditionally based on spectral bands [27,28], has been improved over time by including the time-series analysis of spectral indices [20]. Several metrics have been proposed to quantify environmental seasonality and vegetation phenology based on remotely sensed data [29]. Such metrics commonly focus on describing the cyclic behavior of spectral indices, intended as the combination of spectral reflectance from two or more bands indicating the relative abundance of features of interest, and on identifying key transition dates (e.g., start, peak, end, and length of seasonal periods) [30]. Additionally, other indices, considering the whole array of available bands and used in the past for describing spectral diversity across space, have been recently extended to quantify variations in diversity across time [31,32]. Among them, the Rao's Q index, borrowed from community and landscape ecology [33–35], is able to take into account both the proportion of cells assuming different spectral values and their spectral distance [36]. Consequently, extending Rao's Q rationale to the temporal dimension (e.g., comparing spectral values and distances of a given location in a multi-temporal stack), makes it a good candidate to accurately synthesize ecosystems' seasonality. The calculation of Rao's Q on a temporal stack generates a new layer of temporal diversity, with high or low values indicating seasonal or stable habitats, respectively. Several diversity layers summarizing temporal stacks (bands or spectral indices) can be combined instead of the rough spectral values traditionally used on multi-temporal classification [16,36]. As far as we know, few studies, if any, extending Rao's Q to the spectral temporal diversity for land cover mapping, currently exist in transitional systems.

In this context, the present work sets out to provide and test a land cover classification approach specifically designed for coastal landscapes based on the within-year temporal diversity, computed through Rao's Q index. In particular, the within-year diversity will be calculated upon the seasonal behavior of three spectral indices properly describing the main components of the coastal mosaic (i.e., vegetation, bare surfaces and water surfaces). Such behavior will be then classified to provide a map of natural and seminatural land cover types. While accounting for the dynamic/temporal dimension highly characterizing coastal landscape, we expect that such an approach will provide a high level of classification accuracy.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study was carried out on the Adriatic coast of Central Italy (Molise region, Figure 1). This area, of approximately 9000 km<sup>2</sup> , is mainly composed of sandy beaches, a few river mouths and channels, and one rocky promontory. Dunes, occupying a narrow strip parallel to the seashore, are low and relatively recent (formed in the Holocene period) [37,38].

Along a sea-to-inland gradient, the typical vegetation zonation ranges from embryonic dunes in the seashore, followed by mobile dunes with perennial herbaceous vegetation, fixed dunes covered by evergreen shrub and small sclerophyllous trees and, in the inner sectors, by wooded dunes covered by coniferous forests (Figure 1) [39,40]. The Molise coast hosts several ecosystems of conservation concern in Europe (European Directive 92/43/EEC) [41]. For this reason, such an area is largely included in the European Natura2000 system [41] and is part of the European LTER monitoring network (Long Term Ecological Research network) [40,42].

. **Figure 1.** (**a**) In black, the study area including the coastal dunes of the Molise Region (Italy). Most of the analyzed coastal sectors are included in Sites of European Conservation Concern (SCI, European Directive 92/43/EEC): Foce Trigno-Marina di Petacciato (IT7228221); Foce Biferno-Litorale di Campomarino (IT7222216); Foce Saccione-Bonifica Ramitelli (IT7222217) and belong to the European LTER network [40,42]. (**b**) An example of coastal zonation. Reference system WGS84 UTM32 (epsg: 32632).

#### *2.2. 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: (1) Sentinel-2 imagery selection; (2) spectral indices calculation and temporal variability computation; (3) variables selection and classification; (4) accuracy assessment (Figure 2).

**Figure 2.** Workflow synthesizing the full mapping procedure of coastal dune Semi-natural and natural cover types with temporal MSAVI2, NDWI2, BI2 series and Random Forest classification approach.

#### 2.2.1. 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.

2.2.2. 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.

**Acronym Name Formula Index of References** MSAVI2 Modified Soil-Adjusted Vegetation Index 2 MSAVI2 = 2∗NIR+1− q (2∗NIR+1) 2 −8∗(NIR−RED) 2 photosynthetic biomass [48] NDWI2 Normalized Difference Water Index 2 NDWI2 <sup>=</sup> GREEN−NIR GREEN+NIR presence of water surfaces [51] BI2 Brightness Index 2 BI2 = q (RED∗RED)+(GREEN∗GREEN)+(NIR∗NIR) presence of bare surfaces [52]

3

**Table 1.** Spectral indices selected for analyzing the temporal diversity, serving as proxies of seasonality in vegetation biomass, presence of water surfaces, and bare surfaces.

For each spectral index, we built an annual stack containing the 12 month values (i.e., MSAVI22019, NDWI22019, BI22019), 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., MSAVI22019, NDWI22019, and BI22019) according to the following formula (Equation (1)) [57,58]:

$$Q\_{\text{index\\_}n} = \sum\_{i=1}^{\mathbb{C}-1} \sum\_{j=i+1}^{F} d\_{-}n\_{ij} \* p\_i \* p\_j \tag{1}$$

where:

*Qindex\_n* = Rao's *Q* quantifying the within-year variability of a spectral index (MSAVI22019, NDWI2<sup>2019</sup> 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].

**Figure 3.** Schematic representation of temporal Rao's Q diversity calculation implemented on year stacks of spectral indices MSAVI2, NDWI2 and BI2 to summarize the seasonal behavior of vegetation biomass, water and bare soil surfaces on coastal dunes.

In order to summarize the range of annual variation for each index, we also calculated their mean values (MBI2, MMSAVI2, MNDWI2), along with their 10th and 90th percentiles (10th BI2, 10th MSAVI2, 10th NDWI2, 90th BI2, 90th MSAVI2, 90th NDWI2, Figure 2 box 2 c).

#### 2.2.3. 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 classification, we used 12 variables, which describe the temporal heterogeneity (QMSAVI2, QNDWI2, QBI2), the central tendency (MMSAVI2, MNDWI2, MBI2), and the variation range (10th MSAVI2, 10th NDWI2, 10th BI2, 90th MSAVI2, 90th NDWI2, 90th BI2) of the three spectral indices (i.e., MSAVI2, NDWI2, BI2).

#### (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))

$$\text{BC}\_{\text{pq}} = \frac{\sum\_{i=1}^{n} |\mathbf{x}\_{\text{pi}} - \mathbf{x}\_{\text{qi}}|}{\sum\_{i=1}^{n} (\mathbf{x}\_{\text{pi}} - \mathbf{x}\_{\text{qi}})} \tag{2}$$

where BCpq is the Bray–Curtis dissimilarity between pixels (p, q), xpi are the values of the n selected remotely sensed variables (e.g., QMSAVI2, MNDWI2, 90th BI2) on pixel p and xqi the variables value on pixel q. BCpq 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–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–78], while testing different combinations of the number of variables randomly selected at each node (*Mtry* parameter) and split rules [79–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).

#### 2.2.4. 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, 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–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].

#### **3. Results**

#### *3.1. 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): Mobile Dune Herbaceous Vegetation (MDHV), Fixed Dune Herbaceous Vegetation with Sparse Shrub (FHVSS), Evergreen Woody Vegetation (EWV), and Deciduous and Humid Herbaceous Vegetation (DHHV, 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.

**Table 2.** Description of classification cycles (1◦ , 2◦ ,3◦ cycles) in terms of the optimal number of classes (Op. num. classes) according to Silhouette (S), Calinski–Harabasz (CH) and Davies–Bouldin (DB) indexes; the MESS values (%); the selected variables according to VIF selection, the variables' importance in the definition of classes and the obtained cover classes (Classes). QNDWI2: temporal Rao of NDWI2, QBI2: temporal Rao of BI2, QMSAVI2: temporal Rao of MSAVI2, 10th NDWI2: 10th percentile of NDWI2, 10th BI2: 10th percentile of BI2, 90th BI2: 90th percentile of BI2, 90th MSAVI2: 90th percentile of MSAVI2, MMSAVI2: mean of MSAVI2.


**Figure 4.** Scheme of the obtained natural and semi-natural cover classes (in the top) along with the respective boxplots showing variables' values. Reported variables are those selected by VIF analysis and used for classification on each cycle. Classes are: Water (W), Sand (S), Vegetation (V), Water Edge (WE), Open Sand (OS), Mobile Dune Herbaceous Vegetation (MDHV), Evergreen Woody Vegetation (EWV), Deciduous and Humid Herbaceous Vegetation (DHHV), Fixed Dune Herbaceous Vegetation with Sparse Shrub (FHVSS). Variables are QNDWI2: temporal Rao of NDWI2, QBI2: temporal Rao of BI2, QMSAVI2: temporal Rao of MSAVI2, 10th NDWI2: 10th percentile of NDWI2*,* 10th BI2: 10th percentile of BI2, 90th BI2: 90th percentile of BI2, 90th MSAVI2: 90th percentile of MSAVI2, MMSAVI2: mean of MSAVI2.

#### *3.2. 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). Indeed, the Water Edge class evidences an almost perfect agreement only for the Producer's accuracy, while the agreement of User's accuracy and Matthews

Correlation Coefficient shows a substantial agreement (Producer's accuracy: 96.67, User's accuracy: 72.22, Matthews Correlation Coefficient: 0.66, 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).


**Table 3.** The percentages of overall accuracy assessment values for three classification cycles, in particular Overall accuracy (O ACC), Kappa statistic (K), overall Matthews Correlation Coefficient (O MCC). In the Random Forest accuracy assessment are indicated the mean values and standard deviations for the overall performance metrics.

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).

#### *3.3. Variables Importance*

Among the most important variables in the first RF cycle, the 10th percentile of NDWI2 (10th NDWI2) showed a 63% importance, followed by the 10th percentile of BI2 (10th BI2) with 22%, the temporal Rao of NDWI2 (QNDWI2, 11%), and the temporal Rao of BI2 (QBI2,4%). Water class (W) is characterized by high values of 10th NDWI2 and QNDWI2 (Figure 4) and includes the sea, rivers, channels and wetland. In this class, the 10th BI2 shows overall low values and variability, while QBI2 and QMSAVI2 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%, QMSAVI2) followed by temporal Rao of BI2 (26%, QBI2) and the 90th percentile of BI2 (11%, 90th 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 (90th 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 (90th MSAVI2, Figure 4).

The most important variables in the third RF cycle are the lower values of bare surfaces (10th BI2, 33%) and the temporal variability of biomass (QMSAVI2, 32%) and of bare surfaces (QBI2, 26%, Table 2). Mobile Dune Herbaceous Vegetation (MDHV, Figures 4 and 5) is characterized by high values of bare surfaces (high 10th BI2 and low Q BI2) and moderate values of biomass (low MMSAVI2 and low QMSAVI2). The MSAVI2 index profile depicts a mosaic of bare sand with sparse open perennial and annual vegetation that corresponds to coastal herbaceous vegetation referable to mobile dunes (Figure 5). Fixed Dune Herbaceous Vegetation with Shrubs (FDHVS, Figures 4 and 5) is characterized by a strong seasonality (intermediate MMSAVI2 and very high QMSAVI2), intermediate biomass values and moderate presence of bare surfaces (moderate 10th BI2 and QBI2) and water (moderate QNDWI2). The FDHVS corresponds to herbaceous vegetation with shrubs growing on fixed dunes and dune slacks, more inland with respect to the previous class of ruderal vegetation. Evergreen Woody Vegetation (EWV, Figures 4 and 5) is characterized by high biomass values, a low seasonality (intermediate MMSAVI2 and very low QMSAVI2) and almost no bare surfaces (very low 10th BI2 and QBI2). EWV includes the Mediterranean maquis, pine woods and shrubs and woody evergreen formations. The Deciduous and Humid Herbaceous Vegetation class (DHHV, Figures 4 and 5) is characterized by high biomass values and pronounced phenology (very high MMSAVI2 and very high QMSAVI2), as well as moderate bare surfaces occurrence (moderate 10th BI2 and QBI2). The DHHV includes herbaceous and woody vegetation of humid areas close to river mouths, riparian vegetation and residual patches of lowland woods.

**Figure 5.** An example of the obtained land cover map reporting classes projected on Google Earth View (on the top) along with monthly average values of NDWI2, BI2, and MSAVI2 ± standard deviation. Water (W), Sand (S), Vegetation (V), Water Edge (WE), Open Sand (OS), Mobile Dune Herbaceous Vegetation (MDHV), Fixed Dune Herbaceous Vegetation with Shrub (FDHVS), Evergreen Woody Vegetation (EWV), Deciduous and Humid Herbaceous Vegetation (DHHV).

#### **4. 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 (10th NDWI2) and its temporal diversity (QNDWI2) 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 (QMSAVI2) 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 (QMSAVI2) along with the annual fluctuations on bare soil (QBI2) 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 (QMSAVI2) and the low percentages of bare soil (10thBI2) as the most important predictors for classification. In addition, low temporal variation on bare soil cover (low QBI2) 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 (QMSAVI2) and bare surfaces (QBI2), 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.

#### **5. 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: The percentages of accuracy assessment values by land cover classes for the third classification cycle: Dense Herbaceous Vegetation (DHV), Evergreen Woody Vegetation (EWV), Fixed Dune Herbaceous Vegetation with Shrubs (FDHVS), Deciduous and Humid Herbaceous Vegetation (DHHV). 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.

**Author Contributions:** All authors contributed substantially to the work: F.M., M.M. and M.L.C. conceived and designed the study; F.M. collected the data; F.M., M.D.F., M.L.C. analyzed the data; F.M., M.M., S.G., M.L.C. led the writing of the manuscript; A.T.R.A. and M.L.C. supervised the research. All authors contributed critically to the drafts and gave final approval for publication. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This study was carried out with the partial support of PON-AIM (Italian program for research and innovation 2014-2020—AIM1897595-2) in collaboration with the bilateral program Italy–Israel DERESEMII (Developing state-of-the-art remote sensing tools for monitoring the impact of invasive plant species in coastal ecosystems in Israel and Italy) and INTERREG Italy-Croazia CASCADE (CoAStal and marine waters integrated monitoring systems for ecosystems proteCtion AnD management - Project ID 10255941). The authors acknowledge the Principal Investigator(s) of the Sentinel-2 mission for providing datasets in the archive and the developers of SNAP software used for data analysis. Sentinel-2 images are freely downloadable from the Copernicus Open Access Hub https://scihub.copernicus.eu/). Furthermore, the authors acknowledge the entire team of EnvixLab, in particular Ludovico Frate, Gabriella Sferra, Marco Varricchione, Luca Francesco Russo, for the help given. The authors also thank the editor and the anonymous reviewers for their comments that helped us to improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
