Classiﬁcation of Atlantic Coastal Sand Dune Vegetation Using In Situ, UAV, and Airborne Hyperspectral Data

: Mapping coastal dune vegetation is critical to understand dune mobility and resilience in the context of climate change, sea level rise, and increased anthropogenic pressure. However, the identiﬁcation of plant species from remotely sensed data is tedious and limited to broad vegetation communities, while such environments are dominated by fragmented and small-scale landscape patterns. In June 2019, a comprehensive multi-scale survey including unmanned aerial vehicle (UAV), hyperspectral ground, and airborne data was conducted along approximately 20 km of a coastal dune system in southwest France. The objective was to generate an accurate mapping of the main sediment and plant species ground cover types in order to characterize the spatial distribution of coastal dune stability patterns. Field and UAV data were used to assess the quality of airborne data and generate a robust end-member spectral library. Next, a two-step classiﬁcation approach, based on the normalized di ﬀ erence vegetation index and Random Forest classiﬁer, was developed. Results show high performances with an overall accuracy of 100% and 92.5% for sand and vegetation ground cover types, respectively. Finally, a coastal dune stability index was computed across the entire study site. Di ﬀ erent stability patterns were clearly identiﬁed along the coast, highlighting for the ﬁrst time the high potential of this methodology to support coastal dune management.


Introduction
The formation and evolution of coastal sand dunes, which back sandy beaches, result from complex interactions between marine and aeolian sediment transport, plant ecosystem-engineering effects, coastal and beach topography, and storms [1][2][3][4][5][6][7]. When not destroyed by human activity, foredunes are formed by an accumulation of wind-blown sand, which is trapped by plant species tolerant to sand burial. Indeed, once pioneer species that are tolerant to salinity and sand burial colonize a shoreline, a positive feedback occurs between the trapped sand settling and the pioneer

Study Area
The study site is composed of beach-dune systems extending approximately 18 km alongshore in Aquitaine, southwest France, from the Horizon beach to Grand Crohot beach (Figure 1a). This preserved area, which is characterized by the absence of coastal resort and hard structures, represents an ideal natural laboratory that has been used for assessing beach-dune evolution [35,36], plant-plant interactions, and community composition [37][38][39][40]. Beaches are meso-macrotidal with a mean annual spring tide of 3.7 m [41]. The coast is exposed to high-energy waves generated in the North Atlantic Ocean. Winter wave energy shows large interannual variability, which can lead to dramatic coastal dune erosion during severe winters [41,42]. Contrary to further north along the coast, where chronic erosion can reach 5 m/y, the shoreline at the study area has been relatively stable over the past 65 years [43,44]. The large aeolian sand dunes have been formed through a mixture of natural Aeolian and biological processes and large-scale anthropogenic works, primarily [45]: (1) erection of a foredune by progressive fencing at the end of the 19th century and (2) extensive dune reprofiling with large-scale mechanical aid and marram planting [46,47]. Since then, the coastal dune system has been softly managed by marram grass planting and branching, particularly to prevent the development of blowouts. Current coastal dune height reaches 20 to 25 m above mean sea level, with a width of approximatively 250 m, and is composed, from the sea to the inland, of (1) the incipient foredune formed by Aeolian sand deposition; (2) the established foredune, a shore-parallel dune ridge formed on the backshore by Aeolian sand deposition and vegetation; (3) the transition dune; and (4) the grey dune, occupying the most inland position closest to the forest. Vegetation zonation follows these four morphological units in relation to plant species tolerance to physical disturbance (erosion and accretion due to marine and wind activity) ( Figure 1b) [4,37].
The incipient foredune (the most disturbed unit subject to marine and Aeolian action) is mainly colonized by Elymus farctus (Ely, Figure 1c). This pioneer species can rapidly grow in this habitat because of its high tolerance to sand deposition and the increase in sediment fertility provided by organic matter brought by the sea and the decrease in atmospheric stress due to the water spray [37]. The associated pioneer species (Euphorbia paralias and Eryngium maritimum, Figure 1c) mainly form diffuse patches with low vegetation cover (<30%) [37].
The established foredune (a less disturbed unit mainly subject to Aeolian action) is largely dominated by Ammophila arenaria (Amm, Figure 1c), a species which was used for dune building due to its long leaves, which trap the windblown sand and its dense root network promoting dune stabilization. The plantation of Ammophila arenaria is still used nowadays by coastal dune stakeholders for promoting dune recovery after severe winter storms. Other species are Otanthus maritimus (Ota, Figure 1c), generally pure low-density patches on the established foredune ocean face, and Galium arenarium (Gal, Figure 1c), forming dense patches of several tens of meters on the established foredune inland face when fully developed at the end of spring. The incipient foredune (IFS) and established foredune (EFS) sands are both characterized by a highly mobile surface sand (for mobile dune sand (MDS)), with a median grain size of 414 µm [40].
Behind the established foredune is the transition dune (more stable and less subject to sand burial), colonized by perennial species, primarily Helychrisum stoechas (Hel, Figure 1c) and Artemisia campestris (Art, Figure 1c), the latter developing dense patches extending a few meters horizontally. Here, the grains of sand are slightly finer (median grain size of 389 µm [40]) and are mixed with plant debris (grey dune sand (GDS)).
Phleum arenarium, Cerastium glomeratum, Vulpia membranacea, and Senecio vulgaris). These species are interspersed with each other in very dense features. The eastern boundary of the study site is defined by the iso-contour of topographic elevation z = 12 m, which corresponds to the middle of the grey dune inland face. The western boundary is defined as the beach-incipient foredune limit monitored on 21st May 2019 using an all-terrain vehicle (ATV), equipped with a differential global positioning system (DGPS). Finally, the grey dune (highly stable and mostly sheltered from sand sprinkling during winter storms) is at the junction between the herbaceous and forest environments. It is characterized by the presence of mosses, lichens, the perennial grass Corynephorus canescens, many small annuals (Phleum arenarium, Cerastium glomeratum, Vulpia membranacea, and Senecio vulgaris), and woody species such as Pinus maritima (Pin, Figure 1c). In addition, it is an environment of high floristic interest, Plant species and natural habitats previously described allow us to define 11 plant communities and sand cover classes (Amm, Gal, Hel, Art, Ely, Ota, MDS, GDS, sparse vegetation (SPV), grey dune vegetation (GDV), Pin, see Table 1), which, for the most part, play a major role in the morphological and ecological processes of the coastal dunes in southwest France [37]. SPV (i.e., SParse vegetation, Figure 1c) gathers some species of the incipient foredune (i.e., Euphorbia paralias and Eryngium maritimum), which cannot be radiometrically discriminated between them with a 1 m 2 resolution sensor due to the very low coverage and very small scale of features. For the same reasons, grey dune vegetation (GDV; Figure 1c) gathers some species of the grey dune (i.e., Corynephorus canescens, Phleum arenarium, Cerastium glomeratum, Vulpia membranacea, and Senecio vulgaris). These species are interspersed with each other in very dense features. The eastern boundary of the study site is defined by the iso-contour of topographic elevation z = 12 m, which corresponds to the middle of the grey dune inland face. The western boundary is defined as the beach-incipient foredune limit monitored on 21 May 2019 using an all-terrain vehicle (ATV), equipped with a differential global positioning system (DGPS).

