Automated Coastline Extraction from Satellite Images using Google Earth Engine

By ujaval | January 18, 2023November 22, 2023

In this article, I will outline a method for extracting shoreline from satellite images in Google Earth Engine. This method is scalable and automatically extracts the coastline as a vector polyline. The full code link is available at the end of the post.

The method involves the following steps

Create a Cloud-free Composite Image
Extract All Waterbodies
Remove Inland Water and Small Islands
Convert Raster to Vector
Simplify and Extract Coastline

We will go through the details of each step and review the Google Earth Engine API code required to achieve the results.

Create a Cloud-free Composite Image

To be able to detect the shoreline correctly, we need cloud free images. The recommended way to do this is to create a composite image. You can select a suitable date range (monthly/yearly), apply a cloud mask to remove cloudy pixels and then create a median composite to achieve a cloud-free image.

Tip: Composites for Very Cloudy Regions
If you are working in a region that is very cloudy, a median composite may still have clouds. In such cases, try creating a percentile composite (i.e. a 25-percentile composite) which will give you better results. Learn how to create percentile composites in Chapter F4.1 of EEFA Book.

You can use any multispectral sensor to extract coastline. It is useful to have at least a NIR band and preferably a SWIR band to detect water reliably. Landsat and Sentinel-2 are the most widely used datasets for this task.

In this example, we start with Landsat Level 2 (Surface Reflectance) collection, filter to images collected during the chosen year over the selected region, apply cloud-mask and then create a median composite.

// Define parameters
var startDate = ee.Date('2020-01-01');
var endDate = startDate.advance(1, 'year');

var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

// Apply cloud mask
var maskL8sr = function(image) {
  var qaMask = image.select('QA_PIXEL')
    .bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.')
    .multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*')
    .multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones
  // and apply the masks.
  return image.addBands(opticalBands, null, true)
      .addBands(thermalBands, null, true)
      .updateMask(qaMask)
      .updateMask(saturationMask);
};
var collection = collection.map(maskL8sr);

// Select and Rename Bands
var collection = collection.select(
  ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'],
  ['blue',  'green', 'red',   'nir',   'swir1', 'swir2']
);

// Create a median composite and clip
var composite = collection.median().clip(geometry);


var rgbVis = {
  bands: ['red', 'green', 'blue'],
  min: 0.0, 
  max: 0.3
};
Map.addLayer(composite, rgbVis, 'Composite Image');

Warning: Dealing with Tidal Heights
A median composite of all images from a year will represent the average coastline position from all observations – including high and low tides. Such composites are not suitable for comparing changes in coastline over multiple years. To compare changes over time – you must ensure that the images are in the same phase of the tide. You can read this paper for an approach to mask high and low tide pixels.

Another approach is to use tidal estimates and filter the images which are captured in the same phase of the tide. The source will depend on the country and area of interest.

– HYCOM surface_elevation dataset might be useful in deriving tide information

If you know of any other datasets, please let me know in the comments.

Extract All Waterbodies

There are several methods for detecting and extracting surface water from satellite images. We want an automated method without needing to identify a threshold manually – so we will use Otsu Dynamic Thresholding. We will use the implementation from Gennadii Donchyts’s package users/gena/packages:thresholding.

Here we compute Automated Water Extraction Index (AWEI) and use Otsu Thresholding to find optimal threshold for water pixels in the given image. You may also try using other indices such as NDWI or MNDWI which may work better in your region of interest. If this method does not give satisfactory results, I recommend you look at Unsupervised Clustering method which combines different indices and uses machine learning approach to automatically identify water pixels.

// Code to compute AWEI and detect water 
// using Otsu thresholding
var detectWater = function(image) {
  var awei = image.expression(
    '4 * (GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)', {
      'GREEN': image.select('green'),
      'NIR': image.select('nir'),
      'SWIR1': image.select('swir1'),
      'SWIR2': image.select('swir2'),
  }).rename('awei');
  
  // Otsu Thresholding
  var thresholding = require(
    'users/gena/packages:thresholding');
  var scale = 100;
  var bounds = geometry;
  var cannyThreshold = 0.7;
  var cannySigma = 1;
  var minValue = -0.2;
  var th = thresholding.computeThresholdUsingOtsu(
    awei, scale, bounds, 
    cannyThreshold, cannySigma, minValue);
  // Create a Land-Water Image using Otsu Threshold
  // You can replace th with a manual threshold if
  // Otsu results are not satisfactory
  var water = awei.gt(th).rename('water');
  
  return water;
};
var water = detectWater(composite);
var waterVis = {min:0, max:1, palette: ['white', 'blue']};
Map.addLayer(water, waterVis, 'All Water');

Remove Inland Water and Small Islands

Our primary interest is the Land-Water Boundary, so we remove all inland waterbodies and small islands in the ocean. This is done by removing all clusters of water pixels smaller than a chosen size.

