Assessing the Impacts of Rising Sea Level on Coastal Morpho-Dynamics with Automated High-Frequency Shoreline Mapping Using Multi-Sensor Optical Satellites

: Rising sea level is generally assumed and widely reported to be the signiﬁcant driver of coastal erosion of most low-lying sandy beaches globally. However, there is limited data-driven evidence of this relationship due to the challenges in quantifying shoreline dynamics at the same temporal scale as sea-level records. Using a Google Earth Engine (GEE)-enabled Python toolkit, this study conducted shoreline dynamic analysis using high-frequency data sampling to analyze the impact of sea-level rise on the Malaysian coastline between 1993 and 2019. Instantaneous shorelines were extracted from a test site on Teluk Nipah Island and 21 tide gauge sites from the combined Landsat 5–8 and Sentinel 2 images using an automated shoreline-detection method, which was based on supervised image classiﬁcation and sub-pixel border segmentation. The results indicated that rising sea level is contributing to shoreline erosion in the study area, but is not the only driver of shoreline displacement. The impacts of high population density, anthropogenic activities, and longshore sediment transportation on shoreline displacement were observed in some of the beaches. The conclusions of this study highlight that the synergistic use of multi-sensor remote-sensing data improves temporal resolution of shoreline detection, removes short-term variability, and reduces uncertainties in satellite-derived shoreline analysis compared to the low-frequency sampling approach. To evaluate the automatic shoreline-extraction techniques employed in this study, we compared the results with shorelines extracted from selected satellite imagery acquired from USGS and Sci-hub using manual and semi-automatic techniques. The semi-automatic extraction was conducted in ArcGIS Pro, while the manual on-screen digitizing of shorelines was done in ArcGIS Desktop. Figure 8 shows the shoreline extracted from Landsat 7 with the three techniques. The manually extracted shoreline was taken as the actual shoreline for the validation of the automatic and semi-automatic extracted shorelines. The mean deviations of the transect points from the manually extracted shorelines were computed. To evaluate the automatic shoreline-extraction techniques employed in this study, we compared the results with shorelines extracted from selected satellite imagery acquired from USGS and Sci-hub using manual and semi-automatic techniques. The semi-automatic extraction was conducted in ArcGIS Pro, while the manual on-screen digitizing of shorelines was done in ArcGIS Desktop. Figure 8 shows the shoreline extracted from Landsat 7 with the three techniques. The manually extracted shoreline was taken as the actual shoreline for the validation of the automatic and semi-automatic extracted shorelines.


Introduction
Coastal cities are spread around different regions of the world, with more than two-thirds of the global population living within 100 km of the sea [1,2]. Coastal areas are essential regions of socio-economic and environmental benefits [3], densely populated, and highly urbanized, and are usually highly morpho-dynamic [4]. The shoreline, an intersection of the marine environment and the mobile coastal system, is one of the most dynamic features of coasts [3]. About 31% of the 1.63 million km of world coastlines are sandy coasts that are highly susceptible to rapid advances/retreats and long-term accretion/erosion due to changing environmental conditions and coastal processes [2]. Techniques that couple physical landscape processes associated with climate change with human behavior have been employed in many coastal susceptible and vulnerability assessment studies [5]. Anfuso et al. [6] presented a comprehensive overview of these techniques, outlining data requirements and vulnerability indices used in adaptation planning. Generally, assessing the causes of coastal vulnerability, especially shoreline erosion and coastal flooding, is a difficult task due to the complexity of the processes acting at varying spatial and temporal

Satellite Derived Shorelines
Different techniques have been employed for extracting instantaneous shoreline positions from satellite images. The manual method involves visual inspection and digitizing of shorelines from imagery [36]. The semi-automatic method includes single band, band ratio [36], image classification [37,38], and principal component [39] analyses. The Normalized Difference Water Index (NDWI) [40], Modified Normalized Difference Water Index (MNDWI) [41], and Automated Water Extraction Index (AWEI) [42] are the commonly used indices for shoreline delineation.
With the advancements in machine learning, supervised and unsupervised classification techniques have been also used for shoreline detection. Iterative Self-Organizing Data Analysis (ISODATA) is a common unsupervised classification algorithm for extracting shoreline features from satellite imagery [43]. In particular, ISODATA classification is an improved version of K-means [28]. For example, García-Rubio, et al. [44] assessed the various land-sea classification techniques employed in shoreline studies and confirmed the great performance of ISODATA shoreline classification. Among supervised classification methods, support vector machine is a well-known method for land-sea classification. Elnabwy, et al. [45] classified Landsat images using SVM to extract instantaneous shorelines along the coasts of Nile Delta. Studies have also employed neural network classification to discriminate sand and water features [34]. Feed forward networks such as the multi-layer perceptron (MLP), radial basis function (RBF), probabilistic neural network (PNN), and generalized regression neural network (GRNN) are the commonly used ANN models in remote-sensing applications [46]. This is because neural networks have a strong capability to handle complex phenomena and can improve land-sea segmentation accuracy [47,48]. However, developing a more robust technique for precise and accurate extraction of shorelines from the available satellite data remains a challenging task [38,49]. Moreover, shoreline detection is a complex phenomenon that demands the use of several techniques, rather than a single image-processing technique [28]. As for data management, the file size of Landsat 8 images can be as high as 1 GB per scene, meaning that high-frequency shoreline sampling for long-term studies requires high data storage capabilities and processing power that is difficult to achieve with a conventional Geographic Information System (GIS) approach. Based on the foregoing, it is evident that to understand shoreline dynamics within the context of sea-level rise, there is a need for an innovative approach to manage the large volume of data, a robust detection technique for accurate shoreline delineation, and appropriate analysis techniques to accurately quantify the spatiotemporal evolution of a beach.

