Development of Methods for Satellite Shoreline Detection and Monitoring of Megacusp Undulations

Citation: Angelini, R.; Angelats, E.; Luzi, G.; Masiero, A.; Simarro, G.; Ribas, F. Development of Methods for Satellite Shoreline Detection and Monitoring of Megacusp Undulations. Remote Sens. 2024,16, 4553. https:// doi.org/10.3390/rs16234553 Academic Editor: Chung-Ru Ho Received: 14 October 2024 Revised: 14 November 2024 Accepted: 2 December 2024 Published: 4 December 2024 Copyright: © 2024 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 (https:// creativecommons.org/licenses/by/ 4.0/). Article Development of Methods for Satellite Shoreline Detection and Monitoring of Megacusp Undulations Riccardo Angelini 1, Eduard Angelats 2, Guido Luzi 2, Andrea Masiero 3,* , Gonzalo Simarro 4 and Francesca Ribas 5 1Department of Civil and Environmental Engineering, University of Florence, via di Santa Marta 3, 50139 Florence, Italy; riccar[email protected] 2Geomatics Research Unit, Centre Tecnològic de Telecomunicacions de Catalunya (CTTC/CERCA), Av. Carl Friedrich Gauss 7, 08860 Castelldefels, Spain; [email protected] (E.A.); [email protected] (G.L.) 3Interdepartmental Research Center of Geomatics (CIRGEO), University of Padova, via dell’Università 16, 35020 Legnaro, Italy 4Department of Marine Geosciences, Institut de Ciències del Mar, Passeig Marítim de la Barceloneta 37-49, 08003 Barcelona, Spain; simarr[email protected] 5Physics Department, Universitat Politècnica de Catalunya, Carrer Jordi Girona 1-3, 08980 Barcelona, Spain; [email protected] *Correspondence: andrea.masier[email protected] Abstract: Coastal zones, particularly sandy beaches, are highly dynamic environments subject to a variety of natural and anthropogenic forcings. Instantaneous shoreline is a widely used indicator of beach changes in image-based applications, and it can display undulations at different spatial and temporal scales. Megacusps, periodic seaward and landward shoreline perturbations, are an example of such undulations that can significantly modify beach width and impact its usability. Traditionally, the study of these phenomena relied on video monitoring systems, which provide high-frequency imagery but limited spatial coverage. Instead, this study explored the potential of employing multispectral satellite-derived shorelines, specifically from Sentinel-2 (S2) and PlanetScope (PLN) platforms, for characterizing and monitoring megacusps’ formation and their dynamics over time. First, a tool was developed and validated to guarantee accurate shoreline detection, based on a combination of spectral indices, along with both thresholding and unsupervised clustering techniques. Validation of this shoreline detection phase was performed on three micro-tidal Mediterranean beaches, comparing with high-resolution orthomosaics and in-situ GNSS data, obtaining a good subpixel accuracy (with a mean absolute deviation of 1.5–5.5 m depending on the satellite type). Second, a tool for megacusp characterization was implemented and subsequent validation with reference data proved that satellite-derived shorelines could be used to robustly and accurately describe megacusps. The methodology could not only capture their amplitude and wavelength (of the order of 10 and 100 m, respectively) but also monitor their weekly–daily evolution using different potential metrics, thanks to combining S2 and PLN imagery. Our findings demonstrate that multispectral satellite imagery provides a viable and scalable solution for monitoring shoreline megacusp undulations, enhancing our understanding and offering an interesting option for coastal management. Keywords: shoreline extraction; megacusps; Sentinel-2; PlanetScope; multispectral imagery 1. Introduction Coastal zones are highly complex and dynamic from a geophysical point of view. Simultaneously, their management is challenging because they are intensively used by humans (with strong economic interests) and, at the same time, provide invaluable habitats of high biodiversity. Sandy beaches are subject to direct anthropogenic pressure through the construction of coastal infrastructures, which are already producing chronic erosion [ 1 ]. In Remote Sens. 2024,16, 4553. https://doi.org/10.3390/rs16234553 https://www.mdpi.com/journal/remotesensing Remote Sens. 2024,16, 4553 2 of 26 addition, land subsidence and the decrease in sediment supply by dam presence accelerate the sand loss processes in most delta areas [ 2 ] and climate change and the induced mean sea level rise will worsen these trends [3]. For monitoring and managing coastal changes over time, an essential indicator is the shoreline position, i.e., the interface or physical boundary between land and water [4]. Shorelines along sandy beaches may exhibit diverse morphological configurations at different spatiotemporal scales due to variations in hydrodynamic processes and sediment movement [ 5 , 6 ]. One of these features is megacusps, characterized by undulations along the shoreline often linked to the presence of submerged rhythmic sand bars and channels [ 7 ]. Megacusp horns are the seaward protrusions of the shoreline, whilst their embayments are the landward indentations. Typically, these features show alongshore wavelengths of several hundred meters and cross-shore amplitudes of a few tens of meters [ 8 ]. The formation of megacusps during a few days of low-energy wave conditions is frequently associated with the onshore migration of crescentic bars (also called rip channel systems) and an overall accretion of the shoreline [ 9 ]. Megacusps also occur during storm events and can lead to considerable erosion of beaches and dunes at the locations of megacusp embayments [ 10 ]. Alongshore migration up to tens of meters per day has also been reported. Understanding the dynamics of the formation and disappearance of such coastal patterns is crucial for coastal zone management because their embayments can become erosional hot spots that worsen existing (or future) sand losses. Moreover, potential human interventions would strongly benefit from systematic monitoring of the phenomena. Such monitoring is also important for beach safety because megacusps can be linked to the presence of strong offshore directed currents that might be dangerous for swimmers [7]. Only a few studies have delved into monitoring megacusp dynamics. In situ Global Navigation Satellite System (GNSS) surveys produce detailed data sets (see a list with existing observations in de Swart et al. [8] ), but they typically cover small geographic areas and do not have enough temporal resolution. Airborne imagery collection requires specifically commissioned flights that can be expensive and thereby also produce scarce data. High-time-resolution megacusp monitoring during long time periods has only been conducted so far in a couple of sites with video systems [ 8 , 11 ]. They have the advantage of providing continuous observations, but they cover a limited stretch of the beach (about 1 km). Exceptionally, megacusp dynamics were also monitored using weekly LiDAR data [12], but covering only a 2 km beach sector and a 1-year duration. In the last decade, the use of satellites that offer free medium-resolution images like Sentinel-2 (S2) from the European Space Agency Copernicus program, as well as the Landsat family, has been widely adopted for shoreline monitoring of large areas [ 1 , 13 ]. One disadvantage of the used satellite products to date is attributable to the spatial resolution, between 10 and 30 m. Recently, PlanetScope (PLN) (CubeSat satellites) provided free 3 m spatial resolution images for research purposes [ 14 ]. Parallel to this, the development of “data cubes” like Google Earth Engine (GEE) [ 15 ] has made these data even more accessible and manageable for long-term monitoring purposes. The presence of megacusps can be seen in satellite images [ 16 ] and this opens the door to using them for megacusp monitoring. However, no quantification of megacusp characteristics from satellite imagery has been performed yet. Developing a methodology to achieve this would complement the existing video monitoring studies of megacusps because satellites enable regional-scale monitoring, offering a broader perspective on coastal morphology changes. Moreover, the combination of S2 and PLN images opens the door to capturing the weekly–daily development of megacusps. This methodological study aims to test the potential of satellite multispectral images to automatically detect and characterize megacusps, including their dynamics at a weekly time scale. For this, two methodologies were developed and tested, one for instantaneous shoreline detection, and another for megacusp monitoring. The developed shoreline extraction methodology enabled the use of images from different satellites such as S2 and PLN. Subsequently, different metrics for detecting the presence of megacusps and quantifying Remote Sens. 2024,16, 4553 3 of 26 their wavelength and amplitude were explored. The methodologies were tested on three micro-tidal Mediterranean beaches, where numerous types of in situ validation data (GNSS data and aerial orthomosaics) were available. The novelty of this methodological study lies in investigating the use of satellite-derived shorelines for characterizing coastal morphologies such as megacusps. This includes not only capturing them in a single moment but also monitoring their temporal evolution. This might allow studying in future works the relationship between megacusp dynamics and beach forcing, as well as implementing beach management strategies. Section 2of the document reviews a collection of recent studies related to shoreline extraction from satellite imagery. Section 3outlines the study sites, the data used for analysis and validation, as well as the methodologies developed for shoreline extraction and megacusp characterization. Section 4presents the results, including validation of the two methodologies and an example of their application to characterize two events of megacusps. Section 5discusses the results, compares them with previous studies, and highlights the potential error sources as well as the strengths and limitations of the proposed methodologies. Finally, the most important findings are summarized in Section 6. 2. Background on Satellite Shoreline Extraction Recently, with the increasing availability of data, tools for shoreline extraction from multispectral satellite images, such as CoastSat [ 17 ], SHOREX [ 18 ], CASSIE [ 19 ], and SAET [ 20 ], have emerged. These tools, which can be either open-source or closed, are capable of achieving subpixel accuracy in delineating micro-tidal beach shorelines. During the image preprocessing step, which may include classification, a range of pixel and object-based techniques are used to label the satellite image pixels as water, land, or other classes [ 21 ]. The most prevalent methods in existing studies for distinguishing between water and non-water pixels involve the use of water indices [ 22 ]. These indices execute pixel-level multiband mathematical operations, leveraging the significant drop in surface reflectance from visible to near-infrared wavelengths as the primary differentiation criterion. Various methodologies can be employed to classify between land and water. An example is single thresholding with either fixed or variable thresholds, which the user can set manually or determine automatically [ 21 ]. In recent decades, the advent of machine learning (ML) techniques has supported the classification step, allowing the automated extraction of shorelines at large scales and the identification of patterns and trends that are difficult to discern with conventional methods [23]. Within the domain of water indices, the Normalized Difference Water Index (NDWI), a combination of green and near-infrared (NIR) bands, stands as one of the most used among coastal researchers [ 24 ]. The CASSIE tool [ 19 ] firstly applies the NDWI to a multispectral image, and then, in case the histogram of NDWI is not bimodal, it applies a multi-threshold Otsu algorithm [ 25 ] to separate land, water, and an intertidal zone, followed by a fixed threshold to distinguish the combined intertidal and water zones from the land. However, the use of the NDWI can fail in the presence of foam caused by wave action [ 26 ]. To address this issue, the CoastSat tool [ 17 ] computes the Modified Normalized Difference Water Index (MNDWI), incorporating the short-wave infrared (SWIR − 1) band, which is less impacted by this phenomenon. Then, it performs a supervised classification to differentiate between land, water, and whitewater (foam). Finally, Otsu thresholding [ 25 ] is applied to differentiate between land and water but only using the histogram of the pixels classified as land and water, respectively. The recently published SAET tool [ 20 ] is designed to use various indices, such as the AWEI, to avoid the effects of shadows and/or other dark surfaces, and MNDWI. It also allows using several methods for land/water masking, including K-means unsupervised classification and the Otsu method (both bimodal and multimodal). Additionally, it leverages the previous SHOREX tool’s [ 18 ] advantages during the shoreline extraction phase by also analyzing the gradient of the pixel value in the SWIR − 1 band. Pucino et al. [27] explored various previously mentioned spectral indices, also incorporating the Water Index Remote Sens. 2024,16, 4553 4 of 26 (WI), which showed enhanced capability in delineating shorelines under various challenging conditions. In addition to traditional thresholding methods, they also employed a Convolutional Neural Network (CNN) Unet+++ with Deep Supervision architecture. A common conclusion of many of the previous studies is that no overall shoreline detection solution exists, so the methodologies must be adapted and tested at each type of coast [ 22 , 27 – 29 ]. Notice also that existing software is designed to monitor large-scale shoreline changes whilst our goal is to use it for detecting the fast and small-scale shoreline changes related to megacusps. Finally, the only existing open-source software that can handle both S2 and PLN images is CoastSat (in a recent development). However, this tool is based on only one index and a thresholding-based method, and might not be the best combination for our sites. Therefore, we developed an in-house fully automatic shoreline detection algorithm that can be used in both S2 and PLN imagery and that is versatile (i.e., able to apply multiple combinations of indices, classification methods, and shoreline smoothing levels), in order to maximize shoreline detection accuracy before megacusp characterization. For this, several spectral indices were derived from the combination of bands distributed by the two satellite platforms. Various automatic thresholding methods, both global and local, as well as two unsupervised classification methods, deterministic Kmeans and the probabilistic Gaussian Mixture Model (GMM), were also applied. Moreover, the CoastSat tool was also applied to all satellite images to have a reference of the results that would be obtained from a benchmark methodology. 3. Materials and Methods 3.1. Study Sites This research was conducted on three sandy coastal areas within the Mediterranean Sea: one in Spain, an area located in the southern half of the Llobregat River Delta (SLD), in front of the Castelldefels and Gavà cities (Figure 1), and two in Italy, the embayed Feniglia (FNG) beach and a sandy area covering the northern section of the Ombrone River Delta (NOD), facing the town of Marina di Grosseto. These sandy coasts share some characteristics and differ in others. SLD and FNG beaches were selected as the main study sites, involving also megacusp identification and analysis, because they routinely host megacusp phenomena [ 8 , 30 ]. The NOD was used just as an additional site to evaluate the performance of the proposed shoreline extraction methodology. The wave and tidal conditions are similar in the three sites, as they all belong to the Western Mediterranean basin. They exhibit weak tidal fluctuations of the order of 20 cm. Long periods of calm waves are alternated with short but strong storms with limited fetches. The study area of the Southern Llobregat Delta (SLD) includes a 10 km portion of this sector, facing the towns of Castelldefels and Gavà (Figure 1a). It can be considered a semi-anthropized coast, with variable beach widths and a hard structure behind it. The beach is mainly composed of sand with a median grain size of 270 µ m. Waves generally come from two dominant directions (east–southeast and south–southwest) with average mean wave directions of 100° and 176° (to the north), respectively [31]. Feniglia (FNG) beach (Figure 1b) is located in one of the two spits enclosing the Orbetello lagoon and has a length of 6 km. It consists of a fully natural system of beach/dune ridges, rooted landwards at the Ansedonia promontory. The beach material is fine sand but there is no information about the grain size. FNG is directly exposed to waves from the southern directions [ 32 ]. While those from the SE and S tend to affect the entire coastline, waves from the SW predominantly impact the eastern stretch of the coast due to the obstacle posed by Monte Argentario. The third study area, the Northern Ombrone Delta (NOD, Figure 1b), starts at the port of Marina di Grosseto and extends for about 5 km towards the mouth of the Ombrone River. The grain size is 700 µ m and the prevailing wave direction in this area is from the south [32]. Remote Sens. 2024,16, 4553 5 of 26 Figure 1. Study areas: (a) Southern Llobregat Delta (SLD) coast (Spain), (b) Northern Ombrone Delta (NOD) coast (Italy) and Feniglia beach (FNG) (Italy). The position of wave buoys and tide gauges is also shown. The coordinate reference system is WGS84. 3.2. Shoreline Reference Data To validate the analysis, manually digitized instantaneous shorelines from orthomosaics of the three study sites were used. Products of FNG and NOD beaches were downloaded from the Tuscany region service, GEOscopio WMS. The two orthomosaics available during the operational period of the S2 and PLN missions, dating back to 2019 and 2021, were selected (Table 1). These open access products are distributed in four bands: NIR, red, green, and blue, with a resolution of 20 cm per pixel. FNG and NOD beaches fall into the same flight, and thus the dates coincide. Regarding the SLD coast, orthophotos from the Institut Cartogràfic i Geològic de Catalunya from 2017, 2019, 2020, and 2021 were used (Table 1). The products are distributed across three bands, with a spatial resolution of 25 cm. The manual digitization of the instantaneous shoreline was carried out using Remote Sens. 2024,16, 4553 6 of 26 QGIS 3.22 software, and the points indicating the shoreline position were put in the middle of the swash zone, so in the (approximate) intermediate position between wave runup and rundown. Moreover, GNSS surveys were conducted on a 1.2–1.4 km portion of the SLD coast on four dates (between 2017 and 2018) using an Ashtech Pro.Mark2 system from Thales Navigation, El Camino Real, Santa Clara, CA, USA (Table 1). In addition, a survey of FNG beach conducted in April 2023 with an Emlid Reach RS2 (multifrequency GNSS) was also included (Table 1). Shoreline position was tracked on the two beaches using the same Differential Global Navigation Satellite System (dGNSS). The base receiver was established as a permanent station in a known location, whereas the second was carried by an operator in RTK (Real-Time Kinematic) mode, recording a point every second in the SLD and every 0.2 s in FNG. The operators were experienced researchers who systematically walked in the middle of the swash zone, to make the shorelines consistent with the orthomosaic-derived shorelines. The planimetric and elevation errors associated with these survey campaigns were both on the order of 10 cm. Table 1. Acquisition date of the reference products, aerial orthophotos (*) or GNSS (**), used to validate satellite-derived shorelines, and the corresponding Sentinel-2 and PlanetScope imagery dates and relative time gaps in the three sites (SLD, FNG, and NOD, see the Abbreviations list). Some orthomosaics in the SLD are also used for the validation of the methodology for megacusp characterization (°). Reference Shoreline Date S2 Imagery Date Time Gap (Days) PLN Imagery Date Time Gap (Days) SLD 25 May 2017 *° 23 May 2017 2 27 May 2017 2 31 May 2017 ** 2 June 2017 2 31 May 2017 0 20 November 2017 ** 19 November 2017 1 20 November 2017 0 18 January 2018 ** 18 January 2018 0 18 January 2018 0 14 March 2018 ** 14 March 2018 0 16 March 2018 2 23 May 2019 *° 23 May 2019 0 23 May 2019 0 9 June 2021 *° 11 June 2021 2 8 June 2021 1 FNG 19 July 2019 * 23 July 2019 4 19 July 2019 0 20 July 2021 * 22 July 2021 2 20 July 2021 0 5 April 2023 ** 3 April 2023 2 6 April 2023 3 NOD 19 July 2019 * 23 July 2019 4 19 July 2019 0 20 July 2021 * 22 July 2021 2 20 July 2021 0 3.3. Satellite Images 3.3.1. Image Sources Data from Sentinel-2 (S2) and PlanetScope (PLN) platforms were used in this work. S2 is a European mission from the Copernicus program. This mission encompasses a pair of wide-swath, high-resolution, multispectral imaging satellites, orbiting in tandem along the same path but staggered 180° apart, aiming to achieve a rapid revisit rate of 5 days at the Equator. Equipped with Multispectral Optical Instruments (MSIs), S2 captures imagery across thirteen spectral bands: four bands with a spatial resolution of 10 m, six bands at 20 m, and three at 60 m. The coverage width of its orbital swath extends to 290 km [ 33 ]. S2 images are provided in different levels of processing. In this study, L2A (orthorectified and atmospherically corrected) images were used. The PLN constellation captures visible to near-infrared (NIR) surface reflectance with a spatial resolution of 3 m, achieving a global median average revisit time of 30.3 h. The technology behind the PLN satellites has developed across three generations: Dove-Classic Remote Sens. 2024,16, 4553 7 of 26 (PlanetScope-0), Dove-R (PlanetScope-1), and SuperDove (PlanetScope-2), each featuring unique specifications. The PlanetScope-0 and PlanetScope-1 sensors capture imagery across four spectral bands: blue, green, red, and near-infrared (NIR). The PlanetScope-2 sensor, meanwhile, offers bands with characteristics akin to those of the S2 satellite, also including the red-edge band [ 14 ]. Imagery from all of these satellites was used in our research, depending on the date of the validation data. The megacusp events occurred in the latest available generation of PLN products. PLN imagery is not freely available, but it can be accessed for free for research and academic purposes by applying to the Planet Education and Research Program. In this study, satellite images served multiple purposes. In the validation phase of the satellite-derived shorelines, the satellite images with acquisition times closest to those of the reference shorelines in the three beaches were selected (Table 1). Some SLD reference shorelines were also used to assess the feasibility of detecting megacusps, due to the presence of these phenomena in some of the high-resolution orthomosaics close to a satellite image date. Subsequently, tens of images for studying two megacusp events, one at the SLD coast and another at FNG beach, were also downloaded. The data set for the event on FNG beach comprised nineteen S2 L2A images and sixteen PLN images, between February and June 2022. The event on the SLD coast was tracked with twenty S2 images and twenty PLN images, between March and October 2023. 3.3.2. Satellite Image Preprocessing S2 products were downloaded using the GEE platform. Initially, the three areas of interest were selected, and a filter was applied to limit the number of clouds present in the images (0.60). Then, individual bands of interest were downloaded at different resolutions: 10 m for the blue (B02), green (B03), red (B04) and Near-Infrared 1 (NIR 1, B08) bands and 20 m for the Red-Edge 1 (B05), Short-Wave Infrared 1 (SWIR 1, B11) and ShortWave Infrared 2 (SWIR 2, B12) bands. The bands were subsequently resampled to a final resolution of 10 m, employing a nearest neighbor approach. Additionally, a median filter with a 3 × 3 pixel kernel size was applied to minimize image noise while still maintaining edge detection integrity. Finally, several spectral indices resulting from the combination of the resampled S2 bands were calculated (Table 2): the Normalized Difference Water Index (NDWI) (example provided in Figure 2a), Modified Normalized Difference Water Index (MNDWI), Sentinel Water Index (SWI), Water Index (WI), Sentinel Water Map (SWM), and Automated Water Extraction Index (AWEI) in its two formulations suitable for situations with thin clouds (AWEIsh) and without (AWEInosh). The PLN data were downloaded through the dedicated portal, where it is possible to filter products by cloud coverage and select the area of interest. All the provided bands have a resolution of 3 m, so resampling was unnecessary. The limited number of available bands allowed only the calculation of the NDWI index and the use of the NIR band (Table 2). 3.4. Offshore Wave and Sea Level Data The position of the shoreline has a strong intrinsic variability that depends on tide and wave conditions. The latter produce shoreline movements toward land or sea, due to various phenomena such as wave setup or storm surge. The analyzed coastal stretches (Figure 1), although in a micro-tidal Mediterranean environment, share a low beach slope. Moreover, there is sometimes a time gap between the reference data and the satellite images (Table 1). Therefore, a quantitative evaluation of the potential role of changes in wave and tide conditions is carried out in Section 5. Remote Sens. 2024,16, 4553 8 of 26 Table 2. Formulations of the set of indices derived from different combinations of spectral bands of the Sentinel-2 and PlanetScope platforms (see the Abbreviations list). Sentinel-2 Index Formula NDWI GREEN −NIR GREEN +NIR MNDWI GREEN −SWIR1 (GREEN +SWIR1 SWM BLUE +GREEN (NIR +SWIR1 SWI RED EDGE −SWIR1 RED EDGE +SWIR1 WI (1.7204 +171 ×GREEN +3×RED −70 ×(NIR −45 ×SWIR1 −71 × SWIR2) AWEINOsh (4×(BLUE −SWIR1)−(0.25 ×NIR) + 2.75 ×SWIR2) AWEIsh (BLUE +2.5 ×GREEN −1.5 ×(NIR +SWIR1)−0.25 ×SWIR2) PlanetScope Index Formula NDWI GREEN −NIR GREEN +NIR NIR NIR Figure 2. Exampleof the shoreline extraction method with S2 data: (a) raster file of the spectral index (NDWI), (b) binarization of the image with K-means, (c) contour extraction, (d) comparison between the reference shoreline and the detected one, and (e) validation by using the baseline and transect method. Remote Sens. 2024,16, 4553 9 of 26 The sea level data ztide , were obtained from the tidal gauges of Barcelona harbor (Puertos del Estado, Spain) for the SLD coast, Civitavecchia harbor (Italy) for FNG beach, and Marina di Campo harbor for the NOD coast. The wave conditions were obtained from the Barcelona Buoy II (Puertos del Estado) for the SLD coast, Giannutri buoy (Servizio Idrologico Regionale, Tuscany region, Italy) for FNG beach, and Castiglione della Pescaia buoy (Servizio Idrologico Regionale, Tuscany region) for the NOD coastal stretch. To quantify the setup, the formulation proposed by Stockdon et al. [34] was used: zsetup =0.35 ·tan(β)·pHs·L0(1) where Hs is the deep-water significant wave height, L0=g T2/( 2 π) is the deep-water wavelength (with g being gravity and T the wave period), and tan β is the mean bed slope of the swash zone. Finally, to evaluate the potential contribution to the advancement or retreat of the shoreline by changes in setup and mean sea level, the difference between the z -values obtained at the times of the satellite image and the reference data was divided by tan β. Available LiDAR-based Digital Surface Models (DSMs) were used to compute the beach swash slope β . A DSM with 1 m of spatial resolution, provided by the Tuscany Region, was used to gather the elevation data and assess the slope of the NOD and FNG beaches, matching the same frames as the corresponding orthomosaics. Instead, for SLD beach, LiDAR data were downloaded from the Geoportal de Cartografia of the Area Metropolitana de Barcelona (AMB), which provides a grid with 1 m spatial resolution. The obtained average swash slopes were 0.08 in the SLD, 0.07 in FNG, and 0.06 in the NOD. 3.5. Shoreline Detection Algorithms and Validation The core of the methodology was developed in two separate tools (programmed in Python), one for shoreline extraction and the other one for megacusp characterization (Figure 3). The shoreline extraction procedure was constituted by the selection of four methods to automatically separate the pixels into land and sea classes using the spectral index raster. Two methods were based on automatic threshold definition and the other two on machine learning, in particular, on unsupervised clustering. The procedures described below were repeated for all the considered spectral indices (Table 2). The first thresholding method was the Otsu one [ 25 ]. By analyzing the input raster histogram, the Otsu method determines a threshold that maximizes the interclass variance between water and non-water pixels, enabling the binarization of the input image. An adaptive thresholding technique called Niblack [ 35 ] was also tested. This method calculates a threshold value for each pixel in an image, based on the mean and standard deviation of intensity values in the surrounding pixels. This allows the threshold to be locally adapted according to variations in brightness and contrast around each pixel. Considering the unsupervised clustering methods, K-means and the GMM were employed. The K-means algorithm [ 36 ] employs a deterministic approach to divide the image into a set of classes (clusters) based on the proximity of pixel values. Then, the pixels of one cluster are labeled as water pixels while the remaining ones are labeled as nonwater pixels (example in Figure 2b). Conversely, the GMM [ 37 ] operates on a probabilistic basis. For image clustering, it estimates the parameters of Gaussian distributions that represent different pixel intensities. Each pixel is then assigned a probability of belonging to each distribution, effectively weighting its classification. The image is segmented into classes based on these probabilities, using an iterative process similar to K-means for refinement. These methodologies offer a significant advantage in the binarization of raster files compared to conventional thresholding techniques by adapting to the intricate and varied distributions of pixel intensities and providing a more robust approach to image segmentation [ 38 ]. In our study, the number of clusters was set to 2 for both classification methods. Remote Sens. 2024,16, 4553 16 of 26 Figure 5. First segment of May 2017 used for the validation phase. (a) On top, the four lines correspond to the shorelines of the reference case (green), the best method–index combination in S2 data (GMM–NDWI, orange), the best method–index combination in PLN data (K-means–NIR, grey), and the CoastSat tool (blue). (b) At the bottom, with the same color, the detrended lines show the automatic peaks (red square) and valleys (green dot) for each detected megacusp. The numbers refer to the megacusp embayments that are visible in the orthomosaic. The x-axis is set to zero at the starting point of the segment. Figure 6. Second segment of May 2019 used for the validation phase. (a) On top, the four lines present the shorelines of the reference case (green), the best method–index combination in S2 data (GMM–NDWI, orange), the best method–index combination in PLN data (K-means–NIR, grey), and the CoastSat tool (blue). (b) On the bottom, with the same color, the detrended lines show the automatically detected peaks (red square) and valleys (green dot) for each megacusp. The numbers enumerate the megacusp embayments that are visible in the orthomosaic. The x-axis is set to zero at the starting point of the segment. 4.3. Monitoring of Megacusp Development Events The implemented methodology was finally applied to monitor the time evolution of two megacusp events, in the SLD in 2023 and in FNG in 2022. Given the successful Remote Sens. 2024,16, 4553 17 of 26 identification of megacusps using both satellite sources, S2 data were used for a general analysis of the events and PLN data were added to enrich the analysis during peak periods of these events. The combination of methods used to detect the shorelines was the same as in the megacusp validation phase (Section 4.2). The images with clouds were discarded, as well as those with high waves ( Hs> 2 m) to avoid beach inundation moments that hid megacusps. The 2023 megacusp event at the SLD coast spanned a significant portion of the beach length. Two contiguous segments were analyzed separately, due to their different orientations relative to the north. To characterize the event, a total of twenty S2 images between March and October 2023 and twenty PLN images between May and June 2023 were used. Figure 7shows key moments in the temporal evolution of the SLD event from March 2023 to October 2023 in the first segment. To numerically characterize the megacusps, all the parameters described in Section 4.2 were calculated for each date, as shown in Figure 8. The wave conditions (significant height, peak period, and mean direction) over the study period are also shown. Figure 7. Time series of the megacusp event in the SLD coast in 2023. On the left, Sentinel-2 images in the period between March and October 2023. On the right, the time series enriched by adding PlanetScope images during the peak of the event (May–June 2023). Megacusps of 8–10 m amplitude and 150–200 m wavelength (with s = 1.0075–1.010 and σs = 3–4 m) were present during the whole study period (Figure 8). At the end of May 2023, their amplitude increased up to 15 m (without modifying the wavelength), with s> 1.010 and σs > 4 m, after a 2-day period of 1 m high waves coming from a rather constant SEE direction. The amplitude decreased again in the subsequent weeks of smaller wave height and variable directions. At the end of June 2023, a new period of Hs slightly higher than 1 m and constant SEE direction induced a new increase in amplitude up to 15 m, this time with an increase in wavelength up to 400 m (Figure 8). The 2022 FNG beach megacusp event spanned from February to June. For this study, nineteen S2 images and sixteen PLN images (March 2022) were used. Megacusps primarily formed in the central part of FNG beach, so the analysis focused on a 2 km section in this area. Figure 9showcases the highlights of the event, providing a visual representation of the megacusps. The obtained quantitative megacusp characteristics (Figure 10) show that small amplitude features were already detectable at the beginning of the study period and remained relatively constant from February to the middle of March ( 1.010 <s<1.020, σs∼= 4–5 m , a∼ 10–15 m, and λ∼ 150–250 m). At the end of March, the megacusp amplitude increased up to 20 m, with a more regular wavelength of 160–220 m, during a series of minor storms ( Hs up to 2 m) from the S. At the beginning of April, a couple of Remote Sens. 2024,16, 4553 18 of 26 storms with Hs up to 3–4 m coming from the SW had a destructive effect on the megacusps, again diminishing the amplitude to 10 m with the same wavelength. The megacusps maintained these characteristics until the end of the study period, slowly diminishing the amplitude. Figure 8. Results of the 2023 megacusp event in the SLD coast with the corresponding wave conditions. Time series of (a) significant wave height ( Hs ), (b) peak wave period ( Tp ), (c) direction of wave incidence with respect to the north ( θ ), (d) sinuosity ( s ), (e) shoreline standard deviation ( σs ), (f) mean megacusp amplitude (a), and (g) mean megacusp wavelength (λ) are shown. Remote Sens. 2024,16, 4553 19 of 26 Figure 9. Time series of the megacusp event in FNG beach in 2022. On the left, Sentinel-2 images in the period between February and June 2022. On the right, the time series enriched by adding PlanetScope images at the peak of the event (March 2022). Figure 10. Cont. Remote Sens. 2024,16, 4553 20 of 26 Figure 10. Results of the 2022 megacusp event in FNG beach with the corresponding wave conditions. Time series of (a) significant wave height ( Hs ), (b) peak wave period ( Tp ), (c) direction of wave incidence to north ( θ ), (d) sinuosity ( s ), (e) shoreline standard deviation ( σs ), (f) mean megacusp amplitude (a), and (g) mean megacusp wavelength (λ) are shown. 5. Discussion 5.1. Shoreline Detection Tool In this study, various combinations of spectral indices and shoreline extraction algorithms were compared to detect the shoreline in S2 and PLN images, using in situ GNSS RTK measurements and high-resolution orthomosaics as validation data (Table 1). To prevent misinterpretations of the results, various metrics were utilized, with an emphasis on MAD and RMSD , which highlight the error in absolute terms. The best performance in terms of MAD when using S2 images (10 m of spatial resolution) was about 4 m at the SLD coast and FNG beach and 5 m at the NOD coast (Table 3). For the PLN data (3 m of spatial resolution), the best performance was a MAD of about 1.5 m at the SLD coast and of about 2 m at FGN beach and the NOD coast (Table 4). These results are consistent with those reported in previous benchmark studies for medium-resolution satellites [ 41 ] and in previous PLN results [ 42 ]. The two satellite-derived shorelines exhibited subpixel accuracy, but in PLN this was slightly less pronounced than in S2 data, in line with Bishop-Taylor et al. [28]. Previous studies have demonstrated that no single spectral index or methodology performs optimally in all circumstances [ 22 , 27 – 29 ]. This was also observed in our results, both for S2 and PLN data. However, certain trends could be identified. The WI index proved to be highly effective in detecting shorelines in SLD and FNG S2 images (Table 3). The benefits of using this index were previously identified by Pucino et al. [27] . PlanetScope images in these sites showed maximum accuracy in shoreline detection using the NIR channel directly. The NOD results needed a separate analysis. The stretch of the coast close to the delta suffers from high erosion rates, with the beach almost disappearing and vegetation approaching the water. In such a dynamic scenario, spectral indices that are less sensitive to changes in scene radiance such as the AWEI [ 28 , 43 ] obtained better results. Finally, the best index for PlanetScope in the NOD was the NDWI. Regarding the methodology, combinations of index–method that involve using two unsupervised classification methods, K-means and GMM, generally performed better than combinations that employed thresholding. Unsupervised classification eliminated the need for a training phase, which is required, for example, in the partially supervised classification employed by CoastSat [ 17 ]. This unsupervised approach also minimized the operator intervention. 5.2. Megacusp Monitoring Tool The initial identification of segments with megacusps in Section 3.6 and the identification of the two events analyzed in Section 4.3 took place through visual analysis by experienced users, as in previous studies [ 8 ]. The analysis of the results (Tables 5and 6) indicated the possibility of automatically identifying the presence of megacusps by using a threshold in the sinuosity parameter of s= 1.010. In all the segments of the validation phase, this proposed sinuosity threshold correctly identified segments with and without Remote Sens. 2024,16, 4553 21 of 26 megacusps. On the other hand, the standard deviation parameter σs showed greater variability, and it was impossible to identify a clear threshold, a value of 2 m working well in most but not all the cases. Note also that the obtained values for average wavelength and amplitude were coherent with those measured using a video system by a previous study in the same site but for different years [8]. Regarding the event-scale analysis of the megacusp evolution presented in Section 4.3, the sinuosity parameter performed well in capturing the development of the event in the SLD, but the threshold established during the validation phase seems slightly too selective in comparison with visual analysis. The coastline exhibits undulations even with s< 1.010. Conversely, at FNG beach, s= 1.010 effectively allowed us to identify the start and end of the event. As for the standard deviation parameter σs , the values followed the same trend as s for both sites, but the potential threshold of 2 m established during validation was almost always exceeded during the periods under examination. This event-scale analysis could represent a significant step forward for coastal planning. The embayments created during the development of large-amplitude megacusps constitute erosive hotspots. Traditionally, the long-term monitoring of megacusp phenomena has been conducted through video monitoring systems [ 8 ], which require operator maintenance and have a limited field of view (approximately 1 km). LiDAR, GNSS, and drone surveys have also been used [ 12 , 44 ], but the required repeated measurement campaigns make them complicated to apply. Only one study [ 16 ] observed these phenomena from satellites, but it was based on visual analysis and it did not provide a numerical quantification of megacusp characteristics. Therefore, our methodology, by utilizing both medium-spatial-resolution (10 m) and medium-temporal-resolution (5 days) satellites like S2 and high spatial (3 m) and temporal (about 1 day) resolution satellites like PLN, expands the possibilities for monitoring these phenomena. Finally, a comparative analysis with the CoastSat tool [ 17 ] was also conducted in Section 3.6. Our results suggest that the CoastSat tool appears constrained in its ability to capture megacusps along all the shorelines analyzed. This is probably due to its method of applying thresholding to the index values of pixels near the shoreline. Specifically, CoastSat defines shorelines solely based on pixels with index values close to the threshold level. Consequently, this method utilizes only a subset of shoreline pixels for evaluation, incorporating a type of internal “alongshore smoothing”. This approach turned out to limit its capability to detect alongshore oscillations at the spatial scale of megacusps. 5.3. Error Sources The main sources of errors affecting shoreline extraction algorithms include georeferencing and the resolution of satellite images, and issues related to water level variability due to wave and tide action. Regarding georeferencing, an accuracy of 11 m has been reported for S2 images until 2021, and an improvement up to 6 m for more recent images [ 45 ]. Regarding PLN data, the geolocation accuracy varies between 8.5 and 11.7 m up to 2021 [ 46 ] and between 3.7 and 7.6 m from then on [ 47 ]. Without prior information on the direction of this error, its impact on the cross-shore direction and thereby on the total MAD and Bias remains unknown, but this could be an explanation for part of the (small) remaining errors in shoreline detection of the present study. In fact, among the tools available for shoreline extraction, only SHOREX provides image coregistration, seeking to enhance the absolute geolocation accuracy by fitting all images to a high-resolution orthomosaic. In a comparison study, Vos et al. [41] indeed reported an enhanced accuracy due to SHOREX image coregistration, improving the accuracy of shoreline time series derived from S2. The resolution of satellite images affects the size of objects recognizable in the scene. Many previous studies using medium-resolution images (10 m) have already demonstrated that it is possible to achieve subpixel resolutions [ 19 , 20 , 27 ]. The proliferation of highresolution satellite images (3 m) has not yet allowed for achieving the same level of subpixel accuracy as medium-resolution images [ 28 ], as occurred in our study. This discrepancy is likely due to the increased contribution of other error sources. Remote Sens. 2024,16, 4553 22 of 26 To perform accurate shoreline detection, even in micro-tidal environments, there is a need to combine satellite-derived shorelines with in situ wave and tide data at the time of image acquisition. There is no one-size-fits-all solution and more research is needed to identify how to optimally apply wave corrections across different coastal environments and beach morphologies [ 41 ]. In our study, the wave and tide effects were checked ex post. For both contributions, the differences between the water level values at the time of the satellite images and those at the time of the references (Table 1) were first evaluated. Subsequently, the apparent total shoreline distances due to waves and tides were calculated by considering the swash slope in each site (details in Section 3.4). The results showed maximum apparent distance values of approximately 2 m in only two cases, the rest being smaller than 1.5 m (Table 7). The fact that the sites have small tidal ranges ( < 0.5 m) and low-energy waves most of the time explains why their effect plays a minor role, being smaller than the average deviations (e.g., MAD in Tables 3and 4). However, they are among the factors contributing to the obtained errors. Table 7. Planimetric shoreline displacements due to both tidal and wave setup contributions. The apparent total distances are between the S2and PLN-derived shorelines and the reference ones. The positions of the tide gauges/buoys used for the analyses are reported in Section 3.4 for the three coastal areas. The difference in days between PLN and S2 and the reference are shown in Table 1. Reference Shoreline Date Total Distance (m) S2–Reference Total Distance (m) PLN–Reference SLD 25 May 2017 0.4 1.4 31 May 2017 −1.2 −0.5 20 November 2017 0.5 0.2 18 January 2018 −0.8 −0.2 14 March 2018 −0.2 1.9 23 May 2019 0.4 0.5 9 June 2021 0.1 0.3 FNG 19 July 2019 −2.2 0.1 20 July 2021 0.6 0.6 5 April 2023 1.5 1.3 NOD 19 July 2019 0.3 1.5 20 July 2021 1.0 0.2 It is important to highlight that the megacusp monitoring tool is not affected significantly by any of the potential sources of error mentioned above for shoreline detection. The georeferencing errors do not affect megacusp characteristics because, at this spatial scale, they only imply a positional shift and not a distortion. The same applies to the effect of tides and waves, which produce changes only in the direction orthogonal to the shoreline, so without any impact on the undulation parameters. Finally, it is noticeable that the increase in the resolution, from 10 m of S2 to 3 m for PLN, does not have a major impact on the accuracy in megacusp characterization. 5.4. Strengths and Limitations of the Applied Methodologies The shoreline extraction methodology developed has provided excellent results, in line with the most recent tools [ 20 , 41 ]. However, it is important to highlight some weaknesses and other areas that have not yet been explored. First, this study considered micro-tidal sandy areas, the scenario where shoreline extraction tools perform best. In meso-tidal environments, there can be more difficulties due to the identification of the transition area [ 19 ] and the large number of indices offered by our tool could partially solve this problem. Remote Sens. 2024,16, 4553 23 of 26 Other issues arose in the analysis of images with wave conditions that created foam. This problem can be solved in the way employed by the CoastSat tool [ 17 ], i.e., by performing a previous supervised selection. Depending on the objectives and the level of automation desired, it is possible to set a wave height threshold and avoid analyzing specific images. This would be an additional restriction, similar to what already occurs during the preprocessing phase of satellite images when a threshold for cloud cover in the image is applied. However, this would come at the cost of reducing the number of images available for analysis. The CASSIE tool [ 19 ] attempts to eliminate the effect of the intertidal zone by performing a multi-threshold Otsu and excluding this area in the final division between land and water. In the present work, an attempt to mitigate the effect of foam is made by exploring the use of unsupervised classification techniques. However, a remaining open issue is applying an automatic methodology to determine the optimal number of clusters into which the image should be divided. Further improvements could include incorporating tidal correction and wave motion adjustments into the workflow for the extracted shorelines. Several studies [ 27 ] have investigated the differences in using offshore buoys or modeling, and others have successfully integrated them into their tools [ 41 ]. Similarly, the inclusion of image coregistration, which can be implemented through high-resolution orthomosaics [ 18 ], could enhance the accuracy and usability of the results. However, these latter potential corrections do not affect megacusp characterization, so they were not considered necessary in the present study. Regarding the megacusp tool, the results obtained suggest the difficulty in finding fixed thresholds in the parameters to detect megacusp presence. Different coastal contexts (FNG and SLD) in terms of orientation, wave exposure, and coastal morphology lead to slightly different success rates when applying the threshold values selected during the validation phase. Future developments of this work should involve attempting to find a more robust indicator of megacusp presence suitable for different scenarios and events. 6. Conclusions In this study, a methodology for shoreline detection in Sentinel-2 and PlanetScope satellite images was developed, achieving a good global accuracy compared with reference shorelines, analogous to the current state-of-the-art methods. It was found that a mean absolute error of approximately half a pixel could be consistently achieved in the three tested Mediterranean micro-tidal coasts. This detection accuracy was achievable with both types of satellite images, requiring only minor adjustments in the preprocessing steps. Specifically, mean absolute distances of 4–5 m for S2 and 1.5–2 m for PLN were achieved. The research revealed that the efficacy of index–method combinations changed across the different beach environments and occasionally varied from day to day on the same beach. This fact highlighted the dynamic nature of coastal environments and the requirement for adaptable monitoring techniques. The shorelines extracted with these methodologies were the starting point for developing a semi-automatic approach for identifying and characterizing megacusps, which are shoreline undulations with wavelengths of a few hundred meters and amplitudes of several meters. This novel approach achieved excellent results in terms of megacusp characterization in the satellite images, compared with those in the reference shorelines. Including the ability to modify the level of smoothing, the shoreline extraction method was found to better catch the alongshore variations as opposed to other tools such as CoastSat, which tends to smooth the alongshore undulations. Notably, the methodology was able to determine megacusp presence (using sinuosity as an indicator) and capture the spatial characteristics of the megacusps (wavelength and amplitude). Their temporal evolution, at scales ranging from days to weeks, was also captured, enabling the possibility of establishing correlations with wave forcing and enhancing the understanding of megacusp dynamics. To achieve this level of temporal resolution, the integration of S2 images with PLN images proved essential. Remote Sens. 2024,16, 4553 24 of 26 The use of satellite imagery in the monitoring of megacusps demonstrated the potential for characterizing large geographical areas, presenting an advantage over traditional video monitoring systems (which otherwise remain fundamental in ensuring a higher temporal resolution and avoiding the issue of cloudiness during storms). This study emphasized the importance of using shoreline extraction tools capable of monitoring not just long-term erosion but also rapid changes in beach width, such as those induced by megacusps. This capability is deemed critical for the advancement of beach management strategies. Author Contributions: Conceptualization, R.A., E.A., G.L., A.M. and F.R.; methodology, R.A., E.A. and F.R.; software, R.A.; validation, R.A., E.A., G.L., A.M. and F.R.; formal analysis, R.A., E.A., G.L., A.M. and F.R.; investigation and resources, R.A., E.A., G.L., A.M., G.S. and F.R.; data curation, R.A.; writing—original draft preparation, R.A., E.A. and F.R.; writing—review and editing, R.A., E.A., G.L., A.M. and F.R.; visualization, R.A., E.A. and A.M.; supervision, R.A., E.A., G.L., A.M., G.S. and F.R.; project administration, G.L., A.M. and F.R.; funding acquisition, F.R. All authors have read and agreed to the published version of the manuscript. Funding: The work of F. Ribas was funded by the Spanish Ministry of Science, Innovation and Universities—National Research Agency and EU “NextGenerationEU/PRTR” grant numbers PID2021124272OB-C22 (MOLLY-MOD) and TED2021-130321B-I00 (SOLDEMOR). Data Availability Statement: The codes and data supporting all results shown in the manuscript are available from the corresponding authors upon request. Acknowledgments: The authors are grateful to Planet for access to and use of PlanetScope images through the educational and research standard plan (id. 736919). Conflicts of Interest: The author declares no conflicts of interest. Abbreviations The following abbreviations are used in this manuscript: AWEI Automated Water Extraction Index CNN Convolutional Neural Network dGNSS Differential Global Navigation Satellite System DSM Digital Surface Model FNG Feniglia beach GEE Google Earth Engine GMM Gaussian Mixture Model GNSS Global Navigation Satellite System MAD Mean absolute deviation ML Machine Learning MNDWI Modified Normalized Water Index MSI Multispectral Instrument NDWI Normalized Difference Water Index NIR Near-infrared NOD Northern Ombrone Delta beach PLN PlanetScope RMSD Root Mean Square Deviation RTK Real-Time Kinematic SLD Southern Llobregat Delta beach SWI Sentinel Water Index SWIR Short-Wave Infrared SWM Sentinel Water Map S2 Sentinel-2 WI Water Index Remote Sens. 2024,16, 4553 25 of 26 References 1. Luijendijk, A.P.; Hagenaars, G.; Ranasinghe, R.; Baart, F.; Donchyts, G.; Aarninkhof, S. The State of the World’s Beaches. Sci. Rep. 2018,8, 6641. [CrossRef] [PubMed] 2. Syvitski, J.; Anthony, E.; Saito, Y.; Z˘ainescu, F.; Day, J.; Bhattacharya, J.P.; Giosan, L. Large deltas, small deltas: Toward a more rigorous understanding of coastal marine deltas. Global Planet. Change 2022,218, 103958. [CrossRef] 3. Vousdoukas, M.I.; Ranasinghe, R.; Mentaschi, L.; Plomaritis, T.A.; Athanasiou, P.; Luijendijk, A.; Feyen, L. Sandy coastlines under threat of erosion. Nat. Clim. Change 2020,10, 260–263. [CrossRef] 4. Boak, E.H.; Turner, I.L. Shoreline Definition and Detection: A Review. J. Coast. Res. 2005,2005, 688–703. [CrossRef] 5. Chen, C.; Fu, J.; Zhang, S.; Zhao, X. Coastline information extraction based on the tasseled cap transformation of Landsat-8 OLI images. Estuar. Coast. Shelf Sci. 2019,217, 281–291. [CrossRef] 6. Fogarin, S.; Zanetti, M.; Dal Barco, M.; Zennaro, F.; Furlan, E.; Torresan, S.; Pham, H.; Critto, A. Combining remote sensing analysis with machine learning to evaluate short-term coastal evolution trend in the shoreline of Venice. Sci. Total Environ. 2023, 859, 160293. [CrossRef] 7. Orzech, M.D.; Reniers, A.J.H.M.; Thornton, E.B.; MacMahan, J.H. Megacusps on rip channel bathymetry: Observations and modeling. Coast. Eng. 2011,58, 890–907. [CrossRef] 8. de Swart, R.L.; Ribas, F.; Calvete, D.; Simarro, G.; Guillén, J. Observations of megacusp dynamics and their coupling with crescentic bars at an open, fetch-limited beach. Earth Surf. Process. Landf. 2022,47, 3180–3198. [CrossRef] 9. Segura, L.; Hansen, J.; Lowe, R.; Symonds, G.; Contardo, S. Shoreline variability at a low-energy beach: Contributions of storms, megacusps and sea-breeze cycles. Mar. Geol. 2018,400, 94–106. [CrossRef] 10. Castelle, B.; Marieu, V.; Bujan, S.; Splinter, K.D.; Robinet, A.; Sénéchal, N.; Ferreira, S. Impact of the winter 2013–2014 series of severe Western Europe storms on a double-barred sandy coast: Beach and dune erosion and megacusp embayments. Geomorphology 2015,238, 135–148. [CrossRef] 11. van de Lageweg, W.; Bryan, K.; Coco, G.; Ruessink, B. Observations of shoreline–sandbar coupling on an embayed beach. Mar. Geol. 2013,344, 101–114. [CrossRef] 12. Matsumoto, H.; Young, A.P.; Guza, R.T. Cusp and Mega Cusp Observations on a Mixed Sediment Beach. Earth Space Sci. 2020, 7, e2020EA001366. [CrossRef] 13. Qiao, G.; Mi, H.; Wang, W.; Tong, X.; Li, Z.; Li, T.; Liu, S.; Hong, Y. 55-year (1960–2015) spatiotemporal shoreline change analysis using historical DISP and Landsat time series data in Shanghai. Int. J. Appl. Earth Obs. Geoinf. 2018,68, 238–251. [CrossRef] 14. Team, P. PlanetScope Product Specifications; Planet Team: San Francisco, CA, USA, 2023. Available online: https://assets.planet.com (accessed on 2 July 2024). 15. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017,202, 18–27. [CrossRef] 16. Leuci, R.; Wiles, E.; Thackeray, Z.; Vella, G. Trends in sandy beach variability EThekwini Municipality, South Africa. J. Sea Res. 2022,179, 102149. [CrossRef] 17. Vos, K.; Splinter, K.; Harley, M.; Simmons, J.; Turner, I. CoastSat: A Google Earth Engine-enabled Python toolkit to extract shorelines from publicly available satellite imagery. Environ. Model. Softw. 2019,122, 104528. [CrossRef] 18. Sánchez-García, E.; Palomar-Vázquez, J.; Pardo-Pascual, J.; Almonacid-Caballer, J.; Cabezas-Rabadán, C.; Gómez-Pujol, L. An efficient protocol for accurate and massive shoreline definition from mid-resolution satellite imagery. Coast. Eng. 2020,160, 103732. [CrossRef] 19. Almeida, L.P.; Efraim de Oliveira, I.; Lyra, R.; Scaranto Dazzi, R.L.; Martins, V.G.; Henrique da Fontoura Klein, A. Coastal Analyst System from Space Imagery Engine (CASSIE): Shoreline management module. Environ. Modell. Softw. 2021,140, 105033. [CrossRef] 20. Palomar-Vázquez, J.; Pardo-Pascual, J.E.; Almonacid-Caballer, J.; Cabezas-Rabadán, C. Shoreline Analysis and Extraction Tool (SAET): A New Tool for the Automatic Extraction of Satellite-Derived Shorelines with Subpixel Accuracy. Remote Sens. 2023, 15, 3198. [CrossRef] 21. Toure, S.; Diop, O.; Kpalma, K.; Maiga, A.S. Shoreline Detection using Optical Remote Sensing: A Review. ISPRS Int. J. Geo-Inf. 2019,8, 75. [CrossRef] 22. Kelly, J.T.; Gontz, A.M. Using GPS-surveyed intertidal zones to determine the validity of shorelines automatically mapped by Landsat water indices. Int. J. Appl. Earth Obs. Geoinf. 2018,65, 92–104. [CrossRef] 23. Tsiakos, C.A.D.; Chalkias, C. Use of Machine Learning and Remote Sensing Techniques for Shoreline Monitoring: A Review of Recent Literature. Appl. Sci. 2023,13, 3268. [CrossRef] 24. Apostolopoulos, D.; Nikolakopoulos, K. A review and meta-analysis of remote sensing data, GIS methods, materials and indices used for monitoring the coastline evolution over the last twenty years. Eur. J. Remote Sens 2021,54, 240–265. [CrossRef] 25. Liao, P.S.; Chen, T.S.; Chung, P.C. A fast algorithm for multilevel thresholding. J. Inf. Sci. Eng. 2001,17, 713–727. [CrossRef] 26. Pardo-Pascual, J.; Sánchez-García, E.; Almonacid-Caballer, J.; Palomar-Vázquez, J.; de los Santos, E.; Fernández-Sarría, A.; Balaguer-Beser, A. Assessing the accuracy of automatically extracted shorelines on microtidal beaches from Landsat 7, Landsat 8 and Sentinel-2 imagery. Remote Sens. 2018,10, 326. [CrossRef] 27. Pucino, N.; Kennedy, D.M.; Young, M.; Ierodiaconou, D. Assessing the accuracy of Sentinel-2 instantaneous subpixel shorelines using synchronous UAV ground truth surveys. Remote Sens. Environ. 2022,282, 113293. [CrossRef]