- Email: [email protected]
Coastline extraction from repeat high resolution satellite imagery
Remote Sensing of Environment 229 (2019) 260–270
Contents lists available at ScienceDirect
Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse
Coastline extraction from repeat high resolution satellite imagery a,⁎
a
b
c
Chunli Dai , Ian M. Howat , Eric Larour , Erik Husby a b c
T
Byrd Polar and Climate Research Center, The Ohio State University, Columbus, OH, USA Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA Polar Geospatial Center, University of Minnesota, St. Paul, MN, USA
ARTICLE INFO
ABSTRACT
Keywords: Satellite images Coastline NDWI Water classification
This paper presents a new coastline extraction method that improves water classification accuracy by benefitting from an ever-increasing volume of repeated measurements from commercial satellite missions. The widely-used Normalized Difference Water Index (NDWI) method is tested on a sample of around 12,600 satellite images for statistical analysis. The core of the new water classification method is the use of a water probability algorithm based on the stacking of repeat measurements, which can mitigate the effects of translational offsets of images and the classification errors caused by clouds and cloud shadows. By integrating QuickBird, WorldView-2 and WorldView-3 multispectral images, the final data product provides a 2 m resolution coastline, as well as a 2 m water probability map and a repeat-count measurement map. Improvements on the existing coastline (GSHHSthe Global Self-consistent, Hierarchical, High-resolution Shoreline Database, 50 m–5000 m) in terms of resolution (2 m) is substantial, thanks to the combination of multiple data sources.
1. Introduction Coastal areas are ecologically and economically important and exceptionally dynamic. Accurate shoreline mapping is essential for navigation, coastal monitoring and modeling, as well as coastal resource management and planning (Liu and Jezek, 2004; Liu et al., 2007). Changes in coastlines occur due to variations in tides and weather, estuarine and ocean circulation, geomorphological processes, as well as anthropogenic activities (urbanization and tourism) (Baiocchi et al., 2012; Palazzo et al., 2012). In particular, the coastal zone is sensitive to global sea-level rise (Vafeidis et al., 2008). Given sea-level rise is expected to accelerate in the coming decades, accurate detection and monitoring of coastlines will become essential to provide needed information for mitigation and adaptation. The most commonly used coastline datasets are the Global Selfconsistent, Hierarchical, High-resolution Shoreline Database (GSHHS) provided by Wessel and Smith (1996) and the 30-m resolution Shuttle Radar Topography Mission (SRTM) water body dataset (SWBD) (NASA/ NGA, 2003; Slater et al., 2006; Farr et al., 2007) ranging from 56° south to 60° north latitude. The GSHHS data was constructed from two military data sets (Wessel and Smith, 1996): 1) the Central Intelligence Agency's World Data Bank II, which was created by U.S. government in the 1980s with a resolution of around 500–5000 m and 2) the World Vector Shoreline (WVS), which was generated with better spatial
⁎
resolution (50–500 m) between 1985 and 1989 by the Defense Mapping Agency, now the National Geospatial-Intelligence Agency (NGA) (Soluri and Woodson, 1990; Wessel and Smith, 1996). Historically, coastlines were manually identified and traced from high-resolution aerial images by experts in cartographic applications (e.g. Liu and Jezek, 2004; Dellepiane et al., 2004). Since manual delineation is labor intensive and often subjective, several automatic coastline extraction methods have been proposed, including the edge tracing algorithm for Synthetic Aperture Radar (SAR) imagery (Lee and Jurkevich, 1990), an image segmentation method for radar and optical satellite imagery (Liu and Jezek, 2004), a coherence thresholding method for interferometric SAR (InSAR) (Dellepiane et al., 2004), a combined method integrating image segmentation, region growing, and edge detection for multispectral imagery (Zhang et al., 2013), an active contour method for Polarimetric SAR images (Liu et al., 2017), the normalized difference vegetation index (NDVI) method (Hsiao et al., 2016) and the Normalized Difference Water Index (NDWI) index method for multispectral satellite images (McFeeters, 1996; Maglione et al., 2014). In addition to the above shoreline detection from image processing, methods based on digital elevation models have also been explored (e.g. Liu et al., 2007). The ever-increasing volume of remote sensing data has a significant impact on studies of Earth surface processes and surface water changes. High resolution optical (visible to near infrared band) satellite images
Corresponding author. E-mail address: [email protected] (C. Dai).
https://doi.org/10.1016/j.rse.2019.04.010 Received 2 November 2018; Received in revised form 5 April 2019; Accepted 9 April 2019 Available online 26 April 2019 0034-4257/ © 2019 Elsevier Inc. All rights reserved.
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 1. Data coverage and test area. The colored map shows the data coverage of ArcticDEM, which includes all land area north of 60°, as well as all territory of Greenland, the State of Alaska in entirety, and the Kamchatka Peninsula of the Russian Federation. Here the surface elevation is plotted using the low resolution (1 km) ArcticDEM data. Our test areas include the white rectangle area of 200 km by 150 km in the Aleutian Islands and the entirety Iceland (white quadrilateral).
have been used successfully for analyzing and extracting surface water boundaries (Feyisa et al., 2014). For multispectral images, the NDWI is the most widely used method initially developed for Landsat images (McFeeters, 1996). Modified NDWI indices were further proposed by Xu (2006) and Feyisa et al. (2014) using shortwave infrared band in Landsat 5 Thematic Mapper (TM). Wolf (2012) modified the index specifically for WorldView-2 and 3 multispectral imageries through the use of coastal and near infrared -2 bands. Based on the widely-used NDWI and modified NDWI formulas, we propose a new water classification and shoreline mapping method that utilizes the large and growing volume of repeat, high resolution multispectral satellite imagery now available. We apply this new algorithm to the much higher resolution data provided by commercial multispectral imagers in order to achieve the most precisely mapped coastline currently available from space.
97 min. The revisit time is about 1 to 3 days with an imaging swath width of 13 to 17 km. WorldView-2 and 3 are capable of collecting around 1 million km2 of imagery per day, which is about seven per mill of the total land surface area of Earth per day. The spatial resolution of the panchromatic band is finer than that of the multispectral band, e.g. WorldView-2 image's resolution is 46 cm in panchromatic band and 1.85 m using 8 multispectral bands. Image band wavelengths are slightly different for each satellite (https://www.digitalglobe.com). We test our algorithm acquired over a portion of Alaska's Aleutian Islands, as well as all imagery of Iceland from the archive of Arctic imagery held by the Polar Geospatial Center (Fig. 1), which is the largest single holding of these imagery available to researchers. The algorithm, however, can be applied over any region where sufficient data is available. The radiometrically corrected level 1 imagery is orthorectified into a polar stereographic projection using the sensor model Rational Polynomial Coefficients (RPCs). We orthorectify the imagery using the 2.5 arc-minute Earth Gravitational Model 2008 (EGM2008) (Pavlis et al., 2012) for terrain correction. The geoid is used rather than a DEM because DEMs are either too low in resolution to provide an accurate surface at our target resolution, or tend to be noisy around the coastline, such as with ArcticDEM (Porter et al., 2018), which increases the error in our results. The imagery is then converted from digital number
2. Data sources QuickBird (2001–2015), GeoEye-1 (since 2008), WorldView-2 (since 2009) and WorldView-3 (since 2014) are commercial satellites operated by DigitalGlobe, Inc. (Anderson and Marchisio, 2012), the operator and owner of a constellation of high-resolution imaging satellites. Their orbits are sun-synchronous, with orbital period around 261
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
0.25
a)
QuickBird Ocean Land
0.2 0.15 0.1 0.05 0 -1.5
-1
-0.5
0
0.5
1
1.5
2
1
1.5
2
1
1.5
2
NDWI 0.25
b)
GeoEye-1 Ocean Land
0.2 0.15 0.1 0.05 Fig. 2. Example of water classification using the NDWI. The colormap is the NDWI calculated from a WorldView-2 image acquired on 25 September 2014 at Nakalilok Bay, Alaska (within the white rectangle in Fig. 1). The transparent background is the shaded relief image of the 5-m ArcticDEM mosaic. For the shaded image, the illumination is from the sun's azimuth direction at the actual acquisition time. The red line is the GSHHS coastline. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
0 -1.5
-1
-0.5
0
0.5
NDWI 0.25
c)
0.2
WorldView Ocean Land
0.15
to top-of-atmosphere reflectance using the absolute radiometric calibration parameters provided by DigitalGlobe (Kuester, 2017; Jawak and Luis, 2013).
0.1 0.05
3. Methods
0 -1.5
3.1. Water classification and adaptive thresholding
NIR2)/(Coastal + NIR2)
NIR)/(Green + NIR)
0
0.5
Fig. 3. Histograms of the probability distribution (the count of observations in each bin normalized by the total count) of mean NDWI from (a) 1319 QuickBird images, (b) 272 GeoEye-1 images, and (c) 11,045 WorldView-2 and 3 images. Blue is for NDWI over ocean only, and red for NDWI over land area only. The mean NDWI over ocean is 1.0, 0.4, 0.98 for all QuickBird, GeoEye-1, WorldView images, respectively; and the mean NDWI over land is −0.1, −0.2 and 0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
(1)
NDWI values as over water. This may result in misclassification (Fig. 4) due to the similarity of shadows and water in reflectance patterns (Feyisa et al., 2014). Clouds also tend to have low NDWI values, similar to land, which may result in false islands. In order to determine a water classification threshold, we evaluate the statistics of NDWI over known areas of water and land using a large number of test images. The water or land area is specified a priori based on the GSHHS coastline, which has a spatial resolution of better than 5 km (Wessel and Smith, 1996). By eroding or dilating the land area encompassed by the a priori coastline (red line in Fig. 2) with a width comparable to its uncertainty (e.g. 2 km), the land or ocean mask can be retrieved with high confidence by removing the pixels in the transitional buffer zone (Liu and Jezek, 2004). The mean and standard deviation of NDWI values within the known land and water areas (outside a 2-km buffer of the GSHHS coastline) are then calculated. The statistical analysis of NDWI is carried out using 10,816 images in Iceland and 1820 images in the 200 km by 150 km test region (Fig. 1)
where, Coastal is the top-of-atmosphere reflectance of the coastal band (band 1) and NIR2 is that of the near infrared-2 band (band 8). For the 4-band GeoEye-1, QuickBird, Ikonos, for which the Coastal and NIR2 bands are not available, the NDWI index is calculated as (McFeeters, 1996):
NDWI = (Green
-0.5
NDWI
Based on the comparison between different algorithms such as NDVI, NDWI, Spectral Angle Mapper (SAM), and Automated Water Extraction Index (AWEI) provided by Baiocchi et al. (2012) and (Fisher et al., 2016), we follow Maglione et al. (2014) and adopt the NDWI index for water classification. The NDWI index is calculated from the 8band WorldView-2 and 3 multispectral, orthorectified images (Wolf, 2012; Maglione et al., 2014) as:
NDWI = (Coastal
-1
(2)
where, Green is the top-of-atmosphere reflectance of the green band, and NIR is that of near infrared band. Water bodies are expected to yield a positive NDWI, while land surfaces should be negative (McFeeters, 1996). As shown in Fig. 2, the ocean surface has an NDWI of above 1, while land surfaces have a value averaging −0.4. The inland water bodies (lakes and rivers) also produce positive NDWI values. Water classification methods have difficulty in distinguishing water from land surfaces with low albedos, such as over asphalt, mountain shadows, and the shadows of buildings and clouds (Feyisa et al., 2014). Fig. 2 gives an example of mountain shadows that result in similar 262
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 4. Example of water masking and strip mosaicking. The water mask is generated from the NDWI map (Fig. 2 for scene 1) with a threshold of 0.3. This strip is composed of two WorldView-2 scenes acquired on 25 September 2014, where the scene alignment offset (d2 in Eq. (5)) of scenes 1 (red) and 2 (magenta) is (0, 0) m and (−0.1567 0.0685) m, respectively, in the polar stereographic coordinate system. After the correction of d2, the two scenes are merged together seamlessly. The transparent background is based on the shaded relief image of the 5-m ArcticDEM mosaic tile. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
of the Aleutian Islands, Alaska, from QuickBird, GeoEye-1, WorldView2 and 3. First, the means and standard deviations of NDWI over the known land and ocean areas are calculated for each image. Then, the probability distribution histograms of these mean values are generated. Fig. 3 shows that, generally, the land has a lower NDWI and the ocean larger NDWI. Furthermore, Fig. 3 demonstrates that NDWI can unambiguously delineate land and water for QuickBird and WorldView-2 and 3 images, where land has NDWI values around 0 and ocean has NDWI values averaging around 1.0. Overlap between values is likely due to inland lakes and/or clouds over ocean. For GeoEye-1 images, however, the separation of land and water is not clear, with substantial ocean areas having an NDWI of less than −1. In addition, the average of all NDWI standard deviations from GeoEye-1 images is much larger (3.6) than that (0.2) from QuickBird and WorldView. Thus, we chose not to not apply GeoEye-1 imagery to water classification. Based on the NDWI statistics of our test images, we use a simple adaptive thresholding method for water classification. For each individual image, the threshold, T, is calculated as:
T = (meanW + meanL)/2
2014) (Fig. S1). Due to the fact that NDWI values are very uniform over the ocean, its probability density function is almost a delta function with near zero variance, which brings the optimal threshold to be very close to the mean NDWI of ocean. Our test shows that our simple method is more practical in the application of water classification using NDWI values. Empirically, we establish the following criteria (Eqs. (4a) and (4b)) to exclude bad quality images, such as ones that have clouds or haze.
< minimum {0.5, (meanW
meanW
meanL > 0.6
meanL)/2}
(4a) (4b)
where σ is the NDWI standard deviation for each image over land or water. Any image that does not meet these criteria will be excluded. Fig. 4 shows an example of the produced high-resolution water mask, where water is denoted in blue and land in green. This example also provides a good illustration of expected misclassifications in mountain shadows.
(3)
3.2. Coregistration
where, meanW and meanL are the mean NDWI values within the known water and land areas, respectively. We compared this simple thresholding method to a widely-discussed optimal adaptive thresholding method (Gonzalez and Woods, 1992; Liu and Jezek, 2004; Wang et al.,
The images we use for coastline mapping using the NDWI index will have systematic errors in geolocation (i.e. biases in the map coordinates) of up to several meters due to errors in the sensor model used 263
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
we find that d3 follows a normal distribution (Fig. 5 (a)), with a nearzero mean of (1.5, 1.0) m and standard deviation of (3.7, 3.7) m, which is close to the sensor model geolocation uncertainty of 4 m as discussed in (Noh and Howat, 2015). The cumulative density function of the offsets, shown in Fig. 5 demonstrates how the accumulation of images can reduce the effect of individual translational offsets. As shown in Fig. 5 (b), when a group of images cover the same coastline, only a very small percentile of the images has large negative translational offsets; and the median (50th percentile) of the translational offsets is zero. This characteristic is used later for selection of the water probability threshold. 3.3. Water probability through stacking of images As shown above, the NDWI threshold may misclassify shadowed surfaces and/or clouds. Averaging of multiple, overlapping pixels using repeat imagery (i.e. stacking), however, can mitigate these errors. For any pixel, we can therefore define water probability as the ratio of the number of measurements that classify a pixel as water over the total number of measurements, i.e.:
Pw = n w / n × 100%
where Pw is the water probability, nw is the number of measurements that classify the pixel as water, and n is the total number of all measurements. Fig. 6 shows the stacking of three individual water masks (Fig. 6 (a)–(c)), and the calculated water probability (Fig. 6 (d)) based on Eq. (6). The final water classification (Fig. 7 (a)) is determined from the water probability. Given a threshold of water probability, e.g. 50%, a pixel is classified as water if Pw is equal to or larger than the threshold, and land if less than the threshold. Fig. 6 also shows how water probability can mitigate misclassification due to clouds and shadows. For the simple NDWI thresholding of an individual image, the cloud shadows introduce false lakes inland and clouds introduce false islands in ocean. When stacking all three water masks together, the cloud shadows only slightly increase the water probability, and the clouds slightly reduce the water probability (Fig. 6 (d)). As long as there are fewer cloudy than clear images, and with the proper thresholding of water probability, the merged water mask (Fig. 7 (a)) may not be contaminated by the clouds and cloud shadows. To save computation resources, the GSHHS coastline is used as a priori information to provide an approximate constraint with a ~10 km buffer, so that only pixels within the buffer are processed. Also, pixels with < 3 repeat measurements are voided. The code for retrieving coastlines from multispectral imageries is available on Github (https:// github.com/Chunli-Dai/SatelliteImageCoastline). Note that for the calculation of NDWI statistics of each individual image (Section 3.1), we only process data outside a 2-km buffer of the GSHHS coastline considering that the typical size of an image is only 13–17 km. There are two main advantages to the proposed water probability method. The first is the mitigation of both translational geolocation errors (bias) in the imagery as well as temporal variability due to coastal dynamics and local sea level. Considering each of these factors separately, as shown in Fig. 5, translational errors follow a normal distribution. Hence, the 50th percentile (median) of the translational errors is near the mean of zero. Assuming a coastline is static, using the 50% water probability threshold is effective at mitigating the translational errors, especially if the total number of measurements is large. On the other hand, if the variation of coastline is only caused by local sea level change and coastal dynamics, such as erosion, the 50th percentile reflects the median of coastline locations. Combining both the effects from translational errors and coastline dynamic variations, the 50th percentile will still likely yield the median of coastline locations (Fig. S2). One benefit of using the median coastline is that occasional coastline changes caused by extreme events, such as flooding and storm surge, generally will not affect the predicted coastline.
Fig. 5. Statistics of sensor model bias (d3 in Eq. (5)) for individual images obtained from DEM coregistration. (a) The histogram probability density distribution of d3, where px, py are the x, y components. (b) The histogram cumulative density function. The median (50%) is shown in black, which corresponds to the horizontal shift at 0 m. This evaluation is based on 249 WorldView-2 and WorldView-3 strip files.
to convert between image and map coordinates (Noh and Howat, 2014, 2015). These errors may be removed by registration to identifiable ground control points, such as surveyed landmarks, or through coregistration to other, precisely geolocated, reference imagery through automated cross correlation methods (e.g. Jeong et al., 2017; Leprince et al., 2007). The former is labor intensive, as control points must be visually located in the imagery, and the latter is often not available. If, on the other hand, the sensor model errors tend to be normally distributed, the stacking of multiple images should partly mitigate such errors in the resulting coastline map. Here we evaluate these errors, and the effect of image stacking, using the registration available from ArcticDEM. For images that are one of a stereoscopic pair used for digital elevation model (DEM) extraction, we coregister that DEM to the registered reference DEM (the 5-m ArcticDEM mosaic which is registered to ICESat satellite altimetry) to determine the three-dimensional offsets. The coregistration is achieved by aligning static surfaces between two DEMs (Noh and Howat, 2015; Dai et al., 2018; Barr et al., 2018). For alignment, we use the iterative least-squares method of Nuth and Kääb (2011). To convert the map coordinate of a point in the orthorectified image, X0, to the same point on the ArcticDEM tile, X1, we have:
X1 = X 0 + d1 + d2 + d3
(6)
(5)
where d1 is the sensor model bias correction shift of each image of the stereo pair with respect to the corresponding DEM (Noh and Howat, 2018), which is zero for the first (left) image of the stereo pair. The term d2 is the scene alignment shift applied for merging scenes into strips (Fig. 4). This term, obtained from the ArcticDEM strip metadata file, is small, with a mean value of zero and standard deviation of 0.2 m based on 40, 000 test scenes. d3 is the strip alignment offset for the 5-m ArcticDEM tile mosaicking (offset applied to align the strip DEM with the mosaic DEM tile) which is obtained through coregistration of the strip DEM to the tile DEM as described above. Based on 249 strip files over a test region (white rectangle in Fig. 1), 264
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 6. Example of water probability mapping. (a)–(c) Water masks generated based on NDWI thresholding of multispectral WorldView-2 and 3 images. (d) Calculated water probability from the three individual water masks. The transparent background is the shaded relief image of the 5-m ArcticDEM mosaic tile.
The second advantage of the water probability approach is the reduction in misclassification of water caused by clouds and shadows. Stacking reduces the impact of these errors because they are temporary anomalies, and thus are reduced or removed as outliers. Nevertheless, there may be consistent commission errors, such as mountain shadows, where images are acquired at the same time of the day and, therefore, with a similar sun angle. These commission errors will not be filtered out by the water probability method.
3.4. Cloud detection Cloud detection is an essential step in classification (Zhu et al., 2015), because clouds and cloud shadows cause errors (Champion, 2012). Specifically, for water classification from individual images, clouds over ocean may result in misclassification of water as islands; or cloud shadows across the shorelines will cause misclassification of coastlines along the shadow edges. Using multi-temporal satellite images (Champion, 2012), we analyze the absolute difference, dM, between each water mask, Mi, and the merged water mask (as a reference), Mm, to detect cloudy images. 265
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 7. Cloud Detection. (a) The merged water mask based on water probability (Fig. 6 (d)). (b) The absolute difference between the merged water mask and the individual water mask as in Fig. 6 (a). (c) The difference after removing smaller clusters. (d) Classify each cluster of difference to identify whether it is surrounded by water (blue edges) or land (green edges). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
dM = |Mi
Mm|
(7)
and discarded. A problem with this approach arises when the number of cloudy images is equivalent to the number of clear images, so that cloudy and clear classifications cannot be determined unambiguously. This approach to cloud detection, however, has the benefit of favoring no classification over false classification. As discussed in the previous section, occasional cloud contamination would not affect the final coastline mapping. Nevertheless, the elimination of cloudy images can improve the water probability map, which can be used later for applications in detecting low-laying areas or estimating coastal dynamic/ tidal range.
As shown in Fig. 7 (b), the absolute difference between the merged water mask, dM, (Fig. 7 (a)) and the individual water mask observed on 26 March 2017 (Fig. 6 (a)) highlights clouds and cloud shadows. Considering differences between images caused by seasonal variations in river width and tides, any clusters of pixels, where dM equals 1, in the difference map that have long and narrow (100 m width) spatial pattern may be removed (Fig. 7 (c)), along with clusters smaller than 1000 m2. Each cluster is classified as clouds on water or extra water on land (due to either flood inundation or false identification by shadows). Since our focus is on coastline mapping, extra water inland (clusters with green edges in Fig. 7 (d)) are discarded. If the total area of clouds above water is > 10,000 m2, the whole image will be treated as cloudy 266
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 8. Modeled ocean tidal heights derived from the 1/6° resolution TPXO-9.1 global inverse tide model. (a) Modeled peak tidal height (defined as the maximum absolute surface tidal heights over a year) in meters on 1/6° grids around Nakalilok Bay, Alaska. Black circles represent the centers of all available images in the vicinity (~20 km) of a tide model grid (white triangle). The cyan square is the boundary of a standard ArcticDEM tile size of 50 km. The colored points show the estimated tidal heights along coastlines within this tile. The tide model is downloaded from http://volkov.oce.orst.edu/tides/global.html (last accessed February 2019). (b) Modeled tidal height time series as hours of the day (Greenwich Mean Time (GMT)) at a grid point (white triangle in (a)). Blue dots are the hourly time series from November 2004 to September 2017, which is the time span of all available images in this area. The blue lines show the mean tidal heights at each hour averaged over all days. The magenta and green lines represent the tidal heights variation for the date of 14 April 2013 and 11 March 2013, respectively. Black circles are the modeled tidal heights at the time and location of all acquired satellite images as in (a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
fjord. Further, the impact of bias in sea level due to tides or other variability on the coastline estimate will depend the slope of the coast. Correcting for this bias requires projecting the coastline onto a DEM and adjusting the coastline horizontally depending on the estimated bias and surface slope. We caution that the relatively coarse resolution and/or noise in existing, global coastal DEMs, however, yields errors that are similar or greater than the tidal bias.
3.5. Tidal variation Any mapped coastline will be biased by the stage of tides at the time (s) of observation. The tide causes the sea level to oscillate around mean sea level at a set of discrete periods, predominantly diurnal (daily) or semidiurnal (twice daily) periods (Stammer et al., 2014; Padman et al., 2018). The maximum range of tide heights is around several meters and can exceed 10 m in a few coastal areas (Padman et al., 2018). Using the 1/6° resolution TPXO-9.1 global inverse tide model (an update to Egbert and Erofeeva (2002)), the tidal height at the location and acquisition time of each individual image can be calculated. As shown in Fig. 8 (a), the peak tidal heights (defined as the maximum absolute surface tidal heights over a year) in the area near Nakalilok Bay are around 2.6 m. Due to the sun-synchronous orbit configuration, DigitalGlobe satellites pass a given location on Earth at about the same time (within 2 h) of a day (black circles in Fig. 8 (b)). The modeled tidal height time series (Fig. 8 (b)) shows a large variation in tidal heights at a given time of a day over the days with observations. For the example shown in Fig. 8 (b), the image acquisition times are close to 22:00, but the tidal heights vary between −1.8 m to 1.6 m. Nevertheless, the mean tidal height on each hour (blue curve with triangles in Fig. 8 (b)) averaged over 13 years has smaller amplitude (around 0.4 m here) than the overall amplitude (~2.6 m). As discussed in Section 3.3, the 50th percentile of water probability gives the median coastline, which is corresponding to the median tidal height if ignoring the random image coregistration errors and other driving factors of dynamic coastline. Here, we adopt the median tidal height of all acquisitions at a location as an estimate of tidal bias in the estimated coastline. Fig. 8 (a) (colored points along coastline) shows the estimated tidal height of the derived coastline, which is within 1 m and is less than the peak tidal heights in the area. It is worth mentioning that when sample size is large (e.g. > 100) the median value reflects the mean value of all observed tidal variations. Since satellites pass at around the same time of the day, the mean tidal height may still be slightly biased (Fig. 8 (b)). We note that the true local variability in tidal height may not be represented at the resolution (1/6°) of the global tide model, such as tidal variations in a deep estuary or narrow
4. Results and discussion By mosaicking multiple satellite images with different areas of coverage, the final product is a 2 m coastline polyline, of similar format to other coastline datasets, such as GSHHS, and that can be readily incorporated into a Geographical Information System (GIS) database. In addition, we produce a 2-m resolution water probability map, as well as a map of the total count of repeat images used in the mapping. As shown in Fig. 9, the final results are produced on a tile-by-tile basis, with a size of 50 km by 50 km. The computation area is expanded by 2 km to avoid edge discontinuation problems. The total count of repeated measurements (Fig. 9 (b)) are calculated after removing cloudy images, where the larger the number of measurements, the more reliable the water probability map and the final coastlines. After retrieving the water probability mask, the coastlines (blue line in Fig. 9 (c)) of islands and inland lakes smaller than 10,000 m2 are discarded as noise. We plot the GSHHS and SRTM coastlines for comparison (Fig. 9 (c)). Within the study area, our coastline seems to be closer to the SRTM coastline with differences within 100 m (e.g. 2–15 m at the middle part of the saddle shape area); while it differs from the GSHHS coastline by up to 200 m. However, our coastline more clearly agrees with the coastline visible in the orthoimage acquired by WorldView-2 on 25 September 2014. The water probability map can be used to identify coastal low-lying areas, which are very sensitive to sea level variations. For example, the pink color in Fig. 10 (a) shows the lower water probability area at the coast of the Nakalilok Bay, Alaska, which implies that this area has a large variability in water coverage. The possible causes of sea level variation include ocean, pole and solid Earth tides, atmospheric 267
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 9. Final results of water probability and coastlines. (a) Water probability over a tile of 50 km by 50 km generated from the stacking of individual water masks. The probability of 0 is set to be transparent for better viewing. The green line is the GSHHS coastline. The cyan square is the boundary of a 50 km by 50 km tile (same as Fig. 8), and the black rectangle marks the region for enlargement as in c. (b) The total number of measurements. Background of (a) and (b) is the hillshade of SRTM 1 arcsec DEM. (c) The comparison of different coastlines. The yellow line is the SRTM water body dataset, and blue line is the coastline from this study. The white rectangle marks the region of Fig. 10 (d). The enlargement of (b) at the black rectangle area. Background of (c) and (d) is the 2 m panchromatic orthoimage acquired by WorldView-2 on 25 September 2014. Imagery© 2014 DigitalGlobe, Inc. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
pressure loading, steric effects introduced by temperature and salinity variations within oceans, and water mass redistribution cycles (Chen et al., 2000; Chambers, 2006). Considering the fact that the near-infrared radiation would be strongly absorbed by water and reflected by land area (McFeeters, 1996), we use the near-infrared band of orthoimages to visually demonstrate the large water coverage variation of this area. As shown in Fig. 10 (b), an over 200 m wide of land area is exposed on 14 April 2013, 21:44 UTC, which is probably caused by the low tide at the observation time (modeled tidal height of −1.1 m (Fig. 8 (b))). And Fig. 10 (c) shows the coastal inundation observed on 11 March 2013, 22:35 UTC, probably due to the high tide (model tidal height of 1.6 m) at the observation time. Thus, our water probability mapping is able to reveal areas of large variability along the coast that
are likely to be most vulnerable to changes in sea level. Furthermore, in areas where coastal dynamics are more dominant than coregistration errors, by classifying water with a low threshold and a high threshold (e.g. 5% and 95%), the horizontal coastal dynamic range can be retrieved. The horizontal precision of the coastline is about 4 m (Noh and Howat, 2015) according to the analysis of translational offsets, becoming more precise with the number of repeat measurements. We note that, due to the dynamic nature of the coastline as a result of tidal fluctuations and the other mechanisms defined above, the coastline changes with time (Baiocchi et al., 2012). Nevertheless, this study is not focusing on the temporal evolution of coastline (Palazzo et al., 2012), and the time dependency is mitigated through our approach of stacking 268
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 10. Water Probability Map Reveals Coastal Low-lying Area. (a) The enlargement of Fig. 9 (c) (white rectangle), which highlights a low-lying area with lower water probability. Markings are the same as in Fig. 9 (c). (b) The near infrared band of a WorldView-2 orthoimage acquired on 14 April 2013, 21:44 UTC. The colorbar is the top-of-atmosphere reflectance of NIR band, where 0 is water, and land otherwise. (c) Same as (b) but observed on 11 March 2013, 22:35 UTC. Imagery© 2013 DigitalGlobe, Inc.
multi-temporal images. In addition to the challenges addressed in the previous section, including image geolocation errors, atmospheric contamination (e.g. clouds, cloud shadows, haze), shadows and tidal dynamics, other factors may complicate the global applicability of this algorithm. For example, the spectral properties of water pixels at different wavelengths may vary according to submerged or aquatic vegetation, turbid sediment plumes, colored dissolved organic matter load, shallow bottoms material, as well as observation conditions (e.g. sun-target-sensor geometry) (Verpoorter et al., 2012, 2014; Pekel et al., 2016). Recent studies have used different techniques to address these issues in order to produce global surface water maps using Landsat satellite imagery, at a much coarser spatial resolution (30 m or 14.25 m) than our application (2 m). For example, Pekel et al. (2016) adopts machine learning techniques, first building a large spectral library containing samples of a wide range of water surface conditions through visual interpretation, and then analyzing the spectral library through visual analytics. We plan to process all images in the domain of ArcticDEM (Fig. 1) and will compare our adaptive thresholding method with Pekel et al.'s (2016) method for different water conditions, eventually improving our algorithm for global application by incorporating suitable techniques. The final data products including 2-m resolution coastline and water probability maps and supporting data will be released in future after validation.
dynamics. These capabilities will continue to improve as the archive of imagery grows. Acknowledgments This work was supported by U.S. National Science Foundation Office of Polar Programs grant A005265701 via a subcontract to the University of Minnesota and also by the NASA Cryosphere program. Geospatial support for this work provided by the Polar Geospatial Center under NSF-OPP awards 1043681 and 1559691. DEMs provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736. The imagery used in this paper was provided by the United States National Geospatial Intelligence Agency and the Polar Geospatial Center at the University of Minnesota through the NextView licensing program. The 30 m SRTM DEMs are from NASA's Earth Observing System Data and Information System (EOSIDS), and SRTM Water Body Data from U.S. Geological Survey (USGS). Some figures in this paper were generated using the Generic Mapping Tools (GMT) (Wessel and Smith, 1991). The code developed in this study is available at https://github.com/Chunli-Dai/SatelliteImageCoastline. We thank Editor Menghua Wang, Associate Editor Tiit Kutser and three anonymous reviewers for their constructive comments and suggestions. We also thank Laurie Padman (Earth & Space Research) for his contribution to tide modeling and helpful discussions about the impact of tide model's coarse spatial resolution.
5. Conclusion
Appendix A. Supplementary data
We present a coastline extraction method that improves the accuracy and resolution of water surface classification by taking advantage of a large volume of repeated, high-resolution, multispectral satellite images. By stacking repeat spectral classifications, we determine the probability of water at a given point on the coast. We employ an adaptive thresholding based on existing, lower resolution shorelines, that improves classification. Stacking of multiple classifications mitigates errors due to clouds and cloud shadows, reduces errors in geolocation and accounts for actual spatial and temporal variability in coastline, although terrain shadowing remains among the main sources of misclassification error. The improvement of the spatial resolution (2 m) and precision (< 4 m) of our derived coastline is substantial compared to the existing coastlines in our study area. Our approach enables efficient, accurate mapping of coastlines in both time and space over large areas, opening the door to a wide range of studies in coastal
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.rse.2019.04.010. References Anderson, N.T., Marchisio, G.B., 2012. WorldView-2 and the evolution of the DigitalGlobe remote sensing satellite constellation: introductory paper for the special session on WorldView-2. In: Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVIII. 8390. International Society for Optics and Photonics, pp. 83900L. Baiocchi, V., Brigante, R., Dominici, D., Radicioni, F., 2012. Coastline detection using high resolution multispectral satellite images. In: Proceedings of FIG Working Week, (Rome, Italy). Barr, I., Dokukin, M.D., Kougkoulos, I., Livingstone, S.J., Lovell, H., Małecki, J., Muraviev, A.Y., 2018. Using ArcticDEM to analyse the dimensions and dynamics of debris-covered glaciers in Kamchatka, Russia. Geosciences 8 (6).
269
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
for surface elevation change measurement using geometric constraints. IEEE Trans. Geosci. Remote Sens. 52 (4), 2247–2260. Noh, M.J., Howat, I.M., 2015. Automated stereo-photogrammetric DEM generation at high latitudes: Surface Extraction with TIN-based Search-space Minimization (SETSM) validation and demonstration over glaciated regions. GISci. Remote Sens. 52 (2), 198–217. Noh, M.J., Howat, I.M., 2018. Automatic relative RPC image model bias compensation through hierarchical image matching for improving DEM quality. ISPRS J. Photogramm. Remote Sens. 136, 120–133. Nuth, C., Kääb, A., 2011. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. The Cryosphere 5 (1), 271–290. Padman, L., Siegfried, M.R., Fricker, H.A., 2018. Ocean tide influences on the Antarctic and Greenland ice sheets. Rev. Geophys. 56 (1), 142–184. Palazzo, F., Latini, D., Baiocchi, V., Del Frate, F., Giannone, F., Dominici, D., Remondiere, S., 2012. An application of COSMO-sky med to coastal erosion studies. Eur. J. Remote Sens. 45 (1), 361–370. Pavlis, N.K., Holmes, S.A., Kenyon, S.C., Factor, J.K., 2012. The development and evaluation of the earth gravitational model 2008 (EGM2008). J. Geophys. Res. Solid Earth 117 (B4). Pekel, J.F., Cottam, A., Gorelick, N., Belward, A.S., 2016. High-resolution mapping of global surface water and its long-term changes. Nature 540 (7633), 418. Porter, C., Morin, P., Howat, I., Noh, M.-J., Bates, B., Peterman, K., Keesey, S., Schlenk, M., Gardiner, J., Tomko, K., Willis, M., Kelleher, C., Cloutier, M., Husby, E., Foga, S., Nakamura, H., Platson, M., Wethington Jr., M., Williamson, C., Bauer, G., Enos, J., Arnold, G., Kramer, W., Becker, P., Doshi, A., D’Souza, C., Cummens, P., Laurier, F., Bojesen, M., 2018. ArcticDEM. Harvard Dataverse, V1https://doi.org/10.7910/DVN/ OHHUKH , Accessed date: 30 April 2019. Slater, J.A., Garvey, G., Johnston, C., Haase, J., Heady, B., Kroenung, G., Little, J., 2006. The SRTM data “finishing” process and products. Photogramm. Eng. Remote Sens. 72 (3), 237–247. Soluri, E.A., Woodson, V.A., 1990. World vector shoreline. Int. Hydrogr. Rev. 67 (1). Stammer, D., Ray, R.D., Andersen, O.B., Arbic, B.K., Bosch, W., Carrère, L., Cheng, Y., Chinn, D.S., Dushaw, B.D., Egbert, G.D., Erofeeva, S.Y., 2014. Accuracy assessment of global barotropic ocean tide models. Rev. Geophys. 52 (3), 243–282. Vafeidis, A.T., Nicholls, R.J., McFadden, L., Tol, R.S., Hinkel, J., Spencer, T., Grashoff, P.S., Boot, G., Klein, R.J., 2008. A new global coastal database for impact and vulnerability analysis to sea-level rise. J. Coast. Res. 24 (4), 917–924. Verpoorter, C., Kutser, T., Tranvik, L., 2012. Automated mapping of water bodies using Landsat multispectral data. Limnol. Oceanogr. Methods 10 (12), 1037–1050. Verpoorter, C., Kutser, T., Seekell, D.A., Tranvik, L.J., 2014. A global inventory of lakes based on high-resolution satellite imagery. Geophys. Res. Lett. 41 (18), 6396–6402. Wang, J., Sheng, Y., Tong, T.S.D., 2014. Monitoring decadal lake dynamics across the Yangtze Basin downstream of Three Gorges Dam. Remote Sens. Environ. 152, 251–269. Wessel, P., Smith, W.H., 1991. Free software helps map and display data. Eos, Trans. Am. Geophys. Union 72 (41), 441–446. Wessel, P., Smith, W.H.F., 1996. A global self-consistent hierarchical high-resolution shoreline database. J. Geophys. Res. 101, 8741–8743. Wolf, A.F., 2012. Using WorldView-2 Vis-NIR multispectral imagery to support land mapping and feature extraction using normalized difference index ratios. In: Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVIII. vol. 8390. International Society for Optics and Photonics, pp. 83900N. Xu, H., 2006. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens. 27 (14), 3025–3033. Zhang, T., Yang, X., Hu, S., Su, F., 2013. Extraction of coastline in aquaculture coast from multispectral remote sensing images: object-based region growing integrating edge detection. Remote Sens. 5 (9), 4470–4487. Zhu, Z., Wang, S., Woodcock, C.E., 2015. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images. Remote Sens. Environ. 159, 269–277.
Chambers, D.P., 2006. Observing seasonal steric sea level variations with GRACE and satellite altimetry. J. Geophys. Res. Oceans 111 (C3), C03010. Champion, N., 2012. Automatic cloud detection from multi-temporal satellite images: towards the use of pléiades time series. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 39, B3 (Melbourne, Australia). Chen, J.L., Shum, C.K., Wilson, C.R., Chambers, D.P., Tapley, B.D., 2000. Seasonal sea level change from TOPEX/Poseidon observation and thermal contribution. J. Geod. 73 (12), 638–647. Dai, C., Durand, M., Howat, I.M., Altenau, E.H., Pavelsky, T.M., 2018. Estimating river surface elevation from ArcticDEM. Geophys. Res. Lett. 45 (7), 3107–3114. Dellepiane, S., De Laurentiis, R., Giordano, F., 2004. Coastline extraction from SAR images and a method for the evaluation of the coastline precision. Pattern Recogn. Lett. 25 (13), 1461–1470. Egbert, G.D., Erofeeva, S.Y., 2002. Efficient inverse modeling of barotropic ocean tides. J. Atmos. Ocean. Technol. 19 (2), 183–204. Farr, T.G., Rosen, P.A., Caro, E., Crippen, R., Duren, R., Hensley, S., Alsdorf, D., 2007. The shuttle radar topography mission. Rev. Geophys. 45, RG2004. https://doi.org/10. 1029/2005RG000183. Feyisa, G.L., Meilby, H., Fensholt, R., Proud, S.R., 2014. Automated water extraction index: a new technique for surface water mapping using Landsat imagery. Remote Sens. Environ. 140, 23–35. Fisher, A., Flood, N., Danaher, T., 2016. Comparing Landsat water index methods for automated water classification in eastern Australia. Remote Sens. Environ. 175, 167–182. Gonzalez, R.C., Woods, R.E., 1992. Digital Image Processing. Prentice Hall, New Jersey, pp. 598–605. Hsiao, Y.S., Hwang, C., Cheng, Y.S., Chen, L.C., Hsu, H.J., Tsai, J.H., Liu, C.L., Wang, C.C., Liu, Y.C., Kao, Y.C., 2016. High-resolution depth and coastline over major atolls of South China Sea from satellite altimetry and imagery. Remote Sens. Environ. 176, 69–83. Jawak, S.D., Luis, A.J., 2013. Improved land cover mapping using high resolution multiangle 8-band WorldView-2 satellite remote sensing data. J. Appl. Remote. Sens. 7 (1), 073573. Jeong, S., Howat, I.M., Ahn, Y., 2017. Improved multiple matching method for observing glacier motion with repeat image feature tracking. IEEE Trans. Geosci. Remote Sens. 55 (4), 2431–2441. Kuester, M., 2017. Absolute radiometric correction of DigitalGlobe products, technical note. DigitalGlobehttps://dg-cms-uploads-production.s3.amazonaws.com/uploads/ document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf, Accessed date: May 2018. Lee, J.S., Jurkevich, I., 1990. Coastline detection and tracing in SAR images. IEEE Trans. Geosci. Remote Sens. 28 (4), 662–668. Leprince, S., Barbot, S., Ayoub, F., Avouac, J.P., 2007. Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Trans. Geosci. Remote Sens. 45 (6), 1529–1558. Liu, H., Jezek, K.C., 2004. Automated extraction of coastline from satellite imagery by integrating canny edge detection and locally adaptive thresholding methods. Int. J. Remote Sens. 25 (5), 937–958. https://doi.org/10.1080/0143116031000139890. Liu, H., Sherman, D., Gu, S., 2007. Automated extraction of shorelines from airborne light detection and ranging data and accuracy assessment based on Monte Carlo simulation. J. Coast. Res. 23 (6), 1359–1369. Liu, C., Xiao, Y., Yang, J., 2017. A coastline detection method in polarimetric SAR images mixing the region-based and edge-based active contour models. IEEE Trans. Geosci. Remote Sens. 55 (7), 3735–3747. Maglione, P., Parente, C., Vallario, A., 2014. Coastline extraction using high resolution WorldView-2 satellite imagery. Eur. J. Remote Sens. 47 (1), 685–699. McFeeters, S.K., 1996. The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int. J. Remote Sens. 17 (7), 1425–1432. https:// doi.org/10.1080/01431169608948714. NASA/NGA, 2003. SRTM Water Body Data Product Specific Guidance, Version 2.0. Noh, M.J., Howat, I.M., 2014. Automated coregistration of repeat digital elevation models
270
Contents lists available at ScienceDirect
Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse
Coastline extraction from repeat high resolution satellite imagery a,⁎
a
b
c
Chunli Dai , Ian M. Howat , Eric Larour , Erik Husby a b c
T
Byrd Polar and Climate Research Center, The Ohio State University, Columbus, OH, USA Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA Polar Geospatial Center, University of Minnesota, St. Paul, MN, USA
ARTICLE INFO
ABSTRACT
Keywords: Satellite images Coastline NDWI Water classification
This paper presents a new coastline extraction method that improves water classification accuracy by benefitting from an ever-increasing volume of repeated measurements from commercial satellite missions. The widely-used Normalized Difference Water Index (NDWI) method is tested on a sample of around 12,600 satellite images for statistical analysis. The core of the new water classification method is the use of a water probability algorithm based on the stacking of repeat measurements, which can mitigate the effects of translational offsets of images and the classification errors caused by clouds and cloud shadows. By integrating QuickBird, WorldView-2 and WorldView-3 multispectral images, the final data product provides a 2 m resolution coastline, as well as a 2 m water probability map and a repeat-count measurement map. Improvements on the existing coastline (GSHHSthe Global Self-consistent, Hierarchical, High-resolution Shoreline Database, 50 m–5000 m) in terms of resolution (2 m) is substantial, thanks to the combination of multiple data sources.
1. Introduction Coastal areas are ecologically and economically important and exceptionally dynamic. Accurate shoreline mapping is essential for navigation, coastal monitoring and modeling, as well as coastal resource management and planning (Liu and Jezek, 2004; Liu et al., 2007). Changes in coastlines occur due to variations in tides and weather, estuarine and ocean circulation, geomorphological processes, as well as anthropogenic activities (urbanization and tourism) (Baiocchi et al., 2012; Palazzo et al., 2012). In particular, the coastal zone is sensitive to global sea-level rise (Vafeidis et al., 2008). Given sea-level rise is expected to accelerate in the coming decades, accurate detection and monitoring of coastlines will become essential to provide needed information for mitigation and adaptation. The most commonly used coastline datasets are the Global Selfconsistent, Hierarchical, High-resolution Shoreline Database (GSHHS) provided by Wessel and Smith (1996) and the 30-m resolution Shuttle Radar Topography Mission (SRTM) water body dataset (SWBD) (NASA/ NGA, 2003; Slater et al., 2006; Farr et al., 2007) ranging from 56° south to 60° north latitude. The GSHHS data was constructed from two military data sets (Wessel and Smith, 1996): 1) the Central Intelligence Agency's World Data Bank II, which was created by U.S. government in the 1980s with a resolution of around 500–5000 m and 2) the World Vector Shoreline (WVS), which was generated with better spatial
⁎
resolution (50–500 m) between 1985 and 1989 by the Defense Mapping Agency, now the National Geospatial-Intelligence Agency (NGA) (Soluri and Woodson, 1990; Wessel and Smith, 1996). Historically, coastlines were manually identified and traced from high-resolution aerial images by experts in cartographic applications (e.g. Liu and Jezek, 2004; Dellepiane et al., 2004). Since manual delineation is labor intensive and often subjective, several automatic coastline extraction methods have been proposed, including the edge tracing algorithm for Synthetic Aperture Radar (SAR) imagery (Lee and Jurkevich, 1990), an image segmentation method for radar and optical satellite imagery (Liu and Jezek, 2004), a coherence thresholding method for interferometric SAR (InSAR) (Dellepiane et al., 2004), a combined method integrating image segmentation, region growing, and edge detection for multispectral imagery (Zhang et al., 2013), an active contour method for Polarimetric SAR images (Liu et al., 2017), the normalized difference vegetation index (NDVI) method (Hsiao et al., 2016) and the Normalized Difference Water Index (NDWI) index method for multispectral satellite images (McFeeters, 1996; Maglione et al., 2014). In addition to the above shoreline detection from image processing, methods based on digital elevation models have also been explored (e.g. Liu et al., 2007). The ever-increasing volume of remote sensing data has a significant impact on studies of Earth surface processes and surface water changes. High resolution optical (visible to near infrared band) satellite images
Corresponding author. E-mail address: [email protected] (C. Dai).
https://doi.org/10.1016/j.rse.2019.04.010 Received 2 November 2018; Received in revised form 5 April 2019; Accepted 9 April 2019 Available online 26 April 2019 0034-4257/ © 2019 Elsevier Inc. All rights reserved.
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 1. Data coverage and test area. The colored map shows the data coverage of ArcticDEM, which includes all land area north of 60°, as well as all territory of Greenland, the State of Alaska in entirety, and the Kamchatka Peninsula of the Russian Federation. Here the surface elevation is plotted using the low resolution (1 km) ArcticDEM data. Our test areas include the white rectangle area of 200 km by 150 km in the Aleutian Islands and the entirety Iceland (white quadrilateral).
have been used successfully for analyzing and extracting surface water boundaries (Feyisa et al., 2014). For multispectral images, the NDWI is the most widely used method initially developed for Landsat images (McFeeters, 1996). Modified NDWI indices were further proposed by Xu (2006) and Feyisa et al. (2014) using shortwave infrared band in Landsat 5 Thematic Mapper (TM). Wolf (2012) modified the index specifically for WorldView-2 and 3 multispectral imageries through the use of coastal and near infrared -2 bands. Based on the widely-used NDWI and modified NDWI formulas, we propose a new water classification and shoreline mapping method that utilizes the large and growing volume of repeat, high resolution multispectral satellite imagery now available. We apply this new algorithm to the much higher resolution data provided by commercial multispectral imagers in order to achieve the most precisely mapped coastline currently available from space.
97 min. The revisit time is about 1 to 3 days with an imaging swath width of 13 to 17 km. WorldView-2 and 3 are capable of collecting around 1 million km2 of imagery per day, which is about seven per mill of the total land surface area of Earth per day. The spatial resolution of the panchromatic band is finer than that of the multispectral band, e.g. WorldView-2 image's resolution is 46 cm in panchromatic band and 1.85 m using 8 multispectral bands. Image band wavelengths are slightly different for each satellite (https://www.digitalglobe.com). We test our algorithm acquired over a portion of Alaska's Aleutian Islands, as well as all imagery of Iceland from the archive of Arctic imagery held by the Polar Geospatial Center (Fig. 1), which is the largest single holding of these imagery available to researchers. The algorithm, however, can be applied over any region where sufficient data is available. The radiometrically corrected level 1 imagery is orthorectified into a polar stereographic projection using the sensor model Rational Polynomial Coefficients (RPCs). We orthorectify the imagery using the 2.5 arc-minute Earth Gravitational Model 2008 (EGM2008) (Pavlis et al., 2012) for terrain correction. The geoid is used rather than a DEM because DEMs are either too low in resolution to provide an accurate surface at our target resolution, or tend to be noisy around the coastline, such as with ArcticDEM (Porter et al., 2018), which increases the error in our results. The imagery is then converted from digital number
2. Data sources QuickBird (2001–2015), GeoEye-1 (since 2008), WorldView-2 (since 2009) and WorldView-3 (since 2014) are commercial satellites operated by DigitalGlobe, Inc. (Anderson and Marchisio, 2012), the operator and owner of a constellation of high-resolution imaging satellites. Their orbits are sun-synchronous, with orbital period around 261
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
0.25
a)
QuickBird Ocean Land
0.2 0.15 0.1 0.05 0 -1.5
-1
-0.5
0
0.5
1
1.5
2
1
1.5
2
1
1.5
2
NDWI 0.25
b)
GeoEye-1 Ocean Land
0.2 0.15 0.1 0.05 Fig. 2. Example of water classification using the NDWI. The colormap is the NDWI calculated from a WorldView-2 image acquired on 25 September 2014 at Nakalilok Bay, Alaska (within the white rectangle in Fig. 1). The transparent background is the shaded relief image of the 5-m ArcticDEM mosaic. For the shaded image, the illumination is from the sun's azimuth direction at the actual acquisition time. The red line is the GSHHS coastline. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
0 -1.5
-1
-0.5
0
0.5
NDWI 0.25
c)
0.2
WorldView Ocean Land
0.15
to top-of-atmosphere reflectance using the absolute radiometric calibration parameters provided by DigitalGlobe (Kuester, 2017; Jawak and Luis, 2013).
0.1 0.05
3. Methods
0 -1.5
3.1. Water classification and adaptive thresholding
NIR2)/(Coastal + NIR2)
NIR)/(Green + NIR)
0
0.5
Fig. 3. Histograms of the probability distribution (the count of observations in each bin normalized by the total count) of mean NDWI from (a) 1319 QuickBird images, (b) 272 GeoEye-1 images, and (c) 11,045 WorldView-2 and 3 images. Blue is for NDWI over ocean only, and red for NDWI over land area only. The mean NDWI over ocean is 1.0, 0.4, 0.98 for all QuickBird, GeoEye-1, WorldView images, respectively; and the mean NDWI over land is −0.1, −0.2 and 0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
(1)
NDWI values as over water. This may result in misclassification (Fig. 4) due to the similarity of shadows and water in reflectance patterns (Feyisa et al., 2014). Clouds also tend to have low NDWI values, similar to land, which may result in false islands. In order to determine a water classification threshold, we evaluate the statistics of NDWI over known areas of water and land using a large number of test images. The water or land area is specified a priori based on the GSHHS coastline, which has a spatial resolution of better than 5 km (Wessel and Smith, 1996). By eroding or dilating the land area encompassed by the a priori coastline (red line in Fig. 2) with a width comparable to its uncertainty (e.g. 2 km), the land or ocean mask can be retrieved with high confidence by removing the pixels in the transitional buffer zone (Liu and Jezek, 2004). The mean and standard deviation of NDWI values within the known land and water areas (outside a 2-km buffer of the GSHHS coastline) are then calculated. The statistical analysis of NDWI is carried out using 10,816 images in Iceland and 1820 images in the 200 km by 150 km test region (Fig. 1)
where, Coastal is the top-of-atmosphere reflectance of the coastal band (band 1) and NIR2 is that of the near infrared-2 band (band 8). For the 4-band GeoEye-1, QuickBird, Ikonos, for which the Coastal and NIR2 bands are not available, the NDWI index is calculated as (McFeeters, 1996):
NDWI = (Green
-0.5
NDWI
Based on the comparison between different algorithms such as NDVI, NDWI, Spectral Angle Mapper (SAM), and Automated Water Extraction Index (AWEI) provided by Baiocchi et al. (2012) and (Fisher et al., 2016), we follow Maglione et al. (2014) and adopt the NDWI index for water classification. The NDWI index is calculated from the 8band WorldView-2 and 3 multispectral, orthorectified images (Wolf, 2012; Maglione et al., 2014) as:
NDWI = (Coastal
-1
(2)
where, Green is the top-of-atmosphere reflectance of the green band, and NIR is that of near infrared band. Water bodies are expected to yield a positive NDWI, while land surfaces should be negative (McFeeters, 1996). As shown in Fig. 2, the ocean surface has an NDWI of above 1, while land surfaces have a value averaging −0.4. The inland water bodies (lakes and rivers) also produce positive NDWI values. Water classification methods have difficulty in distinguishing water from land surfaces with low albedos, such as over asphalt, mountain shadows, and the shadows of buildings and clouds (Feyisa et al., 2014). Fig. 2 gives an example of mountain shadows that result in similar 262
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 4. Example of water masking and strip mosaicking. The water mask is generated from the NDWI map (Fig. 2 for scene 1) with a threshold of 0.3. This strip is composed of two WorldView-2 scenes acquired on 25 September 2014, where the scene alignment offset (d2 in Eq. (5)) of scenes 1 (red) and 2 (magenta) is (0, 0) m and (−0.1567 0.0685) m, respectively, in the polar stereographic coordinate system. After the correction of d2, the two scenes are merged together seamlessly. The transparent background is based on the shaded relief image of the 5-m ArcticDEM mosaic tile. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
of the Aleutian Islands, Alaska, from QuickBird, GeoEye-1, WorldView2 and 3. First, the means and standard deviations of NDWI over the known land and ocean areas are calculated for each image. Then, the probability distribution histograms of these mean values are generated. Fig. 3 shows that, generally, the land has a lower NDWI and the ocean larger NDWI. Furthermore, Fig. 3 demonstrates that NDWI can unambiguously delineate land and water for QuickBird and WorldView-2 and 3 images, where land has NDWI values around 0 and ocean has NDWI values averaging around 1.0. Overlap between values is likely due to inland lakes and/or clouds over ocean. For GeoEye-1 images, however, the separation of land and water is not clear, with substantial ocean areas having an NDWI of less than −1. In addition, the average of all NDWI standard deviations from GeoEye-1 images is much larger (3.6) than that (0.2) from QuickBird and WorldView. Thus, we chose not to not apply GeoEye-1 imagery to water classification. Based on the NDWI statistics of our test images, we use a simple adaptive thresholding method for water classification. For each individual image, the threshold, T, is calculated as:
T = (meanW + meanL)/2
2014) (Fig. S1). Due to the fact that NDWI values are very uniform over the ocean, its probability density function is almost a delta function with near zero variance, which brings the optimal threshold to be very close to the mean NDWI of ocean. Our test shows that our simple method is more practical in the application of water classification using NDWI values. Empirically, we establish the following criteria (Eqs. (4a) and (4b)) to exclude bad quality images, such as ones that have clouds or haze.
< minimum {0.5, (meanW
meanW
meanL > 0.6
meanL)/2}
(4a) (4b)
where σ is the NDWI standard deviation for each image over land or water. Any image that does not meet these criteria will be excluded. Fig. 4 shows an example of the produced high-resolution water mask, where water is denoted in blue and land in green. This example also provides a good illustration of expected misclassifications in mountain shadows.
(3)
3.2. Coregistration
where, meanW and meanL are the mean NDWI values within the known water and land areas, respectively. We compared this simple thresholding method to a widely-discussed optimal adaptive thresholding method (Gonzalez and Woods, 1992; Liu and Jezek, 2004; Wang et al.,
The images we use for coastline mapping using the NDWI index will have systematic errors in geolocation (i.e. biases in the map coordinates) of up to several meters due to errors in the sensor model used 263
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
we find that d3 follows a normal distribution (Fig. 5 (a)), with a nearzero mean of (1.5, 1.0) m and standard deviation of (3.7, 3.7) m, which is close to the sensor model geolocation uncertainty of 4 m as discussed in (Noh and Howat, 2015). The cumulative density function of the offsets, shown in Fig. 5 demonstrates how the accumulation of images can reduce the effect of individual translational offsets. As shown in Fig. 5 (b), when a group of images cover the same coastline, only a very small percentile of the images has large negative translational offsets; and the median (50th percentile) of the translational offsets is zero. This characteristic is used later for selection of the water probability threshold. 3.3. Water probability through stacking of images As shown above, the NDWI threshold may misclassify shadowed surfaces and/or clouds. Averaging of multiple, overlapping pixels using repeat imagery (i.e. stacking), however, can mitigate these errors. For any pixel, we can therefore define water probability as the ratio of the number of measurements that classify a pixel as water over the total number of measurements, i.e.:
Pw = n w / n × 100%
where Pw is the water probability, nw is the number of measurements that classify the pixel as water, and n is the total number of all measurements. Fig. 6 shows the stacking of three individual water masks (Fig. 6 (a)–(c)), and the calculated water probability (Fig. 6 (d)) based on Eq. (6). The final water classification (Fig. 7 (a)) is determined from the water probability. Given a threshold of water probability, e.g. 50%, a pixel is classified as water if Pw is equal to or larger than the threshold, and land if less than the threshold. Fig. 6 also shows how water probability can mitigate misclassification due to clouds and shadows. For the simple NDWI thresholding of an individual image, the cloud shadows introduce false lakes inland and clouds introduce false islands in ocean. When stacking all three water masks together, the cloud shadows only slightly increase the water probability, and the clouds slightly reduce the water probability (Fig. 6 (d)). As long as there are fewer cloudy than clear images, and with the proper thresholding of water probability, the merged water mask (Fig. 7 (a)) may not be contaminated by the clouds and cloud shadows. To save computation resources, the GSHHS coastline is used as a priori information to provide an approximate constraint with a ~10 km buffer, so that only pixels within the buffer are processed. Also, pixels with < 3 repeat measurements are voided. The code for retrieving coastlines from multispectral imageries is available on Github (https:// github.com/Chunli-Dai/SatelliteImageCoastline). Note that for the calculation of NDWI statistics of each individual image (Section 3.1), we only process data outside a 2-km buffer of the GSHHS coastline considering that the typical size of an image is only 13–17 km. There are two main advantages to the proposed water probability method. The first is the mitigation of both translational geolocation errors (bias) in the imagery as well as temporal variability due to coastal dynamics and local sea level. Considering each of these factors separately, as shown in Fig. 5, translational errors follow a normal distribution. Hence, the 50th percentile (median) of the translational errors is near the mean of zero. Assuming a coastline is static, using the 50% water probability threshold is effective at mitigating the translational errors, especially if the total number of measurements is large. On the other hand, if the variation of coastline is only caused by local sea level change and coastal dynamics, such as erosion, the 50th percentile reflects the median of coastline locations. Combining both the effects from translational errors and coastline dynamic variations, the 50th percentile will still likely yield the median of coastline locations (Fig. S2). One benefit of using the median coastline is that occasional coastline changes caused by extreme events, such as flooding and storm surge, generally will not affect the predicted coastline.
Fig. 5. Statistics of sensor model bias (d3 in Eq. (5)) for individual images obtained from DEM coregistration. (a) The histogram probability density distribution of d3, where px, py are the x, y components. (b) The histogram cumulative density function. The median (50%) is shown in black, which corresponds to the horizontal shift at 0 m. This evaluation is based on 249 WorldView-2 and WorldView-3 strip files.
to convert between image and map coordinates (Noh and Howat, 2014, 2015). These errors may be removed by registration to identifiable ground control points, such as surveyed landmarks, or through coregistration to other, precisely geolocated, reference imagery through automated cross correlation methods (e.g. Jeong et al., 2017; Leprince et al., 2007). The former is labor intensive, as control points must be visually located in the imagery, and the latter is often not available. If, on the other hand, the sensor model errors tend to be normally distributed, the stacking of multiple images should partly mitigate such errors in the resulting coastline map. Here we evaluate these errors, and the effect of image stacking, using the registration available from ArcticDEM. For images that are one of a stereoscopic pair used for digital elevation model (DEM) extraction, we coregister that DEM to the registered reference DEM (the 5-m ArcticDEM mosaic which is registered to ICESat satellite altimetry) to determine the three-dimensional offsets. The coregistration is achieved by aligning static surfaces between two DEMs (Noh and Howat, 2015; Dai et al., 2018; Barr et al., 2018). For alignment, we use the iterative least-squares method of Nuth and Kääb (2011). To convert the map coordinate of a point in the orthorectified image, X0, to the same point on the ArcticDEM tile, X1, we have:
X1 = X 0 + d1 + d2 + d3
(6)
(5)
where d1 is the sensor model bias correction shift of each image of the stereo pair with respect to the corresponding DEM (Noh and Howat, 2018), which is zero for the first (left) image of the stereo pair. The term d2 is the scene alignment shift applied for merging scenes into strips (Fig. 4). This term, obtained from the ArcticDEM strip metadata file, is small, with a mean value of zero and standard deviation of 0.2 m based on 40, 000 test scenes. d3 is the strip alignment offset for the 5-m ArcticDEM tile mosaicking (offset applied to align the strip DEM with the mosaic DEM tile) which is obtained through coregistration of the strip DEM to the tile DEM as described above. Based on 249 strip files over a test region (white rectangle in Fig. 1), 264
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 6. Example of water probability mapping. (a)–(c) Water masks generated based on NDWI thresholding of multispectral WorldView-2 and 3 images. (d) Calculated water probability from the three individual water masks. The transparent background is the shaded relief image of the 5-m ArcticDEM mosaic tile.
The second advantage of the water probability approach is the reduction in misclassification of water caused by clouds and shadows. Stacking reduces the impact of these errors because they are temporary anomalies, and thus are reduced or removed as outliers. Nevertheless, there may be consistent commission errors, such as mountain shadows, where images are acquired at the same time of the day and, therefore, with a similar sun angle. These commission errors will not be filtered out by the water probability method.
3.4. Cloud detection Cloud detection is an essential step in classification (Zhu et al., 2015), because clouds and cloud shadows cause errors (Champion, 2012). Specifically, for water classification from individual images, clouds over ocean may result in misclassification of water as islands; or cloud shadows across the shorelines will cause misclassification of coastlines along the shadow edges. Using multi-temporal satellite images (Champion, 2012), we analyze the absolute difference, dM, between each water mask, Mi, and the merged water mask (as a reference), Mm, to detect cloudy images. 265
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 7. Cloud Detection. (a) The merged water mask based on water probability (Fig. 6 (d)). (b) The absolute difference between the merged water mask and the individual water mask as in Fig. 6 (a). (c) The difference after removing smaller clusters. (d) Classify each cluster of difference to identify whether it is surrounded by water (blue edges) or land (green edges). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
dM = |Mi
Mm|
(7)
and discarded. A problem with this approach arises when the number of cloudy images is equivalent to the number of clear images, so that cloudy and clear classifications cannot be determined unambiguously. This approach to cloud detection, however, has the benefit of favoring no classification over false classification. As discussed in the previous section, occasional cloud contamination would not affect the final coastline mapping. Nevertheless, the elimination of cloudy images can improve the water probability map, which can be used later for applications in detecting low-laying areas or estimating coastal dynamic/ tidal range.
As shown in Fig. 7 (b), the absolute difference between the merged water mask, dM, (Fig. 7 (a)) and the individual water mask observed on 26 March 2017 (Fig. 6 (a)) highlights clouds and cloud shadows. Considering differences between images caused by seasonal variations in river width and tides, any clusters of pixels, where dM equals 1, in the difference map that have long and narrow (100 m width) spatial pattern may be removed (Fig. 7 (c)), along with clusters smaller than 1000 m2. Each cluster is classified as clouds on water or extra water on land (due to either flood inundation or false identification by shadows). Since our focus is on coastline mapping, extra water inland (clusters with green edges in Fig. 7 (d)) are discarded. If the total area of clouds above water is > 10,000 m2, the whole image will be treated as cloudy 266
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 8. Modeled ocean tidal heights derived from the 1/6° resolution TPXO-9.1 global inverse tide model. (a) Modeled peak tidal height (defined as the maximum absolute surface tidal heights over a year) in meters on 1/6° grids around Nakalilok Bay, Alaska. Black circles represent the centers of all available images in the vicinity (~20 km) of a tide model grid (white triangle). The cyan square is the boundary of a standard ArcticDEM tile size of 50 km. The colored points show the estimated tidal heights along coastlines within this tile. The tide model is downloaded from http://volkov.oce.orst.edu/tides/global.html (last accessed February 2019). (b) Modeled tidal height time series as hours of the day (Greenwich Mean Time (GMT)) at a grid point (white triangle in (a)). Blue dots are the hourly time series from November 2004 to September 2017, which is the time span of all available images in this area. The blue lines show the mean tidal heights at each hour averaged over all days. The magenta and green lines represent the tidal heights variation for the date of 14 April 2013 and 11 March 2013, respectively. Black circles are the modeled tidal heights at the time and location of all acquired satellite images as in (a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
fjord. Further, the impact of bias in sea level due to tides or other variability on the coastline estimate will depend the slope of the coast. Correcting for this bias requires projecting the coastline onto a DEM and adjusting the coastline horizontally depending on the estimated bias and surface slope. We caution that the relatively coarse resolution and/or noise in existing, global coastal DEMs, however, yields errors that are similar or greater than the tidal bias.
3.5. Tidal variation Any mapped coastline will be biased by the stage of tides at the time (s) of observation. The tide causes the sea level to oscillate around mean sea level at a set of discrete periods, predominantly diurnal (daily) or semidiurnal (twice daily) periods (Stammer et al., 2014; Padman et al., 2018). The maximum range of tide heights is around several meters and can exceed 10 m in a few coastal areas (Padman et al., 2018). Using the 1/6° resolution TPXO-9.1 global inverse tide model (an update to Egbert and Erofeeva (2002)), the tidal height at the location and acquisition time of each individual image can be calculated. As shown in Fig. 8 (a), the peak tidal heights (defined as the maximum absolute surface tidal heights over a year) in the area near Nakalilok Bay are around 2.6 m. Due to the sun-synchronous orbit configuration, DigitalGlobe satellites pass a given location on Earth at about the same time (within 2 h) of a day (black circles in Fig. 8 (b)). The modeled tidal height time series (Fig. 8 (b)) shows a large variation in tidal heights at a given time of a day over the days with observations. For the example shown in Fig. 8 (b), the image acquisition times are close to 22:00, but the tidal heights vary between −1.8 m to 1.6 m. Nevertheless, the mean tidal height on each hour (blue curve with triangles in Fig. 8 (b)) averaged over 13 years has smaller amplitude (around 0.4 m here) than the overall amplitude (~2.6 m). As discussed in Section 3.3, the 50th percentile of water probability gives the median coastline, which is corresponding to the median tidal height if ignoring the random image coregistration errors and other driving factors of dynamic coastline. Here, we adopt the median tidal height of all acquisitions at a location as an estimate of tidal bias in the estimated coastline. Fig. 8 (a) (colored points along coastline) shows the estimated tidal height of the derived coastline, which is within 1 m and is less than the peak tidal heights in the area. It is worth mentioning that when sample size is large (e.g. > 100) the median value reflects the mean value of all observed tidal variations. Since satellites pass at around the same time of the day, the mean tidal height may still be slightly biased (Fig. 8 (b)). We note that the true local variability in tidal height may not be represented at the resolution (1/6°) of the global tide model, such as tidal variations in a deep estuary or narrow
4. Results and discussion By mosaicking multiple satellite images with different areas of coverage, the final product is a 2 m coastline polyline, of similar format to other coastline datasets, such as GSHHS, and that can be readily incorporated into a Geographical Information System (GIS) database. In addition, we produce a 2-m resolution water probability map, as well as a map of the total count of repeat images used in the mapping. As shown in Fig. 9, the final results are produced on a tile-by-tile basis, with a size of 50 km by 50 km. The computation area is expanded by 2 km to avoid edge discontinuation problems. The total count of repeated measurements (Fig. 9 (b)) are calculated after removing cloudy images, where the larger the number of measurements, the more reliable the water probability map and the final coastlines. After retrieving the water probability mask, the coastlines (blue line in Fig. 9 (c)) of islands and inland lakes smaller than 10,000 m2 are discarded as noise. We plot the GSHHS and SRTM coastlines for comparison (Fig. 9 (c)). Within the study area, our coastline seems to be closer to the SRTM coastline with differences within 100 m (e.g. 2–15 m at the middle part of the saddle shape area); while it differs from the GSHHS coastline by up to 200 m. However, our coastline more clearly agrees with the coastline visible in the orthoimage acquired by WorldView-2 on 25 September 2014. The water probability map can be used to identify coastal low-lying areas, which are very sensitive to sea level variations. For example, the pink color in Fig. 10 (a) shows the lower water probability area at the coast of the Nakalilok Bay, Alaska, which implies that this area has a large variability in water coverage. The possible causes of sea level variation include ocean, pole and solid Earth tides, atmospheric 267
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 9. Final results of water probability and coastlines. (a) Water probability over a tile of 50 km by 50 km generated from the stacking of individual water masks. The probability of 0 is set to be transparent for better viewing. The green line is the GSHHS coastline. The cyan square is the boundary of a 50 km by 50 km tile (same as Fig. 8), and the black rectangle marks the region for enlargement as in c. (b) The total number of measurements. Background of (a) and (b) is the hillshade of SRTM 1 arcsec DEM. (c) The comparison of different coastlines. The yellow line is the SRTM water body dataset, and blue line is the coastline from this study. The white rectangle marks the region of Fig. 10 (d). The enlargement of (b) at the black rectangle area. Background of (c) and (d) is the 2 m panchromatic orthoimage acquired by WorldView-2 on 25 September 2014. Imagery© 2014 DigitalGlobe, Inc. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
pressure loading, steric effects introduced by temperature and salinity variations within oceans, and water mass redistribution cycles (Chen et al., 2000; Chambers, 2006). Considering the fact that the near-infrared radiation would be strongly absorbed by water and reflected by land area (McFeeters, 1996), we use the near-infrared band of orthoimages to visually demonstrate the large water coverage variation of this area. As shown in Fig. 10 (b), an over 200 m wide of land area is exposed on 14 April 2013, 21:44 UTC, which is probably caused by the low tide at the observation time (modeled tidal height of −1.1 m (Fig. 8 (b))). And Fig. 10 (c) shows the coastal inundation observed on 11 March 2013, 22:35 UTC, probably due to the high tide (model tidal height of 1.6 m) at the observation time. Thus, our water probability mapping is able to reveal areas of large variability along the coast that
are likely to be most vulnerable to changes in sea level. Furthermore, in areas where coastal dynamics are more dominant than coregistration errors, by classifying water with a low threshold and a high threshold (e.g. 5% and 95%), the horizontal coastal dynamic range can be retrieved. The horizontal precision of the coastline is about 4 m (Noh and Howat, 2015) according to the analysis of translational offsets, becoming more precise with the number of repeat measurements. We note that, due to the dynamic nature of the coastline as a result of tidal fluctuations and the other mechanisms defined above, the coastline changes with time (Baiocchi et al., 2012). Nevertheless, this study is not focusing on the temporal evolution of coastline (Palazzo et al., 2012), and the time dependency is mitigated through our approach of stacking 268
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
Fig. 10. Water Probability Map Reveals Coastal Low-lying Area. (a) The enlargement of Fig. 9 (c) (white rectangle), which highlights a low-lying area with lower water probability. Markings are the same as in Fig. 9 (c). (b) The near infrared band of a WorldView-2 orthoimage acquired on 14 April 2013, 21:44 UTC. The colorbar is the top-of-atmosphere reflectance of NIR band, where 0 is water, and land otherwise. (c) Same as (b) but observed on 11 March 2013, 22:35 UTC. Imagery© 2013 DigitalGlobe, Inc.
multi-temporal images. In addition to the challenges addressed in the previous section, including image geolocation errors, atmospheric contamination (e.g. clouds, cloud shadows, haze), shadows and tidal dynamics, other factors may complicate the global applicability of this algorithm. For example, the spectral properties of water pixels at different wavelengths may vary according to submerged or aquatic vegetation, turbid sediment plumes, colored dissolved organic matter load, shallow bottoms material, as well as observation conditions (e.g. sun-target-sensor geometry) (Verpoorter et al., 2012, 2014; Pekel et al., 2016). Recent studies have used different techniques to address these issues in order to produce global surface water maps using Landsat satellite imagery, at a much coarser spatial resolution (30 m or 14.25 m) than our application (2 m). For example, Pekel et al. (2016) adopts machine learning techniques, first building a large spectral library containing samples of a wide range of water surface conditions through visual interpretation, and then analyzing the spectral library through visual analytics. We plan to process all images in the domain of ArcticDEM (Fig. 1) and will compare our adaptive thresholding method with Pekel et al.'s (2016) method for different water conditions, eventually improving our algorithm for global application by incorporating suitable techniques. The final data products including 2-m resolution coastline and water probability maps and supporting data will be released in future after validation.
dynamics. These capabilities will continue to improve as the archive of imagery grows. Acknowledgments This work was supported by U.S. National Science Foundation Office of Polar Programs grant A005265701 via a subcontract to the University of Minnesota and also by the NASA Cryosphere program. Geospatial support for this work provided by the Polar Geospatial Center under NSF-OPP awards 1043681 and 1559691. DEMs provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736. The imagery used in this paper was provided by the United States National Geospatial Intelligence Agency and the Polar Geospatial Center at the University of Minnesota through the NextView licensing program. The 30 m SRTM DEMs are from NASA's Earth Observing System Data and Information System (EOSIDS), and SRTM Water Body Data from U.S. Geological Survey (USGS). Some figures in this paper were generated using the Generic Mapping Tools (GMT) (Wessel and Smith, 1991). The code developed in this study is available at https://github.com/Chunli-Dai/SatelliteImageCoastline. We thank Editor Menghua Wang, Associate Editor Tiit Kutser and three anonymous reviewers for their constructive comments and suggestions. We also thank Laurie Padman (Earth & Space Research) for his contribution to tide modeling and helpful discussions about the impact of tide model's coarse spatial resolution.
5. Conclusion
Appendix A. Supplementary data
We present a coastline extraction method that improves the accuracy and resolution of water surface classification by taking advantage of a large volume of repeated, high-resolution, multispectral satellite images. By stacking repeat spectral classifications, we determine the probability of water at a given point on the coast. We employ an adaptive thresholding based on existing, lower resolution shorelines, that improves classification. Stacking of multiple classifications mitigates errors due to clouds and cloud shadows, reduces errors in geolocation and accounts for actual spatial and temporal variability in coastline, although terrain shadowing remains among the main sources of misclassification error. The improvement of the spatial resolution (2 m) and precision (< 4 m) of our derived coastline is substantial compared to the existing coastlines in our study area. Our approach enables efficient, accurate mapping of coastlines in both time and space over large areas, opening the door to a wide range of studies in coastal
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.rse.2019.04.010. References Anderson, N.T., Marchisio, G.B., 2012. WorldView-2 and the evolution of the DigitalGlobe remote sensing satellite constellation: introductory paper for the special session on WorldView-2. In: Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVIII. 8390. International Society for Optics and Photonics, pp. 83900L. Baiocchi, V., Brigante, R., Dominici, D., Radicioni, F., 2012. Coastline detection using high resolution multispectral satellite images. In: Proceedings of FIG Working Week, (Rome, Italy). Barr, I., Dokukin, M.D., Kougkoulos, I., Livingstone, S.J., Lovell, H., Małecki, J., Muraviev, A.Y., 2018. Using ArcticDEM to analyse the dimensions and dynamics of debris-covered glaciers in Kamchatka, Russia. Geosciences 8 (6).
269
Remote Sensing of Environment 229 (2019) 260–270
C. Dai, et al.
for surface elevation change measurement using geometric constraints. IEEE Trans. Geosci. Remote Sens. 52 (4), 2247–2260. Noh, M.J., Howat, I.M., 2015. Automated stereo-photogrammetric DEM generation at high latitudes: Surface Extraction with TIN-based Search-space Minimization (SETSM) validation and demonstration over glaciated regions. GISci. Remote Sens. 52 (2), 198–217. Noh, M.J., Howat, I.M., 2018. Automatic relative RPC image model bias compensation through hierarchical image matching for improving DEM quality. ISPRS J. Photogramm. Remote Sens. 136, 120–133. Nuth, C., Kääb, A., 2011. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. The Cryosphere 5 (1), 271–290. Padman, L., Siegfried, M.R., Fricker, H.A., 2018. Ocean tide influences on the Antarctic and Greenland ice sheets. Rev. Geophys. 56 (1), 142–184. Palazzo, F., Latini, D., Baiocchi, V., Del Frate, F., Giannone, F., Dominici, D., Remondiere, S., 2012. An application of COSMO-sky med to coastal erosion studies. Eur. J. Remote Sens. 45 (1), 361–370. Pavlis, N.K., Holmes, S.A., Kenyon, S.C., Factor, J.K., 2012. The development and evaluation of the earth gravitational model 2008 (EGM2008). J. Geophys. Res. Solid Earth 117 (B4). Pekel, J.F., Cottam, A., Gorelick, N., Belward, A.S., 2016. High-resolution mapping of global surface water and its long-term changes. Nature 540 (7633), 418. Porter, C., Morin, P., Howat, I., Noh, M.-J., Bates, B., Peterman, K., Keesey, S., Schlenk, M., Gardiner, J., Tomko, K., Willis, M., Kelleher, C., Cloutier, M., Husby, E., Foga, S., Nakamura, H., Platson, M., Wethington Jr., M., Williamson, C., Bauer, G., Enos, J., Arnold, G., Kramer, W., Becker, P., Doshi, A., D’Souza, C., Cummens, P., Laurier, F., Bojesen, M., 2018. ArcticDEM. Harvard Dataverse, V1https://doi.org/10.7910/DVN/ OHHUKH , Accessed date: 30 April 2019. Slater, J.A., Garvey, G., Johnston, C., Haase, J., Heady, B., Kroenung, G., Little, J., 2006. The SRTM data “finishing” process and products. Photogramm. Eng. Remote Sens. 72 (3), 237–247. Soluri, E.A., Woodson, V.A., 1990. World vector shoreline. Int. Hydrogr. Rev. 67 (1). Stammer, D., Ray, R.D., Andersen, O.B., Arbic, B.K., Bosch, W., Carrère, L., Cheng, Y., Chinn, D.S., Dushaw, B.D., Egbert, G.D., Erofeeva, S.Y., 2014. Accuracy assessment of global barotropic ocean tide models. Rev. Geophys. 52 (3), 243–282. Vafeidis, A.T., Nicholls, R.J., McFadden, L., Tol, R.S., Hinkel, J., Spencer, T., Grashoff, P.S., Boot, G., Klein, R.J., 2008. A new global coastal database for impact and vulnerability analysis to sea-level rise. J. Coast. Res. 24 (4), 917–924. Verpoorter, C., Kutser, T., Tranvik, L., 2012. Automated mapping of water bodies using Landsat multispectral data. Limnol. Oceanogr. Methods 10 (12), 1037–1050. Verpoorter, C., Kutser, T., Seekell, D.A., Tranvik, L.J., 2014. A global inventory of lakes based on high-resolution satellite imagery. Geophys. Res. Lett. 41 (18), 6396–6402. Wang, J., Sheng, Y., Tong, T.S.D., 2014. Monitoring decadal lake dynamics across the Yangtze Basin downstream of Three Gorges Dam. Remote Sens. Environ. 152, 251–269. Wessel, P., Smith, W.H., 1991. Free software helps map and display data. Eos, Trans. Am. Geophys. Union 72 (41), 441–446. Wessel, P., Smith, W.H.F., 1996. A global self-consistent hierarchical high-resolution shoreline database. J. Geophys. Res. 101, 8741–8743. Wolf, A.F., 2012. Using WorldView-2 Vis-NIR multispectral imagery to support land mapping and feature extraction using normalized difference index ratios. In: Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVIII. vol. 8390. International Society for Optics and Photonics, pp. 83900N. Xu, H., 2006. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens. 27 (14), 3025–3033. Zhang, T., Yang, X., Hu, S., Su, F., 2013. Extraction of coastline in aquaculture coast from multispectral remote sensing images: object-based region growing integrating edge detection. Remote Sens. 5 (9), 4470–4487. Zhu, Z., Wang, S., Woodcock, C.E., 2015. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images. Remote Sens. Environ. 159, 269–277.
Chambers, D.P., 2006. Observing seasonal steric sea level variations with GRACE and satellite altimetry. J. Geophys. Res. Oceans 111 (C3), C03010. Champion, N., 2012. Automatic cloud detection from multi-temporal satellite images: towards the use of pléiades time series. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 39, B3 (Melbourne, Australia). Chen, J.L., Shum, C.K., Wilson, C.R., Chambers, D.P., Tapley, B.D., 2000. Seasonal sea level change from TOPEX/Poseidon observation and thermal contribution. J. Geod. 73 (12), 638–647. Dai, C., Durand, M., Howat, I.M., Altenau, E.H., Pavelsky, T.M., 2018. Estimating river surface elevation from ArcticDEM. Geophys. Res. Lett. 45 (7), 3107–3114. Dellepiane, S., De Laurentiis, R., Giordano, F., 2004. Coastline extraction from SAR images and a method for the evaluation of the coastline precision. Pattern Recogn. Lett. 25 (13), 1461–1470. Egbert, G.D., Erofeeva, S.Y., 2002. Efficient inverse modeling of barotropic ocean tides. J. Atmos. Ocean. Technol. 19 (2), 183–204. Farr, T.G., Rosen, P.A., Caro, E., Crippen, R., Duren, R., Hensley, S., Alsdorf, D., 2007. The shuttle radar topography mission. Rev. Geophys. 45, RG2004. https://doi.org/10. 1029/2005RG000183. Feyisa, G.L., Meilby, H., Fensholt, R., Proud, S.R., 2014. Automated water extraction index: a new technique for surface water mapping using Landsat imagery. Remote Sens. Environ. 140, 23–35. Fisher, A., Flood, N., Danaher, T., 2016. Comparing Landsat water index methods for automated water classification in eastern Australia. Remote Sens. Environ. 175, 167–182. Gonzalez, R.C., Woods, R.E., 1992. Digital Image Processing. Prentice Hall, New Jersey, pp. 598–605. Hsiao, Y.S., Hwang, C., Cheng, Y.S., Chen, L.C., Hsu, H.J., Tsai, J.H., Liu, C.L., Wang, C.C., Liu, Y.C., Kao, Y.C., 2016. High-resolution depth and coastline over major atolls of South China Sea from satellite altimetry and imagery. Remote Sens. Environ. 176, 69–83. Jawak, S.D., Luis, A.J., 2013. Improved land cover mapping using high resolution multiangle 8-band WorldView-2 satellite remote sensing data. J. Appl. Remote. Sens. 7 (1), 073573. Jeong, S., Howat, I.M., Ahn, Y., 2017. Improved multiple matching method for observing glacier motion with repeat image feature tracking. IEEE Trans. Geosci. Remote Sens. 55 (4), 2431–2441. Kuester, M., 2017. Absolute radiometric correction of DigitalGlobe products, technical note. DigitalGlobehttps://dg-cms-uploads-production.s3.amazonaws.com/uploads/ document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf, Accessed date: May 2018. Lee, J.S., Jurkevich, I., 1990. Coastline detection and tracing in SAR images. IEEE Trans. Geosci. Remote Sens. 28 (4), 662–668. Leprince, S., Barbot, S., Ayoub, F., Avouac, J.P., 2007. Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Trans. Geosci. Remote Sens. 45 (6), 1529–1558. Liu, H., Jezek, K.C., 2004. Automated extraction of coastline from satellite imagery by integrating canny edge detection and locally adaptive thresholding methods. Int. J. Remote Sens. 25 (5), 937–958. https://doi.org/10.1080/0143116031000139890. Liu, H., Sherman, D., Gu, S., 2007. Automated extraction of shorelines from airborne light detection and ranging data and accuracy assessment based on Monte Carlo simulation. J. Coast. Res. 23 (6), 1359–1369. Liu, C., Xiao, Y., Yang, J., 2017. A coastline detection method in polarimetric SAR images mixing the region-based and edge-based active contour models. IEEE Trans. Geosci. Remote Sens. 55 (7), 3735–3747. Maglione, P., Parente, C., Vallario, A., 2014. Coastline extraction using high resolution WorldView-2 satellite imagery. Eur. J. Remote Sens. 47 (1), 685–699. McFeeters, S.K., 1996. The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int. J. Remote Sens. 17 (7), 1425–1432. https:// doi.org/10.1080/01431169608948714. NASA/NGA, 2003. SRTM Water Body Data Product Specific Guidance, Version 2.0. Noh, M.J., Howat, I.M., 2014. Automated coregistration of repeat digital elevation models
270