Malaysia: Sea-Level Rise and Coastal Vulnerability
A rising sea level has been reported along the Malaysian coastline [50,51], and a corresponding increase in coastal hazards such as flooding and erosion have been observed in many coastal areas [52]. The coasts of Selangor and Batu Pahat have experienced severe coastal erosion, leading to the loss of 18.785 km 2 and 4.155 km 2 of land, respectively [53]. Similarly, coastal flooding in Johor has caused an economic loss of an estimated  [53]. Existing studies of sea levels are fragmental in terms of space and time coverage. In addition, the lack of accurate local sea-level information, in addition to high uncertainties in shoreline monitoring, makes it difficult to quantify the impact of sea-level rise on shoreline dynamics. A nationwide study of sea-level change and shoreline dynamic is lacking. To address this gap, this study leveraged advances in geospatial technology to conduct a holistic evaluation of long-term shoreline dynamics and sea-level change along the Malaysian coastline. For an effective management of the shoreline data, Google Earth Engine (GEE) was used to acquire and process multi-mission satellite images of Landsat series 4-8 and Sentinel 2. In particular, GEE is a free-to-use, cloud-based geospatial platform for large-scale analysis of big EO data [54]. It provides a petabyte of satellite data, high-level API, and machine-learning algorithms for processing and analyzing big earth observation data, thereby overcoming the limitations of the conventional digital image-processing workflows [55]. For shoreline extraction, a robust shoreline detection method based on a supervised classification and sub-pixel resolution segmentation technique developed by Vos, Harley, Splinter, Simmons and Turner [34] was employed to improve the detection accuracy of the sand and land interface.
The rest of this paper is arranged as follows: descriptions of the study area and data used are provided in Section 2; the methodology employed to estimate sea level and shoreline displacement is provided in Section 3; and Section 4 presents the results of the study. A discussion of the results is then presented in Section 5, followed by the conclusions of the study in Section 6.

Study Area and Data Used
This work was performed around the Malaysian coastline, which measures about 4675 km, making it the 29th-longest coastline in the world. The country has a large population inhabiting lowlands in coastal areas; more than 13% of its land area is within 5 km of the coast, and thereby vulnerable to the harmful impacts of the sea-level rise [56]. The aesthetic natural environment of most islands in Malaysia draws millions of tourists to the country every year. Therefore, the development of tourist attractions in the natural-settings area results in a large impact on the coastal zone.
The data employed to accomplish the aim of this study were acquired from different sources. These data had different formats and dimensions, as presented in Table 1 below. The table also provides information on the source of the data and their specific use in this study. To quantify sea-level change over the study area, monthly averaged sealevel data were acquired from the Permanent Service of Mean Sea Level (PSMSL) and Copernicus Marine Environment Service (CMEMS), respectively. PSMSL provides highfrequency sea-level data measured by tide gauges, while CMEMS provides multi-mission altimeter satellite gridded sea-level anomalies (SLAs). A total of 21 tide-gauge stations cover the Malaysian coastline-12 of these stations are along the coast of Peninsular Malaysia (West Malaysia), and 9 tidal stations are along the coast of Sabah and Sarawak (East Malaysia) ( Figure 1).
Optical satellites were employed to evaluate shoreline dynamics at regions of extreme SLR. These datasets, including Landsat (4,5,7,8) and Sentinel-2 imagery, were acquired from GEE. Landsat imagery of the free and open satellite data were organized in two tiers. Tier 1 images had well-characterized geometry and valid geometric accuracy, and were thus suitable for time series analysis, while Tier 2 images were not too fit for this purpose [57]. Landsat Tier 1 images have three products, namely Raw scenes, Top-of-Atmosphere (TOA), and Surface Reflectance. Raw scenes contain just the digital number (DN) as measured by the sensor; the DN values from the sensor are converted to reflectance for Top-of-Atmosphere (TOA) reflectance images [58]; while the Surface Reflectance image is atmospherically corrected. TOA provides a standardized comparison between images acquired on different dates. The Level-1C product of Sentinel-2 has similar standardized TOA reflectance images that are appropriate for time-series analysis. The Shuttle Radar Topography Mission (SRTM) digital elevation dataset and the 2014 Finite Element Solution (FES2014) ocean tide model, which are required for tidal correction of instantaneous shorelines, were acquired from GEE and AVISO, respectively. The next section describes the methodology employed to process and analyze the data.  For estimating tide level required for tidal correction Optical satellites were employed to evaluate shoreline dynamics at regions of extreme SLR. These datasets, including Landsat (4,5,7,8) and Sentinel-2 imagery, were acquired from GEE. Landsat imagery of the free and open satellite data were organized in two tiers. Tier 1 images had well-characterized geometry and valid geometric accuracy, and were thus suitable for time series analysis, while Tier 2 images were not too fit for this purpose [57]. Landsat Tier 1 images have three products, namely Raw scenes, Top-of-Atmosphere (TOA), and Surface Reflectance. Raw scenes contain just the digital number

