Three decades of shoreline evolution were captured from a multi-sensor longitudinal dataset consisting of four temporal snapshots from 1995, 2005, 2015 and 2025. These decadal epochs were selected to capture long-term shoreline evolution, while maintaining temporal consistency across the multi-sensor dataset. Data were accessed and processed in the Google Earth Engine (GEE) cloud-computing environment. The 1995 and 2005 epochs were constructed using Landsat 5 Thematic Mapper (TM) Surface Reflectance (SR) data, while the 2015 epoch used Landsat 8 Operational Land Imager (OLI) SR data, and the 2025 snapshot employed Sentinel-2 Multi-Spectral Instrument (MSI) Level-2A Bottom-of-Atmosphere products (Table 1). The use of decadal snapshots reflects a trade-off between temporal frequency and the need for consistent image quality, shoreline visibility and sufficient spatial detail in a tropical coastal environment where cloud cover, haze, vegetation and tidal-stage variability can complicate shoreline interpretation vegetation ((Anderson et al., 2021); (Chakma & Akter, 2021); (Dewi & Bijker, 2020); (Jayanthi et al., 2023); (Makris et al., 2023)). Intermediate years between 2015 and 2025 were not included in the main analysis because the study was designed to assess net decadal shoreline change rather than annual or event-scale variability. Therefore, the 2015–2025 period should be interpreted as one of net change, and shorter-term fluctuations within that decade may not be fully captured.

To minimize cloud contamination, candidate images for each target year were filtered in Google Earth Engine using image-level cloud-cover metadata. Landsat 5 TM and Landsat 8 OLI Surface Reflectance scenes were filtered using the CLOUD_COVER metadata property, while Sentinel-2 MSI Level-2A scenes were filtered using the CLOUDY_PIXEL_PERCENTAGE metadata property. For all sensors, only images with reported cloud cover of ≤10% were retained. This procedure was used as a scene-level screening step rather than a pixel-level quality-assessment masking procedure. All the retained images were visually inspected to ensure that the shoreline and adjacent salt-farming areas were not obscured by residual cloud, cloud shadow, haze or other atmospheric artefacts. Images with visible obstruction over the coastline or salt-pan areas were excluded from the final composite.

For each target year (see Table 1), the retained images were composited using a median reducer. The median composite was used to reduce the influence of short-term atmospheric noise, residual cloud contamination and spectral outliers. This approach produced cloud-minimized annual composites for shoreline interpretation, while maintaining a transparent and reproducible image-selection procedure across all target years. This multi-platform factor allowed the study to take advantage of both optical and radar remote sensing for consistent monitoring, which is challenging in a tropical climate such as that of the Cirebon region ((Abd‐Elhamid et al., 2023); (Gómez‐Pazo et al., 2021)). In addition, a median-reducer algorithm was implemented for the collection of images for each target year in order to yield a cloud-free composite, effectively minimizing transient atmosphere noise and outliers. To enable planimetric accuracy and ensure spatial integrity during the quantitative analysis of the region, all the spatial data were projected into the Universal Transverse Mercator (UTM) Zone 49S coordinate system.

Quantitative analysis of shoreline evolution was performed using the QGIS Shoreline Change Analysis Tool (QSCAT), which enables systematic calculation of shoreline displacement and shoreline-change statistics from multi-temporal shoreline vectors ((Facun et al., 2025); (Galan et al., 2026); (Manurung et al., 2025)). An offshore baseline was constructed following the general orientation of the Cirebon coastline, and orthogonal transects were generated at fixed 50 m intervals. In total, 668 transects were used in the analysis, with transect numbering starting from the westernmost section and increasing eastward from Transect 1 to Transect 668. For spatial interpretation, the transects were divided into three coastal segments: west coast, transects 1–223; central coast, transects 224–445; and east coast, transects 446–668. The transects were generated with a length of 3,000 m to ensure intersection with all shoreline positions across the 1995, 2005, 2015 and 2025. The baseline was smoothed using a 500 m window to reduce geometric artefacts caused by localized shoreline irregularities and to improve transect-orientation consistency. Net Shoreline Movement (NSM) and End Point Rate (EPR) were calculated for each transect to quantify the magnitude and rate of shoreline displacement over the study period.

