Change Vector Analysis to Monitor the Changes in Fuzzy Shorelines

: Mapping of shorelines and monitoring of their changes is challenging due to the large variation in shoreline position related to seasonal and tidal patterns. This study focused on a ﬂood-prone area in the north of Java. We show the possibility of using fuzzy-crisp objects to derive shoreline positions as the transition zone between the classes water and non-water . Fuzzy c-means classiﬁcation (FCM) was used to estimate the membership of pixels to these classes. A transition zone between the classes represents the shoreline, and its spatial extent was estimated using fuzzy-crisp objects. In change vector analysis (CVA) applied to water membership of successive shorelines , a change category was deﬁned if the change magnitude between two years, T 1 and T 2 , differed from zero, while zero magnitude corresponded to no-change category. Over several years, overall change magnitude and change directions of the shoreline allowed us to identify the trend of the ﬂuctuating shoreline and the uncertainty distribution. The fuzzy error matrix (FERM) showed overall accuracies between 0.84 and 0.91. Multi-year patterns of water membership changes could indicate coastal processes such as: (a) high change direction and high change magnitude with a consistent positive direction probably corresponding to land subsidence and coastal inundation, while a consistent negative direction probably indicates a success in a shoreline protection scheme; (b) low change direction and high change magnitude indicating an abrupt change which may result from spring tides, extreme waves and winds; (c) high change direction and low change magnitude which could be due to cyclical tides and coastal processes; and (d) low change direction and low change magnitude probably indicating an undisturbed environment, such as changes in water turbidity or changes in soil moisture. The proposed method provided a way to analyze changes of shorelines as fuzzy objects and could be well-suited to apply to coastal areas around the globe.


