Fuzzy Classification for Shoreline Change Monitoring in a Part of the Northern Coastal Area of Java, Indonesia

This study presents an unsupervised fuzzy c-means classification (FCM) to observe the shoreline positions. We combined crisp and fuzzy methods for change detection. We addressed two perspectives of uncertainty: (1) uncertainty that is inherent to shoreline positions as observed from remote sensing images due to its continuous variation over time; and (2) the uncertainty of the change results propagating from object extraction and implementation of shoreline change detection method. Unsupervised FCM achieved the highest kappa (κ) value when threshold (t) was set at 0.5. The highest κ values were 0.96 for the 1994 image. For images in 2013, 2014 and 2015, the κ values were 0.95. Further, images in 2003, 2002 and 2000 obtained 0.93, 0.90 and 0.86, respectively. Gradual and abrupt changes were observed, as well as a measure of change uncertainty for the observed objects at the pixel level. These could be associated with inundations from 1994 to 2015 at the northern coastal area of Java, Indonesia. The largest coastal inundations in terms of area occurred between 1994 and 2000, when 739 ha changed from non-water and shoreline to water and in 2003–2013 for 200 ha. Changes from water and shoreline to non-water occurred between 2000 and 2002 (186 ha) and in 2013–2014 (65 ha). Urban development in flood-prone areas resulted in an increase of flood hazards including inundation and erosion leading to the changes of shoreline position. The proposed methods provided an effective way to present shoreline as a line and as a margin with fuzzy boundary and its associated change uncertainty. Shoreline mapping and monitoring is crucial to understand the spatial distribution of coastal inundation including its trend.


Introduction
Shoreline location and its change in position are critical for understanding coastal structures [1], safe navigation [2,3], sustainable coastal resource management [4], and for flood protection and other risk management [5,6]. Furthermore, shoreline mapping and monitoring can help to understand the spatial distribution of coastal inundation including its trend over time.
In the literature, shoreline is defined as an intersection of coastal land and water surface indicating water edge movements as the tides rise and fall [7][8][9]. Even though shoreline position can be defined as the waterline at various stages of the tides, e.g., high tide, mid tide, and low tide, the shoreline is largely associated with the sea level [10]. Ideally, the shoreline is the physical interface of land and water with its position changing through time (Doland et al. (1980) in Boak and Turner [8]). Its change  The area is characterized by a flat topography ranging between 0 and 5 m above mean sea level (AMSL). As a lowland coastal area, it is dominated by alluvial clay and sand sedimentation with more than 10 rivers running across the area and carrying sediments from the upstream areas [32,33]. The poor drainage system leads to a regular occurrence of floods during the rainy season. Meanwhile, the areas adjacent to the sea are prone to the impact of tidal flood which occurs daily in line with the Remote Sens. 2016, 8,190 3 of 24