To account for positional and processing uncertainty, an uncertainty budget was incorporated into the QSCAT workflow. Because field-based shoreline error measurements were unavailable for the historical epochs, the uncertainty values were adopted as pragmatic processing parameters rather than site-specific field-validated error estimates. A shoreline uncertainty value of 25 m and a QSCAT-specific processing uncertainty, qs_uncer, of 15 m were applied consistently across all shoreline epochs and transects. The 25 m shoreline uncertainty was considered reasonable for this multi-sensor analysis because it is broadly comparable to the spatial scale of the Landsat-based epochs and accounts for uncertainty related to sensor resolution, possible co-registration differences, mixed land-water pixels, high-water-line interpretation, and tidal-stage variability. These sources of uncertainty are commonly recognized in optical shoreline extraction and multi-temporal shoreline-change analysis ((Apostolopoulos & Nikolakopoulos, 2020); (Dewidar & Bayoumi, 2021); (Sunny et al., 2022)). The 15 m qs_uncer parameter was used to represent additional uncertainty associated with shoreline delineation and transect-intersection processes during QSCAT calculation. These values should not be interpreted as field-validated positional errors, but as conservative and consistent uncertainty settings used to support cautious interpretation of shoreline-change statistics, particularly for low-magnitude shoreline changes that may be affected by image resolution, shoreline-proxy ambiguity and tidal-stage variability.

The geostatistical outputs were then grouped according to their morphodynamical trends and relevant risks to assess the implications of shoreline evolution for salt-production feasibility. A classification framework of four regions is proposed based on global vector models of Shoreline Change Envelope (SCE), Net Shoreline Movement (NSM) and End Point Rate (EPR), adjusted accordingly. NSM is signed (negative = retreat; positive = progradation) and SCE is unsigned in this case. Therefore, analyses linking SCE with cumulative displacement employed |NSM| in comparisons of magnitude and signed NSM in classifications of erosion–accretion. This transformation was made only to compare magnitudes with SCE; signed NSM was preserved to facilitate directional interpretation across results. Moreover, landward accretionary trends were grouped into three classes: stable accretion rate, moderate accretion rate, and high chance for the amelioration of a high chance of improved progress. Using these categories, we assessed the “coastal distance” challenge of excessive land growth that potentially decouples salt-pan infrastructure from seawater intake over time in a timely fashion. In contrast, zones showing negative long-term displacement or landward migration were categorized as at risk of inundation. In these high-vulnerability areas, Weighted Linear Regression (WLR) was selected as the key indicator to derive trends for at-risk land submergence over the long term, and for permanent loss of land. All rates were assessed at a 95% confidence interval, meaning that resultant retreat/accretion rates across the Cirebon salt-production landscape are significant and represent true geomorphological trends. Notably, this coupled assessment enabled the linkage of changes to physical shorelines with threats to operational stability in the Cirebon salt-farming sector and thus provided a metric for determining which regions may require urgent adaptive action. The overall research workflow is summarized in Figure 2.

Shoreline-change assessment in Cirebon is based on the standard shoreline-cover metrics generated from QGIS-QSCAT transect export (n = 668 transects): End Point Rate (EPR, m/year), Linear Regression Rate (LRR, m/year), Net Shoreline Movement (NSM, m) and Shoreline Change Envelope (SCE, m). Overall spatial pattern is directionally mixed (again, accretion dominates erosion by transect count), with highly uneven magnitudes across the coastal expanse; this heterogeneity is reflected in the divergence between mean and median statistics, in which a few extreme erosion sites bring the average down, masking the typical accretionary behavior seen across most of the region. EPR in this research also varies between −81.86 m/year and +59.14 m/year (Figure 3a), while NSM also varies from −2476.28 m to +1789.01 m (shown in Table 3 and Figure 3b). The median ensemble planform retreat (EPR) is positive (+1.83 m/year), as is the median net-shoreline migration (NSM) (+55.33 m), indicating that typical transects undergo seaward displacement, whereas the mean EPR (−0.97 m/year) and mean NSM (−29.38 m) are slightly negative due to a smaller number of extreme erosional segments having disproportionally strong influence. LRR displays a similar distributional structure (median +1.89 m/year; mean −0.78 m/year) (presented in Figure 3c), suggestive of bulk accretionary behavior with concentrated hotspots for ephemeral erosion.

