Capturing Coastal Dune Natural Vegetation Types Using a Phenology-Based Mapping Approach: The Potential of Sentinel-2

: Coastal areas harbor the most threatened ecosystems on Earth, and cost-e ﬀ ective ways to monitor and protect them are urgently needed, but they represent a challenge for habitat mapping and multi-temporal observations. The availability of open access, remotely sensed data with increasing spatial and spectral resolution is promising in this context. Thus, in a sector of the Mediterranean coast (Lazio region, Italy), we tested the strength of a phenology-based vegetation mapping approach and statistically compared results with previous studies, making use of open source products across all the processing chain. We identiﬁed ﬁve accurate land cover classes in three hierarchical levels, with good values of agreement with previous studies for the ﬁrst and the second hierarchical level. The implemented procedure resulted as being e ﬀ ective for mapping a highly fragmented coastal dune system. This is encouraging to take advantage of the earth observation through remote sensing technology in an open source perspective, even at the ﬁne scale of highly fragmented sand dunes landscapes. and Random Forest classification approach.


Introduction
Environmental monitoring is essential to identify and understand the structure, integrity and conservation status of different habitats forming landscape mosaics [1]. Next to traditional field-based techniques [2], Remote Sensing (RS) methods are useful tools for ecosystems monitoring, as they are able to capture a wide range of properties of vegetation in a standard and replicable way [3].
During the last few decades, satellite images have supported vegetation mapping and monitoring of wide landscapes [4], with a continuous improvement in spatial resolution and use of multi/hyperspectral sensors consistently boosting the performance of remotely sensed data for mapping highly fragmented areas [5][6][7]. In particular, for the interpretation of particularly complex or very fine-grained vegetation mosaics such as those commonly encountered on sand dunes, high-resolution data are an essential requirement [8,9]. Indeed, complex landscapes have long represented a challenge for vegetation mapping and multi-temporal monitoring applications, which still need further development [10,11]. In this context, several space agencies (e.g., ESA, NASA/USGS, CBERS, ISRO) deliver free remotely sensed products with several resolutions that represent a reliable support for different applications of ecosystem monitoring and management. Among these, the free access Sentinel-2 mission of ESA (European Space Agency), with a 10 m spatial resolution for bands in the spectral range of the blue (B2-490 nm), green (B3-560 nm), red (B4-665 nm) and infra-red (B8-842 nm) and a revisit rate of approximately five days, represents an important support for environmental mapping in the Mediterranean area [12,13]. The relatively recent release of Sentinel-2 mission (Sentinel-2A launched on 23 June 2016 and Sentinel-2B launched on 7 March 2017) still has several potentialities to be explored [14].
In order to achieve their full potential for accurate mapping of complex vegetation mosaics, RS techniques must be coupled with proper approaches able to capture spatial and temporal vegetation patterns [1]. For instance, the analysis of vegetation phenological properties, describing recurring biological events (e.g., seasonality) [15], is very effective for mapping intense seasonal biomass variations [16] as those characterizing Mediterranean landscapes. Among the remotely sensed data, vegetation indices, depicting ecosystem spectral properties, are very efficient tools to map vegetation and its temporal pattern [17,18]. Among these, Normalized Difference Vegetation Index (NDVI) [19,20] is a good proxy of canopy biomass [21], and its application for environmental monitoring is highly appreciated [22,23]. Furthermore, the monthly variation of NDVI values across an entire year proved to be a sound surrogate of ecosystem phenology [18,24,25], which allows for discriminating contiguous vegetation cover types featuring different seasonality [26][27][28].
Coastal dune landscapes are complex mosaics that develop in the transition zones between terrestrial and marine environments, occupying strips parallel to the seashore [29]. Along the sea-inland gradient, coastal dunes are ruled by a large variety of constraining environmental conditions, such as soil salinity, substrate instability, wind and marine aerosol [30][31][32]. By shaping the biomass levels, this gradient determines the occurrence of a mosaic of highly specialized and diverse plant communities coexisting in a relatively narrow area which represents a hotspot of exclusive biodiversity [32]. In spite of their high biodiversity value and complex ecosystem functioning, coastal dunes are among the most threatened ecosystems worldwide [33,34]. In the Mediterranean areas, the loss and degradation of coastal dune ecosystems have been particularly severe in the last few decades [34][35][36], with the main threats being urban expansion [10,37,38], coastal erosion [39] and invasion by alien species [40][41][42][43]. In order to prevent these and other endangered habitats from further degradation, all European Member States adopted the Council Directive 92/43/EEC (hereafter Habitats Directive, HD). By signing the HD, the States committed themselves to maintain, restore and monitor habitats and species of European conservation concern (listed in dedicated annexes) and to report their conservation status every six years. Each habitat type (mainly identified by plant communities) is characterized by specific biotic and abiotic factors [44]. In this light, innovative and scientifically sound instruments are needed for setting conservation priorities and providing management indications for the coastal dune habitat types [45][46][47][48]. Until now, coastal dunes mapping procedures have mostly been based on the integration of visual interpretation (e.g., photointerpretation of aerial imagery) and floristic data [10,49]. However, the use of these mapping procedures presents some shortcomings in linear and fragmented ecosystems characterized by low biomass such as coastal dunes. For instance, photo-interpretation is a time-consuming procedure and its results vary depending on the subjectivity of the interpreter, his experience and personal knowledge of mapped landscapes [50,51]. Moreover, as small patches (below the minimum mapping area) are neglected, the photo-interpreted maps can be limited for characterizing the coastal fine scale mosaics and related landscape processes [52,53].
In consideration of the above, the present work sets out to explore the potential of Sentinel-2 in capturing coastal dune natural vegetation types using a phenology-based mapping approach. In particular, by a multi-temporal analysis of NDVI images of coastal dunes in central Italy, we focused on two main questions: (i) does NDVI phenological profiles allow for identifying and correctly mapping different vegetation types distributed along a coastal zonation? (ii) does the product of a phenology-based classification agree with existing coastal dunes classification systems (i.e., photo interpreted land cover maps and 92/43/EEC habitat distribution)?