Study Area
This study was carried out in a part of the Sayung sub-district in the northern coastal area of Central Java Province extending 5 km from east to west (Figure 1). The area has faced a severe impact of coastal inundation which has led to a significant change of shoreline position over the past two decades ( Figure 2). The average tidal range in this location is 60 cm with the highest tidal ranges occurring in December and June during the rainy and dry seasons, and the lowest tidal ranges in March and September during the transitional seasons.  The area is characterized by a flat topography ranging between 0 and 5 m above mean sea level (AMSL). As a lowland coastal area, it is dominated by alluvial clay and sand sedimentation with more than 10 rivers running across the area and carrying sediments from the upstream areas [32,33]. The poor drainage system leads to a regular occurrence of floods during the rainy season. Meanwhile, the areas adjacent to the sea are prone to the impact of tidal flood which occurs daily in line with the tides [34,35]. In this study area, four villages are located amidst extensive fishponds and rice fields. Settlements are found along the riverbanks or adjacent to the shorelines. Since the 1990s, the productive rice fields became inundated and were converted into fishponds, or merely abandoned as swamp areas [32]. The area is characterized by a flat topography ranging between 0 and 5 m above mean sea level (AMSL). As a lowland coastal area, it is dominated by alluvial clay and sand sedimentation with more than 10 rivers running across the area and carrying sediments from the upstream areas [32,33]. The poor drainage system leads to a regular occurrence of floods during the rainy season. Meanwhile, An extensive change of shoreline has been occurring for more than 20 years (Figure 3) caused by natural processes, e.g., coastal inundation, erosion, and sedimentation and by the development of man-made structures, e.g., beach reclamation and extended seaport. Moreover, the coastal inundation accelerating these changes has increased recently in terms of frequency and duration as a result of many factors. These include at the short term: extreme winds, heavy rains, and at the long term: land subsidence, sea level rise, mangrove conversion, groundwater extraction, construction load, and the increase of impervious surface [34,[36][37][38]. Many efforts have been made by the Demak local government to overcome the erosion of 80 km 2 of land [39] and curb the coastal inundation which leads to the dramatic changes of the shoreline. Efforts include embankment, elevated road, breakwater, and mangrove planting.
Pre-processing consisted of histogram minimum adjustment to reduce the effect of atmospheric path radiance [40,41], followed by geo-referencing using >100 ground control points (GCP) that were carefully selected on both Landsat and ortho-rectified WorldView-2 images from road intersections, building corners, wall boundaries, river and other prominent features, and re-sampling to a 30 m pixel size using the nearest neighbour resampling method and third order polynomial transform algorithm. The root mean square error (RMSE) was less than 0.1 pixels. (a-d) The comparison of normal (a,c) and flooded (b,d) situations due to coastal inundation at two locations at Sayung sub-district. Over a longer period, this cyclic flood leads to a permanent inundation.

Satellite Images and Data Pre-Processing
Landsat images from three different sensors were used to monitor the shoreline change between 1994 and 2015. Landsat has a 16-day revisit time, and passes Indonesia at approximately 02.00-03.00 GMT. Six spectral bands of Landsat TM and ETM and seven spectral bands of Landsat OLI/TIRS were used. Spectral bands of Landsat TM and ETM applied in this research cover the blue (0.45-0.515 µm), green (0.525-0.605 µm), red (0.63-0.69 µm), near infrared (0.75-0.90 µm), shortwave infrared 1 (1.55-1.75 µm) and shortwave infrared 2 (2.09-2.35 µm) parts of the electromagnetic spectrum. In addition, the spectral bands of Landsat OLI/TIRS included in FCM consisted of coastal and aerosol (0.43-0.45 µm), blue (0.45-0.51 µm), green (0.53-0.59 µm), red (0.64-0.67 µm), near infrared (0.85-0.88 µm), shortwave infrared 1 (1.57-1.65 µm) and shortwave infrared 2 (2.11-2.29 µm). Table 1 shows the images used in this study supplemented by tidal data. Tidal data in accordance with the time of acquisition of the images was collected from the Indonesia Geospatial Information Agency. Pre-processing consisted of histogram minimum adjustment to reduce the effect of atmospheric path radiance [40,41], followed by geo-referencing using >100 ground control points (GCP) that were carefully selected on both Landsat and ortho-rectified WorldView-2 images from road intersections, building corners, wall boundaries, river and other prominent features, and re-sampling to a 30 m pixel size using the nearest neighbour resampling method and third order polynomial transform algorithm. The root mean square error (RMSE) was less than 0.1 pixels.

Fuzzy C-Means (FCM) Classification and Parameter Estimation
The FCM iterative clustering method developed by Bezdek, et al. [31] was performed on the images to discriminate the land and water classes. FCM separates data clusters with fuzzy means and fuzzy boundaries allowing for partial membership. Let Y " ty 1 , y 2 , . . . , y N u be a sample of the N pixels on an image, with y k R n where n is the number of bands in an image, i.e., n = 6 or n = 7 for Landsat images used here. Let c denote the number of subsets (clusters or partitions) with 2 ď c ď N. In this research, we have c = 2 for the classes water and non-water, respectively, since we considered the boundary between these two classes as the shoreline position. FCM minimizes the following objective function J m [31]: where µ ik is the membership value of k th pixel to class i, m is the fuzzy weight controlling the level of fuzziness, and v i " pv i1 , v i2 , . .., v in q is the mean vector for class i. The membership value µ ik for class i and pixel k satisfies the following constraints: There should be at least one class for which the membership value µ ik of the k th pixel larger than 0. Meanwhile, the sum of all memberships µ ik in a pixel should be equal to 1. The membership values of the classification corresponding to the pixel value Y follow the trapezoidal membership function in Figure 4 and Equation (5): of fuzziness, and = ( , , . . . , ) is the mean vector for class i. The membership value for class i and pixel k satisfies the following constraints: 0 ≤ ≤ 1 ∈ 1, … . , , ∈ 1, … , > 0, ∈ 1, … , There should be at least one class for which the membership value of the pixel larger than 0. Meanwhile, the sum of all memberships in a pixel should be equal to 1. The membership values of the classification corresponding to the pixel value Y follow the trapezoidal membership function in Figure 4 and Equation (5):  If m equals 1, clusters that minimize the objective function are hard clusters and FCM is a hard classifier. An increase of m tends to an increase in fuzziness. Bezdek, et al. [31] further explained that no evidence distinguishes an optimal m, but for most data, 1.5 ď m ď 3.0 give good results. In addition, Foody [42] stated that in most studies, m = 2.0 produces an accurate fuzzy classification. In this work, the values m = 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.5, and 3.0 were used to test the influence of m on the classification results.
FCM was finalized by labelling one of the two membership images resulting from each (unsupervised) FCM as the water membership image. To do so, we used the combination of near infrared (NIR) and shortwave infrared (SWIR) of Landsat bands. Infrared bands exhibit a strong contrast between water and land features, because water absorbs these wavelengths while they are reflected by land [43]. In the visible part of the spectrum, the differences between land and water are less outspoken, especially if the water contains some sediment. Therefore, the water label was given to the class which has the lowest value of the sum of the cluster means in the infrared bands, defined as: where vi c1 is the sum of mean vector for the first class in the infrared bands IR, and vi c2 is the sum of mean vector for the second class in the infrared bands IR. Table 2 shows an example of cluster means of each subset in the infrared bands. The labelling of c1 or c2 as water was decided by using Equation (6). In this example, the water label was given to c1 as it has the lowest value of the sum of the cluster means in the infrared bands.

Deriving Water Class Images
Water membership images, resulting from the FCM classification, show the membership of pixels to the water class. Thresholding was applied to transform the water membership image into several hard classifications. The class C in water class images was defined as: where 1 is water class, 0 is non-water class, and t is threshold value. The possible ranges of the threshold values are between 0 and 1. In this study, we set the value of t between 0.1 and 0.9 to observe the influence of t on the results of classification. The results of thresholding were binary images called water class images.

Accuracy Assessment
Reference data are described in Table 1. For the 2015 image, reference data were derived from fieldwork conducted at the end of March 2015. In this case, 150 points from fieldwork data were selected based on the same tide condition. Furthermore, an image made available via Google Earth 2014 captured during the high tide, a 2013 Pleiades image during low tide, and a 2003 QuickBird image captured during the rising tide were interpreted visually to generate 150 reference points. Further, because of limited availability of high-resolution images and maps, the 2003 QuickBird image was used as well for accuracy assessment of images in 2000 and 2002. Finally, the topographic map published in 2000 and generated from aerial photographs of 1994 was manually digitized on-screen to produce water and non-water classes. For this map, 150 points were randomly selected, and were used as reference data against the classification result of the 1994 image. To evaluate how well the remotely sensed classifications agree with the reference data, error matrices were generated. A kappa (κ) coefficient was computed for each error matrix [43].

Shoreline Generation
Two methods were followed to identify the shoreline. The first method determined shoreline as a single line, as has been widely considered in the previous studies, whereas the second method assumed shoreline as a margin, which reflects the possible locations of shoreline influenced by the membership to the water class in a pixel.

Shoreline as a Single Line
First, the shoreline was derived by generating water class image and set t = 0.5. Two sub-areas were identified, namely water and non-water. Shoreline was located at the boundary of the two sub-areas and obtained by converting the water class image to line features in GIS. For this research, a sub-area was defined as a set of contiguous cells with the same value.

Shoreline as a Margin
Secondly, we considered the shoreline as an area (margin). The shoreline margin was generated by creating a crisp sub-area determined by t = 0.3 and t = 0.7 as the lower and upper thresholds obtained in the parameter estimation. Afterwards, each water class image was converted into polygon feature in GIS. Thus, three sub-areas were identified, namely water if µ ik ě 0.7, non-water if µ ik ă 0.3 and shoreline margin if 0.3 ď µ ik ă 0.7. Over time, the changes of shoreline position are due to the exchanges between the shoreline margin and water or non-water sub-area.

Uncertainty Estimation
Considering the imprecise position of the shoreline in remote sensing images and the uncertainty propagated through the change detection method, shoreline margin, water and non-water sub-areas were associated with values reflecting the uncertainty of pixels belonging to any of these classes. Water membership values were used for this purpose. The certainty of pixel k to belong to any class was assessed using possibility and necessity measures [44][45][46][47]. If C is a subset of a universe of discourse U and x P C, then π x puq is the degree of possibility that x takes value u. The value of π x puq is evaluated by the degree of membership µ C puq. This can be written as: The value of x P C is then estimated by assessing possibility measure: The possibility measure Π C corresponds to the element of C that has the highest possibility degree according to π x . Further, to inform that the event will be realized, the certainty of C is defined as the impossibility of the complement: The N pixels in Y are therefore indicated as C if Π C ą Π C and N C ą N C . Further, the uncertainty of C is defined as:

Shoreline Change Detection
To detect the changes in the positions of the shoreline, shoreline margin, water and non-water areas, results for two dates T 1 and T 2 had to be superimposed in GIS. In order to have a detailed "from-to" change trajectory information, the post-classification comparison approach was used [48]. Topological relation between two sub-areas can be characterized by considering the nine-intersection model of interiors and exteriors [49,50]. Based on this method, a sub-area identified at time T 1 is denoted as R T1 , with the boundary BpR T1 q, interior I pR T1 q, and exterior E pR T1 q. It intersects with another sub-area identified in time T 2 and denoted as R T2 , with the boundary BpR T2 q, interior I pR T2 q, and exterior E pR T2 q. These intersections define the nine-intersection matrix as: with intersections being either empty p∅q or non-empty p ∅q. Figure 5 shows eight topological relationships of two sub-areas for each intersection value in the matrix M including disjoint, meet, overlap, contains, inside, covers, covered by, and equal [49]. Following the aforementioned methods in the shoreline generation, the changes of shoreline and its change uncertainty can be presented as follows: with intersections being either empty (∅) or non-empty (¬∅). Figure 5 shows eight topological relationships of two sub-areas for each intersection value in the matrix M including disjoint, meet, overlap, contains, inside, covers, covered by, and equal [49]. Following the aforementioned methods in the shoreline generation, the changes of shoreline and its change uncertainty can be presented as follows:

Shoreline as a Single Line
In the first method, the changes of this single shoreline may have occurred as a consequence of changes between water and non-water sub-areas. We determined the changes of these sub-areas at times and and analyzed the uncertainty of the changes at pixel level. Water class images from two dates and were superimposed, and abrupt and gradual changes were identified. An abrupt change is defined when a sub-area emerges at without a corresponding sub-area at . Also, a sub-area present at without a corresponding sub-area at indicates an abrupt change. A gradual change specifies an increase or decrease of sub-areas that were both present at and . The process occurs in small stages over a period rather than suddenly. The overlay analysis between images at different epochs permitted the identification of changes categorized as: water to non-water, and non-water to water. Considering the topological relationship as given in Figure 5, disjoint and meet were found where a sub-area emerges or disappears. Meanwhile, any of the other six relationships account for gradual changes.

Shoreline as a Single Line
In the first method, the changes of this single shoreline may have occurred as a consequence of changes between water and non-water sub-areas. We determined the changes of these sub-areas at times T 1 and T 2 and analyzed the uncertainty of the changes at pixel level. Water class images from two dates T 1 and T 2 were superimposed, and abrupt and gradual changes were identified. An abrupt change is defined when a sub-area emerges at T 2 without a corresponding sub-area at T 1 . Also, a sub-area present at T 1 without a corresponding sub-area at T 2 indicates an abrupt change. A gradual change specifies an increase or decrease of sub-areas that were both present at T 1 and T 2 . The process occurs in small stages over a period rather than suddenly. The overlay analysis between images at different epochs permitted the identification of changes categorized as: water to non-water, and non-water to water. Considering the topological relationship as given in Figure 5, disjoint and meet were found where a sub-area emerges or disappears. Meanwhile, any of the other six relationships account for gradual changes.

Shoreline as a Margin
The second method measured the changes of the shoreline margin at different epochs. Change uncertainty was presented at the pixel level. Water class images from times T 1 and T 2 were superimposed in a GIS. We distinguished again abrupt and gradual change. In this method, the overlay analysis between water class images at different epochs, however, permitted more changes to be identified categorized as: (1) shoreline to non-water; (2) water to shoreline; (3) water to non-water; (4) non-water to shoreline; (5) shoreline to water; and (6) non-water to water. The changes of a sub-area where sedimentation has taken place resulted in changes of shoreline to non-water, water to shoreline, and water to non-water. These types of changes were considered positive changes to non-water sub-area (+). Coastal inundation led to the changes of non-water to shoreline, shoreline to water, and non-water to water. These types of changes were considered negative changes to non-water sub-area (´). These changes were identified either as abrupt or gradual changes, using the same criteria on corresponding objects as in the previous section.

Change Uncertainty and Change Area Estimation
Shoreline changes and the related sub-areas were associated with change uncertainty values. Change uncertainty was derived based upon the uncertainty of pixels belonging to the specified sub-areas in T 1 and T 2 [51]: where U c pR T1 q is the uncertainty of pixel k belonging to sub-area R T1 at T 1 , and U c pR T2 q is the uncertainty of pixel k belonging to sub-area R T2 at T 2 . Based upon the results of the change uncertainty estimation, the changed area of a specific change category, e.g., water to non-water, or shoreline to water was defined as: A pChq " P k pChqˆA pkq (15) where P k pChq is number of pixels belonging to the area of a specific change category, A pkq is area of pixel k equal to 30ˆ30 (m 2 ).
2.9.1. Change Area of Shoreline as a Single Line Figure 6 illustrates the procedure to estimate the net change between T 1 and T 2 defined as: where A and B represent the area described in Figure 6, and P k pCh A q and P k pCh B q are the number of pixels belonging to the change area A and B, respectively.
water to water. These types of changes were considered negative changes to non-water sub-area (−). These changes were identified either as abrupt or gradual changes, using the same criteria on corresponding objects as in the previous section.

Change Uncertainty and Change Area Estimation
Shoreline changes and the related sub-areas were associated with change uncertainty values. Change uncertainty was derived based upon the uncertainty of pixels belonging to the specified subareas in and [51]: where ( ) is the uncertainty of pixel k belonging to sub-area at , and ( ) is the uncertainty of pixel k belonging to sub-area at . Based upon the results of the change uncertainty estimation, the changed area of a specific change category, e.g., water to non-water, or shoreline to water was defined as: where ( ℎ) is number of pixels belonging to the area of a specific change category, ( ) is area of pixel k equal to 30 × 30 (m 2 ).
2.9.1. Change Area of Shoreline as a Single Line Figure 6 illustrates the procedure to estimate the net change between and defined as: where A and B represent the area described in Figure 6, and ( ℎ ) and ( ℎ ) are the number of pixels belonging to the change area A and B, respectively. The negative sign (´) shows that the change in B categorized as the change from non-water to water has produced a negative change to the non-water area. Meanwhile, the positive sign (+) represents the change in A categorized as the change from water to non-water which has caused a positive change to the non-water area. Figure 7 illustrates the procedure to estimate the total changed area in the second approach. The net changed area was determined as:

Change Area of the Shoreline as a Margin
where A, B, C, D, E, and F represent the area described in Figure 7, and P k pCh A q , P k pCh B q , P k pCh C q , P k pCh D q , P k pCh E q and P k pCh F q are the number of pixels belonging to changed areas A, . . . ,F, respectively.
represent shoreline at whereas dashed lines refer to shoreline at .
The negative sign (−) shows that the change in B categorized as the change from non-water to water has produced a negative change to the non-water area. Meanwhile, the positive sign (+) represents the change in A categorized as the change from water to non-water which has caused a positive change to the non-water area. Figure 7 illustrates the procedure to estimate the total changed area in the second approach. The net changed area was determined as:

Change Area of the Shoreline as a Margin
where A, B, C, D, E, and F represent the area described in Figure 7, and ( ℎ ), ( ℎ ), ( ℎ ), ( ℎ ), ( ℎ ) and ( ℎ ) are the number of pixels belonging to changed areas A,…,F, respectively. The positive sign (+) indicates that the changes in A, B and C have caused a positive change to the shoreline margin and non-water areas. On the other hand, the negative sign (−) shows that the changes in D, E, and F have induced a negative change to the shoreline margin and non-water areas.  The positive sign (+) indicates that the changes in A, B and C have caused a positive change to the shoreline margin and non-water areas. On the other hand, the negative sign (´) shows that the changes in D, E, and F have induced a negative change to the shoreline margin and non-water areas.          Figure 9 shows FCM results for m = 1.7 presenting the membership to the water class ranging from 0 to 1, together with classified images for t = 0.5. Areas with a high membership to the water class include marine areas; e.g., Figure 9b grid cell A2, fishponds; e.g., Figure 9n grid cell D3, and water-covered agricultural areas; e.g., Figure 9j grid cell E5. Muddy areas are located on the border of water and non-water; e.g., Figure 9k grid cell A4. Further, the shrinking of non-water areas over two decades could also be distinguished. This can be seen by the change of the shape of the non-water class from wide strips to the thin elongated shapes over the series of images in Figure 9a-n; e.g., grid cells C3. On the other hand, non-water sub-areas emerged in several locations, such as when mangroves were planted; e.g., see Figure 9i grid cells C2, and in reclamation areas; e.g., see Figure 9a,c grid cells A5.   Figure 9 shows FCM results for m = 1.7 presenting the membership to the water class ranging from 0 to 1, together with classified images for t = 0.5. Areas with a high membership to the water class include marine areas; e.g., Figure 9b grid cell A2, fishponds; e.g., Figure 9n grid cell D3, and water-covered agricultural areas; e.g., Figure 9j grid cell E5. Muddy areas are located on the border of water and non-water; e.g., Figure 9k grid cell A4. Further, the shrinking of non-water areas over two decades could also be distinguished. This can be seen by the change of the shape of the non-water class from wide strips to the thin elongated shapes over the series of images in Figure 9a-n; e.g., grid cells C3. On the other hand, non-water sub-areas emerged in several locations, such as when mangroves were planted; e.g., see Figure 9i grid cells C2, and in reclamation areas; e.g., see Figure 9a,c grid cells A5.  5 (b,d,f,h,j,l,n). The shrinking of non-water sub-areas over two decades can be identified by the change of the shape of the non-water class from wide strips to the thin elongated shapes over the series of images (see (a-n); e.g., grid cells C3). Whereas non-water sub-areas emerged when mangroves were planted (see (i) grid cells C2), and in coastal reclamation areas (see (a,c) grid cells A5).  5 (b,d,f,h,j,l,n). The shrinking of non-water sub-areas over two decades can be identified by the change of the shape of the non-water class from wide strips to the thin elongated shapes over the series of images (see (a-n); e.g., grid cells C3). Whereas non-water sub-areas emerged when mangroves were planted (see (i) grid cells C2), and in coastal reclamation areas (see (a,c) grid cells A5). Figure 10a presents the shoreline positions derived for t = 0.5. This t value was selected because it yielded the best κ result when applying the threshold to derive the water class images from membership images. Shoreline was thus assessed through the uncertainty value of water and non-water sub-areas (see Figure 10b). The uncertainty values represent the uncertainty that the pixel belongs to water, determined following Equation (12). A dark blue colour indicates a higher uncertainty that the pixels to belong to the water class, whereas a light blue colour denotes pixels having a lower uncertainty to be classified as water. Generally, pixels which are closer to shoreline have a higher uncertainty of belonging to the water class; e.g., see Figure 10d grid cells C2 and D2.  Figure 10a presents the shoreline positions derived for t = 0.5. This t value was selected because it yielded the best κ result when applying the threshold to derive the water class images from membership images. Shoreline was thus assessed through the uncertainty value of water and nonwater sub-areas (see Figure 10b). The uncertainty values represent the uncertainty that the pixel belongs to water, determined following Equation (12). A dark blue colour indicates a higher uncertainty that the pixels to belong to the water class, whereas a light blue colour denotes pixels having a lower uncertainty to be classified as water. Generally, pixels which are closer to shoreline have a higher uncertainty of belonging to the water class; e.g., see Figure 10d grid cells C2 and D2.  Figure 11 presents the second method, an illustration of the shoreline margin generated by setting t = 0.3 and t = 0.7. Shoreline margin are represented by blue polygons and their width is determined by the shoreline condition. A wider margin shows the more gradual transition from water to nonwater occurring for instance in a low-lying muddy area, or near a swamp area; see Figure 11c grid cells C3 and D3. Meanwhile, a narrow margin reflects a more abrupt transition, as for example at a steep coast or at a shoreline with embankment or shoreline protection; see Figure 11c grid cells B2, and B3. Shoreline margin was assessed through different levels of uncertainty (see Figure 11d-g). The uncertainty values represent the uncertainty that a pixel belongs to the water class estimated following Equation (12). A dark blue colour indicates a higher uncertainty of pixels to belong to the water class, whereas a light blue colour denotes pixels having lower uncertainty. Generally, pixels which are closer to water have a higher membership to the water class. Consequently, these pixels may have a higher certainty to be classified as water.  Figure 11 presents the second method, an illustration of the shoreline margin generated by setting t = 0.3 and t = 0.7. Shoreline margin are represented by blue polygons and their width is determined by the shoreline condition. A wider margin shows the more gradual transition from water to non-water occurring for instance in a low-lying muddy area, or near a swamp area; see Figure 11c grid cells C3 and D3. Meanwhile, a narrow margin reflects a more abrupt transition, as for example at a steep coast or at a shoreline with embankment or shoreline protection; see Figure 11c grid cells B2, and B3. Shoreline margin was assessed through different levels of uncertainty (see Figure 11d-g). The uncertainty values represent the uncertainty that a pixel belongs to the water class estimated following Equation (12). A dark blue colour indicates a higher uncertainty of pixels to belong to the water class, whereas a light blue colour denotes pixels having lower uncertainty. Generally, pixels which are closer to water have a higher membership to the water class. Consequently, these pixels may have a higher certainty to be classified as water.

The Results of Change for Shoreline as a Single Line
The results of the first method in shoreline change detection are presented in Figure 12. We distinguished two types of change, i.e., water to non-water (red colour), and non-water to water (blue colour). Figure 12. (a-f) Shoreline change analysis at t = 0.5. Two changes were identified, namely non-water to water and water to non-water. Large areas changed from non-water to water such as due to inundation

The Results of Change for Shoreline as a Single Line
The results of the first method in shoreline change detection are presented in Figure 12. We distinguished two types of change, i.e., water to non-water (red colour), and non-water to water (blue colour).

The Results of Change for Shoreline as a Single Line
The results of the first method in shoreline change detection are presented in Figure 12. We distinguished two types of change, i.e., water to non-water (red colour), and non-water to water (blue colour). Figure 12. (a-f) Shoreline change analysis at t = 0.5. Two changes were identified, namely non-water to water and water to non-water. Large areas changed from non-water to water such as due to inundation Figure 12. (a-f) Shoreline change analysis at t = 0.5. Two changes were identified, namely non-water to water and water to non-water. Large areas changed from non-water to water such as due to inundation and erosion which were indicated between 1994 and 2000 (a). Whereas large areas changed from water to non-water and were distinguished between 2000 and 2002 (b).
The maps in Figure 13 demonstrate shorelines with their associated change uncertainties derived from Equation (14). Two categories of change uncertainty were identified: (1) change uncertainty to water (shades of blue); and (2) change uncertainty to non-water (shades of red). For both colours, darker shades indicate a higher change uncertainty than lighter shades. Table 4 summarizes the changes in a number of pixels between water and non-water at different levels of uncertainty for site in the yellow rectangle in Figure 13a The maps in Figure 13 demonstrate shorelines with their associated change uncertainties derived from Equation (14). Two categories of change uncertainty were identified: (1) change uncertainty to water (shades of blue); and (2) change uncertainty to non-water (shades of red). For both colours, darker shades indicate a higher change uncertainty than lighter shades. Table 4 summarizes the changes in a number of pixels between water and non-water at different levels of uncertainty for site in the yellow rectangle in Figure 13a. Number of pixels increase with the increase of change uncertainty values.  The trend of shoreline changes was thus assessed in the period 1994-2015 ( Figure 14 and Table  5). The changes in water and non-water sub-areas have been observed by comparing values for two consecutive years. The results for changed areas are reported in Table 5. Table 5 shows that in the period 1994-2000 in which the largest inundation occurred, non-water areas were inundated on one side (−670.1 ha), while a small change to non-water can be found on the other side (e.g., +20.0 ha). The Figure 13. (a) Shoreline change uncertainty at t = 0.5; (b-f) Change uncertainty is highlighted at different levels for the period 1994-2000 for the yellow rectangle site. The number of red pixels indicates that the change uncertainty from water to non-water increase with the increase of uncertainty values, as also can be seen for the blue pixels. Table 4. Changed area (in number of pixels) between water and non-water at different change uncertainty levels (see yellow rectangle site in Figure 13a). The number of pixels increases with the increase of change uncertainty values. Obvious changes were observed by a change uncertainty value ď0.1. Water to non-water +9 +12 +15 +20 +27 Non-water to water´190´219´235´241´250

Change Area
Note: + gain of non-water,´loss of non-water.
The trend of shoreline changes was thus assessed in the period 1994-2015 ( Figure 14 and Table 5). The changes in water and non-water sub-areas have been observed by comparing values for two consecutive years. The results for changed areas are reported in Table 5. Table 5 shows that in the period 1994-2000 in which the largest inundation occurred, non-water areas were inundated on one side (´670.1 ha), while a small change to non-water can be found on the other side (e.g., +20.0 ha). The net change was inundation (´650.2 ha), as shown in Figure 14a, e.g., grid cells C3, and D3. Another large change to water (´186.8 ha) was identified from 2002 to 2003 (see Figure 13c; e.g., grid cells C3, and D3). Whereas extensive change to non-water (+165.5 ha) occurred over the period 2000-2002 (Figure 14b; e.g., grid cells C3 and D3). Table 5. Changed area (in ha) between water and non-water at t = 0.5 and CU ≤ 0.1 during the period 1994-2015. Inundation has been distinguished during four periods (1994-2000, 2002-2003, 2003-2013 and 2014-2015), while change to non-water has been identified for two periods (2000-2002 and 2013-2014 Note: + gain of non-water, − loss of non-water.  Figure 15 provides change maps of the shoreline margin and related sub-areas in the period 1994-2015, in which we identified six changes. A wider change area from non-water to water (blue colour) can be seen in e.g., Figure 15a grid cells C3 and D3. On the other hand, narrow change areas from shoreline to water are present in e.g., Figure 15f grid cells C2. Between 2000 and 2002, large areas changed from water and shoreline to non-water, e.g., Figure 15b grid cells C2, C3 and D3. Some of those areas changed again to water between 2002 and 2003. Those changes could be due to a different growing phase of crops since this location has an extensive agricultural area such as paddy field [32]. Meanwhile, some areas changed from water and shoreline to non-water in the period 2003-2013 (see Figure 15c grid cell B2) which was caused by a successful mangrove planting program in Bedono village.  Table 5. Changed area (in ha) between water and non-water at t = 0.5 and CU ď 0.1 during the period 1994-2015. Inundation has been distinguished during four periods (1994-2000, 2002-2003, 2003-2013 and 2014-2015), while change to non-water has been identified for two periods (2000-2002 and 2013-2014). Water to non-water +20.0 +197.5 +23.2 +51.4 +64.5 +21.7 Non-water to water´670.1´32.0´210.1´182.8´20.3´26. 8 Net change´650.2 +165.5´186.8´131.4 +44.3´5.1

The Results of Change for Shoreline as a Margin
Note: + gain of non-water,´loss of non-water. Figure 15 provides change maps of the shoreline margin and related sub-areas in the period 1994-2015, in which we identified six changes. A wider change area from non-water to water (blue colour) can be seen in e.g., Figure 15a grid cells C3 and D3. On the other hand, narrow change areas from shoreline to water are present in e.g., Figure 15f grid cells C2. Between 2000 and 2002, large areas changed from water and shoreline to non-water, e.g., Figure 15b grid cells C2, C3 and D3. Some of those areas changed again to water between 2002 and 2003. Those changes could be due to a different growing phase of crops since this location has an extensive agricultural area such as paddy field [32]. Meanwhile, some areas changed from water and shoreline to non-water in the period 2003-2013 (see Figure 15c grid cell B2) which was caused by a successful mangrove planting program in Bedono village.

The Results of Change for Shoreline as a Margin
Change uncertainty of shoreline margin, water and non-water are presented in Figure 16. Meanwhile, Table 5 shows the changes in the number of pixels between shoreline margin, water and non-water at different levels of uncertainty for the yellow rectangle site in Figure 16a. The number of pixels in the change area decreases with a decrease in the level of uncertainty. Obvious changes from shoreline margin to water due to inundation were 88 pixels (see Table 6 column 1 and Figure 16a). Obvious changes from shoreline margin to non-water due to reclamation or deposition were indicated for one pixel (see Table 6 column 1 and Figure 16a).   Change uncertainty of shoreline margin, water and non-water are presented in Figure 16. Meanwhile, Table 5 shows the changes in the number of pixels between shoreline margin, water and non-water at different levels of uncertainty for the yellow rectangle site in Figure 16a. The number of pixels in the change area decreases with a decrease in the level of uncertainty. Obvious changes from shoreline margin to water due to inundation were 88 pixels (see Table 6 column 1 and Figure 16a). Obvious changes from shoreline margin to non-water due to reclamation or deposition were indicated for one pixel (see Table 6 column 1 and Figure 16a). Table 6. Changed area (in the number of pixels) between shoreline margin, water and non-water at different uncertainty levels (see yellow rectangle site in Figure 16a). Obvious changes were observed for a level of uncertainty ď0.1. Shoreline to non-water +1 +2 +3 +6 +10 Water to shoreline +17 +14 +12 +21 +23 Non-water to shoreline´11´20´21´26´34

Change Area
Shoreline to water´88´94´101´101´103 Non-water to water´149´175´189´189´190 Water to non-water +6 +8 +10 +11 +11 Note: + gain of non-water,´loss of non-water.  The changes of shoreline margin, water and non-water sub-area were analyzed over the period 1994-2015. Figure 17 and Table 7 present the results of shoreline change analysis at the uncertainty level ď 0.1. Table 7 shows that the largest coastal inundations in terms of area occurred between 1994 and 2000, when 739.4 ha changed from non-water and shoreline to water (see Figure 17a, e.g., grid cells C3, and D3), and in 2003-2013 for 199.7 ha (see Figure 17d, e.g., grid cells B3 and C3). Between 2000 and 2002, 186.3 ha changed from water and shoreline to non-water (see Figure 17b, e.g., grid cells C3 and D3), and another 64.6 ha in 2013-2014 (Figure 17e, e.g., grid cell C1). In general, we can conclude that there was a considerable change of shoreline margin and its associated sub-areas over these two decades. Table 7. Changed area (in ha) at CU ď 0.1 in the period 1994-2015. The largest change from non-water to water due to coastal inundation occurred in the period 1994-2000, and the largest change from water to non-water due to different planting periods in an agricultural area occurred in the period 2000-2002.  Figure 17 and Table 7 present the results of shoreline change analysis at the uncertainty level ≤ 0.1. Table 7 shows that the largest coastal inundations in terms of area occurred between 1994 and 2000, when 739.4 ha changed from non-water and shoreline to water (see Figure 17a, e.g., grid cells C3, and D3), and in 2003-2013 for 199.7 ha (see Figure 17d, e.g., grid cells B3 and C3). Between 2000 and 2002, 186.3 ha changed from water and shoreline to non-water (see Figure 17b, e.g., grid cells C3 and D3), and another 64.6 ha in 2013-2014 (Figure 17e, e.g., grid cell C1). In general, we can conclude that there was a considerable change of shoreline margin and its associated sub-areas over these two decades.

Discussion
In this paper, we have demonstrated two change detection methods for shoreline considering its complex and fuzzy nature including its uncertainty. Both methods were successful in identifying abrupt and gradual changes of shoreline at an object level and in estimating the spatial distribution of uncertainty at the pixel level. The first method derived the shoreline by applying a threshold t = 0.5 on the water membership images. Shoreline was then considered as a single line; its position was influenced by the spatial extent of its associated sub-areas. The uncertainty in the shoreline position could be assessed by means of the uncertainty of its associated sub-areas. The second method derived shoreline margin as an area in which water moves to and fro as tides rise and fall. This margin was presented as a crisp object with a boundary determined by t values resulting from parameter estimation. Comparing these showed that the second method allowed us to assess the spatial extent of the shoreline and measure its change in uncertainty at different levels. This method provided more insight into the spatial distribution of changes and their uncertainty and more spatial detail of the process of change from non-water via shoreline margin to water and vice versa. In estimating the net change category, both methods gave similar results. Net changes identified by the second method covered a larger area both for negative changes (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2002)(2003) Tables 5  and 7). The different results were mainly due to differences in threshold values when generating the shoreline. Shoreline and its changes have been presented as crisp sub-areas. The changed areas were thus associated with the distribution of change uncertainty.
In deriving shorelines as a fuzzy object, we used FCM to calculate the membership. Membership values obtained by applying FCM were used to deal with uncertain information on the position of objects, such as the location of fuzzy shorelines. In this research, we implemented FCM for two classes, since shoreline is the set of locations where water and non-water have physical interactions. In addition, the difference between water and non-water provided the largest spectral differences in images, as was confirmed in [52]. For similar situations, a suitable number of clusters need to be specified either by users based on their a priori knowledge or estimated from the images.
In the literature, several methods exist to measure the cluster validity index for finding a suitable number of clusters, as for example exponential cluster validity index, non-fuzziness index, fuzziness performance index, and entropy measure [53][54][55][56].
Reference data were collected in the field in 2015. In this research, images were selected based on the same low tide condition, and therefore, reference points were also selected under similar low tide conditions. Due to limited availability of high resolution images, however, the 2003 QuickBird image was used for accuracy assessment of images in 2000, 2002 and 2003, even though it was collected during rising tide. Differences in tides between images and these reference data have led to a low accuracy of classified images for 2000, 2002 and 2003, as compared to 2013 and 2015. Moreover, because of limited availability of Landsat images and severe cloud cover problems in the study area, other criteria such as availability of images with the same seasonal condition have not been applied yet. For images captured during rainy seasons, the rain increased the wetness of soil surfaces, which led to a false impression of more flooded conditions in the classified images, and increased the uncertainty in the shoreline position. Landsat images captured in 2000 and 2002 were recorded during the rainy seasons, while the 2003 QuickBird image, used as reference data, was captured during the dry season. Therefore, differences in seasons and tides could have contributed to the misclassification of water and non-water leading to lower classification accuracy. Considering the larger availability of new satellite images in the future, we recommend to select images recorded under similar, preferably dry, seasonal conditions to reduce the uncertainty due to seasonal influence.
The shoreline change detection methods developed and applied in this research involved the generation of a shoreline and its uncertainty using selected m and t values, and produced a hard classification by means of thresholding. A variation in m values had little influence on the results of the classification. An important reason was the small changes in membership values obtained when m values were set between 1.5 and 3.0. Within that range, the change in membership value was generally less than 0.1. However, t variation strongly influenced the results. Small variations in t causing large variation in area indicate the presence of gradual boundaries, as for example in a muddy area. Meanwhile, sharp or sudden boundaries caused little change in area with varying t, such as at a steep coast or at the shoreline with embankment.
The change uncertainty value expresses how sure we are that a change really occurred. Therefore, the uncertainty addressed in this research corresponds to the existential uncertainty of the identified changes. A complete model of uncertainty, however, should also include an extensional uncertainty which considers the spatial extent of the change [57]. This would involve modelling shoreline as a fuzzy object as proposed by Cheng [51]. In fact, the membership and t values of the water class in this work may be used to represent shoreline as a conditional boundary since spatial extents of water (and non-water) depend on the choice of t value. In this case, points where µ " 1 are the interior of fuzzy sub-areas (cores), whereas points with 0 ă µ ă 1 would belong to transition zones and points with µ " 0 belong to non-water objects.
Analysis of shoreline changes in the study area in Sayung sub-district revealed an extensive change of shoreline with the largest change between 1994 and 2000. Land subsidence was indicated as one of the causes of this permanent submergence. Chaussard, et al. [58] mentioned that the rate of the land subsidence in this location reaches approximately 6.0 cm¨year´1. Subsidence was accelerated by an extreme ground water extraction for industrial purposes and as the consequence of a rapid population growth, classified as an anthropogenic subsidence. Putranto and Rüde [59] cited the Directorate of Environmental Geology and Mining Regions of Indonesia, stating the number of registered deep wells in Semarang Demak in the early 1900s was only 16. The number of deep wells increased significantly up to 1194 in 2002. Besides, land subsidence was also triggered by a natural compaction of clayey sediment and settlement loading [38,58,60]. Furthermore, the threat of flooding comes from sea level rise as well. As mentioned by Sofian [61], the rise of sea level in Indonesia is approximately 0.2-1.0 cm¨year´1 with an average of 0.6 cm¨year´1. Different planting periods in agricultural areas could be a major cause of the changes from water and shoreline to non-water in the period 2000-2002 followed by a large reverse change in the period 2002-2003 (see Tables 5 and 7). In addition, we identified coastal land reclamation activities to install an industrial area and a settlement in Sriwulan village. After the period 2003-2013, more agricultural areas were converted into fishponds as a result of expanding saltwater infiltration. On the other hand, some changes from water to non-water occurred as a result of a successful mangrove planting program which started in 2003.
The study area was located in a river delta formed from deposits carried by many rivers discharging into the Java Sea. The remarkable coastal inundation and erosion, which Sayung sub-district has faced for more than two decades, was probably also due to a change in sediment-carrying capacity of the longshore current [62,63]. This longshore current, generated by waves, exceeded the quantity of sediment supplied to the beach. Further, land use change in upstream areas resulted in an increase of erosion and water discharge in particular during rainy seasons, yielding more sediment to the downstream area. In fact, not all of this sediment could be discharged into the sea; some was deposited along the river bottom, irrigation canals, estuaries and other water bodies in its path. This led not only to the narrowing and to silting up of the canals and rivers, but also to the reduction of sediment supply to the littoral zone inducing coastal erosion. Moreover, massive coastal reclamation in a neighbouring area such as Terboyo industrial complex extended seaward after 1994, and Tanjung Emas Harbour first developed in 1985 could have changed currents and material transport along the coast in the study area as well. To prove that there were some influences of harbour development and beach reclamation to the severe inundation impact, however, is beyond the scope of this study.

Conclusions
This research presents two methods to identify shoreline positions: as a line and as a margin, including a measure of change uncertainty at different epochs. Both methods used FCM classification to determine partial membership of water and non-water. While shoreline changes can be detected by both methods, the shoreline as a margin provides a more detailed estimation of change area than the shoreline as a line. Moreover, by having shoreline as a margin, we can assess its spatial extent and measure its change uncertainty at different levels. Abrupt and gradual changes of shoreline were identified at an object level and the spatial distribution of uncertainty estimated at a pixel level.
Some challenges for the improvement of the method consist of: (1) including an extensional uncertainty and represent shoreline as a conditional boundary; (2) including seasonal condition in the selection of images; and (3) integrating the results with a digital elevation model. In addition, the integration with other datasets such as land use and land cover pattern, hydrological data and other information associated with the study area is important as well to have a better understanding of processes which influence the shoreline change.
Both methods have been successfully implemented in a coastal area in the Sayung sub-district in Java using a series of Landsat images. The change area estimation and its change uncertainty may support local government and other stakeholders in monitoring shoreline changes. Integrating the results with the distribution of elements at risk such as settlements and other important facilities can help to analyze which location needs to be prioritized in disaster response. For example, priority could be given to a change area from shoreline to water which has low change uncertainty value as it was obvious that shoreline has changed due to inundation or erosion. Priority could also be given, however, to the area with higher change uncertainty, for example if the location is a densely populated area. Continuous monitoring of shoreline changes in order to understand risks and to anticipate the potential impacts is highly important.