Airborne Hyperspectral Data
On 18 June 2019, the 18 km coastline of the study site was flown over by a Chieftain PA31-350 aircraft from PIPER (Vero Beach, FL, USA) operated by PIXAIR company, equipped with a Visible Near InfraRed (VNIR) hyperspectral HySpex VNIR 1600 camera from Norsk Elektro Optikk, Skedsmokorset, Norway, operated by the Laboratory of Planetology and Geodynamics (LPG) and measuring 160 spectral channels (from 409.23 nm to 987.08 nm in steps of 3.61 nm, with a 4.5 nm spectral resolution) and an angle of aperture of 17 • . The record of the aircraft attitude, with the fast inertia measurement unit (IMU), coupled to the accurate global positioning system (GPS) of POS AV AP50 OEM (IMU-8) from Applanix (Richmond Hill, Ontario, Canada) mounted in the dual wavelength Ligth Detection and Ranging (LiDAR) Titan of Teledyne Optech Incorporated (Vaughan, Ontario, Canada) was provide by the LiDAR Nantes-Rennes university platform, and the trajectory was calculated by GEOFIT Expert Company (Nantes, France). The images acquired by the VNIR camera were transformed into luminance from a sensor factory calibration and were georeferenced using the airborne trajectory with Airborne Optical Scanner Data (PARGE) v3.4 software (ReSe Applications LLC, Wil, Switzerland) [48]. The overlap between each image was 45% to avoid gaps between irregular flight lines, but only the pixels with the best angle of incidence were kept to minimize the directional effect, as discussed in [48]. Moreover, as discussed in [49], this also minimizes the strong deviations induced by the different spatial resolutions and the distances to the target between field measurements and remote acquisitions [50,51]. Atmospheric corrections were then applied with Atmospheric and Topographic Correction (ATCOR-4) v7.3 software (ReSe Applications LLC, Wil, Switzerland) [48] in order to obtain a VNIR image with a resolution of 1 × 1 m. Using the ASD FieldSpec3 FR spectrometer (Malvern Panalytical Ltd., Royston, UK) of the LPG laboratory, a final spectral intensity adjustment was applied to insure the best matching between airborne and a single ground reflectance of dry sand, removing the weak anomalous spectral undulation that remained after the atmospheric compensation [52].

