# Mapping Recent Shoreline Changes Spanning the Lateral Collapse of Anak Krakatau Volcano, Indonesia
## 3 Materials and Methods
Following [PERSON] (1996) [37], the NDWI has been calculated in GEE by using the 'image.normalizedDifference' function as follow:  
\[\mathrm{NDWI}=\frac{\mathrm{GREEN-NIR}}{\mathrm{GREEN+NIR}}, \tag{1}\]  
where GREEN is the TOA spectral reflectance acquired in the green waveband (central wavelength of 560 nm in Sentinel-2, band 3) and NIR is the TOA spectral reflectance acquired in the near-infrared region (central wavelength of 842 nm for Sentinel-2, band 8). The resulting NDWI map contains pixel values ranging between \(-\)1.0 and \(+\)1.0 (Figure 4b), where, in our imagery, values \(\leq\)0.2 correspond to land surfaces (vegetated and not) and turbid waters and values \(\geq\)0.2 are associated with the presence of clouds, cloud shadows and clear waters. The subsequent differentiation between land and sea is straightforward because the NDWI values have a bimodal distribution due to the distinct spectral characteristics of the two surfaces types (Figure 4b). In Figure 4b crater and crater walls have higher NDWI than the outer slopes to the north, east and south, which might reflect hydrothermal alteration of the deposits in the crater walls and floor or be a shadowing effect. Figure 3: Workflow diagram showing the step-by-step procedure to extract shorelines. Next, a two-class segmentation is computed to discern land from water by finding the optimum threshold that separates the two classes through an automatic data-driven method, called Otsu's method [38]. The threshold is found by whatever partition of the data that maximizes the between-class (i.e., inter-class) variance by considering the following formula:  
\[BSS=\sum_{k=1}^{p}\left(\overline{NDWI}_{k}-\overline{NDWI}\right)^{2} \tag{2}\]  
where \(BSS\) is the between-sum-of-squares, \(p\) is the number of classes (2 in our case), \(\overline{NDWI}_{k}\) indicates the mean NDWI pixel value in class \(k\) and \(\overline{NDWI}\) is the mean NDWI pixel value of the entire dataset. Class \(k\) is defined by every NDWI less than some threshold and the between-class variance is defined as \(BSS\)/\(p\). [PERSON]'s method looks at every possible partition of the input data based on different \(\overline{NDWI}\) values for class \(k\) (Figure 4c) and requires only a single pass over the data. At each bin of the histogram,  
Figure 4: Example of the application of the proposed methodology. **(a)** Red-Green-Blue (RGB) color composite image of the Top of Atmosphere Sentinel-2A acquisition on 25 April 2019. **(b)** NDWI map and relative histogram. **(c)** BSS distribution for the bins of NDWI values with the threshold (red vertical solid line) that maximizes the between-class variance. **(d)** Binary map within which land and water have been discriminated. the approach defines class \(k\) as all the pixels in that bin and lower while class \(k+1\) represents everything else. The threshold is set where the _BSS_ calculated with (2) is highest. Once found, the threshold is then applied to the NDWI map to generate a binary raster map where areas of sea and land have been discriminated (Figure 4d). The [PERSON]'s method works better when a region in the image that has water and land in approximately equal proportion is strategically chosen so we subset the Sentinel-2 image to an area of interest extended \(\sim\)10 km\({}^{2}\). The edge detection stage consists of converting the boundary raster into polygons and then extracting the land-sea boundary as the coastline. Additionally, for the period June 2018-October 2018, we have the possibility to compare the shoreline changes with the deposition of the volcanic material thanks to the thermal emission recorded by 12 Sentinel-2 images (Table A2) by using the SWIR spectrum (see Figure 2) who is sensible to the temperature of hot emitting surfaces [39].