A Simple, Fully Automated Shoreline Detection Algorithm for High-Resolution Multi-Spectral Imagery

This paper develops and validates a new fully automated procedure for shoreline delineation from high-resolution multispectral satellite images. The model is based on a new water–land index, the Direct Difference Water Index (DDWI). A new technique based on the buffer overlay method is also presented to determine the shoreline changes from different satellite images and obtain a time series for the shoreline changes. The shoreline detection model was applied to imagery from multiple satellites and validated to have sub-pixel accuracy using beach survey data that were collected from the Lake Michigan (USA) shoreline using a novel backpack-based LiDAR system. The model was also applied to 132 satellite images of a Lake Michigan beach over a three-year period and detected the shoreline accurately, with a >99% success rate. The model out-performed other existing shoreline detection algorithms based on different water indices and clustering techniques. The resolution shoreline position timeseries is the first satellite image-extracted dataset of its kind in terms of its high spatial and temporal resolution, and paves the road to obtaining other high-temporal-resolution datasets to refine models of beaches worldwide.


Introduction
Coastal areas are some of the most densely populated areas in the world, with average population densities three times higher than the average world density [1]. Most of the world's megacities are located in the low-level coastal areas [2], and migration to coastal areas continues [2]. However, coastal areas are also increasingly exposed to many hazards that are enhanced by sea level rises, such as flooding, shoreline erosion, and coastal storms [3], all of which can cause dramatic changes to the shoreline position, in turn threatening infrastructure and communities. The ability to detect and quantify these changes is essential for the development and refinement of models that can predict such shoreline changes, in turn guiding the design and management of coastal areas that are more resilient to ever-increasing coastal hazards.
Remote-sensing approaches have long been utilized to determine shoreline positions, allowing for the accurate mapping of shorelines over large geographical areas. Shoreline mapping was first carried out using aircraft-derived black and white photographs, but many types of aerial imagery, ranging from black and white images to hyperspectral images, can be used to quantify the shoreline location once the images are georeferenced [4,5]. In addition to aerial imagery, other remote-sensing modalities, e.g., airborne LiDAR, can be used for shoreline detection. LiDAR technology, in terms of performance, size, and cost, has seen tremendous developments over the last two decades [6][7][8][9][10]. In related recent work, we presented a LiDAR-equipped unmanned aerial vehicle (UAV) that was utilized for high-resolution shoreline and beach topography quantification for Lake Michigan [8,9]. However, the drawback of UAV-and aircraft-derived imagery for use in shoreline mapping is their high instrumentation and logistical costs as well as their coverage limitation, which can preclude their application and limit the temporal frequency of image acquisition.
A key step in shoreline detection from multispectral imagery is the combination of multiple spectral bands to create an index that distinguishes land from water following the assignment of a threshold value. Numerous indices have been developed and refined, including the Normalized Difference Water Index (NDWI), Normalized Difference Vegetation Index (NDVI), Modified Normalized Difference Water Index (MNDWI), Coastal Water Index (CWI), Tasseled Cap Wetness (TCW), Coal Mine Index (CMI), and Automated Water Extraction Index (AWEI) Table 1. The main purpose of these indices is to increase the contrast between the water and the land pixels in the image. Furthermore, new mathematical methods, indices or a combination of several indices have been examined to increase the contrast between the water and the land to improve the shoreline detection [12,13,18,20]. The appropriate index for an application depends largely on the spectral bands that are available from the satellite imagery. Many indices were developed specifically for the available bands in the Landsat images (e.g., CWI, TWC, CMI, AWEI), which limits their applications. However, other indices, such as the NDWI and NDVI, can be applied to any satellite with RGB and near-infrared (NIR) bands, and that is why they are the most widely used indices. Many of these indices were developed as normalized indices in the hope of finding a constant thresholding value between the water and the land pixels. However, this turned out to be inapplicable or not good for accurate prediction of the shorelines [14,21,22]. Landsat, Sentinel [14,15] Coastal Water Index CW I = Max(SW IR2) * B Max(B) * SW IR2 Landsat [25] Tasseled Cap Wetness Landsat [12] Automated Water Extraction Index [19,20] * ρ1, ρ2, ρ3, . . . indicates the number of the band in Landsat image.
A second key step in shoreline detection from multispectral imagery is the application of an automated or semi-automated technique to delineate the shoreline from an image in the classifier index (NDWI, etc.). Unsupervised classification techniques (e.g., K-means, mean shift, and Otsu) are the most common methods for shoreline delineation and were used by many researchers [11,12,15,17,25]. However, since these classifiers are centroidbased, they are greatly affected by the outliers [26,27]. Therefore, the optimum thresholding values obtained from these methods do not have to be the best values to distinguish between the water and the land. This will be further illustrated in Sections 4 and 5. Another Remote Sens. 2022, 14, 557 3 of 18 equally common method is using a constant threshold after applying one or more of the water indices, or after applying some operations to the image's bands [16,19,23,24,28]. Nevertheless, it was found to be impossible to obtain the constant thresholding that is applicable to different beaches' imageries from different satellites, with different topologies, and seasons [14,21]. Some other models [15,29] use supervised classifiers (supported vector machine (SVM), maximum likelihood (ML), and neural networks (NN)). However, to make a robust supervised classifier that can work for different satellite missions and different beach topologies, many images from different satellites and environments are needed to train the classifier. Therefore, since the current efforts are motivated by the lack of data for shoreline positions, a comprehensive supervised shoreline delineator cannot be achieved.
Despite the extensive refinement of shoreline detection algorithms, challenges still exist in their widespread application. This is particularly true as satellite imagery becomes more widely available [30]. The present work was motivated by the failure of multiple existing algorithms when applied to multispectral images, from several different satellites, of the Lake Michigan shoreline, for reasons detailed later in the manuscript. In brief, it was found that there is no robust adaptive method that can be used to determine the threshold between the water and the land pixels, which, in turn, results in wrong or inaccurate detection of the shoreline position.
In this paper, we present and directly validate, using a newly-developed LiDAR system, a novel method for the automated detection of shorelines from multispectral imagery. The model is able to quantify the shoreline advancing and retreating between two input satellite images for the same place. The model is applied to satellite imagery of the Lake Michigan (USA) shoreline, which has seen extensive erosion between 2013 and present in response to high water levels [8].