Methodology
The approach to accomplish the goals of this study, including estimation of sea-level trend, extraction of shorelines, and the analysis of displacement pattern is presented in this section. Figure 2 below summarizes the study's methodology.
image is atmospherically corrected. TOA provides a standardized comparison between images acquired on different dates. The Level-1C product of Sentinel2 has similar standardized TOA reflectance images that are appropriate for time-series analysis. The Shuttle Radar Topography Mission (SRTM) digital elevation dataset and the 2014 Finite Element Solution (FES2014) ocean tide model, which are required for tidal correction of instantaneous shorelines, were acquired from GEE and AVISO, respectively. The next section describes the methodology employed to process and analyze the data.

Methodology
The approach to accomplish the goals of this study, including estimation of sea-level trend, extraction of shorelines, and the analysis of displacement pattern is presented in this section. Figure 2 below summarizes the study's methodology.

Estimation of Sea-Level Trend
The monthly mean sea-level record of the 21 tide gauge stations acquired from PSMSL was preprocessed, and the mean sea level was subtracted to derive the relative sea-level anomaly. Sea-level anomaly values were also extracted from gridded altimetry satellites at the locations of the tide-gauge stations to estimate the corresponding absolute sea-level trend. The trend in absolute and relative sea levels was determined for all stations. The trend was estimated using robust fit regression techniques, since robust fit analysis optimizes solution determination and outlier detection [59]. This approach was used to fit a linear trend at each station with an iteratively re-weighted least squares (IRLS) technique, wherein weights of observations were adjusted according to the deviations from the trendline to re-fit the trendline repeatedly until the solution converged [60].

Shoreline Extraction
Mapping of the coastline at sites of extreme sea-level change was conducted using CoastSat [4], a GEE-enabled Python package. The toolkit exploits the capabilities of Google Earth Engine, machine learning, and other image-processing packages, such as scikit-learn and scikit-image, to facilitate access, preprocessing, detection, and automatic extraction of instantaneous shoreline positions from multi-spectral optical imagery. The

Estimation of Sea-Level Trend
The monthly mean sea-level record of the 21 tide gauge stations acquired from PSMSL was preprocessed, and the mean sea level was subtracted to derive the relative sea-level anomaly. Sea-level anomaly values were also extracted from gridded altimetry satellites at the locations of the tide-gauge stations to estimate the corresponding absolute sealevel trend. The trend in absolute and relative sea levels was determined for all stations. The trend was estimated using robust fit regression techniques, since robust fit analysis optimizes solution determination and outlier detection [59]. This approach was used to fit a linear trend at each station with an iteratively re-weighted least squares (IRLS) technique, wherein weights of observations were adjusted according to the deviations from the trendline to re-fit the trendline repeatedly until the solution converged [60].

Shoreline Extraction
Mapping of the coastline at sites of extreme sea-level change was conducted using CoastSat [4], a GEE-enabled Python package. The toolkit exploits the capabilities of Google Earth Engine, machine learning, and other image-processing packages, such as scikit-learn and scikit-image, to facilitate access, preprocessing, detection, and automatic extraction of instantaneous shoreline positions from multi-spectral optical imagery. The steps are further described below; to support open reproducible science, the code written for this process is available as a supplemental document.