Study Area
The study was carried out on a representative tract of the Mediterranean coast, placed in the Tyrrhenian seashore of Central Italy (Lazio Region). The study area includes ca. 250 km of sandy coast, mainly formed by recent (Holocenic) dunes ( Figure 1) [45]. These dunes are relatively simple in structure, with a single low dune ridge (lower than 10 m height) occupying a narrow strip (usually no more than 500 m with) along the seashore [32]. In this area, coastal plant communities mostly range from pioneer vegetation near the shoreline to Mediterranean shrubs on landward fixed dunes [33,45]. Previous studies discriminate eight different habitats occurring along this coastal zonation, of which two have been listed as priority (Table 1) [38,53]. In the analyzed littoral zone, human activities in the analyzed littoral zone have intensified in the course of the 20th century [10,54]. Specifically, during the last 60 years, the Tyrrhenian coast faced consistent processes of fragmentation [55], simplification [45] and biodiversity loss [56]. Nevertheless, the Tyrrhenian coast still hosts a good number of plant communities of conservation concern in Europe (92/43/EEC Habitats Directive; EEC, 1992) for which monitoring and conservation strategies must be improved.

Methodology
We performed a phenology-based classification of coastal dune ecosystem following a sequence of steps ( Figure 2): (1) multitemporal dataset collection and image preprocessing, (2) NDVI calculation and data processing, (3) data classification, (4) accuracy assessment, and (5) comparison of phenology-based classes vs. previous studies.

Methodology
We performed a phenology-based classification of coastal dune ecosystem following a sequence of steps ( Figure 2): (1) multitemporal dataset collection and image preprocessing, (2) NDVI calculation and data processing, (3) data classification, (4) accuracy assessment, and (5) comparison of phenology-based classes vs. previous studies.