Improved Procedure for Fully Automated Shoreline Detection from Multispectral Imagery
In this section, an improved method is presented for the fully automated detection of the shoreline position from multispectral imagery. This method introduces a new index to distinguish land from water and an automated procedure for the determination of a threshold value. The detection algorithm flow chart is shown in Figure 1. The shoreline detection algorithm was coded in MATLAB [31] and is available from the authors.
As input data, the model utilizes the green and the NIR band images. The first step is the application of a 2D median filter with a 3 × 3 stencil size to reduce the noise in the green and NIR images. The model then calculates the direct (non-normalized) difference water index (DDW I) as the difference between the Green and NIR bands Figure 1c: The direct difference water index is introduced here in lieu of numerous normalized indices used in previous studies Table 1 because it was seen to heighten the contrast between the magnitude of water and land pixel values, in turn allowing for the straightforward, automated determination of an image-specific index threshold that distinguishes water from land.
To determine an image-specific value of the DDW I threshold, the DDW I probability density function (pdf) for the image is calculated and then smoothed with an adaptive normal kernel, as in Figure 1d,e. As highlighted by [11,14], the probability density functions of existing normalized indices generally have bimodal distributions. This is also the case for the DDW I used in this study. However, DDW I has a higher intensity for water than that of the other indices. Moreover, unlike the other indexes mentioned by reference [14], the DDW I has been found to have always one main peak for the water, which indicates the ability of index to classify water pixels as a class. The threshold value is selected as the minimum value between the two peaks in the DDW I distribution, and this is found to provide a very robust, image-specific threshold value to separate water and land in the DDW I image. Additionally, to improve the reliability of the automated threshold determination, a hundred-class frequency distribution of the pixel values of the DDW I image is used to avoid the local peaks in the original probability distribution. We further discuss this improvement in the Section 4 of the paper.
Once the DDW I threshold value is determined, it is applied to the DDW I image to create a binary image, Figure 1f. Then, the model fills gaps such as small water bodies and misclassified pixels in the water and the land masks by using open, close, and fill mathematical morphologic operations [25,28]. At that point, if successful, the model will have two separate solid masks for the water and the land. The model will then calculate the boundary of the water mask and remove the edges of the image, if any, from that boundary. The resultant boundary is the shoreline detected from this image. As input data, the model utilizes the green and the NIR band images. The first step is the application of a 2D median filter with a 3 × 3 stencil size to reduce the noise in the green and NIR images. The model then calculates the direct (non-normalized) difference water index ( ) as the difference between the Green and NIR bands Figure 1c: The direct difference water index is introduced here in lieu of numerous normalized indices used in previous studies Table 1because it was seen to heighten the contrast be-