Data Acquisition
Using CoastSat, Top-of-Atmosphere reflectance images from the Landsat 5, 7, and 8 Tier 1 collection, as well as Level-1C products of Sentinel-2 at regions of interest (ROI), were retrieved from GEE. The four satellite missions were chosen to maximize the frequency of sampling the coastline. To reduce the size of files to download, the toolkit cropped the image to the ROI and downloaded only the bands necessary for shoreline detection. Details Remote Sens. 2021, 13, 3587 7 of 22 of satellite data retrieved by CoastSat, such as pixel size of required bands, revisit period of mission, and the GEE collections for the image storage, are provided in Table 2 below.

Data Preprocessing
The acquired images were orthorectified by GEE. Sharpening of Landsat 7 and Landsat 8 images with higher resolution panchromatic bands was performed by Coast-Sat to enhance the spatial resolution of the images from 30 m to 15 m, with the aim of achieving an optimal shoreline detection. Landsat 5 scenes were down-sampled from 30 m to 15 m, as they do not include a panchromatic band. Sentinel-2 images have a higher spatial resolution: 10 m visible and NIR bands, and 20 m SWIR1. The SWIR band was also down-sampled to 10 m using bilinear interpolation to maintain equal resolution among the S2 bands. For cloud-masking, we chose a threshold of 0.2 to discard all images with more than 20% cloud cover in our ROI.

Shoreline Detection
The instantaneous boundary between water and sand captured at the time of image acquisition was taken as the indicator of the shoreline here. The CoastSat toolkit applied a robust shoreline detection algorithm to the pre-processed satellite images [4]. The technique involves two major steps-image classification and sub-pixel segmentation. A multi-layer perceptron, a neural-network-supervised classifier algorithm, was pretrained, with 99% accuracy based on cross-validation, to categorize the pixels in the image into four classes-sand, water, white water, and other features (containing features such as vegetation, rock, buildings, and coastal-defense structures). A sample classification of a Landsat 5 image for Kukup Beach is shown in Figure 3 below. The sand/water boundary for each classified image was then delineated using the Modified Normalized Difference Water Index (MNDWI).

Tidal Correction
The satellite images were acquired at different stages of tides. So, a linear tidal correction [61] was applied to the extracted shorelines using Equation (1) below to remove the displacement due to tidal fluctuation.
where ∆x is the cross-shore horizontal shift, Zre f is the reference tidal datum, Zlocal is the measured (or modelled) tide level at the time of image acquisition, and tanβ is a characteristic beach face slope of the site. The slope of the beach at each site was derived from the SRTM elevation dataset, and the tide level was acquired from the FES2014 ocean tide model [62]. Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 22

Tidal Correction
The satellite images were acquired at different stages of tides. So, a linear tidal correction [61] was applied to the extracted shorelines using Equation (1) below to remove the displacement due to tidal fluctuation.
where ∆ is the cross-shore horizontal shift, is the reference tidal datum, is the measured (or modelled) tide level at the time of image acquisition, and is a characteristic beach face slope of the site. The slope of the beach at each site was derived from the SRTM elevation dataset, and the tide level was acquired from the FES2014 ocean tide model [62].

Shoreline Analysis
After the extraction of the shorelines, baseline and shore-normal transects were created. Using the Digital Shoreline Analysis System (DSAS) [63], significant shoreline change statistics such as the shoreline change envelope (SCE), net shoreline movement (NSM), end point rate (EPR), linear regression rate (LRR), and weighted linear regression (WLR) were calculated. The principle behind the calculation was based on the measured differences in shoreline position over time. SCE is the range of all the shorelines that intersect a given transect, and NSM is the distance between the oldest and the youngest shorelines for each transect. The end point rate (EPR) was calculated by dividing the distance of shoreline movement by the time elapsed between the oldest and the most recent shoreline. The major advantages of the EPR are the ease of computation and the minimal requirement of only two shoreline dates. A linear regression rate-of-change statistic can be determined by fitting a least-squares regression line to all shoreline points for a transect. The regression line is placed so that the sum of the squared residuals (determined by squaring the offset distance of each data point from the regression line and adding the squared residuals together) is minimized. The linear regression rate is the slope of the line. In a weighted linear regression, the more reliable data are given greater emphasis or