UAV Data
On 19 June 2019, a UAV survey covering 4 km of beach dune system at Truc Vert beach was conducted with a DJI Phantom 4 Pro quadricopter, equipped with a 20 MPix camera. Photogrammetry algorithms from the Agisoft Metashape software v1.5, together with the presence of 36 adequately distributed permanent ground control points, allowed generating of a 0.02 m resolution orthomosaic, with a planimetric error of 0.019 m and 0.026 m for the x axis and y axis, respectively [36]. Based on the primary authors expertise on photointerpretation of existing cover types from UAV imagery, the georeferenced orthomosaic was used to generate a suite of masks that allowed extracting control points from airborne VNIR images. A total of 1488 control points (ranging from 63 to 203 per class) were used to generate validation points in order to assess the performances of supervised classification. For only three classes (GDV, Pin, and SPV), control points also allowed the extraction of random training points, which are independent of validation points. First, dense and extensive areas associated with the 11 ground cover types were visually located and digitized on the orthomosaic. These masks were then used to extract the control points from the airborne VNIR image ( Figure 2).

Ground Hyperspectral Data
On 25 June 2019, a radiometric ground survey was carried out at Truc Vert beach in order to collect the in situ remote sensing reflectance spectra of different coastal dune cover types. Hyperspectral measurements were performed in the 400-900 nm spectral range (every 3 nanometers) from two radiometrically-calibrated TriOS-RAMSES sensors [53]. The radiance sensor was pointed downward to measure the upward signal of the substrate (Lu(λ), W·m −2 ·sr −1 ·nm −1 ), while the irradiance sensor was pointed upward to measure the downward irradiance (Ed (λ), W·m −2 ·nm −1 ). The radiance sensor had a field of view of 7°. It was mounted on a tripod 1.2 m above the ground ( Figure  3a) and provided field samples with a unit area of 0.017 m 2 . Radiance and irradiance measurements were performed simultaneously, as in [54]. The remote sensing reflectance (Rrs(λ), sr −1 ) was then calculated as the ratio of Lu(λ) to Ed(λ). Note that, for brevity, hereafter, we did not include the reference to the spectral dependency.
In complex and very fine-scale fragmented vegetation mosaics, due to a high intra-patch variability of density and vegetation cover, ground Rrs spectra are the result of a linear spectral mixture between plant species and soil contributions. Sampling strategies and the number of collected points, summarized in Table 1, were designed to collect the ground reference Rrs spectra (RG-REF) of each plant species and sand types with the lowest mixing percentage. Ground measurements were performed in 1 × 1 m 2 representative field sites, dominated by one ground cover type. For plant species with patches equal to or larger than 1 m 2 (Amm, Gal, Hel, Art), a 1 × 1 m 2 quadra was used to collect Rrs at nine points homogeneously distributed over features with dense cover (Figure 3b). For plant species characterized by isolated individuals or patches with sparse

Ground Hyperspectral Data
On 25 June 2019, a radiometric ground survey was carried out at Truc Vert beach in order to collect the in situ remote sensing reflectance spectra of different coastal dune cover types. Hyperspectral measurements were performed in the 400-900 nm spectral range (every 3 nanometers) from two radiometrically-calibrated TriOS-RAMSES sensors [53]. The radiance sensor was pointed downward to measure the upward signal of the substrate (L u (λ), W·m −2 ·sr −1 ·nm −1 ), while the irradiance sensor was pointed upward to measure the downward irradiance (E d (λ), W·m −2 ·nm −1 ). The radiance sensor had a field of view of 7 • . It was mounted on a tripod 1.2 m above the ground ( Figure 3a) and provided field samples with a unit area of 0.017 m 2 . Radiance and irradiance measurements were performed simultaneously, as in [54]. The remote sensing reflectance (R rs (λ), sr −1 ) was then calculated as the ratio of L u (λ) to E d (λ). Note that, for brevity, hereafter, we did not include the reference to the spectral dependency.

End-Member Spectra
End-member spectra, associated with the 11 vegetation and sand cover classes, were generated from airborne VNIR training points. For SPV, GDV, and Pin habitats, 50 Rrs spectra per ground cover type were randomly extracted from the control points associated with their respective masks. For other types of ground cover, methodology was based on the extraction of 3 × 3 sub-images centered on field measurement locations. These sub-images allowed us to compute for each ground cover type the mean Rrs spectra, so called here, the airborne reference Rrs spectra (RA-REF), and the airborne reference standard deviation. The analysis of the spectra extracted from control points shows that the spectral variability of Rrs for MDS, Ely, Amm, Gal, and Art follows a normal distribution for any In complex and very fine-scale fragmented vegetation mosaics, due to a high intra-patch variability of density and vegetation cover, ground R rs spectra are the result of a linear spectral mixture between plant species and soil contributions. Sampling strategies and the number of collected points, summarized in Table 1, were designed to collect the ground reference R rs spectra (R G-REF ) of each plant species and sand types with the lowest mixing percentage. Ground measurements were performed in 1 × 1 m 2 representative field sites, dominated by one ground cover type. For plant species with patches equal to or larger than 1 m 2 (Amm, Gal, Hel, Art), a 1 × 1 m 2 quadra was used to collect R rs at nine points homogeneously distributed over features with dense cover (Figure 3b). For plant species characterized by isolated individuals or patches with sparse vegetation covers (Ota and Ely) and for sand types (IFS, EFS, and GDS), three replicates of R rs measurements were performed on individuals, showing the highest vegetation cover (Figure 2c), and on spatially homogeneous sites (Figure 3d), respectively. Then, for each ground cover type, R G-REF was computed by averaging all R rs spectra.
It is important to note that no field measurements were performed for SPV, GDV, or Pin. SPV and GDV habitats were characterized by a high spatial variability of intra-patch diversity, which was difficult to record. In addition, the instruments were not suitable to measure Pin species reflectance, due to their height (>2 m).
Associated with radiometric measurements, geo-positioning data, using a differential global navigation satellite system (DGNSS), were collected in order to generate matchups for comparison between ground and airborne hyperspectral data. Finally, a linear interpolation was applied to field data in order to make them compatible with the spectral bands of airborne data.