Sentinel-2 Imagery and Multitemporal Dataset
Sentinel-2 is a European wide-swath, multi-spectral imaging mission, constituted by a twosatellite platform: Sentinel-2A and Sentinel-2B [57]. The Multi Spectral Instrument (MSI) on-board Sentinel-2 can provide images with a temporal resolution of five days at the equator, and a 12-bit radiometric resolution from 492 nm to 1377 nm, which includes the Visible (VIS), Near Infra-Red (NIR), and Short Wave Infra-Red (SWIR) spectra (13 bands). Sentinel-2 images are freely downloadable from the Copernicus Open Access Hub [58]. In this study, we used the red band (R, around 665 nm in the VIS spectrum) and the NIR band (around 833 nm), at 10 meters of resolution [12]. The study area is included in two Sentinel-2 dataset tiles (T33TTG, T33TUF).
We extracted monthly NDVI images recorded in the year 2017 for the two Sentinel-2 tiles describing the analyzed coast and we built a multi-temporal dataset (stack) ( Table S1). For the analysis, we selected only those images with low cloud coverage (< 15%), excluding March, due to

Sentinel-2 Imagery and Multitemporal Dataset
Sentinel-2 is a European wide-swath, multi-spectral imaging mission, constituted by a two-satellite platform: Sentinel-2A and Sentinel-2B [57]. The Multi Spectral Instrument (MSI) on-board Sentinel-2 can provide images with a temporal resolution of five days at the equator, and a 12-bit radiometric resolution from 492 nm to 1377 nm, which includes the Visible (VIS), Near Infra-Red (NIR), and Short Wave Infra-Red (SWIR) spectra (13 bands). Sentinel-2 images are freely downloadable from the Copernicus Open Access Hub [58]. In this study, we used the red band (R, around 665 nm in the VIS spectrum) and the NIR band (around 833 nm), at 10 meters of resolution [12]. The study area is included in two Sentinel-2 dataset tiles (T33TTG, T33TUF).
We extracted monthly NDVI images recorded in the year 2017 for the two Sentinel-2 tiles describing the analyzed coast and we built a multi-temporal dataset (stack) (Table S1). For the analysis, we selected only those images with low cloud coverage (< 15%), excluding March, due to the excess in cloud coverage. The clean stack for the phenology-based classification was composed of 22 images (one per tile and month, excluding March). Only part of the downloaded data was already corrected from atmospheric noise, while the other part was rough. We corrected these others using Sen2Cor version 2.5.5 [59,60].

NDVI Calculation and Masking
For the entire stack, we calculated NDVI (Equation (1)) as follows: NDVI ranges from -1 to 1, with increasing values related with growing photosynthetic biomass [24]. The NDVI stack was masked using a fine scale land cover map (1:5000; SITR-Sistema Informativo Territoriale Regionale Lazio) in order to exclude from the classification, at least partially, urban areas, agriculture fields and water bodies.

Data Classification
We classified the resulting NDVI stack by implementing a Random Forest algorithm (RF) with a hierarchical logic, using ESA's Sentinel-2 toolbox-ESA Sentinel Application Platform 6.0 (SNAP). RF is a machine learning classification method that operates by constructing a multitude of decision trees [61,62]. RF algorithm is widely used for classification because of its speed, stability and ability to discriminate differences [63][64][65].
We cyclically performed different RF analyses, defining for each cycle two parameters: the training set of pixels (Mtry, see the next paragraph) and the number of trees (Ntree, in our case 100 runs per cycle). In order to maximize the efficacy of RF, in the training set, we used a number of pixels that was higher than the square root of the total of the pixels [61,62].
In each cycle, the identification of the training set was supported by a Principal Component Analysis (PCA) of the monthly NDVI values (pixels x monthly NDVI values matrix). We projected all the pixels in the PCA1 and PCA2 ordination space [66]. As objects that are close in the ordination space (similar component values) describe similar cover types [67], we therefore selected as the training set for classification the two furthest away clouds of pixels with maximum variance between them [68]. For each tree (run), the training set was split through a bootstrapping procedure in two groups: seed pixels (inbag), to build the classification tree, and validation pixels (out-of-bag), to estimate the classification performance. Then, all pixels were compared with the inbag set using the Gini inequality index [69] and assigned to a class based on the Lorenz Curve. The latter ranges from 0 (perfectly equality) to 1 (perfectly inequality) [69]. At the end of each run, RF conferred to each pixel an ordinal vote (the minimum value of Lorenz Curve). After the entire cycle of 100 runs, each pixel was definitively assigned to the more frequently attributed class [61]. In each cycle, the performance of classification was assessed by repeating the classification procedure on the basis of the out-of-bag data and comparing both classifications [61,70].
In each RF cycle, we classified the stack in two vegetation classes, and we carried on all the cycles where the out-of-bag error was <50%.