Shoreline Movement Quantification
A new method for the quantification of shoreline movement is also presented here. The traditional technique for calculating shoreline movement is the tracking of the shoreline position along shore-perpendicular transects (e.g., [15,32]). However, the results of transect-Remote Sens. 2022, 14, 557 5 of 18 based methods depend on the choice of transects (spacing and orientation), which, in turn, can lead to method-dependent results associated with shoreline change studies.
The new shoreline movement technique presented here is based on the buffer-overlay method proposed by [33] as an accuracy measure for linear features. This method has the advantage of being fully automated, resulting in its consistent application between shoreline studies. The shoreline evolution algorithm was developed and coded in Matlab [31].  The algorithm defines the first shoreline as a reference. Then, seaward and landward buffers are created within specific increments from the reference shoreline. Consequently, the number of pixels of the target shoreline that intersect the buffers are counted. Since the shoreline is composed of pixels from the input satellite image, the buffer increment must be a multiple of the pixel size. After including all the pixels in the second shoreline, the differences between the number of pixels in each buffer are calculated. Finally, the number of pixels in each class is converted to distance by multiplying by the effective pixel size. Since the shoreline can move in a horizontal or diagonal direction, the effective pixel size is taken as the average of the horizontal and the diagonal distance of the pixel.
As a result, the model produces a distribution for the different values of movements in the shoreline small segments (each segment is one pixel size). Moreover, the model returns a timeseries of the average shoreline movements for all the images.

Study Areas
The primary study location for this work is the Indiana coastline of Lake Michigan (USA) Figure 3. Lake Michigan is considered the third largest Great Lake by area and the second largest by volume [34]. Lake Michigan is not influenced by tides, and experiences the largest wind waves in the fall and winter [35]. The southern part of Lake Michigan along the Indiana coastline has recently experienced widespread erosion and shoreline damage in response to the nearly 2 m water level increase between 2013 and 2020 [8]. The coastline examined here is primarily sandy beaches backed by vegetated dunes, and includes the Indiana Dunes National Park. Validation of the shoreline detection method with in situ data was based on LiDAR surveys carried out on beaches in Beverly Shores and Dune Acres (described in this section). This coastline can experience large storm-generated waves during fall and winter, with significant wave heights of 5-6 m [36]. The application of the shoreline detection method to historical imagery was conducted on the Jeorse Park Beach, a 1.6 km-long closed littoral cell located east of Burns Harbor Figure 3. The algorithm defines the first shoreline as a reference. Then, seaward and landward buffers are created within specific increments from the reference shoreline. Consequently, the number of pixels of the target shoreline that intersect the buffers are counted. Since the shoreline is composed of pixels from the input satellite image, the buffer increment must be a multiple of the pixel size. After including all the pixels in the second shoreline, the differences between the number of pixels in each buffer are calculated. Finally, the number of pixels in each class is converted to distance by multiplying by the effective pixel size. Since the shoreline can move in a horizontal or diagonal direction, the effective pixel size is taken as the average of the horizontal and the diagonal distance of the pixel.
As a result, the model produces a distribution for the different values of movements in the shoreline small segments (each segment is one pixel size). Moreover, the model returns a timeseries of the average shoreline movements for all the images.

Study Areas
The primary study location for this work is the Indiana coastline of Lake Michigan (USA) Figure 3. Lake Michigan is considered the third largest Great Lake by area and the second largest by volume [34]. Lake Michigan is not influenced by tides, and experiences the largest wind waves in the fall and winter [35]. The southern part of Lake Michigan along the Indiana coastline has recently experienced widespread erosion and shoreline damage in response to the nearly 2 m water level increase between 2013 and 2020 [8]. The coastline examined here is primarily sandy beaches backed by vegetated dunes, and includes the Indiana Dunes National Park. Validation of the shoreline detection method with in situ data was based on LiDAR surveys carried out on beaches in Beverly Shores and Dune Acres (described in this section). This coastline can experience large storm-generated waves during fall and winter, with significant wave heights of 5-6 m [36]. The application of the shoreline detection method to historical imagery was conducted on the Jeorse Park Beach, a 1.6 km-long closed littoral cell located east of Burns Harbor Figure 3.