Introduction
The study of changing shorelines is essential to assist in the design of effective coastal protection [1,2], verifying numerical models [3,4], developing hazard maps [5], formulating policies regarding coastal development [6], and for coastal research and monitoring [7]. A shoreline is defined as the intersection of coastal land and water surface indicating the water edge movements of which the position is changing through time due to different water levels during high tide and low tide [8][9][10]. Oertel [11] referred to a shoreline as the line associated with sea level rather than with high and low tides. When considering only the tide, many shorelines are due to the shifting of water with tidal differences. Tidal differences vary and are influenced by the changes in the magnitude of gravitational attractions on the water body of the Earth, winds and waves. Furthermore, shorelines have changed their dynamics at varying rates as a response to coastal processes such as sediment erosion, transportation and deposition along the shore. Rapid changes occur during an extreme event such as storms, whereas gradual changes occur during an intervening period [12]. Shoreline changes can be estimated over various time scales and result into long-term, cyclic and local random variation. Long-term variation includes variation due to the land subsidence, relative sea level rise and sediment storage. Cyclic variation is related to the tide cycles or seasons, whereas waves and storms cause random variation of a local character.
As the shoreline positions vary over time, shoreline indicators are used as proxies to represent shoreline positions, including: (a) distinguishable coastal features, for example, a previous high-tide line; (b) the intersection of coastal profile with specific vertical water elevation, e.g., the mean sea level; and (c) shoreline features observable from remote sensing images. An example of the latter is the boundary between water and non-water pixels [13][14][15][16]. The shoreline boundary between coastal land and water is fuzzy since there is a gradual transition from coastal land to water. Given the nature of the fuzzy shoreline and its changing position, detection of shoreline requires dealing with uncertainty. Fisher [17] mentioned three types of uncertainty: (a) errors: if a shoreline is clearly identified, the uncertainty may arise from error, for example in data processing, spatial generalization, and measurement; (b) vagueness: if it is not possible to define the spatial extent of coastal land, water, and the transition zone (Williamson (1994) after Fisher [17]); and (c) ambiguity: relating to the confusion of land and water definition considering a different classification system or a different perception.
Previous studies have proposed several ways of generating shoreline positions. Shoreline survey and photogrammetry have been primary technology for shoreline mapping, yet both methods are time consuming and expensive [18]. Therefore, image classification is used widely nowadays to detect shoreline positions. Most studies regarding shoreline detection have used hard classification such as thresholding, water indices, iterative self-organizing data analysis (ISODATA), binary slicing, maximum likelihood classification (MLC) and manual digitizing [7,15,[19][20][21], whereas only a few applied soft classifications [16,22,23].
Due to the fuzziness of shoreline positions, using hard classification for shoreline mapping could produce errors on the classification results, since hard classification assigns a single label to a pixel, based on its highest membership. To overcome this limitation, this paper explores fuzzy classification to detect shoreline positions from a remote sensing image. In our previous work, we proposed two procedures to derive fuzzy shorelines: (a) we derived shorelines by applying a threshold equal to 0.5 to the membership and depicted shorelines as a single line; and (b) we derived shorelines as a margin determined by the choice of thresholds on the membership function [16]. In the current paper, we proposed a third procedure to distinguish shoreline proxies from digital images. A shoreline is represented as the transition zone between water and land. In this case, pixels at which the membership value (µ) exceed 0.99 are the core of a class, whereas pixels with 0.01 < µ < 0.99 belong to transition zones and pixels with µ < 0.01 do not belong to objects. In this way, we can account for the gradual transition between water and land (vagueness of the boundary). Moreover, in change detection, use of transition zones instead of crisp shorelines allows us to account for the influence of ambiguity resulting from comparing images recorded under different circumstances, such as weather, and have a more detailed description, of not only the magnitude and direction of the changes, but also of the related uncertainty.
Various change detection techniques have been developed. They can be divided into two groups, namely bi-temporal change detection and temporal trajectory analysis [24,25]. The former measures changes based on two separate time periods, for example image differencing and post classification comparison. Image differencing does not provide a detailed change matrix while post classification comparison does not allow the detection of subtle changes within a class. The latter, temporal trajectory analysis, is based on the temporal development curve or trajectory for successive times. It focuses both on what has changed between dates, as well as on the trend of the change over the period [24][25][26][27]. Change detection in this research utilizes the second method. To measure the change of the fuzzy shoreline, change vector analysis (CVA) based upon pixel-wise comparison was used to estimate the changes of successive shorelines. CVA identifies changes of features which were acquired at different times. In previous studies, CVA was applied to the brightness and greenness indices [28,29], normalized difference vegetation index [25,30], near infrared band and vegetation index [31], wetness and bare soil index [32], and spectral bands and textural images [33,34]. In this study, CVA was applied to the water membership values of shoreline images. Furthermore, in the earlier studies CVA has been applied in the multi-spectral space [29,35], and then extended to be applied in multi-temporal observation vectors of an indicator variable measured at different times [25]. CVA provides an overall change magnitude and change direction showing the trend of the fluctuating shoreline.
The objective of this study was to develop a method that is useful for monitoring the changes of a fuzzy shoreline. The method is based on fuzzy classification and CVA. A series of Landsat images is used to detect shoreline positions as a transition zone while taking tides into account. For this study, the uncertainty of shoreline positions was estimated by means of confusion indices. We focus on inherent uncertainty caused by continuous variation of a shoreline over time, and on uncertainty as it propagates from extraction and implementation of the shoreline change detection method. The method is applied to an area in Java, in the northern coastal area of the Central Java Province, Indonesia, where extensive shoreline changes associated with coastal inundation have increased in term of frequency and duration.

Study Area
The study area is located at the northern part of the Central Java Province in Indonesia. It is characterized by a low-land landscape ( Figure 1) with an elevation less than 5 m above mean sea level (AMSL). The extent of the study area is approximately 7.5 km from east to west, and 6.5 km from north to south. The central point of the area is at UTM coordinates 444,243 • E and 9,234,731 • N zone 49S or geographic coordinates 6 • 55 S and 110 • 29 E. Alluvial and sand sedimentation dominate its soil type [36,37]. As a coastal area, this area has a mixed semi-diurnal tide with two high tides and two low tides each day. These two highs and lows differ in height, whereas the average tidal range is 0.6 m. The highest tidal ranges occur in December and June during the rainy and dry seasons, while the lowest tidal ranges occur in March and September during the transitional seasons. Two types of flood regularly occur: (a) floods caused by a tidal flood occurring daily in line with tidal cycles [38,39]; and (b) floods due to poor drainage systems during rainy seasons.
Remote Sens. 2017, 9,147 3 of 28 used to estimate the changes of successive shorelines. CVA identifies changes of features which were acquired at different times. In previous studies, CVA was applied to the brightness and greenness indices [28,29], normalized difference vegetation index [25,30], near infrared band and vegetation index [31], wetness and bare soil index [32], and spectral bands and textural images [33,34]. In this study, CVA was applied to the water membership values of shoreline images. Furthermore, in the earlier studies CVA has been applied in the multi-spectral space [29,35], and then extended to be applied in multi-temporal observation vectors of an indicator variable measured at different times [25]. CVA provides an overall change magnitude and change direction showing the trend of the fluctuating shoreline. The objective of this study was to develop a method that is useful for monitoring the changes of a fuzzy shoreline. The method is based on fuzzy classification and CVA. A series of Landsat images is used to detect shoreline positions as a transition zone while taking tides into account. For this study, the uncertainty of shoreline positions was estimated by means of confusion indices. We focus on inherent uncertainty caused by continuous variation of a shoreline over time, and on uncertainty as it propagates from extraction and implementation of the shoreline change detection method. The method is applied to an area in Java, in the northern coastal area of the Central Java Province, Indonesia, where extensive shoreline changes associated with coastal inundation have increased in term of frequency and duration.

Study Area
The study area is located at the northern part of the Central Java Province in Indonesia. It is characterized by a low-land landscape ( Figure 1) with an elevation less than 5 m above mean sea level (AMSL). The extent of the study area is approximately 7.5 km from east to west, and 6.5 km from north to south. The central point of the area is at UTM coordinates 444,243°E and 9,234,731°N zone 49S or geographic coordinates 6°55′S and 110°29′E. Alluvial and sand sedimentation dominate its soil type [36,37]. As a coastal area, this area has a mixed semi-diurnal tide with two high tides and two low tides each day. These two highs and lows differ in height, whereas the average tidal range is 0.6 m. The highest tidal ranges occur in December and June during the rainy and dry seasons, while the lowest tidal ranges occur in March and September during the transitional seasons. Two types of flood regularly occur: (a) floods caused by a tidal flood occurring daily in line with tidal cycles [38,39]; and (b) floods due to poor drainage systems during rainy seasons.  Extensive fishponds and rice fields are covering the study area. Settlements are found bordering the sea and along the riverbanks, which are threatened if tidal floods become higher. This area has suffered from a changing shoreline position leading to severe coastal inundation and erosion. The coastal inundation has increased recently both in terms of frequency and duration. Some factors such as extreme winds and waves contribute to this increase. Furthermore, in the long run, other causes such as land subsidence, sea level rise, mangrove conversion, beach reclamation and a seaport extension are potential causes for increased coastal inundation [38,[40][41][42].
2.2. Satellite Images, Data Pre-Processing and Reference Data Generation

Satellite Images and Data Pre-Processing
Multi-temporal images from the Landsat 8 OLI/TIRS (Operational Land Imager/Thermal Infrared Sensor) with 30 m spatial resolution were used to monitor the shoreline change between 2013 and 2015 ( Table 1). We obtained terrain corrected Landsat images (L1T product) from USGS EarthExplorer [43]. Those images were acquired at the low tide. Tidal data relating to the time of acquisition of the images were collected from the Indonesian Geospatial Information Agency. Pre-processing of Landsat 8 OLI/TIRS comprises two steps: (a) histogram minimum adjustment; it was applied to remove the influence of atmospheric path radiance [44,45]; and (b) geo-referencing; it was implemented using >100 ground control points (GCP) collected from road intersections, rivers and other prominent features. The root mean square error (RMSE) values were less than 0.1 pixels. Geo-registration of Landsat images was conducted using geometrically corrected reference images: (1) a Pleiades image at a 2 m spatial resolution; (2) a SPOT 6 (Satellite Pour l'Observation de la Terre) image at a 6 m spatial resolution; and (3) a Sentinel 2 image at a 10 m spatial resolution. The spectral band information for each reference image including Landsat 8 OLI/TIRS is available in Table 2.

Reference Data Generation
To evaluate the accuracy of a fuzzy classification, it is necessary to use soft reference data [46,47]. We generated soft reference data from available fine resolution datasets [48,49]. These datasets (Pleiades, Spot 6 and Sentinel 2) were rectified using a 2015 orthoimage. To reduce the variance of the Pleiades image, smoothing was performed using the average filter applied to a 3 × 3 window size. Afterwards, we applied fuzzy c-means (FCM) with the number of classes c = 2 and the fuzzy weight m = 1.7 [16]. Further, membership images generated using FCM classification from these high resolution datasets were used as reference images.
For accuracy assessment purpose, the pixel size of Spot 6 image was resampled to 10 m using nearest neighbour resampling, so that the spatial resolution of Pleiades, Spot 6, Sentinel 2 and Landsat images were in the ratio 15:3:3:1. Hence, 225 pixels (15 × 15) of Pleiades, 9 pixels (3× 3) of Spot 6, and 9 pixels (3 × 3) of Sentinel 2 were combined (pixel values averaged) to achieve the pixel dimension of Landsat images. Furthermore, an effective comparison could be made between images of different resolutions.
For the alternative methods, MLC and hardened classification, we visually interpreted Pleiades and Spot 6 images as hard reference data for the year 2013 and 2014 respectively, whereas ground data were used as the 2015 reference data.

FCM Classification
To discriminate water classes from non-water, we applied a fuzzy c-means (FCM) classification [50]. FCM iteratively separates data clusters with fuzzy means and fuzzy boundaries and the results assign each pixel to a partial membership of land cover classes. The membership values (µ) range from 0 to 1, and add up to 1 for each pixel. In this work, the membership values of the classification follow the trapezoidal membership function. In the literature, there are two possible ways of generating the membership function: Similarity Relation Model (SRM) and Semantic Import Model (SIM) [51][52][53]. The former derives the membership function using classifiers like for example fuzzy k-means, fuzzy c-means, and neural networks [53,54]. The first two are data driven, partitioning the observations based on multivariate attributes. The latter, SIM, generates membership based on expert knowledge [51,55].
The FCM results assign each pixel to membership of the two classes. The value of m determines the level of fuzziness in FCM classification. If m = 1, FCM is hard classifier. FCM was carried out by labeling two membership images resulting from each FCM classification as the water and non-water images. To do so, the combination of near infrared (NIR) and shortwave infrared (SWIR) of Landsat bands were used. The water label was given to the class which has the minimum value of the sum of the cluster means in the infrared bands. Detailed descriptions regarding the FCM algorithm are available in Bezdek et al. [50], whereas detailed explanations regarding membership function, pixel labeling, and parameter estimation for FCM classification can be found in Dewi et al. [16].

Validation
To quantify the accuracy of the FCM classifier, a conventional error matrix cannot be used. In this study, we used a fuzzy error matrix which has non-negative real numbers [48,56,57], since pixels have a partial membership to two classes.
For accuracy assessment, soft reference images were generated by applying an FCM classification to Pleiades, Spot 6, and Sentinel 2 images which were all captured during low tides. Let the value of µ ik and µ jl represent membership values of the kth pixel for class i in the classified image and lth pixel for class j in the reference images. It was assumed that the rows of the matrix are classes of the classified image and the columns are classes of the reference image. The fuzzy error matrix (FERM) is obtained using minimum operator showing the maximum possible overlap between the classified and reference images and indicating the agreement between classes in both images [48,49,56]: To calculate the agreement in FERM, a group of 225 Pleiades pixels (15 × 15), 9 Spot 6 pixels (3 × 3) and 9 Sentinel 2 pixels (3 × 3) were averaged to achieve pixel dimension of Landsat images. Using this reference data, membership of 200 pixels randomly selected from both classified and reference images were computed to obtain the overall accuracies (OA) of the FCM classifications: where k s represents the number of pixels used to generate the FERM. Moreover, we compared the quality of the FCM results with respect to alternative pixel-based classification methods. Firstly, we classified the multi-spectral bands of Landsat using the MLC classifier being the most commonly used supervised classification technique for remote sensing images [58]. Secondly, we classified the multi-spectral bands of Landsat images using FCM and then labelled each pixel to the class to which it has the highest membership. It was assumed that hard output is the highest membership value which is actually computed from the soft output [59][60][61][62]. We called this the hardened classification. After classification, post classification comparisons were applied to detect the changes of the shorelines by superimposing the classification results in GIS.

Deriving Fuzzy Shoreline
FCM classification derives two raster layers, namely: water and non-water membership images. Each layer consists of fuzzy regions with fuzzy boundaries. Estimation of the spatial extent of objects i.e., water, non-water and shoreline, and their representations is related to the interpretation of the fuzziness of objects [52]. To derive shorelines at the locations where water and non-water objects meet, we modified the fuzzy-crisp object model based upon Cheng [52]. The two classes (water and non-water) are spatially disjoint, but their boundary is vaguely defined, whereas their interiors are crisp. Given this concept, we consider the boundary between water and non-water as fuzzy and form a transition zone that we call shoreline. To determine the spatial extent of water, non-water and shoreline, it is necessary to combine class objects from different layers into a single layer. The decision function d wk assigns pixel k with water membership value µ wk to a sub-area of water class based upon the following conditions: which means that the pixels belong to sub-areas water. Threshold 0.99 was set to represent the highest water membership values indicating the core of water.
This equation classifies pixels as shoreline.
Pixels not belonging to water or shoreline constitute non-water. A threshold of 0.01 represents the lowest water membership values. This indicates pixels with membership below that threshold not belong to water or shoreline areas. The results after deriving fuzzy shorelines by applying Equations (3)-(5) were called as shoreline images.

Uncertainty Estimation
The uncertainty in class assignment was estimated by a measure of the confusion index CI for each pixel resulting from FCM classification as follows [63][64][65]: If CI approaches 1 then the difference in membership values between the first and the second highest membership values are small meaning that both membership values are almost equal. Thus, it is more likely that the pixel defines a fuzzy boundary and the uncertainty of the pixel to belong to the class with the largest membership is high. If CI approaches 0, however, then the difference in membership values between the first and the second highest membership values are high and the uncertainty of the pixel to belong to the class with the largest membership is low.

Shoreline Change Detection
For establishing the changes over time, shoreline images obtained using Equations (3)-(5) of the same year were stacked and compared with the stack of shoreline images of the next year with corresponding seasons. If membership values to water (µ wk ) of shoreline images within year T 1 and T 2 are given by G = (g 1 , g 2 , .. . . . , g z ) T 1 and H = (h 1 , h 2 , .. . . . , h z ) T 2 , respectively, and z is the number of shoreline images, a change vector is defined as: Here, ∆CV includes all the change information between two years for a given pixel. The final result of CVA is an image of vector changes. The shoreline change is defined as the vector difference between successive time periods and is represented by a vector in a multi-dimensional space. The length of the change vector indicates the magnitude of change and its direction indicates the nature of the change [25,66].

Change Magnitude
The change magnitude ∆CV was derived by determining the Euclidean distance between shoreline images as: ∆CV represents the total membership differences between two years and measures the intensity of the shoreline change. Two categories of change were identified, namely change and no-change. A change category was defined when the water membership difference between T 1 and T 2 is larger than zero, whereas a no-change category is related to a magnitude equal to zero. A higher change magnitude corresponds with a large water membership difference between shoreline images in T 1 and T 2 . When the change magnitude is low, the water membership difference between shoreline images in T 1 and T 2 is small.

Change Direction
For all pixels classified as change, we estimated the change directions. Change direction was determined by evaluating the water membership difference between shoreline images in two successive years. It quantifies the variation of water membership in each pixel and shows how frequent the changes have occurred. Change direction estimation started by calculating the number of change combinations (CC) as: where d refers to the types of change direction which can be distinguished when comparing the stack of shoreline images from both years for corresponding seasons and p refers to the number of shoreline image pairs. We identified three types of change direction to water: positive change direction (or in short positive direction), negative change direction (negative direction) and unclear change direction (unclear direction).
The change vector (CV) showing water membership difference between a pair of shoreline images in T 1 and T 2 from corresponding seasons needs to be estimated: (a) if the water membership difference between pair of shoreline images within years T 1 and T 2 is less than zero then CV = −1 showing a decrease of water membership in T 2 ; (b) if the water membership difference is larger than zero then CV = +1 showing an increase of water membership in T 2 ; (c) if the water membership difference is equal to zero then CV = 0 showing that the water membership in T 1 and T 2 were the same. The total change vector (TCV) values are defined as: and CVz refers to (h z − g z ). Finally, the change direction (Chg.dir) categories showing the degree of change direction to water membership were obtained by grouping the direction values: (a) TCV values from +1 up to +z were grouped as positive direction; (b) TCV values from −1 up to −z were grouped as negative direction; (c) TCV values equal to 0 showing unclear change directions were classified as unclear direction; and (d) TCV values equal to 0 having water membership differences equal to 0 at all time periods were classified as no-change. Table 3 shows the procedure to determine the change direction categories by using four pairs of image used in this study.
Based upon these results, the change area of a specific change direction category (positive direction, negative direction, and unclear direction) and the no-change area were defined as: where Pix k (Chg) is the number of pixels belonging to the area of change and no-change, and Ar(k) is area of pixel k (30 × 30 m 2 ).

Change Uncertainty
Based upon the change detection results, the change uncertainty of related areas was estimated by the confusion index CI. If CI of two images for T 1 and T 2 are given by Q = (q 1 , q 2 , . . . .., q z ) T 1 and R = (r 1 , r 2 , . . . ., r z ) T 2 , respectively, then the change confusion is derived as: A high ∆CU value is related to a large difference of confusion indices between images for T 1 and T 2 , whereas a low change confusion corresponds to a small difference of confusion indices between images for T 1 and T 2 .
Notes: CV: Change vector (based on Equation (7); CC: Change combinations number (based on Equation (9)); TCV: Total change vector (based on Equation (10)); and Chg.Dir: Change direction. Table 4 presents the accuracy assessment of classification results using FCM and alternative classification methods. The FCM classifier outperformed MLC and the accuracy values of FCM are generally higher than the hardened classification.  Figure 2 presents an example of FCM outputs, together with MLC and hardened classification. In the image, MLC overestimated the non-water area shown by the larger area of non-water (see Figure 2d-f, e.g., grid cells B1 and B2), whereas the hardened classification underestimated the non-water area (see Figure 2g-i, e.g., grid cells B1 and B2). Both methods failed to distinguish the gradual transition between water and non-water.

FCM Classification and Accuracy Assessment
Remote Sens. 2017, 9, x FOR PEER REVIEW 12 of 27 Table 4 presents the accuracy assessment of classification results using FCM and alternative classification methods. The FCM classifier outperformed MLC and the accuracy values of FCM are generally higher than the hardened classification.  Figure 2 presents an example of FCM outputs, together with MLC and hardened classification. In the image, MLC overestimated the non-water area shown by the larger area of non-water (see Figure 2d-f, e.g., grid cells B1 and B2), whereas the hardened classification underestimated the non-water area (see Figure 2g-i, e.g., grid cells B1 and B2). Both methods failed to distinguish the gradual transition between water and non-water.  The results of FCM classification are presented in Figure 3 with the values ranging from 0 to 1 for both membership images of water (Figure 3a-c) and non-water (Figure 3d-f). Areas with higher water membership values were located for example in marine areas, fishponds, and water-covered agricultural areas (Figure 3a, e.g., grid cell A2). In Figure 3d, higher non-water membership pixels are located near settlements adjacent to the shorelines, and mangrove forests (Figure 3d, e.g., grid cells B2 and C3).

FCM Classification and Accuracy Assessment
Remote Sens. 2017, 9,147 13 of 28 The results of FCM classification are presented in Figure 3 with the values ranging from 0 to 1 for both membership images of water (Figure 3a-c) and non-water (Figure 3d-f). Areas with higher water membership values were located for example in marine areas, fishponds, and water-covered agricultural areas (Figure 3a, e.g., grid cell A2). In Figure 3d, higher non-water membership pixels are located near settlements adjacent to the shorelines, and mangrove forests (Figure 3d, e.g., grid cells B2 and C3). To derive shoreline position, we combined both membership images using fuzzy-crisp object model (g-i). Blue pixels indicate core of water, orange pixels represent the core of non-water and shoreline is represented by light green pixels.   The width of these shorelines is determined by natural conditions of the coastal areas, for example, a wider shoreline is more likely to be found in a muddy coastal area or at a gently sloping beach, whereas a narrow shoreline is usually found along a steeper slope beach and coastal area with embankment and other man-made structures. To derive shoreline position, we combined both membership images using fuzzy-crisp object model (g-i). Blue pixels indicate core of water, orange pixels represent the core of non-water and shoreline is represented by light green pixels.   The width of these shorelines is determined by natural conditions of the coastal areas, for example, a wider shoreline is more likely to be found in a muddy coastal area or at a gently sloping beach, whereas a narrow shoreline is usually found along a steeper slope beach and coastal area with embankment and other man-made structures.

Change Magnitude and Change Uncertainty
The change magnitude and the change and no-change categories of shoreline are displayed in Figure 5. Low change magnitude values correspond to a small water membership difference between shoreline images at and . They cover marine areas (Figure 5a, e.g., grid cells A2 and B3) and a relatively undisturbed coastal land (Figure 5c e.g., grid cells D1 and D2). In addition, high change magnitude values correspond to a large water membership difference. Those pixels cover muddy areas (Figure 5a,c, e.g., grid cells C1 and D1) and coastal land which was highly-influenced by tidal floods (Figure 5b,d, e.g., grid cells B2 and B3).

Change Magnitude and Change Uncertainty
The change magnitude and the change and no-change categories of shoreline are displayed in Figure 5. Low change magnitude values correspond to a small water membership difference between shoreline images at T 1 and T 2 . They cover marine areas (Figure 5a, e.g., grid cells A2 and B3) and a relatively undisturbed coastal land (Figure 5c e.g., grid cells D1 and D2). In addition, high change magnitude values correspond to a large water membership difference. Those pixels cover muddy areas (Figure 5a,c, e.g., grid cells C1 and D1) and coastal land which was highly-influenced by tidal floods (Figure 5b,d, e.g., grid cells B2 and B3).

Change Magnitude and Change Uncertainty
The change magnitude and the change and no-change categories of shoreline are displayed in Figure 5. Low change magnitude values correspond to a small water membership difference between shoreline images at and . They cover marine areas (Figure 5a, e.g., grid cells A2 and B3) and a relatively undisturbed coastal land (Figure 5c e.g., grid cells D1 and D2). In addition, high change magnitude values correspond to a large water membership difference. Those pixels cover muddy areas (Figure 5a,c, e.g., grid cells C1 and D1) and coastal land which was highly-influenced by tidal floods (Figure 5b,d, e.g., grid cells B2 and B3).  The fuzziness of shoreline changes is presented in Figure 6. Low change confusion correspond to small CI differences between images in T 1 and T 2 . This indicates a low uncertainty that the changes have occurred as can be seen in Figure 6a,c, e.g., grid cells A1, A2 and B1. High values are associated with large CI differences and indicate a high uncertainty that the changes have occurred (Figure 6b,d, e.g., grid cells B1, B2 and B3).
The fuzziness of shoreline changes is presented in Figure 6. Low change confusion correspond to small differences between images in and . This indicates a low uncertainty that the changes have occurred as can be seen in Figure 6a,c, e.g., grid cells A1, A2 and B1. High values are associated with large differences and indicate a high uncertainty that the changes have occurred (Figure 6b,d, e.g., grid cells B1, B2 and B3).

Change Direction
The representation of change direction of shoreline showing the variation of water membership in each pixel can be seen in Figure 7. Positive directions to water membership correspond to the increase of water membership at time (Figure 7a, e.g., grid cells A2, B1, and C1). On the contrary, negative directions to water membership were associated with the decrease of water membership values at time (see Figure 7d, e.g., grid cells A2 and B1). The no-change category indicates an undisturbed environment (see Figure 7a, e.g., grid cell B1), whereas the unclear direction category indicates an ambiguous condition since the changes occurred without an obvious trend (see Figure 7b, e.g., grid cells A2 and C1).
The change areas for each category are presented in Table 5

Change Direction
The representation of change direction of shoreline showing the variation of water membership in each pixel can be seen in Figure 7. Positive directions to water membership correspond to the increase of water membership at time T 2 (Figure 7a, e.g., grid cells A2, B1, and C1). On the contrary, negative directions to water membership were associated with the decrease of water membership values at time T 2 (see Figure 7d, e.g., grid cells A2 and B1). The no-change category indicates an undisturbed environment (see Figure 7a, e.g., grid cell B1), whereas the unclear direction category indicates an ambiguous condition since the changes occurred without an obvious trend (see Figure 7b, e.g., grid cells A2 and C1).
The change areas for each category are presented in Table 5 Figures 9 and 10 show the change detection of shoreline using post classification comparison of MLC and hardened classification, respectively. Both MLC and hardened classification present shoreline as a single line. The changes of this single shoreline have occurred due to the changes of water and non-water area. Binary images from two dates T 1 and T 2 were superimposed in GIS and four types of change were identified, namely: non-water to water, water to non-water, water to water and non-water to non-water. Figures 9a-c and 10a- (Figures 9d and 10d) and from 2014 to 2015 (Figures 9e and 10e). Figures 9 and 10 show the change detection of shoreline using post classification comparison of MLC and hardened classification, respectively. Both MLC and hardened classification present shoreline as a single line. The changes of this single shoreline have occurred due to the changes of water and non-water area. Binary images from two dates and were superimposed in GIS and four types of change were identified, namely: non-water to water, water to non-water, water to water and non-water to non-water. Figures 9a-c and 10a- (Figures 9d and 10d) and from 2014 to 2015 (Figures 9e and 10e).    Figures 9 and 10 show the change detection of shoreline using post classification comparison of MLC and hardened classification, respectively. Both MLC and hardened classification present shoreline as a single line. The changes of this single shoreline have occurred due to the changes of water and non-water area. Binary images from two dates and were superimposed in GIS and four types of change were identified, namely: non-water to water, water to non-water, water to water and non-water to non-water. Figures 9a-c and 10a- (Figures 9d and 10d) and from 2014 to 2015 (Figures 9e and 10e).    Figure 11 shows the comparison between the proposed method and the alternative method at the selected study area. In this example, both methods agree on the results of change detection as can be seen in Figure 11, e.g., grid cells B3 and C2. From CVA results (Figure 11c,d), the area in yellow polygons shows a negative change to water membership with high change magnitudes as shown in Figure 11e,f. The negative change to water membership means a change towards non-water. Similarly, post classification results also denote that these yellow polygon sites experienced a change from water to non-water without further information on the intensity of the change (Figure 11a,b). Figure 11 shows the comparison between the proposed method and the alternative method at the selected study area. In this example, both methods agree on the results of change detection as can be seen in Figure 11, e.g., grid cells B3 and C2. From CVA results (Figure 11c,d), the area in yellow polygons shows a negative change to water membership with high change magnitudes as shown in Figure 11e,f. The negative change to water membership means a change towards non-water. Similarly, post classification results also denote that these yellow polygon sites experienced a change from water to non-water without further information on the intensity of the change (Figure 11a,b).

Multi-Year Pattern of Water Membership Changes
Each pixel from the resulting change vectors provides information regarding its change direction and magnitude. Each combination represents specific types of change processes that may occur in the field and shows a multi-year pattern of water membership changes over the observation periods. Four combinations of change and their related processes are interpreted as follows:

(a) High change direction and high change magnitude
The areas with high change direction and high change magnitude values are observed for both positive and negative directions. Both conditions indicate a continuous change of an area to a certain direction with a relatively large intensity. A consistency to positive direction indicates a persistence of enhanced water influence as those pixels show an increase of water membership in multi-temporal images (see Figure 12a-d). This probably corresponds to the land subsidence and coastal inundation. As the land subsides and the water level increases, some mangrove trees located closely to the sea are falling down. The RGB 542 of Landsat images in Figure 12e,f depict these changes indicated by the decrease of vegetation cover between 2013 and 2015. Figure 13 presents areas characterized by continuously decreasing water membership in multitemporal images categorized as negative direction (Figure 13a,b) with high change magnitude (Figure 13c,d). This may indicate a success in shoreline protection scheme that caused sediment accretion to occur allowing mangroves to grow as can be seen from the RGB 542 of Landsat images in white-dashed ellipses in Figure 13e-h.

Multi-Year Pattern of Water Membership Changes
Each pixel from the resulting change vectors provides information regarding its change direction and magnitude. Each combination represents specific types of change processes that may occur in the field and shows a multi-year pattern of water membership changes over the observation periods. Four combinations of change and their related processes are interpreted as follows: (a) High change direction and high change magnitude The areas with high change direction and high change magnitude values are observed for both positive and negative directions. Both conditions indicate a continuous change of an area to a certain direction with a relatively large intensity. A consistency to positive direction indicates a persistence of enhanced water influence as those pixels show an increase of water membership in multi-temporal images (see Figure 12a-d). This probably corresponds to the land subsidence and coastal inundation. As the land subsides and the water level increases, some mangrove trees located closely to the sea are falling down. The RGB 542 of Landsat images in Figure 12e,f depict these changes indicated by the decrease of vegetation cover between 2013 and 2015. Figure 13 presents areas characterized by continuously decreasing water membership in multi-temporal images categorized as negative direction (Figure 13a,b) with high change magnitude (Figure 13c,d). This may indicate a success in shoreline protection scheme that caused sediment accretion to occur allowing mangroves to grow as can be seen from the RGB 542 of Landsat images in white-dashed ellipses in Figure 13e-h.

(b) Low change direction and high change magnitude
This category indicates an abrupt change which may be influenced by random events. Figure 14b shows a low positive direction with high change magnitude (Figure 14d) which may result from coastal flooding triggered by spring tides, extreme waves and winds. Since the magnitude of the changes is high and the change is sudden, this type of change may indicate a higher risk. Images made available by Google Earth from 2013 to 2015 in Figure 14e-h show the decrease of mangrove coverage. In fact, mangroves can act as sediment trap and can reduce the energy of the high waves, therefore, when the mangroves disappear, the threat from tidal floods increases. Figure 15a,b shows the embankment which protects settlements from high tide; however, during an extreme event for example when a higher tide combines with an extreme wind, the water level may increase and overflow this embankment.

(b) Low change direction and high change magnitude
This category indicates an abrupt change which may be influenced by random events. Figure  14b shows a low positive direction with high change magnitude (Figure 14d) which may result from coastal flooding triggered by spring tides, extreme waves and winds. Since the magnitude of the changes is high and the change is sudden, this type of change may indicate a higher risk. Images made available by Google Earth from 2013 to 2015 in Figure 14e-h show the decrease of mangrove coverage. In fact, mangroves can act as sediment trap and can reduce the energy of the high waves, therefore, when the mangroves disappear, the threat from tidal floods increases. Figure 15a,b shows the embankment which protects settlements from high tide; however, during an extreme event for example when a higher tide combines with an extreme wind, the water level may increase and overflow this embankment.   This category indicates an abrupt change which may be influenced by random events. Figure  14b shows a low positive direction with high change magnitude (Figure 14d) which may result from coastal flooding triggered by spring tides, extreme waves and winds. Since the magnitude of the changes is high and the change is sudden, this type of change may indicate a higher risk. Images made available by Google Earth from 2013 to 2015 in Figure 14e-h show the decrease of mangrove coverage. In fact, mangroves can act as sediment trap and can reduce the energy of the high waves, therefore, when the mangroves disappear, the threat from tidal floods increases. Figure 15a,b shows the embankment which protects settlements from high tide; however, during an extreme event for example when a higher tide combines with an extreme wind, the water level may increase and overflow this embankment.   (c) High change direction and low change magnitude A gradual, continuous increase of wet conditions was observed by an increase in water membership with low magnitude values. This type of change was categorized as positive direction which could be due to cyclical tides and coastal processes, for example flooded land (Figure 16a,c), and water turbidity (Figure 16b,d). Even though the magnitude of the change is low, the changes occur frequently. Hence, this type of change may give a higher risk. In a longer observation, if the areas persistently become wetter, this location may have a risk of coastal inundation as well.

(c) High change direction and low change magnitude
A gradual, continuous increase of wet conditions was observed by an increase in water membership with low magnitude values. This type of change was categorized as positive direction which could be due to cyclical tides and coastal processes, for example flooded land (Figure 16a,c), and water turbidity (Figure 16b,d). Even though the magnitude of the change is low, the changes occur frequently. Hence, this type of change may give a higher risk. In a longer observation, if the areas persistently become wetter, this location may have a risk of coastal inundation as well.

(d) Low change direction and low change magnitude
This type of change probably indicates an undisturbed environment with a low change magnitude (see black-dashed circle sites in Figure 17c,d). This category mainly occurs in water areas, probably due to the changes in water turbidity (see black-dashed circle sites in Figure 17a,b). In addition, this type of change was observed in small patches of the coastal land probably resulting from changes in soil moisture (see black rectangle sites in Figure 17b,d).

(d) Low change direction and low change magnitude
This type of change probably indicates an undisturbed environment with a low change magnitude (see black-dashed circle sites in Figure 17c,d). This category mainly occurs in water areas, probably due to the changes in water turbidity (see black-dashed circle sites in Figure 17a,b). In addition, this type of change was observed in small patches of the coastal land probably resulting from changes in soil moisture (see black rectangle sites in Figure 17b,d).

Discussion
In this study, the dynamics of fuzzy shorelines have been assessed using fuzzy classification and a raster-based change detection technique. FCM classification was used to discriminate the land and water classes and to estimate their memberships. FCM is a well-known clustering method which is less dependent on the initial state of clustering [62] and capable of describing phenomena such as water and non-water which is changing gradually. Instead of FCM, there are various ways to derive a fuzzy classification for example from fuzzy maximum likelihood classification [59], and artificial neural network fuzzy classification [67]. Membership values obtained by applying FCM are used to deal with the uncertain information on the position of fuzzy shorelines. In FCM classification, we set c = 2 since we were interested in identifying the boundary between water and non-water and because both classes give the largest spectral differences in image [16,68]. Finding the suitable number of clusters in the beginning of the classification could be difficult. A priori knowledge regarding the study area, for example by observing an aerial photo, can be used to define the suitable number of clusters [69,70]. For other situations, by assessing the homogeneity measure using a posteriori indicators, the number of clusters could be determined using entropy and non-fuzziness index [70], exponential cluster validity [71], and spatial fuzzy clustering [72].
Shorelines and their changes were presented as fuzzy areas. The fuzzy-crisp object model in this study was successful in identifying the extent of shoreline positions as the transition zone between water and non-water. Setting the threshold to the highest (0.99) and the lowest (0.01) memberships are intended to find the core of water and non-water, respectively. The uncertainty addressed in this

Discussion
In this study, the dynamics of fuzzy shorelines have been assessed using fuzzy classification and a raster-based change detection technique. FCM classification was used to discriminate the land and water classes and to estimate their memberships. FCM is a well-known clustering method which is less dependent on the initial state of clustering [62] and capable of describing phenomena such as water and non-water which is changing gradually. Instead of FCM, there are various ways to derive a fuzzy classification for example from fuzzy maximum likelihood classification [59], and artificial neural network fuzzy classification [67]. Membership values obtained by applying FCM are used to deal with the uncertain information on the position of fuzzy shorelines. In FCM classification, we set c = 2 since we were interested in identifying the boundary between water and non-water and because both classes give the largest spectral differences in image [16,68]. Finding the suitable number of clusters in the beginning of the classification could be difficult. A priori knowledge regarding the study area, for example by observing an aerial photo, can be used to define the suitable number of clusters [69,70]. For other situations, by assessing the homogeneity measure using a posteriori indicators, the number of clusters could be determined using entropy and non-fuzziness index [70], exponential cluster validity [71], and spatial fuzzy clustering [72].
Shorelines and their changes were presented as fuzzy areas. The fuzzy-crisp object model in this study was successful in identifying the extent of shoreline positions as the transition zone between water and non-water. Setting the threshold to the highest (0.99) and the lowest (0.01) memberships are intended to find the core of water and non-water, respectively. The uncertainty addressed in this research corresponds to the existential and extensional uncertainty of shoreline objects as has been mentioned by Molenaar and Cheng [73] and Cheng [52]. Existential uncertainty expresses the uncertainty of the existence of shoreline in reality. It refers to the possibility of existence of a shoreline to be detected on an image. Extensional uncertainty implies that the area indicated as shoreline can be determined with limited certainty, for example with boundaries that reflect the transition zone between water and non-water. Moreover, when the values of an adjacent grid are very similar, the zones of confusion divide regions indicating the presence of gradual transitions. The extensional uncertainty in shoreline identification includes differences in applied threshold values when defining the core of water and non-water, the applied shoreline definition, tides condition during image acquisition, time series of observation, and the nature of the beaches (such as flat or steep slope beaches, and muddy or rocky beaches). In addition, the changed areas of the fuzzy shoreline are thus associated with the distribution of changed confusion indices. The change uncertainty represented by changed confusion indices shows the degree of uncertainty of the changes that have occurred. It can be seen from the results that a location having a higher change magnitude, has a higher change confusion value as well. It corresponds to the higher differences of both water memberships and confusion indices between corresponding images in T 1 and T 2 . Explicit handling of uncertainty by addressing the shoreline as a transition zone allows decision makers and planners to include this uncertainty in spatial planning. Moreover, it visualizes not only the changes in shoreline, but also the uncertainty of these changes for every location, thereby providing a better base for a debate on the combined effects of land subsidence and sea level rise in this area.
The change of shoreline was explained in terms of change magnitude and change direction using CVA. Information provided by CVA allows us to see the trend of the fluctuating shoreline over time, whereas the change detection results of the alternative method could provide only "from-to" change information and the detail of subtle within-class changes was lost because it only compared images from two dates [24,25]. In our previous study, we used post classification comparison because we were interested in observing the changes of shoreline over a longer period from 1994 up to 2015 [16]. Given different methods have been implemented in monitoring the change position of shoreline, both studies confirmed that shoreline changes associated with coastal inundation have occurred in this study area.
The analysis of information provided by the change magnitude and direction reveals that each change combination represents one specific change process type. The processes could vary depending on the characteristic of the coastal areas. For example, shorelines could change due to floods triggered by land subsidence, and floods caused by seasonal variation, abrupt shoreline changes due to extreme tides and waves, and the changes of water turbidity and soil moisture triggered by daily weather events. These specific type processes were explained on the basis of the analysis of four images for each observation period. In fact, the number of images could easily be extended to more than four images. The seasonal variation of shorelines and other information regarding whether the changes would lead to a permanent coastal inundation could not fully be assessed in this study because only four observations each year were compared. This shows that the change vector analysis is sensitive to the length of the stacked period and the number of images stacked over that period, as is also confirmed by Lambin and Strahler [25].
To have images captured during similar tides and with corresponding seasons for each pair is important, hence shorelines can be compared equally. Images captured during wet seasons gave more uncertainty to the classification results since the rain affected the wetness of soil surfaces, thus the classified images produced false impression of higher flooding. Therefore, to reduce the uncertainty due to seasonal influence, images captured during dry season should be preferred. Furthermore, the uncertainty of the results was also influenced by tide differences. Images captured during high tides and low tides produced different positions of the boundary between water and non-water. Hence, this increased the uncertainty of shoreline position as well, because shoreline by its definition is an intersection of coastal land and water surface. To have images acquired at exactly the same tide level is hardly possible. Therefore, all images were acquired at the low tides with negligible differences. We only considered astronomical tide level assuming it was more influential than meteorological factors, as confirmed by Pugh [74]. Astronomical and meteorological factors have different influences on different slope conditions and the magnitude of the influences may become larger if the slope is gentler. However, if there would be any remaining small influence of meteorological factors, it would be accounted for by the use of fuzzy classification in deriving fuzzy shoreline.
In the accuracy assessment, soft reference data were generated from various higher resolution images. To obtain the required image resolution, resampling and aggregating have been implemented. Resampling of image and aggregation of pixel values were potential sources of error but were ignored in this work since the error was likely to be very small [49,75]. Furthermore, using a soft classifier to generate soft reference data is more likely to reduce the uncertainty due to the vagueness in class definition and mixed pixel problems. Hence, the finer resolution dataset was not assumed to be pure and no information was lost due to the hardening of the soft classification [60,61,76]. This could be an explanation for the higher accuracies obtained by using the proposed classification in comparison with alternative methods, as confirmed by Zhang and Foody [61] and Chawla [75]. Furthermore, although the differences of the accuracy results were only small, the advantage of a fuzzy approach is not only in the improved accuracy of the shoreline, but also in the fact that it makes clear what the margins of uncertainty around the shoreline are, which provides a better basis for decision.
Analysis of the shoreline changes in the northern coastal area of Central Java shows the changes of shoreline positions from 2013 to 2015. This could be related to the processes that shape shorelines determined by the interaction of several factors, including: (a) the change of sea-level; (b) the amount of sediment supplied to the beach by rivers; (c) the movement of the sediment by marine processes; and (d) the role of waves, currents, tides and winds in moving the sediment [77]. Furthermore, sediment transport is not constant, and it is constantly subject to change. The alteration of sediment transport can come from changes in water flow, water level, weather events and human influence. In addition, previous studies have mentioned that this location has suffered from the changing of shorelines for more than 20 years due to coastal inundation accelerated by, for example, land subsidence, sea level rise, seaport development, and ground water extraction [16,20,38]. Many attempts have been made to combat coastal inundation and erosion along the 1.3 km coast in Demak. Elevated roads, raised floor of the houses, breakwater, and mangrove planting have been applied. In 2013, a Dutch-Indonesian consortium agreed to start a pilot project "Building with Nature" building a permeable dam of natural material called "hybrid engineering" [78,79]. This development could be one reason for the increase of negative changes to water membership in the period of 2014-2015. Hybrid engineering is one type of coastal protection combining technical and ecosystem-based solutions referred to as sediment traps [79]. Netherlands-Water-Partnership [78] reported that after a year, new sediment layers were deposited at the surrounding areas.

Conclusions
In this article, we present a method to identify shoreline positions and their changes as a fuzzy area including a measure of change confusion. Shoreline changes could be detected, and the method provided information regarding the change magnitude and the trend of water membership in each pixel. Our results reveal that this information represents specific type of change processes showing multi-year patterns of water membership changes over time. These include: (a) high change direction and high change magnitude with a consistent positive change direction, which probably corresponds to land subsidence and coastal inundation, while a consistent negative direction may indicate a success in shoreline protection scheme; (b) low change direction and high change magnitude indicates an abrupt change which may be influenced by random events, such as flooding triggered by spring tides, extreme waves and winds; (c) high change direction and low change magnitude, which, if in positive direction, could be due to cyclical tides and coastal processes, for example flooded land; and (d) low change direction and low change magnitude probably indicates an undisturbed environment, such as water areas with changes in water turbidity and coastal land with changes in soil moisture. Finally, we conclude that the proposed method can assess changes in a shoreline by taking into account that it is a fuzzy boundary.
The change area estimation, change magnitude and direction of the shorelines may support local government and stakeholders in monitoring the change of fuzzy shorelines. Combining information given by this research with other information such as distribution of population could help to determine priority locations prioritized in disaster preparedness and response. To include digital elevation model in the processing phase is important for further research, because it allows the change analysis to focus more on the area affected by tidal floods. We further realized that astronomical and meteorological factors have different influences on different slope conditions. Observation data regarding the wave run-up and other incident wave conditions in the study area were not available for each observation period. This information may be collected and used in a near-shore wave model. This might be included in future studies as well.