The SCE, as a non-directional measure of variability in the position of the shoreline, ranges from 10.01 m to 2,476.28 m, with a median value of 170.31 m, and the large upper tail indicating that certain sectors deployed wide positional envelopes throughout the observation window, resulting in pronounced morphodynamic variability (see Figure 3d). In order to interpret how this positional envelope relates to net change, SCE was compared on a per-transect basis against the cumulative displacement magnitude (|NSM|), with SCE unsigned and NSM signed. In relation to the three alongshore divisions created (West: transects 1–223; Central: 224–445; East: 446–668) the relationship between SCE and |NSM| remains consistently positive and strong across each sector (r ≈ 0.83 in both West and Central sectors, r ≈ 0.97 in East sector), implying that over displaced transects, larger amounts of shoreline envelopes are also generally representative, with strongest coupling within the eastern segment, where the extreme values of displacement are concentrated (see Figure 4). Another factor raised with respect to interpretation is that of uncertainty handling; the dataset indicates that all transects possess an EPR_unc = 0.70 m/year (for all transects), suggesting that stability may have been defined using some form of uncertainty-based neutral band around zero.

The QSCAT trend classification indicates that shoreline accretion is more widespread than erosion when assessed by transect count. Based on the End Point Rate trend classification (EPR_trend), 421 of the 668 transects, or 63.0%, were classified as accreting, while 192, or 28.7%, were eroding, and 55, or 8.2%, were stable. A similar pattern was obtained from the Net Shoreline Movement trend classification (NSM_trend), in which 430 transects, or 64.4%, were accreting, 204, or 30.5%, were eroding, and 34, or 5.1%, were stable (Table 4). The consistency between the EPR and NSM_trends suggests that the dominance of accreting transects is not merely an artefact of a single shoreline-change metric; rather, both indicators show that much of the Cirebon coastline experienced seaward shoreline movement or relative stability during the observation period. Nevertheless, this count-based dominance should not be read as evidence that erosion is negligible, as transect frequency does not capture the magnitude of shoreline retreat at specific locations.

However, the dominance of accretion by transect count must be interpreted together with the magnitude and spatial distribution of both the EPR and NSM values. Although most transects were classified as accreting, the magnitude distribution was skewed by a smaller number of high-erosion ones. This pattern is reflected in the coexistence of positive median values with negative mean values, particularly for NSM, where the median remained positive at +55.33 m, while the mean was negative at −29.38 m. This indicates that the median better represents the more common condition across the coastline, whereas the mean is strongly influenced by localized retreat segments with extreme values. Therefore, negative mean values should not be interpreted as evidence of uniform erosion along the entire Cirebon coast, but rather as an indication of spatially uneven shoreline dynamics and localized erosion hotspots.

Hotspot identification based on EPR further shows that the highest accretion rates occurred at Transects 657 (+59.14 m/year), 387 (+51.96 m/year), 668 (+50.87 m/year), 658 (+49.49 m/year) and 369 (+47.66 m/year), while the highest erosion rates were found at Transects 651 (−81.86 m/year), 650 (−81.20 m/year), 619 (−78.19 m/year), 640 (−77.92 m/year) and 630 (−76.28 m/year) (as presented in Figure 5). The NSM range, from −2476.28 m to +1789.01 m, confirms that the shoreline-change distribution contains a limited number of high-magnitude retreat and advance values. These extreme transects may represent genuine geomorphic hotspots, especially where they form spatially coherent clusters, but they should also be interpreted cautiously because satellite-derived shoreline proxies can be affected by shoreline-extraction uncertainty, tidal-stage variability, mixed land-water pixels, salt-pan inundation, or localized spectral noise. Therefore, the interpretation emphasizes clustered erosion and accretion patterns shown in the segmentation map rather than isolated transect-level extremes.

QSCAT outputs, including EPR, NSM and SCE, were interpreted using a category-based framework to translate transect-scale shoreline-change metrics into salt-farming exposure classes. Category 1 represents stable or low-change segments, in which EPR falls within the uncertainty-defined neutral band, EPR_unc = 0.70 m/year, and NSM and SCE values remain comparatively low. These segments are interpreted as having relatively high shoreline predictability and lower direct exposure for salt-pan infrastructure. Category 2 represents moderate accretion or seaward displacement beyond the stability threshold, but below the highest observed magnitudes. This condition may indicate gradual shoreline progradation that could require periodic adjustment of seawater-intake alignment, access routes or pond boundaries. Category 3 represents high accretion and/or high-envelope behavior, indicated by strongly positive EPR or NSM values and large SCE magnitudes. Although accretion is often interpreted as land gain, in salt-farming systems excessive accretion can be operationally problematic because it may increase sedimentation, obstruct seawater-intake channels, alter drainage pathways, reduce hydraulic connectivity between ponds and the sea, and require repeated channel maintenance or pond-boundary adjustment.