End-Member Spectra
End-member spectra, associated with the 11 vegetation and sand cover classes, were generated from airborne VNIR training points. For SPV, GDV, and Pin habitats, 50 R rs spectra per ground cover type were randomly extracted from the control points associated with their respective masks. For other types of ground cover, methodology was based on the extraction of 3 × 3 sub-images centered on field measurement locations. These sub-images allowed us to compute for each ground cover type the mean R rs spectra, so called here, the airborne reference R rs spectra (R A-REF ), and the airborne reference standard deviation. The analysis of the spectra extracted from control points shows that the spectral variability of R rs for MDS, Ely, Amm, Gal, and Art follows a normal distribution for any wavelengths (p-value of the Shapiro-Wilk normality test <0.01). For GDS, Ota, and Hel, the normality is rejected only for near-infrared (NIR) bands at the 0.01 level but is accepted at the 0.001 level. Based on these results, 50 R rs spectra per ground cover type were further randomly computed using a normal probability law in order to generate the training dataset. In order to take into account all the cross-shore variability associated with MDS, this class was specifically composed of 25 spectra of IFS and 25 spectra of EFS. In total, the vegetation and sand reference dataset consisted of 550 spectra. This procedure allowed reproducing of a spectral mixing range statistically representative of the study site.
The mixing percentage, α, ranging from 0 to 1, could be derived from the linear spectral mixture model: where R mod is the modeled R rs , R G-REF;Sdi is the ground reference R rs spectrum of sand type i, and R G-REF;Spj isthe ground reference R rs spectrum of plant species j. A root mean square cost function between airborne R rs and R mod was implemented to estimate α values and then its statistical distribution was generated for each plant species.

Normalized Difference Vegetation Index Filter
Prior to applying a supervised classification of the airborne VNIR image, we conducted a first experiment using a spectral index in order to differentiate pixels with and without vegetation covers. This first step was crucial in sparse and very fine-scale fragmented vegetation mosaics. In 1 × 1 m 2 pixels with a low but ecologically significant presence of vegetation, the percentage of sand contribution in the mixing R rs spectrum can be high and may generate misclassified points between vegetation and sand. The normalized difference vegetation index (NDVI) has been proven to be useful for ground cover type characterization [55][56][57], including, more specifically, coastal dune areas [24]. The NDVI was computed as: Remote Sens. 2020, 12, 2222 10 of 25 where the subscripts NIR and R are associated with the near-infrared and red bands, respectively. The differentiating procedure of sand and vegetation covers is based on a NDVI threshold value. The most relevant value was objectively determined using a statistical analysis on the end-member NDVI values.

Random Forest
After having separated the observations into two separate classes (sand or vegetation cover) based on the NDVI index, a Random Forest (RF) classifier [58] was considered for pixel-based classification, as shown in Figure 4. RF is a suitable algorithm for coastal vegetation classification, because of its stability and ability to discriminate ground cover differences [24,32,54,59,60]. When applying the RF classifier, the classification map is accompanied by an accuracy map, which contains the indicator of the confidence degree of the classification. More precisely, the RF classifier is composed of a set of different individual base classifiers (decision trees), each tree giving a decision, and the final decision corresponds to the majority vote. The degree of confidence is obtained as the percentage of individual classifiers that provided a decision equal to the final one.

Stability Index
Based on seasonal surveys of vegetation abundance at this coastal dune system, [40] showed that plant species were spread out according to their tolerance to disturbance, which, in this system, is primarily driven by sand burial. This distribution can be seen in Figure A1 in the correspondence analysis (CA), where plant species have a perfect cross shore distribution along the first axis, representing a disturbance gradient, going from the most disturbed habitat (incipient foredune) to the most stable one (grey dune). Thus, it is possible to extract the CA axis 1 coordinate of each plant species and to associate a value between 0 (highly disturbed) and 1 (highly stable). A value of 0 was associated to the bare sand (MDS and GDS), 1 to Pins and the value of SPV and GDV corresponded to the mean axis value of species representing theses classes (i.e., Euphorbia paralias and Eryngium maritimum for SPV and Phleum arenarium, Cerastium glomeratum, Corynephorus canescens, Vulpia membranacea, and Senecio vulgaris for GDV) ( Table 2), calculating a mean stability index for each cross shore transect.  An accuracy assessment was performed through the confusion matrix, with dimensions equal to the number of classes, from which the overall accuracy was extracted. The agreement between classified points and ground truth was based on a dataset of 550 validation points, i.e., 50 validation points were randomly extracted for each of the 11 cover masks. For SPV, GDV, and Pin vegetation, we removed training points prior to the extraction of validation points.
For SPV, GDV, and Pin vegetation, we removed training points prior to the extraction of validation points. Validation points, training points, and mean ground hyperspectral reflectance spectra are freely available at the Supplementary Materials.