Accuracy Assessment
The accuracy of the phenology classification map was assessed through an error matrix (Table S2) calculating overall accuracy (Equation S1), producer's accuracy (Equation S2) and user's accuracy (Equation S3) and Kappa statistic (0 ≤ Kappa ≤ 100; Equation S4) [71]. We based this assessment on 250 random checkpoints visually inspected on Google Earth images [72][73][74]. The positional errors of objects in Google Earth images are much lower compared to of the minimum spatial resolution of Sentinel-2 images [73] and their spatial resolution (~1m) is high enough to allow clear visual interpretation of the land cover [75]. In our case, the random checkpoints were inspected in Google Earth images for a 5 m radius buffer area (comparable scale of the 10 m pixel of Sentinel-2 images) to limit the scale mismatch errors. We built the error matrix, reporting the assigned phenology-based class in rows, and the Google Earth visual attribution in columns.

Phenology-Based Map vs. Previous Vegetation Studies
The error matrix was further used to compare the phenology-based map with two previous studies conducted in the same area. The first study reports an actual vegetation map produced at 1:5000 scale by visual interpretation of aerial orthophotos [76] and the second one refers to the natural dune vegetation types classified according to the Habitats directive (92/43/EEC; Table 1) and identified by floristic field data ( Figure S1). To compare the phenology-based classes with the photointerpreted map, we used 250 random checkpoints, while the congruence with habitat types was tested using 135 floristic plots extracted from "RanVegDunes", a database of randomly distributed floristic surveys [77]. For testing the congruence of the phenology-based classification with the existing documents, we aggregated and homogenized the classes of the vegetation map (Table S4) and the habitat types (Table S5) on the basis of similarities in physiognomies and ecological conditions [78]. Then, we tested the correspondence among maps and defined their respective levels of agreement through the Kappa statistic (Table S3).