Category 4 identifies the highest-risk shoreline-instability segments, characterized by landward displacement, negative EPR and/or NSM values, and/or very large SCE values. Because high-resolution elevation, tidal-range and field-based inundation data were unavailable, Category 4 is interpreted as a shoreline-instability exposure proxy rather than a direct inundation-risk zone. Nevertheless, this category is considered highly relevant for salt-pan vulnerability because shoreline retreat may reduce production space, damage pond margins, increase stress on bunds and access routes, and expose intake or drainage infrastructure to more frequent disturbance. Sector-level summaries indicate that the eastern sector, Transects 446–668, contains the most extreme displacement magnitudes, including both cumulative retreat and rapid progradation, while the western and central sectors show more positive medians and less severe extremes. Therefore, priority areas for higher-resolution salt-pan mapping are those that spatially intersect or immediately adjoin Category 4 transects, followed by Category 3 transects, where rapid accretion or shoreline reconfiguration may impose repeated operational adjustment.

To make the exposure framework less subjective, the expected severity classes were anchored to measurable shoreline-change indicators and operational proxies rather than qualitative judgement alone. These indicators include EPR direction and magnitude, NSM displacement, SCE envelope size, transect location within west-central-east sectors, and spatial proximity or overlap with salt-farming parcels. In this framework, Category 4 is prioritized when negative EPR/NSM and large SCE indicate shoreline retreat or instability near salt-pan infrastructure, while Category 3 is prioritized when high positive EPR/NSM and large SCE suggest rapid accretion that may affect seawater intake, drainage, sediment management and pond access. Because direct data on maintenance costs, intervention frequency or production loss were not available, the severity ranking should be interpreted as an exposure-based prioritization tool rather than a direct economic damage estimate. Table 5 summarizes the expected impact pathways and measurable indicators used to identify potential vulnerability zones where shoreline dynamics intersect with salt-farming infrastructure.

The results from Cirebon confirm a pattern widely reported for dynamic tropical and Global South coastlines: shoreline behavior is directionally mixed and spatially clustered, where an accretion-dominant transect count can coexist with localized but severe retreat hotspots that strongly influence magnitude-based interpretation. This pattern explains why the mean values may become slightly negative despite positive median values, indicating that a smaller number of high-magnitude erosion transects exert a strong influence on coastline-wide statistics. Similar hotspot-based shoreline dynamics have been reported in multi-decadal studies, with erosion and accretion being shaped by interacting hydro-oceanographic and anthropogenic drivers rather than by uniform shore-wide trends ((Gallina et al., 2020); (Jayanthi et al., 2023); (Mageswaran et al., 2021); (Pantusa et al., 2022)). The alongshore segmentation in this study further shows that the eastern part of the Cirebon coastline contains the highest extremes, supporting the interpretation that shoreline-change “cells” and modified sediment systems can create adjoining zones of contrasting erosion and accretion ((Dewi & Bijker, 2020); (Pantusa et al., 2022); (Zoysa et al., 2023)). When benchmarked against other Indonesian and Southeast Asian coastal studies, this pattern is not unusual. In Quang Nam, Vietnam, DSAS-based analysis found both erosional and accretional transects over the period 1990–2019, with strong localized erosion reaching −42.4 m/year near the Cua Dai estuary (Quang et al., 2021). In East Java, Indonesia, multi-sensor remote-sensing analysis showed that accretion was generally more pronounced than erosion, but some deltaic and aquaculture-mangrove areas experienced very large local shoreline changes, with EPR values ranging from strong erosion to very high accretion (Arjasakusuma et al., 2021). Recent work in the Upper Gulf of Thailand also shows that erosion and accretion can vary strongly across periods and regions due to sediment transport, sea-level rise, monsoonal processes, coastal protection, mangrove restoration and human activity (Chawalit et al., 2025). These comparisons suggest that Cirebon should not be interpreted as a uniformly eroding or accreting coast, but as a spatially heterogeneous coastal system where localized hotspots are especially important for salt-farming exposure.