LiDAR Shoreline Surveys
Shoreline surveys were carried out with a newly developed, LiDAR-equipped backpack Mobile Mapping System (Backpack MMS) Figure 4. The Backpack MMS comprises a Velodyne VLP-16 Hi-Res LiDAR and a Sony α7R II 43.6 MP full-frame camera with a 35 mm lens. A Novatel SPAN-CPT GNSS/INS is used for direct georeferencing of the LiDAR and camera. Table 2. LiDAR manufacturer specifications for Backpack and UAV MMS. provides the manufacturer specifications for the LiDAR unit onboard the Backpack [37]. The Novatel SPAN-CPT GNSS/INS provides a post-processing accuracy of 2-5 cm for position, 0.025 • for roll/pitch angles, and 0.080 • for heading angle [38]. To reconstruct accurate point clouds from the backpack LiDAR, a rigorous system calibration is required to estimate the mounting parameters-lever arm and boresight angles-relating onboard LiDAR sensors with the GNSS/INS unit. In this study, the Backpack MMS was calibrated using the approach proposed by Ravi et al. [39]. The expected accuracy of the point cloud was estimated based on the individual sensor specifications using the LiDAR Error Propagation Calculator developed by Habib et al. [40]. The calculator suggests an accuracy of ±3 cm at a range of 50 m for the Backpack MMS.

LiDAR Shoreline Surveys
Shoreline surveys were carried out with a newly developed, LiDAR-equipped backpack Mobile Mapping System (Backpack MMS) Figure 4. The Backpack MMS comprises a Velodyne VLP-16 Hi-Res LiDAR and a Sony α7R II 43.6 MP full-frame camera with a 35 mm lens. A Novatel SPAN-CPT GNSS/INS is used for direct georeferencing of the LiDAR and camera. Table 2. LiDAR manufacturer specifications for Backpack and UAV MMS. provides the manufacturer specifications for the LiDAR unit onboard the Backpack [37]. The Novatel SPAN-CPT GNSS/INS provides a post-processing accuracy of 2-5 cm for position, 0.025° for roll/pitch angles, and 0.080° for heading angle [38]. To reconstruct accurate point clouds from the backpack LiDAR, a rigorous system calibration is required to estimate the mounting parameters-lever arm and boresight angles-relating onboard LiDAR sensors with the GNSS/INS unit. In this study, the Backpack MMS was calibrated using the approach proposed by Ravi et al. [39]. The expected accuracy of the point cloud was estimated based on the individual sensor specifications using the LiDAR Error Propagation Calculator developed by Habib et al. [40]. The calculator suggests an accuracy of ±3 cm at a range of 50 m for the Backpack MMS.   Walking shoreline surveys with the Backpack MMS were carried out on 1 June 2020, where a total of 4.4 km of the beach was surveyed. As long as there are no obstructions to the GNSS signal, the accuracy of the point cloud derived from Backpack LiDAR is expected to be in the range of 3-5 cm (this considers the accuracy of the GNSS/INS position and orientation system, LiDAR ranges, and system calibration). For shoreline surveys, this level of accuracy can easily be attained. The long range of the LiDAR sensor enabled the mapping of a wide extent across the shoreline, from the land-water intersection to vegetated areas on the landward side. Nonetheless, to ensure that the acquired point data were dense, a data acquisition configuration was designed to cover the shoreline from two tracks separated by a wide lateral distance, well under LiDAR's maximum range (one close to the shoreline and the other near the vegetated area on the landward side). Other than scanning various surface features at a high point density, this data acquisition configuration served to verify the alignment of point clouds from the two tracks, thus validating the system calibration and georeferencing accuracy. Moreover, this double track configuration allows for the acquisition of imagery on both sides of the shoreline (i.e., facing water and inland). In other Remote Sens. 2022, 14, 557 8 of 18 words, since the camera is mounted on the left side of the Backpack, forward and backward tracks along the shoreline would ensure the coverage of both sides by the captured imagery. A sketch of the track configuration of the Backpack system along a hypothetical shoreline is shown in Figure 5.   Walking shoreline surveys with the Backpack MMS were carried out on 1 June 2020, where a total of 4.4 km of the beach was surveyed. As long as there are no obstructions to the GNSS signal, the accuracy of the point cloud derived from Backpack LiDAR is expected to be in the range of 3-5 cm (this considers the accuracy of the GNSS/INS position and orientation system, LiDAR ranges, and system calibration). For shoreline surveys, this level of accuracy can easily be attained. The long range of the LiDAR sensor enabled the mapping of a wide extent across the shoreline, from the land-water intersection to vegetated areas on the landward side. Nonetheless, to ensure that the acquired point data were dense, a data acquisition configuration was designed to cover the shoreline from two tracks separated by a wide lateral distance, well under LiDAR's maximum range (one close to the shoreline and the other near the vegetated area on the landward side). Other than scanning various surface features at a high point density, this data acquisition configuration served to verify the alignment of point clouds from the two tracks, thus validating the system calibration and georeferencing accuracy. Moreover, this double track configuration allows for the acquisition of imagery on both sides of the shoreline (i.e., facing water and inland). In other words, since the camera is mounted on the left side of the Backpack, forward and backward tracks along the shoreline would ensure the coverage of both sides by the captured imagery. A sketch of the track configuration of the Backpack system along a hypothetical shoreline is shown in Figure 5.  The shoreline position was extracted from the Backpack MMS LiDAR point cloud as the lakeward extent of valid points from the LiDAR point cloud. The obtained shoreline position was validated using the georeferenced RGB images. The datapoints associated with floating objects in the water (buoys and boats) were manually removed from the point cloud.