Sentinel-2 NDVI Classification
The phenology-based classification allowed for identifying five vegetation classes organized in three hierarchical levels, each one characterized by a specific phenological pattern ( The Dense Herbaceous Vegetation and Ruderals class (Figure 3b2) is characterized by NDVI values slightly over 0.5 from November to April that decrease during summer. Dense Herbaceous Vegetation and Ruderals include a highly seasonal herbaceous vegetation. This class occurs close to the seashore in sectors exposed to environmental stress and on inner dune sectors characterized by high anthropic pressure.
The Sparse Woody class is composed by sparse evergreen vegetation (Figure 3c1) with high monthly NDVI values (> 0.5) that slightly decrease in summer. This class tends to occur at intermediate distances from the seashore and in the inner sectors of the dune in contact with densely wooded dunes.
Lastly, the Dense Woody Vegetation class presents high monthly NDVI values (> 0.7; Figure 3c2) throughout the year. It occurs in the back-dune zones corresponding to dense shrublands and forests.

Classification Accuracy Assessment
The classification of multi-temporal NDVI data (phenology-based classification) showed a good level of accuracy when compared with the visual interpretation of Google Earth images. The first level of classification (two classes: Vegetation and Open Sand) exhibited very high values of overall accuracy (96%) with both producer's and user's accuracy over 88% and Kappa statistic of 86% (Figure 4a).
Similarly, the overall accuracy of the second hierarchical level (Herbaceous and Woody Vegetation) was high (88%) with a moderate agreement among classes and reference data given by a Kappa statistic of 79% (Figure 4b).
All land cover classes identified at the third level of detail evidenced values of accuracy with producer's accuracy ranging between 69% (Herbaceous Vegetation) and 99% (Woody Vegetation). Similarly, the user's accuracy ranged between 86% (Woody Vegetation) and 90% (Open Sand and Herbaceous Vegetation) (Figure 4c).
The correspondence between the third hierarchical level of classification and the set of Google images resulted as moderate as underlined by both overall accuracy (79%) and Kappa statistic (71%) (Figure 4c). The producer's accuracy was adequate for all classes. In particular, Open Sand and Dense Woody Vegetation had the highest values of producer's accuracy (88% and 96%, respectively). This result defined an elevated precision of the map in these two classes. Moreover, both herbaceous classes showed moderate agreement (75% Sparse Herbaceous Vegetation, 65% Dense Herbaceous Vegetation). Finally, Sparse Woody Vegetation class featured the lowest value of producer's accuracy detected, even though the value showed moderate agreement (56%). The user's accuracy showed the higher values in Open Sand and Dense Herbaceous Vegetation classes (90% and 98%, respectively). Sparse Herbaceous Vegetation and Sparse Woody Vegetation had the lowest values (respectively 63% and 64%). Finally, Dense Woody Vegetation presented user's accuracy (75%).

Classification Accuracy Assessment
The classification of multi-temporal NDVI data (phenology-based classification) showed a good level of accuracy when compared with the visual interpretation of Google Earth images. The first level of classification (two classes: Vegetation and Open Sand) exhibited very high values of overall accuracy (96%) with both producer's and user's accuracy over 88% and Kappa statistic of 86% ( Figure  4a).
Similarly, the overall accuracy of the second hierarchical level (Herbaceous and Woody Vegetation) was high (88%) with a moderate agreement among classes and reference data given by a Kappa statistic of 79% (Figure 4b).
All land cover classes identified at the third level of detail evidenced values of accuracy with producer's accuracy ranging between 69% (Herbaceous Vegetation) and 99% (Woody Vegetation). Similarly, the user's accuracy ranged between 86% (Woody Vegetation) and 90% (Open Sand and Herbaceous Vegetation) (Figure 4c).
The correspondence between the third hierarchical level of classification and the set of Google images resulted as moderate as underlined by both overall accuracy (79%) and Kappa statistic (71%) (Figure 4c). The producer's accuracy was adequate for all classes. In particular, Open Sand and Dense Woody Vegetation had the highest values of producer's accuracy (88% and 96%, respectively). This result defined an elevated precision of the map in these two classes. Moreover, both herbaceous classes showed moderate agreement (75% Sparse Herbaceous Vegetation, 65% Dense Herbaceous Vegetation). Finally, Sparse Woody Vegetation class featured the lowest value of producer's accuracy detected, even though the value showed moderate agreement (56%). The user's accuracy showed the higher values in Open Sand and Dense Herbaceous Vegetation classes (90% and 98%, respectively). Sparse Herbaceous Vegetation and Sparse Woody Vegetation had the lowest values (respectively 63% and 64%). Finally, Dense Woody Vegetation presented user's accuracy (75%).

Harmonization and Agreement Test with Existing Documents
The NDVI classification showed significant values of agreement with the photointerpreted map [76], and floristic field data with the habitat types [53] at both the first and the second hierarchical levels ( Figure 5).
The first hierarchical level showed strong agreement values with the photointerpreted classification map, with 95% of overall accuracy and 83% of Kappa statistic (Table S6), and also the producer's and user's accuracy showed high agreement values (~100%). The agreement of the second classification level resulted in being quite significant for all of the classes, with high overall accuracy (80%) and Kappa statistic (66 %) depicting a moderate congruence between them. The user's accuracy for the three classes was over 77%. The producer's accuracy ranges between 58 % and 95% (Table S7). Finally, at the third level of classification, the agreement test resulted in being moderate, with overall accuracy and Kappa statistic values being approximately 66% and 55%, respectively. However, values of user's accuracy were high only for Open Sand (94%) and Dense Woody Vegetation (85%), while the other classes showed values under 50%. The producer's accuracy ranged between 26% and 81% (Table S8). The NDVI classification showed significant values of agreement with the photointerpreted map [76], and floristic field data with the habitat types [53] at both the first and the second hierarchical levels ( Figure 5).
The first hierarchical level showed strong agreement values with the photointerpreted classification map, with 95% of overall accuracy and 83% of Kappa statistic (Table S6), and also the producer's and user's accuracy showed high agreement values (~100%). The agreement of the second classification level resulted in being quite significant for all of the classes, with high overall accuracy (80%) and Kappa statistic (66 %) depicting a moderate congruence between them. The user's accuracy for the three classes was over 77%. The producer's accuracy ranges between 58 % and 95% (Table S7). Finally, at the third level of classification, the agreement test resulted in being moderate, with overall accuracy and Kappa statistic values being approximately 66% and 55%, respectively. However, values of user's accuracy were high only for Open Sand (94%) and Dense Woody Vegetation (85%), while the other classes showed values under 50%. The producer's accuracy ranged between 26% and 81% (Table S8).  see table 1) and vegetation classes mapped in a previous study [76].
On the other hand, the agreement values of NDVI phenology-based classification and habitat types assigned by floristic data (Habitats directive 92/43/EEC) showed great differences among the hierarchical levels. At the first two levels, agreement values indicated moderate congruence between the classification systems. At the first hierarchical level, the overall accuracy was ~79% and the Kappa statistic denoted moderate agreement value. User's accuracy was relatively high for both classes, Open Sand (81%) and Vegetation (78%). Producer's accuracy ranged between 53% to 93% (Table S9).
The agreement test of the second hierarchical level indicated similar consistency, with overall accuracy ~71% and moderate value for the Kappa statistic (58%). The user's accuracy was relatively high for Open Sand and Woody Vegetation (81% and 90% respectively), and the Herbaceous Vegetation manifested moderate values (62%). Similarly, producer's accuracy ranged between 53% and 84% (Table S10). Finally, the third level of classification exhibited the lowest agreement values.  Table 1) and vegetation classes mapped in a previous study [76].
On the other hand, the agreement values of NDVI phenology-based classification and habitat types assigned by floristic data (Habitats directive 92/43/EEC) showed great differences among the hierarchical levels. At the first two levels, agreement values indicated moderate congruence between the classification systems. At the first hierarchical level, the overall accuracy was~79% and the Kappa statistic denoted moderate agreement value. User's accuracy was relatively high for both classes, Open Sand (81%) and Vegetation (78%). Producer's accuracy ranged between 53% to 93% (Table S9).
The agreement test of the second hierarchical level indicated similar consistency, with overall accuracy~71% and moderate value for the Kappa statistic (58%). The user's accuracy was relatively high for Open Sand and Woody Vegetation (81% and 90% respectively), and the Herbaceous Vegetation manifested moderate values (62%). Similarly, producer's accuracy ranged between 53% and 84% (Table S10). Finally, the third level of classification exhibited the lowest agreement values. The overall accuracy and Kappa statistic were approximately 53% and 43%, respectively. Moreover, user's accuracy was relatively high only for Open Sand (81%), and moderate agreement for Dense Herbaceous Vegetation and Ruderals (61%) and Dense Woody Vegetation (67%), with the other classes featuring values under 50%. Producer's accuracy ranged between 11% and 73% (Table S11).