Stability Index
Based on seasonal surveys of vegetation abundance at this coastal dune system, [40] showed that plant species were spread out according to their tolerance to disturbance, which, in this system, is primarily driven by sand burial. This distribution can be seen in Figure A1 in the correspondence analysis (CA), where plant species have a perfect cross shore distribution along the first axis, representing a disturbance gradient, going from the most disturbed habitat (incipient foredune) to the most stable one (grey dune). Thus, it is possible to extract the CA axis 1 coordinate of each plant species and to associate a value between 0 (highly disturbed) and 1 (highly stable). A value of 0 was associated to the bare sand (MDS and GDS), 1 to Pins and the value of SPV and GDV corresponded to the mean axis value of species representing theses classes (i.e., Euphorbia paralias and Eryngium maritimum for SPV and Phleum arenarium, Cerastium glomeratum, Corynephorus canescens, Vulpia membranacea, and Senecio vulgaris for GDV) ( Table 2), calculating a mean stability index for each cross shore transect.

Analysis of Spectral Signatures of Ground Cover Types
Mean R G-REF displays sharp spectral characteristics, which allows discriminating the different ground cover types (Figure 5a). The spectral signature of MDS and GDS is monotonically increasing from blue to near-infrared (NIR) wavelengths. For any spectral bands, GDS shows lower R rs values than MDS, due to the presence of high concentrations of black organic debris strongly absorbing the incident irradiance. R G-REF spectra of Ely, Gal, and Art show a typical spectral behavior of plant covers with a minimum value at 680 nm and a sharp increase in reflectance values in the red edge region (680-750 nm), due to leaf pigment absorption in red and leaf cellular scattering in NIR spectral bands, respectively [30]. For Ota, Amm, and Hel, R G-REF spectra are the result of the spectral mixture between sand and plant species. The sand percentage is higher for Ota, which has very scattered patches. It is weaker for Hel. The sampled features have very sharp outlines and are very dense. However, their feature surface is less than 1 m 2 , and two of the nine sampled pixels are dominated by GDS cover.
Field and airborne hyperspectral data show good agreement, except for Gal and Art. IFS and EFS, which are characterized by a high optical homogeneity of the substrate into the 3 × 3 box, have the lowest mean relative difference (MRD), with a value of 1.0% and 3.7%, respectively (Figure 5b). Spectral agreement is better for visible wavelengths (0.8% and 1.4% for IFS and EFS, respectively) than for NIR bands (1.3% and 6.9%). GDS, Amm, and Hel, which are characterized by a slight heterogeneity into the 3 × 3 box, have MRD values of 6.3%, 9.9%, and 7.5%, respectively (Figure 5c-e). On the other hand, field and airborne R rs spectra of Gal and Art are significantly different in magnitude, with MRD values of 57.7% and 92.5%, respectively. These differences are due to the non-alignment between airborne grid and field quadrats. While field measurements recorded R rs over dense and homogeneous plant features, airborne measurements recorded R rs of ground covers composed of sand and vegetation. The mixing percentages (see Equation (2)) associated with the Gal and Art R A-REF are estimated to 51% (MRD = 3.5%) and 46% (MRD = 2%), respectively. than for NIR bands (1.3% and 6.9%). GDS, Amm, and Hel, which are characterized by a slight heterogeneity into the 3 × 3 box, have MRD values of 6.3%, 9.9%, and 7.5%, respectively (Figure 5ce). On the other hand, field and airborne Rrs spectra of Gal and Art are significantly different in magnitude, with MRD values of 57.7% and 92.5%, respectively. These differences are due to the nonalignment between airborne grid and field quadrats. While field measurements recorded Rrs over dense and homogeneous plant features, airborne measurements recorded Rrs of ground covers composed of sand and vegetation. The mixing percentages (see Equation (2)) associated with the Gal and Art RA-REF are estimated to 51% (MRD = 3.5%) and 46% (MRD = 2%), respectively. Rrs spectra of sand ground covers (Figure 6a,b) display a high spectral variability, which takes into account the cross-shore variability of sand and the variability of black organic debris concentration for MDS and GDS, respectively. Spectral characteristics of plant ground cover Rrs of established foredunes and transition dunes show a high mixing percentage with sand. High similarity is observed between Ely (Figure 6c  R rs spectra of sand ground covers (Figure 6a,b) display a high spectral variability, which takes into account the cross-shore variability of sand and the variability of black organic debris concentration for MDS and GDS, respectively. Spectral characteristics of plant ground cover R rs of established foredunes and transition dunes show a high mixing percentage with sand. High similarity is observed between Ely (Figure 6c), Ota (Figure 6d), SPV (Figure 6e), and MDS spectra and Amm (Figure 6f), Hel (Figure 6h), and GDS spectra, respectively. On the other hand, R rs spectra of GDV and Pin classes show spectral signatures significantly different than R rs spectra of sand classes. The spectral similarity between some mixed plant species and sand classes can generate misclassified pixels, which impact the performance of supervised classifications. Prior to applying supervised classifications, NDVI was used to differentiate sand and plant ground covers. Several combinations of spectral bands were tested, and the best results were achieved using the red edge bands of 680 nm and 750 nm. Figure 7 shows statistics of NDVI values for the 11 ground cover classes. The first, second, and third quartiles of MDS and GDS were 0.019, 0.025, 0.038, and 0.038, 0.049, 0.061, respectively, while they were 0.066, 0.074, 0.081 for Hel, which displayed the lowest difference. A nonparametric one-way analysis of variance (Kruskal-Wallis test), followed by nonparametric pair-wise comparisons (Wilcoxon-Mann-Whitney test), highlights that NDVI median values between sand and plant species classes are statistically different (p-value ≤ 0.001). Based on the analysis of the NDVI value distribution of GDS and Hel, we selected the threshold value of 0.061 to differentiate sand and plant ground covers. The choice of the third quartile of GDS NDVI values allows to classify sparse vegetation covers as vegetation (see Figure A2 for the NDVI map).
nonparametric one-way analysis of variance (Kruskal-Wallis test), followed by nonparametric pairwise comparisons (Wilcoxon-Mann-Whitney test), highlights that NDVI median values between sand and plant species classes are statistically different (p-value ≤ 0.001). Based on the analysis of the NDVI value distribution of GDS and Hel, we selected the threshold value of 0.061 to differentiate sand and plant ground covers. The choice of the third quartile of GDS NDVI values allows to classify sparse vegetation covers as vegetation (see Figure A2 for the NDVI map).

Assessment of Pixel-Based Supervised Classification Performances
The results of the RF classifier on the 550 validation points are shown in Tables 3 and 4 for sand and vegetation classes, respectively. The discretization between the two sand types is excellent, with an overall accuracy of 100% (Table 3). For vegetation, the overall accuracy is 82.67% (Table 4). Only three classes have a precision lower than 90%; Ely (62%), Ota (10%), and Gal (86%). The poor performance for Ota classification can be attributed to a low variability of its spectral signature in the training dataset (see Figure 6d) and a high similarity of Rrs spectra with the other plant classes of the established foredune. For example, 32% of misclassified pixels are associated with SPV. In coastal

Assessment of Pixel-Based Supervised Classification Performances
The results of the RF classifier on the 550 validation points are shown in Tables 3 and 4 for sand and vegetation classes, respectively. The discretization between the two sand types is excellent, with an overall accuracy of 100% (Table 3). For vegetation, the overall accuracy is 82.67% (Table 4). Only three classes have a precision lower than 90%; Ely (62%), Ota (10%), and Gal (86%). The poor performance for Ota classification can be attributed to a low variability of its spectral signature in the training dataset (see Figure 6d) and a high similarity of R rs spectra with the other plant classes of the established foredune. For example, 32% of misclassified pixels are associated with SPV. In coastal dune environments, Ota population is composed of isolated individuals and patches with sparse vegetation cover. A spatial resolution of 1 m 2 is clearly a major limitation for the identification of scattered plant species with low ground coverage in complex vegetation mosaic systems. Table 3. Confusion matrix of the Random Forest (RF) classification for the two sand classes.  Despite an overall accuracy reaching 82.67%, which already represents a good agreement [62], a second RF model was learned by integrating Otanthus maritimus in the SPV class (Table 5). For this second RF vegetation classification, Ely and Gal are still under 90% accuracy, with 64% and 82%, respectively. According to the confusion matrix, the misclassified Ely and Gal pixels are in SPV. Indeed, in their undeveloped stage, Ely and Gal are characterized by small plants of about 10 cm and are quite isolated on a surface area of 1 m 2 . It is only when they are well developed that they can easily cover this surface and thus be well classified. However, the overall accuracy of the vegetation species classification reaches 92.25%, which represents a strong agreement [62].
The RF classifier also highlights the importance of the spectral channels used in the classification. Thus, it was possible to study the impact of the number of bands used (from the most to the least important) on the overall accuracy of the vegetation classification, which increases non-linearly with the number of channels used (left axis on Figure 8a). In fact, there seems to be a threshold of 19 bands that allows the accuracy to increase from 76.75% (with 18 channels) to 83.5% before asymptotically increasing. This sudden increase in accuracy was obtained by adding the wavelength 675 nm, whereas all the bands used before were between 750 nm and 900 nm, highlighting the importance of NIR wavelengths in the classification of plant species (right axis on Figure 8a). In addition to the spectral resolution, the impact of the spatial resolution on the vegetation classification accuracy was evaluated by progressively degrading the pixel size from 1 m 2 to 20 m 2 . Our classification procedure was then applied to the entire image for the eight vegetation classes. The classification accuracy ranged from 92.75% for 1 m 2 pixel to 73.02% for 2 m 2 pixel then decreased drastically below 70% for larger sizes (Figure 8b). Table 5. Confusion matrix of the second RF classifier for the eight vegetation classes (after merging Ota and SPV).

Ely
Spv The RF classifier also highlights the importance of the spectral channels used in the classification. Thus, it was possible to study the impact of the number of bands used (from the most to the least important) on the overall accuracy of the vegetation classification, which increases non-linearly with the number of channels used (left axis on Figure 8a). In fact, there seems to be a threshold of 19 bands that allows the accuracy to increase from 76.75% (with 18 channels) to 83.5% before asymptotically increasing. This sudden increase in accuracy was obtained by adding the wavelength 675 nm, whereas all the bands used before were between 750 nm and 900 nm, highlighting the importance of NIR wavelengths in the classification of plant species (right axis on Figure 8a). In addition to the spectral resolution, the impact of the spatial resolution on the vegetation classification accuracy was evaluated by progressively degrading the pixel size from 1 m 2 to 20 m 2 . Our classification procedure was then applied to the entire image for the eight vegetation classes. The classification accuracy ranged from 92.75% for 1 m 2 pixel to 73.02% for 2 m 2 pixel then decreased drastically below 70% for larger sizes (Figure 8b).

Spatial Distribution of Coastal Dune Vegetation
The aerial image (Figure 9a) shows a relatively alongshore-uniform coastal dune, which is well vegetated, except at the southern end between −6 and −10 km, characterized by a much less dense and heterogeneous dune vegetation cover and a large bare sand surface. The RF classification map (Figure 9b) shows that each plant species is well classified in its ecosystem (e.g., Ely in incipient foredunes, Amm in established foredunes), with an accuracy generally better than 50% in the transition dune and better than 70% in the other dune units (Figure 9c). The alongshore stability index is in line with visual observation (Figure 9a), with an average value of 0.37 (std = 0.09, se = 7.52 × 10 −4 ) between the northern limit and −6 km alongshore. Some variations are also locally observed due to the presence of access tracks. In the southern sector, between −6 and −10 km, a high alongshore variability is observed, with a mean stability index of 0.33 (std = 0.15, se = 2.5 × 10 −3 ) (Figure 9d).
( Figure 9b) shows that each plant species is well classified in its ecosystem (e.g., Ely in incipient foredunes, Amm in established foredunes), with an accuracy generally better than 50% in the transition dune and better than 70% in the other dune units (Figure 9c). The alongshore stability index is in line with visual observation (Figure 9a), with an average value of 0.37 (std = 0.09, se = 7.52 × 10 −4 ) between the northern limit and −6 km alongshore. Some variations are also locally observed due to the presence of access tracks. In the southern sector, between −6 and −10 km, a high alongshore variability is observed, with a mean stability index of 0.33 (std = 0.15, se = 2.5 × 10 −3 ) (Figure 9d). The RF classification also highlights the percentage of appearance of each class (Table 6) with the predominance of the dominant species, such as Amm representing the established foredune (10.64%), Hel representing the transition dune (17.77%), and GDV (7.40%) representing the grey dune. The presence of MDS is generally higher than GDS because of the greater plant cover in the transition dune and grey dune than in the incipient and established foredunes. Only the central area ( Figure 9j) shows a lower percentage of MDS than GDS. This can be mainly attributed to a narrower width of the established foredune.   The large bare sand areas in the southern area, between −6 km and −10 km, are clearly visible in the 1.1 km alongshore zoom (Figure 9m,n), where they represent 62.83% of the surface. These areas, with a low vegetation cover and quite a low and variable stability index (mean value of 0.23, with std = 0.10, se = 2.9 × 10 −3 ) characterize a less stable and therefore possibly more mobile sand dune ( Figure 9p and Table 6). To limit such dune mobility and prevent blowout development, coastal dune stakeholders have developed various techniques, including the laying of tree branches. This technique, which aims to stop Aeolian erosion and favor sand deposition, is often used to stabilize active dune blowout and/or to protect the inland infrastructures (roads, buildings) or forests from the windblown sand [46,47]. Tree branches within the established foredune are visible in the zoom into the southern zone. Nevertheless, the class associated with this type of cover has not been taken into account in this study due to the low number of control points. As a consequence, tree branch cover is mainly classified as GDV generating misclassified pixels. This class should be taken into account in future research.
The zoom into the northern zone (Figure 9e,f) shows a greater stability, with an average stability index of 0.31 with some alongshore variations (std = 0.06, se = 1.7 × 10 −3 ) (Figure 9h). Indeed, the established foredune width is not uniform alongshore and shows variations by several tenths of meters in width. Such variability is caused by the localized presence of developing blowouts, leading to the presence of sand deposition fields across the transition dune. In turn, this promotes the presence of vegetation tolerant to physical disturbance (i.e., Ely and SPV) over a large cross-shore distance of about 75 m.
Finally, the zoom on the central area (Figure 9i,j) shows greater stability than the two previous sectors, with an almost alongshore-uniform stability index with an average of 0.36 (std = 0.05, se = 1.6 × 10 −3 ) (Figure 9i). Plant species characteristics of disturbed environments are present only over a short cross-shore distance of less than 50 m. Moreover, the classification highlights the recolonization of the vegetation of a sand deposition field, extending approximately 100 m in the cross shore direction across the transition dune that formed during the severe 2017-2018 winter storms [36], covering plant species that are not tolerant to sand burial with about 10 cm of sand.

Discussion
The first objective of this paper was to develop a pixel-based classification approach in order to generate a comprehensive coastal dune cover mapping, including a large number of classes. Such classes are critical for coastal dune managers in the short term, e.g., to address blowout development or the impact of fencing or branching on vegetation dynamics. In the long-term, for example, decade(s), such mapping will provide new insight into the large-scale variability of vegetation ground cover type in relation to natural (e.g., marine erosion) and anthropogenic (e.g., dune reshaping) forcings. However, a major drawback for remotely sensed application in fragmented and fine-scale environments is the high spectral similarity between different vegetation cover types due to spectral mixing with sand.
To address this issue, recent research has focused on phenology-based [24] or spectral library-based [34] approaches. For high-energy eroding sandy coasts, the high mobility of MDS significantly reduces the interest of the former approach.
Spectral library approaches require specific procedures to generate end members [34,63]. To take into account the range of natural variability of spectral mixing into each class, matchups between ground and airborne data are used to compute airborne reference reflectance spectra (R A-REF ) and its standard deviation. The normal probability law is then used to randomly generate end-member spectra of each class. In a pre-analysis, two other procedures were tested using ground reference reflectance spectra (R G-REF ) and R A-REF spectra, respectively. End-member spectra were then modeled using the linear spectral mixture model for different ranges of α values, selected from the statistical analysis of training points. The highest classification performances using RF classifiers were obtained for the procedure described in this paper. The performances of our approach, with an overall accuracy of 92.3% for plant species classes, highlight that the normal probability law procedure is relevant to generate end members representative of the natural variability of spectral mixing into each vegetation cover type.
However, some limitations should be noted. At the dune system scale, the cross-shore distribution of species was in line with the location of their ecological niche. At the community scale, species with dense, large, and well-delimited patches displayed high accuracy (from 96% to 100%), while species with isolated or scattered patches (especially during the early stages of development) showed a lower accuracy (from 64 to 82%). Moreover, if the overall accuracy of SPV and GDV classes reaches 100%, no proportion of the different merged species can be provided. To enhance the monitoring and detection of fragmented and interspersed species, end member spectral libraries should be based on hyperspectral ground data [34,63]. This requires a specific sampling design to match the scale of ground data with that of large-scale observations [63] or to apply up-scaling procedures [61].
The second objective was to characterize the coastal dune stability patterns using sediment and vegetation cover mapping. Results highlight a high spatial variability of the stability index, which can be used to identify areas that may be more vulnerable to the storms. Vulnerability of coastal dunes is often attributed to a lack of vegetation cover, which reduces stability. This was observed in the southern part of the study area. This sector was characterized by a low vegetation cover and wide sand bare surfaces, where the average stability index was only 0.33, whereas it reached an average of 0.37 in the central and northern part, where the vegetation cover was denser, leaving less space for the sand. This methodological approach can become an efficient tool for researchers and coastal dune stakeholders. For instance, the mapping of plant species and sandy areas may allow the monitoring and detection of the different stages of blowout development and detection of the colonization of sandy areas by pioneer species, shrubs, and Pinus maritima colonization in the transition dune. However, for specific monitoring of small and sparse plant species, hybrid hyperspectral pixel-based and high resolution object-oriented approaches should be considered, either with coupling of optical satellite missions or with sensors on board UAVs, depending on the surface area to be monitored.

Conclusions
Recent technological and methodological advances in remote sensing and mapping tools have greatly facilitated the monitoring of spatiotemporal distribution of vegetation at species and habitat scales. In order to improve our understanding of coastal dune evolution, we developed an original pixel-based classification approach to detect a large range of ground cover classes for environments characterized by fragmented and fine-scale landscape patterns. The development of our classification approach was based on a comprehensive multi-scale dataset combining UAV, hyperspectral ground, and airborne data, collected along a 20-km stretch of a coastal dune system in southwest France. The results showed high classification performances with regard to the large number of cover classes and their low dissimilarities due to high spectral mixing levels. The high performances highlight that the methodology used to generate the end member spectral library meets representativeness requirements. Additionally, the two-step approach allows us to reduce misclassified pixels. Finally, the stability index mapping, derived from the sediment and vegetation mapping, provides large scale insight into the vulnerability of coastal dunes. We believe that the methodology developed in this paper can become an efficient tool for supporting coastal dune management.  The pixels at the western position of the beach-incipient foredune limit and the pixels at the eastern position of the iso contour z = 12 m have been filtered beforehand.