In the context of salt farming, the main implication is that operational exposure may arise from a limited number of shoreline sectors, even when the broader coastline appears accretion-dominant by transect count. This is consistent with empirical work from northern Java which showed that coastal hazards and shoreline instability could place operational stress on salt-production landscapes, including pond-boundary disturbance, access disruption and changes in management routines ((Nirwansyah & Braun, 2019), (Nirwansyah & Braun, 2021); (Sutrisno et al., 2021)). However, the impact framework developed here should be interpreted as a remote-sensing-based screening tool rather than direct evidence of confirmed salt-pan damage. The categories identify where shoreline dynamics may intersect with salt-farming infrastructure, but they do not directly measure embankment failure, yield reduction, repair costs, maintenance frequency or impacts on household livelihoods. This distinction is important because field validation, farmer interviews, production records and historical photographs of affected salt pans were not available for the full multi-decadal period. Therefore, the most feasible interpretation is that EPR, NSM and SCE identify potential exposure pathways between shoreline dynamics and salt-pan infrastructure, but not confirmed economic loss or observed operational failure.

Thematically, the study contributes an operational workflow that links cloud-based multi-decadal shoreline extraction and preprocessing with QSCAT transect statistics, establishing a repeatable basis for shoreline-change assessment in settings where consistent long-term monitoring is limited ((Gómez‐Pazo et al., 2021); (Tang et al., 2023); (Yasir et al., 2020)). A second contribution is the use of multiple complementary metrics rather than reliance on a single rate indicator. The combined interpretation of EPR, NSM and SCE allows the analysis to distinguish between shoreline direction, cumulative displacement, and positional variability.

This is important for salt farming because large-scale cumulative displacement and wide shoreline envelopes may indicate locations where pond boundaries, intake channels, drainage routes and access paths are more likely to require adjustment. Similar studies in Southeast Asia also demonstrate the value of using multi-temporal satellite imagery and shoreline-change tools such as DSAS/QSCAT to quantify long-term coastal dynamics and identify hotspots for coastal management ((Abd‐Elhamid et al., 2023); (Chawalit et al., 2025); (Pantusa et al., 2022); (Quang et al., 2021)). However, this study differs from many previous shoreline-monitoring ones by translating physical shoreline-change metrics into salt-farming-relevant exposure categories ((Pantusa et al., 2022); (Tang et al., 2023)). This helps address a gap in the literature, where shoreline change is often mapped as a physical process, but is less frequently connected to livelihood-specific infrastructure such as salt pans, bunds, seawater-intake channels and access routes.

Several limitations and uncertainties should guide interpretation of the study findings. First, because high-resolution elevation, complete tidal-stage information, and field-based inundation data were not available, high-change categories should be interpreted as shoreline-instability exposure proxies rather than direct inundation-risk classes. This restricts conclusions about flooding frequency, overtopping thresholds or depth of inundation in salt pans. Second, as in other multi-decadal optical shoreline studies, the results are sensitive to shoreline-proxy definition, water-level differences between image dates, sensor resolution, georeferencing uncertainty, and mixed land-water pixels. These factors may inflate apparent variability when shoreline changes are small, which is why this study uses uncertainty-aware interpretation and emphasizes spatially coherent patterns rather than isolated transect-level patterns ones ((Abd‐Elhamid et al., 2023); (Gómez‐Pazo et al., 2021); (Jayanthi et al., 2023); (Pantusa et al., 2022)).

Third, the selected epochs of 1995, 2005, 2015 and 2025 were designed to capture net decadal shoreline change, not annual or event-scale variability. The 2015–2025 period should therefore be understood as one of net-change; shorter-term erosion, accretion, storm-driven change or temporary recovery within that decade may not be fully captured. Fourth, while the analysis identifies where shoreline dynamics are strongest, it does not quantify salt-production loss, repair costs, intervention frequency or household-level impacts. Future work should combine category-coded transects with parcel-level salt-pan boundaries, UAV or drone surveys, GPS mapping of embankments and intake channels, participatory mapping with salt farmers, farmer interviews, production records, and maintenance-cost documentation. This would allow related research to move from identifying potential exposure towards measuring actual operational and socioeconomic impacts.