Shoreline Analysis
After the extraction of the shorelines, baseline and shore-normal transects were created. Using the Digital Shoreline Analysis System (DSAS) [63], significant shoreline change statistics such as the shoreline change envelope (SCE), net shoreline movement (NSM), end point rate (EPR), linear regression rate (LRR), and weighted linear regression (WLR) were calculated. The principle behind the calculation was based on the measured differences in shoreline position over time. SCE is the range of all the shorelines that intersect a given transect, and NSM is the distance between the oldest and the youngest shorelines for each transect. The end point rate (EPR) was calculated by dividing the distance of shoreline movement by the time elapsed between the oldest and the most recent shoreline. The major advantages of the EPR are the ease of computation and the minimal requirement of only two shoreline dates. A linear regression rate-of-change statistic can be determined by fitting a least-squares regression line to all shoreline points for a transect. The regression line is placed so that the sum of the squared residuals (determined by squaring the offset distance of each data point from the regression line and adding the squared residuals together) is minimized. The linear regression rate is the slope of the line. In a weighted linear regression, the more reliable data are given greater emphasis or weight towards determining a best-fit line. In the computation of rate-of-change statistics for shorelines, greater emphasis is placed on data points for which the position uncertainty is smaller. The weight (w) is defined as a function of the variance in the uncertainty of the measurement (e): where e is the shoreline uncertainty value. NSM and SCE are distance-measurement statistics calculated in meters, while EPR, LRR, and WLR are rates of change calculated in m/y. The relationship between sea-level change and shoreline dynamics at all tide-gauge sites was determined based on Pearson correlation analysis. Pearson's correlation coefficient evaluates the statistical relationship, or association, between two continuous variables. It is based on the concept of covariance and recognized as an optimal approach for measur- ing the association between variables of interest [64]. The correlation formula is as shown in Equation (3) below: where r is the Pearson correlation coefficient, the value of which ranges between −1 and 1, where +1 indicates a perfect positive relationship, −1 indicates a perfect negative relationship, and 0 denotes the absence of a relationship; x i is the yearly averaged sea-level change values; y i is the yearly average shoreline displacement values; χ is the mean of the yearly average sea-level change values; and γ is the mean of the yearly averaged shoreline change values.

Results
Relative and absolute sea-level trends estimated for the Malaysian coastline are presented in Figures 4 and 5, respectively.
Shorelines were extracted from the multi-sensor satellites at all 21 locations with tide gauges using the shoreline-detection method described above. The application of the automated shoreline-detection method was first validated at a test site prior to its use at other tide gauge sites. Teluk Nipah was selected as the test site because a previous study [52] had been conducted at this location using the conventional approach. The change statistics for the shorelines extracted at the test sites using the robust shorelinedetection algorithm employed in this study are presented here. A total of 428 shorelines were correctly extracted from the 825 images acquired over Teluk Nipah Beach. The inability to extract shorelines from some images was due to excessive cloud cover. The result of shoreline analysis at the test sites presented in Table 3 Figure 6. The mean shoreline change envelope was 59.14 m. Figure 7 shows the map and plot of the shoreline change distance for SCE and NSM. The mean of erosion/accretion values at all transects showed that erosion was more dominant across the coastline of Teluk Nipah Beach at a rate of −0.711 m/y, −0.23 m/y, and −0.30 m/y based on EPR, LRR, and WLR, respectively. The result of a similar analysis conducted at the 21 tide gauge stations and the corresponding sea level change information is summarized in Table 4 below. The relationship between shoreline displacement and sea-level change based on the correlation of mean annual shoreline displacement and sea-level data is also provided in this table. Sea-level data were obtained from tide gauges and satellite altimetry to estimate the relative and absolute sea-level trend, respectively. However, the tide-gauge data had missing values for some timestamps, and some stations had stopped providing records. In addition, some stations were flagged for errors in the PSMSL database. In contrast, absolute sea-level data were not affected by these limitations. So, we used absolute sea-level change for the correlation analysis because it was not affected by missing data and had more recent records than those for relative sea level, which are usually affected by tide-gauge instability and vertical land movement.     Shorelines were extracted from the multi-sensor satellites at all 21 locations with tide gauges using the shoreline-detection method described above. The application of the automated shoreline-detection method was first validated at a test site prior to its use at other tide gauge sites. Teluk Nipah was selected as the test site because a previous study [52] had been conducted at this location using the conventional approach. The change PSMSL database. In contrast, absolute sea-level data were not affected by these limitations. So, we used absolute sea-level change for the correlation analysis because it was not affected by missing data and had more recent records than those for relative sea level, which are usually affected by tide-gauge instability and vertical land movement.      To evaluate the automatic shoreline-extraction techniques employed in this study, we compared the results with shorelines extracted from selected satellite imagery acquired from USGS and Sci-hub using manual and semi-automatic techniques. The semi-automatic extraction was conducted in ArcGIS Pro, while the manual on-screen digitizing of shorelines was done in ArcGIS Desktop. Figure 8 shows the shoreline extracted from Landsat 7 with the three techniques. The manually extracted shoreline was taken as the actual shoreline for the validation of the automatic and semi-automatic extracted shorelines. The mean deviations of the transect points from the manually extracted shorelines were computed. To evaluate the automatic shoreline-extraction techniques employed in this study, we compared the results with shorelines extracted from selected satellite imagery acquired from USGS and Sci-hub using manual and semi-automatic techniques. The semiautomatic extraction was conducted in ArcGIS Pro, while the manual on-screen digitizing of shorelines was done in ArcGIS Desktop. Figure 8 shows the shoreline extracted from Landsat 7 with the three techniques. The manually extracted shoreline was taken as the actual shoreline for the validation of the automatic and semi-automatic extracted shorelines. The mean deviations of the transect points from the manually extracted shorelines were computed. Furthermore, shoreline-change assessment of the manual and semi-automatic extracted shorelines was performed at low-frequency sampling at an interval of approximately a decade. The resulting change statistics for the shoreline are presented in Table 5. Figure 9 below shows the elevation profile and developments along the test sites.  Furthermore, shoreline-change assessment of the manual and semi-automatic extracted shorelines was performed at low-frequency sampling at an interval of approximately a decade. The resulting change statistics for the shoreline are presented in Table 5. Figure 9 below shows the elevation profile and developments along the test sites.