Discussion
Our results suggest the analysis of remote sensed data (Sentinel-2 images) by a phenology-based classification as an effective approach for monitoring natural landscapes. Sentinel-2 images confirmed their high potential for vegetation mapping [12], while the multitemporal analysis of NDVI provided complementary and useful information, proving its convenience even in complex vegetation mosaics, that is to say, beyond their traditional field of application [18,26,28].
The phenology-based classification using a Random Forest algorithm on a Mediterranean wide coastal area allowed for identifying and mapping five vegetation classes organized in three hierarchical levels. Such classes, each one characterized by specific phenological and ecological features, exhibited high levels of accuracy and clearly depicted the coastal ecosystem zonation ranging from Open Sand, occurring near the seashore line, to Dense Woody Vegetation on the inner dunal sectors [45,79,80]. Sparse Herbaceous Vegetation and Sparse Woody Vegetation occurred discontinuously, while Open Sand and Dense Herbaceous Vegetation and Ruderals formed a continuous strip close to the seashore running along all the analyzed coast. Finally, Dense Woody Vegetation formed regular shaped patches, and, as previously observed [32,38], occurred in the back-dune zone.
The classification of multi temporal NDVI images, which was successfully used for land cover mapping [26,80], is extended here for vegetation mapping on Mediterranean coasts. The phenological analysis that allows for depicting vegetation seasonality [26,81,82] enabled to discriminate woody evergreen from herbaceous annual vegetation [26]. Furthermore, by exploring phenological spatial variations occurring in correspondence with biomass transitions, it was possible to distinguish between densely and sparsely vegetation formations and to identify edges [83,84] between vegetation classes occurring in the analyzed complex mosaic.
Our results respond to the scope of Sentinel mission [12,57] and give new evidence of its potential for monitoring and mapping coastal dunes with reduced costs and time efforts. The high temporal resolution that assures a continuous release of new clean and fine resolution images (~each 10 days) postulates Sentinel-2 as one of the most effective supports for phenology-based coastal dunes monitoring [18].
The good agreement between phenology-based classification and the photo-interpreted vegetation map [76] suggests the adequacy of Sentinel 2 multi-temporal NDVI classification for producing new vegetation maps of the coastal dune landscapes. Indeed, the phenology based map is quite consistent with the photo-interpreted vegetation map (scale 1:5000) produced using panchromatic digital aerial ortho-photographs with about one meter of resolution of the year 2008. Furthermore, the hierarchical nature of RF classification offers a good basis for comparing the new remotely sensed classification with existing documents produced with different methodological procedures, data sources and spatial resolution. Linking the new remotely sensed classes with previous maps and mapping supports (as aerial photos) is essential for building a long-term ecological series and for monitoring coastal dune landscapes across time [78,85].
The agreement test of NDVI classes and floristic data referable to EC habitats (92/43/EEC) was significant at the first and second hierarchical level for all classes. At the third level, only Open Sand, Dense Herbaceous Vegetation and Ruderals, and Dense Woody Vegetation showed a relative high congruence with vegetation plots, consequently showing the possibility to discriminate the habitats included in these classes. The lower agreement of NDVI classification and vegetation plots classified in habitat types is probably related with differences in the spatial resolution between the remotely sensed instrument (10x10 m) and the field floristic plots (2x2 m). Furthermore, the Mediterranean coastal dunes, naturally conformed by fine scale ecosystem mosaics [33,45], are often disturbed by fragmentation processes that alter their biodiversity [56]. This probably reduces the possibility of a match between conventional small floristic plots and the 10 m resolution of Sentinel-2 images. This scale mismatch is the principal restriction in this study, and the impossibility to discriminate the single habitat type through this phenology-based approach using 10 m resolution images limits its use to coarse vegetation classes (general physiognomies) [86,87]. In coastal dunes, both subtle variations of environmental factors and human pressure promote the formation of fine scale mosaics of habitats, some of them featuring similar physiognomies but differing in their floristic composition [88][89][90][91][92][93].
To deal with such shortcomings, the integration of sentinel data with finer resolution data and tools are advisable. For instance, the use of multispectral satellite images with higher spatial resolution, or the implementation of other classification methodologies as spectral unmixing algorithms able to quantify the percentage of different cover classes inside the single pixels [94,95] should improve the performance of the proposed classification procedure. In any case, the classification performance of each cycle, elaborated by out-of-bag pixels, estimated the uncertainty of the Random Forest result giving an idea of the presence of mixed pixels.
Overall, the potential of Sentinel-2 data in a phenology-based mapping was accomplished with a relative high degree of accuracy assessment and significant congruence with existing previous classifications. It is very promising in the discrimination of annual, deciduous and evergreen vegetation. Moreover, the integration of remotely sensed maps with field data could contribute to continuous update of coastal dune habitats maps, reducing costs and risks of delaying the periodical reporting requested by the Habitat Directive (92/43/EEC).