Satellite Imagery
In this paper, recent satellite images from Planetscope, RapidEye, Sentinel 2, and Landsat 8 were used to develop and test the automated shoreline detection algorithm. A summary of the satellite specifications used in this study is shown in Table 3. The images were manually downloaded from the Planet Explorer website (https://www.planet.com/ explorer/ (accessed on 30 October 2021)). The website can be used to download commercial imageries (Planetscope, RapidEye) as well as publicly available images (Landsat 8 and sentinel 2). Images with extensive cloud cover and/or clarity problems were discarded; cloud cover obscured approximately 70% of available images. Additionally, images from winter months when the shoreline was covered in ice were not used.

Test 1-Comparison with LiDAR Surveys
For direct comparison with the LiDAR shoreline surveys, the new shoreline detection algorithm was applied to Planetscope, Sentinel 2, and Landsat 8 images of the 4.4-km-long section of the Indiana shoreline of Lake Michigan that was surveyed in Table 4. All three satellite images were acquired within two days of the LiDAR survey. The waves and the water level conditions during and between the image and the LiDAR surveys are shown in Table 4, which shows that the period of time including the image acquisitions and LiDAR survey had consistent water levels and very small waves. The water level change between the survey and the images was less than 0.04 m, as recorded at the Calumet Harbor, with the IL water level gauge located approximately 30 km from the survey location (NOAA Gauge 9087044, [41]). Assuming an average waterline slope of 1:12, as inferred from LiDAR observations [8], this change in the water level would result in a passive horizontal change in the shoreline location of 45 cm, which is substantially less than the imagery pixel resolutions Table 4. Additionally, the significant wave heights between LiDAR data collection and the satellite image acquisition were less than 0.1 m [42], which should generate a maximum vertical runup (R u2% ) of approximately 0.12 m following standard estimation techniques [43]. This run-up would generate a change in the shoreline location equal to 1.4 m, which is also less than the pixel resolution of the Planetscope image (3 m). Thus, while the LiDAR survey and the satellite image acquisitions were not exactly synchronous, the difference in the horizontal shoreline positions for the two estimates should be substantially less than the imagery resolution, indicating that, in general, the observed differences between the LiDAR-detected shoreline and the image-derived shoreline can be attributed to failures in the shoreline detection methodology.

Test 2-Comparison with Other Shoreline Detection Algorithms
To test the ability of the algorithm to work in an automated fashion on many images, satellite images were downloaded for the location of Jeorse Park Beach (IN), a small 1.6-km-long Lake Michigan beach located west of the LiDAR survey Figure 3 This beach was chosen because it is a challenging shoreline to detect due to the presence of numerous shoreline-parallel visual features, as seen in Figure 3. With littoral barriers on both sides, it is considered a close littoral cell [44]. Overall, 132 Planetscope and Rapideye images were downloaded for Jeorse Park beach for image acquisition dates between August 2018 and September 2021. The new shoreline detection algorithm was applied to these images and compared to different existing combinations of different water indexes and thresholding techniques Table 1. The resultant shorelines that were derived from all detection algorithms were then subjected to a visual screening in which the digitized shorelines were qualitatively evaluated relative to the shorelines seen on the images used for the shoreline digitization. The resulting comparison was categorized into three categories: accurate detection; good but not perfect detection; and failed detection.

Validation of Shoreline Detection Algorithm
Two types of tests were applied to the shoreline detection procedure. The first test was a direct comparison between the LiDAR-derived shoreline positions and the shoreline detected from the satellite imagery acquired within several days of the LiDAR surveys. This allows for a direct quantification of the error between the actual and the detected shorelines. The second test is the blind application of the shoreline detection algorithm to a large number of images for a common location.

Test 1-Comparison with LiDAR Surveys
The shorelines extracted from the detection algorithm and from the Planetscope image are plotted with the LiDAR-derived shoreline in Figure 6. The satellite-derived shoreline and the LiDAR shoreline have very good visual agreement, and the previously described buffer-overlay method found that 51% of the LIDAR points fall within the same pixel of the extracted shoreline from the Planetscope satellite image, and more than 99% of the points fall within a buffer of 1 pixel shoreward or landward, as shown in Table 5. Thus, the average error of the shoreline extracted from the Planetscope image is 1.48 m, which is less than the pixel resolution of the image.   The errors from the shorelines extracted from the coarser-resolution Sentinel 2 and Landsat 8 images are shown in Table 5  The errors from the shorelines extracted from the coarser-resolution Sentinel 2 and Landsat 8 images are shown in Table 5

Test 2-Comparison with Other Shoreline Detection Methods
The results of the large-scale test (132 images) comparing the new shoreline detection method with existing methods are shown in Figure 9. Of the 132 tested images, the DDWIbased detection algorithm produced accurate shoreline detection for 122 images, with good detection for nine images and failed detection for one image. In contrast, the OTSU and K-Means techniques failed to detect the shoreline using the NDWI index for 52 and 50 images, respectively. The K-means technique with RGB images was unable to detect the shoreline for any images, which is not surprising given the shoreline section used for the test Figure 10, which has numerous features that confound a clear distinction between land and water. When applied to the NDWI and NIR images, it can be seen that the proposed simpler thresholding technique works substantially better than the more complex OTSU and K-Means techniques.

Test 2-Comparison with Other Shoreline Detection Methods
The results of the large-scale test (132 images) comparing the new shoreline detection method with existing methods are shown in Figure 9. Of the 132 tested images, the DDWIbased detection algorithm produced accurate shoreline detection for 122 images, with good detection for nine images and failed detection for one image. In contrast, the OTSU and K-Means techniques failed to detect the shoreline using the NDWI index for 52 and 50 images, respectively. The K-means technique with RGB images was unable to detect the shoreline for any images, which is not surprising given the shoreline section used for the test Figure 10, which has numerous features that confound a clear distinction between land and water. When applied to the NDWI and NIR images, it can be seen that the proposed simpler thresholding technique works substantially better than the more complex OTSU and K-Means techniques.

Shoreline Evolution Module Results
The shoreline evolution technique was applied to the shorelines detected for the large-scale test to calculate the shoreline change over the imaging period (August 2018-

Shoreline Evolution Module Results
The shoreline evolution technique was applied to the shorelines detected for the large-scale test to calculate the shoreline change over the imaging period (August 2018-September 2021). The model considered the first shoreline, 1 August 2018, as the reference shoreline. Changes along the shoreline were averaged over the entire beach at each time and the resulting timeseries is shown in Figure 11. The shoreline position is seen to change in response to changing water levels and storms, with the general shoreline movement mirroring the water level variation. This is generally consistent with a simple inundation model where the change in the shoreline position ∆y can be calculated as ∆y = −tan(α) * ∆s with tan(α) = 1/20 corresponding to the average slope of the beach face and ∆s the change in the water level. Using this concept, the shoreline retreats when water levels increase and advances when water levels decrease. in response to changing water levels and storms, with the general shoreline move mirroring the water level variation. This is generally consistent with a simple inund model where the change in the shoreline position ∆ can be calculated as − ( ) * ∆ with ( ) = 1/20 corresponding to the average slope of the beach and ∆ the change in the water level. Using this concept, the shoreline retreats when ter levels increase and advances when water levels decrease.

Discussion
The results show the ability of the proposed fully automated shoreline dete method to accurately and reliably delineate shorelines using multispectral imagery method was tested with success on several satellites, including Planetscope, Rapi Sentinel 2, and Landsat 8. The accuracy of the technique was validated to subpixel racy using direct shoreline surveys carried out with a novel, backpack-based LiDAR tem. The detection algorithm was also applied to 132 high-resolution multispectral im and comparisons of the detection shoreline with existing techniques showed that the nique was more accurate than these techniques.
There are several key improvements associated with the new shoreline detectio gorithm that leads to improved accuracy and reliability. First, the method relies on a non-normalized index, the direct difference water index (DDWI = Green − NIR). Re to other existing water indices-both normalized and non-normalized-the DDWI heightens the contrast between land and water pixels. Perhaps most importantly DDWI shows a clear, significantly narrower minimum in its probability density fun

Discussion
The results show the ability of the proposed fully automated shoreline detection method to accurately and reliably delineate shorelines using multispectral imagery. The method was tested with success on several satellites, including Planetscope, Rapideye, Sentinel 2, and Landsat 8. The accuracy of the technique was validated to subpixel accuracy using direct shoreline surveys carried out with a novel, backpack-based LiDAR system. The detection algorithm was also applied to 132 high-resolution multispectral images, and comparisons of the detection shoreline with existing techniques showed that the technique was more accurate than these techniques.
There are several key improvements associated with the new shoreline detection algorithm that leads to improved accuracy and reliability. First, the method relies on a new non-normalized index, the direct difference water index (DDWI = Green − NIR). Relative to other existing water indices-both normalized and non-normalized-the DDWI index heightens the contrast between land and water pixels. Perhaps most importantly, the DDWI shows a clear, significantly narrower minimum in its probability density function in comparison to other indices Figure 12. In contrast to the most commonly applied water indices, the DDWI is not normalized, which precludes the a priori designation of a land-water threshold value that can be applied to an entire set of images (this can be an advantage of normalized indices). However, the narrow minimum value in the DDWI pdf can readily be found for each image, resulting in a robust image-specific threshold that can be used to classify land and water pixels.
Regarding the application of the more sophisticated Otsu and K-Means clustering techniques to the DDWI images, it was found that the simple technique of selecting the DDWI value corresponding to the minimum of its probability density function produced more accurate shoreline detection. Thus, it was determined that the complexity of the Otsu and K-Means clustering techniques was not warranted if the DDWI index is used. Additionally, if the NDWI or NIR indices are used instead of the proposed DDWI, the simple procedure of determining the threshold value as the minimum of the probability density function was found to be more reliable and accurate than the more complex Otsu and K-Means techniques as well Figure 9. This is because these techniques can be greatly affected by the outliers, which can shift the detected threshold away from a suitable value [26,27]. As such, the simple probability density function thresholding procedure may be applied as an improvement to existing shoreline detection algorithms based on other indices, offering a simpler alternative. in comparison to other indices Figure 12. In contrast to the most commonly applied water indices, the DDWI is not normalized, which precludes the a priori designation of a landwater threshold value that can be applied to an entire set of images (this can be an advantage of normalized indices). However, the narrow minimum value in the DDWI pdf can readily be found for each image, resulting in a robust image-specific threshold that can be used to classify land and water pixels.

Figure 12
Probability density functions (pdf) for different water indices calculated from the shoreline image shown in Figure 1.
Regarding the application of the more sophisticated Otsu and K-Means clustering techniques to the DDWI images, it was found that the simple technique of selecting the DDWI value corresponding to the minimum of its probability density function produced more accurate shoreline detection. Thus, it was determined that the complexity of the Otsu and K-Means clustering techniques was not warranted if the DDWI index is used. Additionally, if the NDWI or NIR indices are used instead of the proposed DDWI, the simple procedure of determining the threshold value as the minimum of the probability density function was found to be more reliable and accurate than the more complex Otsu and K-Means techniques as well Figure 9. This is because these techniques can be greatly affected by the outliers, which can shift the detected threshold away from a suitable value [26,27]. As such, the simple probability density function thresholding procedure may be applied as an improvement to existing shoreline detection algorithms based on other indices, offering a simpler alternative.
The new shoreline detection algorithm, applied to newly available high temporal and spatial resolution imagery, allows for unprecedented measurement of shoreline movement in response to changing water levels and storms. The Jeorse Park Beach shoreline position time series shown in Figure 11 is unique in terms of its fine temporal resolution and spatial accuracy, capable of resolving shoreline changes resulting from both seasonal water level fluctuations and individual storms. This dataset, and others derived from this method in the future, will allow for improved shoreline position modeling by providing important calibration and validation data at high temporal and spatial resolutions. The fully automated nature of the detection algorithm makes it ideally suited for high numbers of images with varying image intensities and conditions.
Beyond the Great Lakes, the shoreline detection technique shows promise for producing similar high-frequency, high-accuracy shoreline position timeseries with newly The new shoreline detection algorithm, applied to newly available high temporal and spatial resolution imagery, allows for unprecedented measurement of shoreline movement in response to changing water levels and storms. The Jeorse Park Beach shoreline position time series shown in Figure 11 is unique in terms of its fine temporal resolution and spatial accuracy, capable of resolving shoreline changes resulting from both seasonal water level fluctuations and individual storms. This dataset, and others derived from this method in the future, will allow for improved shoreline position modeling by providing important calibration and validation data at high temporal and spatial resolutions. The fully automated nature of the detection algorithm makes it ideally suited for high numbers of images with varying image intensities and conditions.
Beyond the Great Lakes, the shoreline detection technique shows promise for producing similar high-frequency, high-accuracy shoreline position timeseries with newly available multispectral imagery. Additional large-scale tests carried out on ocean shorelines with a range of shoreline types (rocky, sandy, and artificial) using both high-resolution and large-scale multispectral imagery (not shown) showed similar success in shoreline delineation for calm conditions, with some errors for shorelines with strong wave-breaking in the surf zone (whitecapping), a common difficulty for shoreline detection algorithms [15,23]. The model also fails with imagery with large amounts of cloud cover. Further tests are needed to determine the suitability of the algorithm for other coastal environments where the effects of wave-breaking, wave run-up, and different coastal geologies are present.
The use of a Backpack system in coastal applications has several advantages. As an active sensor, LiDAR can scan poorly lit, hard-to-reach shoreline features such as caves and sinkholes. On the other hand, satellite imagery cannot scan these features. While it is possible to generate image-based 3D surface models, intricate subsurface features are difficult to resolve through distant aerial/spaceborne imagery.

•
Using LiDAR, it is possible to delineate the actual shoreline (land-water intersection). In areas where a clear distinction cannot be established, imagery can be used to verify the ground condition.

•
The point density of a Backpack-based LiDAR is very high compared to UAV-based point clouds due to the proximity of the former to the surface. However, this comes at the price of an increased data acquisition time in the case of the Backpack system. • The use of a Backpack-based LiDAR system may not require extensive regulatory approval, unlike vehicle-based or unmanned aerial platforms (as in Indiana dunes). More specifically, the Backpack system can be used regardless of the nature/restrictions of the airspace at the vicinity of shorelines. • Unlike UAVs, Backpack-based LiDAR can be used in less favorable weather conditions. Additionally, it does not have any requirements for a clear takeoff/landing site, as in the case of UAV systems.
A backpack-based mapping system can thus be a valuable tool to assess coastal conditions together with other established state-of-art technologies.
Additionally, a new technique to quantify shoreline movement, based on the buffer overlay method, was developed and tested on the detected shorelines. When used with digitized shorelines, the buffer-overlay method allows for the unambiguous quantification of the shoreline evolution timeseries, without the need to manually define cross-shore transects. The choice of the transects can be problematic especially with curved shorelines, and the application of the new method for shoreline movement quantification may lead to more consistency between shoreline movement studies.

Conclusions
In this paper, a new, fully automated model was developed to extract the shoreline position from multispectral satellite images. The main advantage of the model is its ability to automatically determine the needed threshold for different images using a novel technique and a modified water index, the "direct difference water index" (DDWI = NIR − Green). The model was applied to Planetscope, Rapideye, Sentinel 2, and Landsat 8 satellite images. The accuracy of the shoreline model was directly validated using LiDAR data collected with a novel backpack-based system that surveyed the Lake Michigan shoreline. The model was able to detect the shoreline from the different satellite images with subpixel accuracy. Additionally, the shoreline detection model was applied to 132 high-resolution multispectral images and qualitatively compared with other shoreline detection algorithms. The comparison showed that the model consistently and reliably detected the shoreline, with a success rate exceeding 99%, which exceeded all other existing models. Moreover, the shoreline detection model is supplemented with a transect-free shoreline evolution module that is able to compare different shoreline positions and construct a shoreline change timeseries. The Lake Michigan time series resulting from the new shoreline detection method's application provides an unprecedented satellite-based quantification of shoreline movements on timescales that were previously unresolved by shoreline detection methods for the Great Lakes.
It is hoped that the Lake Michigan dataset produced here, and other similar highresolution shoreline timeseries produced by the new shoreline detection method and its improvements, will aid in the refinement of shoreline evolution models capable of resolving storm-scale and seasonal-scale effects. Better shoreline evolution models, in turn, may help coastal communities make informed decisions that improve coastal resilience.  Wave data are available from the NOAA National Data Buoy Center (NOAA NDBC) website (https: //www.ndbc.noaa.gov/, accessed on 30 October 2021). The LiDAR and shoreline position data presented in this study are available on request from the corresponding author, and will be publicly available following the completion of the research projects funding this work.