// Export resolution in meters is at which the coastline will
// be vectorized.
// Higher resolution (such as 10) will give smooother results
var exportResolutionMeters = 30;

// We need a threshold in meters to remove inland waterbodies
// and coastal islands. If you have a large waterbody close
// to the coastline, set this higher.
var waterbodySizeMeters = 2000;
var islandSizeMeters = 1000;

// This function takes a binary Land-Water image and
// removes inland water and small islands
function removeInlandWaterAndIslands(waterImage) {
  // reduceConnectedComponents expects an interger image
  waterImage = waterImage.int();
  
  // Define neighborhood based on user parameters
  var connectedPixelsLand = ee.Number(waterbodySizeMeters)
    .divide(exportResolutionMeters).int();
    
  var connectedPixelsWater = ee.Number(islandSizeMeters)
    .divide(exportResolutionMeters).int();

  // Remove inland water
  var landFilled = waterImage.addBands(waterImage)
   .reduceConnectedComponents(
     ee.Reducer.median(), 'water', connectedPixelsLand)
   .unmask(99).eq(99).and(waterImage.neq(0));
  
  // Remove small islands  
  var waterFilled = landFilled.addBands(landFilled)
    .reduceConnectedComponents(
      ee.Reducer.median(), 'water_1', connectedPixelsWater)
    .unmask(99).eq(99).and(landFilled.neq(1));    
  
  // Land-Water Boundary
  return waterFilled;
}
var landWaterBoundary = removeInlandWaterAndIslands(water);
var landWaterBoundaryVis = {
  min:0,
  max:1,
  palette: ['blue', 'white']
};

Map.addLayer(landWaterBoundary, landWaterBoundaryVis,
  'Land-Water Boundary (Raster)');

Convert Raster to Vector

Next, we vectorize the binary land-water image to obtain a polygon outline of the land area. This can also serve as the final output for many types of analysis. So you may export this as a shapefile and use it further.

// Convert the coastline image to vector
var vectors = ee.Image(landWaterBoundary).selfMask()
  .reduceToVectors({
    geometry: geometry,
    scale: exportResolutionMeters,
    eightConnected: true,
    maxPixels: 1e10,
    tileScale: 16
  });

Map.addLayer(vectors, {color: 'blue'},
  'Land-Water Boundary (Vector)');

Simplify and Extract Coastline

Many applications will require the resulting geometry as a line instead of polygon. So we take the resulting polygon, simplify it and extract the coordinates. We then use the coordinates to construct ee.Geometry.MultiLineString() features.

// This function takes vectorized polygons and 
// extracts a polyline
var simplifyAndExtractCoastline = function(vectors){
  // Simplify vectors
  var processedVectors = vectors.map(function(f) {
    var coords = f.geometry()
      .simplify({maxError: exportResolutionMeters})
      .coordinates();
    
    // Buffer the geometry by a pixel to avoid rasterizing
    // the boundary polygon
    var bufferDistance = ee.Number(
      exportResolutionMeters).multiply(-1);
    return f
      .setGeometry(
        ee.Geometry.MultiLineString(coords)
          .intersection(geometry.buffer(bufferDistance)));
  });
  return processedVectors;
};

var coastlineVector = simplifyAndExtractCoastline(
    vectors);
Map.addLayer(coastlineVector, {color: 'red'},
  'Coastline (Vector)');

Export the Results

The code editor environment of GEE is limited to running interactive computations with a limit of 5 minutes. If you are trying to extract coastlines from a large region, you may you may get tile or memory errors. This simply means you can’t visualize the results interactively. To get your results in such cases, you have to Export the results. Exports can run for longer duration and have access to a larger pool of resources. You may export the results to Google Drive or even to your own Asset storage.

// Exports
// Exports can run for longer time and have more resources
// If any of your layers time-out or give tile errors,
// Run the exports instead
Export.table.toDrive({
  collection: coastlineVector,
  description: 'Extracted_Coastline_Vector',
  folder: 'earthengine',
  fileNamePrefix: 'coastline',
  fileFormat: 'SHP'})

Export.image.toDrive({
  image: landWaterBoundary,
  description: 'Extracted_Land_Water_boundary_Raster',
  folder: 'earthengine',
  fileNamePrefix: 'land_water_boundary_raster',
  region: geometry,
  scale: exportResolutionMeters,
  maxPixels: 1e10})
  
Export.image.toDrive({
  image: composite,
  description: 'Composite',
  folder: 'earthengine',
  fileNamePrefix: 'composite',
  region: geometry,
  scale: exportResolutionMeters,
  maxPixels: 1e10})

Here’s the full script for automated extraction of shoreline from Landsat images. You may use and adapt the code for different sensors or different methodologies.

If you are new to Earth Engine and want to learn how to apply it to such problems, check out my course Google Earth Engine for Water Resources Management.