High-Frequency Sampling and Automatic Extraction with GEE
Regular and periodic monitoring of shoreline changes is vital for coastal management and future planning. While the contemporaneous observation of the Earth by multiple satellite missions facilitates regular and long-term monitoring of shoreline changes, big-data management has been a challenge in leveraging this opportunity. Consequently, earlier studies employed low-frequency data sampling for shoreline analysis. A major limitation with low-frequency sampling is that it does not resolve the shoreline dynamics well enough to quantify the impact of sea level, which is usually sampled at hourly, daily, or monthly frequencies. In addition, there are several uncertainties due to short-term variability and tidal effects on instantaneous shorelines. So, a method to automatically extract shorelines from multi-mission optical satellites at high-frequency sampling with improved data management and a robust extraction algorithm was applied in this study. Between January 1993 and October 2019, a total of 845 images were available for shorelinechange analysis at the test site. Consequently, 845 satellite images (266 from Landsat 5, Figure 9. Elevation profile and development along the coast of Teluk Nipah Beach. The profile was drawn over the July 2020 shoreline, along the north direction.

High-Frequency Sampling and Automatic Extraction with GEE
Regular and periodic monitoring of shoreline changes is vital for coastal management and future planning. While the contemporaneous observation of the Earth by multiple satellite missions facilitates regular and long-term monitoring of shoreline changes, bigdata management has been a challenge in leveraging this opportunity. Consequently, earlier studies employed low-frequency data sampling for shoreline analysis. A major limitation with low-frequency sampling is that it does not resolve the shoreline dynamics well enough to quantify the impact of sea level, which is usually sampled at hourly, daily, or monthly frequencies. In addition, there are several uncertainties due to short-term variability and tidal effects on instantaneous shorelines. So, a method to automatically extract shorelines from multi-mission optical satellites at high-frequency sampling with improved data management and a robust extraction algorithm was applied in this study.
Between January 1993 and October 2019, a total of 845 images were available for shorelinechange analysis at the test site. Consequently, 845 satellite images (266 from Landsat 5, 252 from Landsat 7, 135 from Landsat 8, and 192 from Sentinel-2) returned after temporal and spatial filtering were processed with CoastSat, a GEE-enabled toolkit. GEE supports the acquisition of only the required bands of satellites and cropping to the region before downloading, thus significantly reducing the data size. A total of 428 of the 845 satellite images that met the quality requirement of cloud cover of less than 20% were just 20.83 MB The total size of the four satellite images was 1.7 GB, and they were required to be downloaded in their entirety before further processing such as clipping and pan sharpening. One image from Landsat 8 acquired from USGS was 41 times larger than the total size of 428 multi-sensor satellite images retrieved from GEE. Thus, GEE addressed the big-data storage and management issue of SDS change analysis. This facilitated the effective combination of multiple satellite missions at high sampling frequency, which made it possible to offer continuous data throughout the year. By ensuring that shorelines were extracted for almost every month of the study period, the impact of sea-level rise on shoreline displacement could be estimated.
When compared to the manually digitized shoreline, the automatically extracted shoreline appeared to be in closer proximity than the semi-automatic extracted shorelines as shown in Figure 8. The mean deviation of the transect points of the automatic and the semi-automatic extracted shorelines from the actual shoreline were 0.104 m and 0.119 m, respectively. The shoreline extracted manually through digitizing of the land-sea boundary from a false-colour composite of satellite scenes was accurate, but time-consuming and labour-intensive for high-frequency shoreline sampling. Semi-automatic extraction of shorelines using tasseled cap transformation and NDVI also takes time depending on the number of images and user experience. In addition, there is a high propensity for false classification due to clouds and breaking waves. In contrast, automatic shoreline extraction enables faster and easier extraction of multi-shorelines while eliminating errors due to human intervention. It took approximately 10 min to extract 428 shorelines from the multi-mission satellites after training with the robust ML algorithm. The extracted shorelines based on the automated shoreline detection method were as accurate as the manually and carefully digitized shorelines.
In an earlier study, Foo, Teh and Babatunde [52] reported the decadal shoreline change rates of Teluk Nipah in the 1990s, 2000s, and 2010s using EPR and LRR. The shoreline variation for the 1990s was relatively stable, with maximum erosion and accretion rates of −1.64 m/y (−0.74 m/y) and 2.48 m/y (2.95 m/y), respectively, based on EPR (LRR). Contrary to the 1990s, Teluk Nipah experienced severe erosion in the 2010s. The study estimated the maximum erosion rate for the 2010s as −6.83 m/y and −5.68 m/y based on EPR and LRR, respectively, but acknowledged that it is most likely that this indication of critical erosion was an overestimation of the actual result. The result was obtained from the analysis of four instantaneous satellite-derived shorelines extracted using a semiautomatic approach with no tidal correction. The study suggested that the unusually high erosion rate was perhaps due to the tidal effect on the instantaneous shorelines acquired at different times. The use of high-frequency sampling and tidal correction of a time series of shoreline displacement in the present study effectively reduced uncertainties in the satellitederived shoreline analysis due to short-term variability and tidal effects, respectively. It was noteworthy that there was no difference in some change analyses between high-and low-frequency sampling. The NSM and EPR estimates were approximately the same for high-and low-frequency analyses, both indicating that 85.71% and 14.29% of the shoreline experienced erosion and accretion, respectively. This was not unexpected, because these change statistics only considered two endpoints (oldest and earliest date) and ignored other sampled data in between. Other change statistics; however, showed a significant difference between the low-and high-frequency sampling analyses. The linear regression rate (LRR) statistic for the low-frequency sampling analysis indicated that 92.31% of the shoreline experienced erosion at an average rate of −0.79 ± 1.32 m/y, while the highfrequency sampling analysis indicated 50% shoreline erosion at a rate of −0.17 ± 0.07 m/y. The maximum erosion rates of −2.26 m/y and −1.37 m/y were estimated based on LRR statistics of the low-and high-frequency sampling, respectively. With this acceptable lowerintensity erosion rate, this study established the potential of high-frequency-sampled, tidalcorrected, satellite-derived shorelines in reducing uncertainties in instantaneous shorelines.

SLR and Shoreline Dynamics
Sea-level change is usually perceived to be responsible for long-term shoreline displacement, but our study showed that both erosion and accretion occurred along the shoreline, whereas there was a consistent rise in sea level at all stations along the Malaysian coastline. Results of shoreline displacement (Table 4) demonstrated that the larger part of the Malaysian coastline (about 76%) was eroding, but some coastal areas also showed an accretion pattern. An uptrend was observed in the absolute sea-level change at all stations. The highest trend of 4.46 ± 2.13 mm/y was estimated at the Sejinkat station in Sabah and Sarawak (Figure 5b), while Kukup in West Peninsular Malaysia (Figure 5a) had the lowest trend estimate of 2.63 ± 1.22 mm/y. The shoreline displacement pattern at these two tide-gauge stations was accretion (for Kukup) and erosional (for Sejinkat) at an LRR of 0.09 m/y and −0.41 m/y, respectively. LRR statistics showed that Pulau Tioman Beach, with an absolute sea-level trend of 3.72 ± 0.85 mm/y, was the most erosional, at a rate of −0.44 m/y. Cendering Beach, with an absolute sea-level trend of 3.86 ± 0.89 mm/y, was the most accretional, at a rate of 0.31 m/y. The correlation of yearly mean absolute sea-level change and shoreline displacement showed a positive relationship at all stations. Tanjung Sedili, Pulua Pinang, and Miri showed a positive correlation of 0.54, 0.43, and 0.42, respectively. Positive correlation values greater than 0.2 were also observed at the majority of coasts with erosion patterns. These correlation estimates were small but significant in the context of SLR-shoreline dynamics, considering that shoreline erosion is also connected to other energetic processes such as waves and winds. However, positive but low correlation was also found at the coasts with an accretion pattern-Lumut, Kukup, Cendering, Sandaku, and Lahad Datu. It is possible that the input of other factors, such as sediment supply, overrode the contribution of SLR in those areas. This implies that SLR contributed to shoreline erosion, but was not the major driver of shoreline displacement for this study area. Investigation of other processes and potential factors might be required to holistically explain the causes of shoreline dynamics experienced along the coastline.

Other Processes Responsible for Shifting Shorelines
The impact of the high density of human population and activities, such as the development of residential areas, industries, tourist destinations, and commercial projects, as well as construction of harbours and jetties, was observed in some coastal areas. As one of Malaysia's most popular recreational beaches [52], Teluk Nipah, for example, witnesses a high inflow of tourists annually, and is thus a region of high anthropogenic activity. The growing human activities have influenced commercial developments in the area. Figure 9 shows the development and profile of Teluk Nipah Beach. Restaurants, hotels and resorts, homestays, and hawker stalls were distributed along the coastline, but were concentrated at the mid-section, thereby exerting enormous pressure on this region. Inspection of the landscape revealed erosion in this region of intense activities and development. The erosion regions were concentrated in areas showing restaurants, homestays, resorts, hawkers stalls, and commercial shops along the coastline. The beach profile revealed a downward slope from the eroding mid-section to the upper north section of the beach. We believe that this arrangement enabled the eroded sediments from the mid-sections to be transported to the lower elevation upper section. The accretion observed at the end section was perhaps due to this longshore transportation of eroded sediments from the higher elevation mid-sections of the coast. Foo, Teh and Babatunde [52] had previously reported this contribution of sediment supply over the test site. The study conducted an analysis of the grain-size distribution of the soil sample using a sieve analysis, and discovered that the sediment in the northern region was finer than that in the southern region. It was later discovered during site studies that the southern beach's top sand had been covered by a thin layer of finer, brighter sand. The fine sand was reported to have been transported from the beach in the northern region by sediment movement along the shore. Although shoreline erosion is connected to energetic (e.g., sea-level rise, marine regime, and winds) and mass balances (e.g., the inputs/outputs of sediments), this study largely considered the impact of SLR on shoreline dynamics. In future studies, adopting a more holistic approach by considering other pertinent parameters could provide further insights into shoreline dynamics.

Future Impact of Rising Sea Level
From this study, the rate of SLR along the Malaysian coastline estimated from the average sea-level trend at all stations was 3.72 mm/y for the relative sea level and 3.68 mm/y for the absolute sea level. Based on these rates, absolute sea-level rise of about 0.11 m and 0.29 m is expected by 2050 and 2100, respectively. The relative sea level is predicted to rise by 0.12 m and 0.30 m over the same period. This is threatening, considering that 76% of the coastline already shows erosion patterns. For a high global temperature rise scenario due to increased greenhouse gas emissions (RCP 8.5), some studies based on semi-empirical models [65] and the Intergovernmental Panel on Climate Change (IPCC) [13] predicted a global rise of more than 1 m by the end of the 21st century. This estimate would cause a devastating submersion of most of the Malaysian coastline. Pulau Pinang, Port Kelang, Johor Bahru, Tanjung Sedili, Bintulu, Lanbuan 2, and Miri, where a moderate positive correlation (greater than 0.3) was observed between the sea level and shoreline change, might be most affected by this future projection of the sea level.

Conclusions
Satellites provide long-term, continuous, and quasi-global data that are vital for assessing the morpho-dynamics of shorelines; however, there are challenges in data management, processing, and analysis. By employing a Google Earth Engine-enabled Python toolkit and multi-mission optical satellites, we demonstrated an approach to automatically extracting instantaneous shorelines at high-frequency sampling towards quantifying the response of the coastlines to rising sea levels. These techniques were employed to assess 30 years of shoreline displacement at a test site and 21 coastal areas in Malaysia with tide gauges. Compared to semi-automatic and manual approaches, the approach was faster, less laborious, and required minimal storage space, and automatically extracted multiple shorelines from Landsat 5, 7, and 8 and Sentinel-2 images in a short time with high accuracy. By ensuring that shorelines were extracted for almost every month of the study period, the impact of sea-level rise on shoreline displacement could be estimated. Rising absolute sea level was discovered at all stations, but the shoreline displacement pattern showed both erosion and accretion trends, with the former more prevalent. A significant positive correlation existed between SLR and shoreline dynamics at all eroding coasts. A positive but low correlation was also found at the coasts with an accretion pattern. This confirmed that a rising sea level is contributing to shoreline erosion in the study area, but is not the only driver of shoreline displacement. Further investigation revealed that a combination of high population density, anthropogenic activities, and long shore-sediment transport are contributing to shoreline displacement over the study area. Malaysia's coastline might experience a considerable coastal retreat in the face of future estimates of SLR.