Conclusions
The performed phenology-based classification emphasizes the potential of Sentinel-2 images for mapping natural vegetation and extends its field of application to low biomass and highly fragmented systems as coastal dunes. The combined use of NDVI multi-temporal data, machine learning (Random Forest) algorithms, and a pixel-oriented approach allowed for adequately describing with high values of accuracy the complex mosaic of coastal dune vegetation.
The phenology-based classification approach with Sentinel-2 data proposed here is a time saving and more objective approach, complemented with open source earth observation data and implemented through free ESA software, effective and inexpensive instruments for coastal monitoring. Furthermore, the good levels of agreement of phenology-based with previous vegetation maps should allow for building long-term ecological series necessary for exploring and monitoring coastal ecosystems dynamics over time.
Nevertheless, there are several possibilities to improve this phenology-based classification and enhance its potential. For instance, the integration of phenological classification with LiDAR and other remotely sensed data or the implementation of spectral unmixing algorithms could improve the agreement with floristic filed data and should represent new research frontiers to explore.
From an applied perspective, the phenology-based vegetation classification provides 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/11/12/1506/s1, Table S1: The total Sentinel-2 imagery dataset. It shows for each Sentinel-2 image the month, day, platform (Sentinel-2A or Sentinel-2B), the hour of acquisition, the cloud percentage, the processing level (top of the atmosphere-1C or bottom of the atmosphere-2A), and the tiles (T33TTG for Lazio north, T33TUF for Lazio south). Table S2: Example of error matrix. It is a contingency table (k x k array, where k is the number of classes in the classification). Equation S1: Overall accuracy, defined as the total of the correctly classified checkpoints on the total number of the checkpoints where nii indicates the number of checkpoints classified in the same category both in the satellite mapped classes and the Google Earth reference data; in other words, the elements of the major diagonal. Equation S2: Producer's accuracy: fraction of correctly classified checkpoints in all checkpoints of the produced classification. Equation S3: User's accuracy: fraction of correctly classified checkpoints in all checkpoints of the reference data. Equation S4: Kappa statistic (K) computed as follows where n_ij is the number of observation in row i and column j, n_(i+) and n_(+j) are respectively the total number of observation of rows, and the second total number of observation of columns. Table S3: Classes of Kappa statistic interpretation. The Kappa statistic is a measure of the difference between the actual agreement of real objects observable on Google Maps with resulted classes, and an agreement due to chance (where real objects are compared with a random classification). Kappa varies between 0 and 100, where values close to 0 represent a poor agreement, and values close to 100 are indicated as excellent level of agreement. Figure S1: A subset of the photointerpreted vegetation map produced at 1:5000 scale by visual interpretation of aerial ortophotos, and a subset of the floristic field data classified according with the Habitats directive (92/43/EEC). Table S4: Nomenclature homogenization between the produced phenology-based map and the vegetation map. Table S5: Nomenclature homogenization of the EC habitats (92/43/EEC) types: 1210 (Annual vegetation of drift lines), 2110 (Embryonic shifting dunes), 2120 (Shifting dunes along the shoreline with Ammophila arenaria), 2210 (Crucianellion maritimae fixed beach dunes), 2230 (Malcolmietalia dune grasslands), 2250 (Coastal dunes with Juniperus spp.), 2260 (Cisto-Lavanduletalia dune sclerophyllous scrubs), 2270 (Wooded dunes with Pinus pinea and/or P. pinaster). Table S6: Results of the harmonization test (error matrix and Kappa statistic) between phenology-based classes in the first hierarchical level of classification and the photo-interpreted classification map. Table S7: Results of the harmonization test (error matrix and Kappa statistic) between phenology-based classes and the photo-interpreted classification map in the second hierarchical level of classification. Table S8: Results of the harmonization test (error matrix and Kappa statistic) between phenology-based classes in the third hierarchical level of classification and the photo-interpreted classification map. Table S9: Results of the harmonization test (error matrix and Kappa statistic) between phenology-based classes in the first hierarchical level of classification and habitats of conservation concern (Habitats Directive; 92/43/EEC; Table 1) assigned on 2 m floristic plots collected in the field. Table S10: Results of the harmonization test (error matrix and Kappa statistic) of between phenology-based classes in the second hierarchical level of classification and habitats of conservation concern (Habitats Directive; 92/43/EEC; Table 1) assigned on 2 m floristic plots collected in the field. Table S11: Results of the harmonization test (error matrix and Kappa statistic) between phenology-based classes in the second hierarchical level of classification and habitats of conservation concern (Habitats Directive; 92/43/EEC; Table 1) assigned on 2 m floristic plots collected in the field.