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.
> I am currently not aware of any open global tidal estimate datasets that have this information
Sea surface heights datasets show tides. This is available to some extend. Though might lack spacial/time resolution.
https://developers.google.com/earth-engine/datasets/catalog/HYCOM_sea_surface_elevation
The Sentinel-3 mission SSH goal accuracy is 10 cm for Near Real-Time (NRT) products and 3-5 cm for Slow Time Critical (STC) products.
Thank you. The HYCOM product might work. I am going to explore this.
Pls how do I get this as a document to print and use for practice?
Thanks
In Chrome, you can just do Ctrl+P and print as a PDF.
Hello Mr. Ujaval. Thankyou for this article. is there any automated process you recommend for extracting the coastline from SAR images (eg. sentinel-1) ??
Once you are able to detect water, the process is the same for extracting the coastline. Detecting water is much easier with SAR than with optical images. Here’s a starter script to use a static threshold to extract water. You can just as easily use Otsu-thresholding I used in this post instead of a static threshold.
https://code.earthengine.google.co.in/7507e30acc01095f8fac34d989aefa90
Thanks for the suggestion. I shall try the same.
Lucid as well as Brilliant article! Thank you Sir❤️.
Merci, j ai une erorr waterbodySizeMeters is not defined
The complete script is linked at the bottom. https://code.earthengine.google.co.in/08ac8b9630010644ac4783f0f0645b0a
what if I want to compare the value of shoreline changes over time?
Put everything in a function and map() it over a collection of images. If you are new to map/reduce, read the Section F4 of the EEFA Book https://www.eefabook.org/
Thanks for the code Mr. Ujaval, its impressive, i’m new in GEE, and i was trying to use this script with landsat 4, 5 and 7 images but i dont know how to make the necesary variations in order for it to perform correctly, how can i make it work? thanks in advance for your attention
There are GEE collections for different Landsat satellites. You need to replace the Landsat 8 collection with the appropriate collection from another satellite. You also need to change the cloud mask function that is appropriate for the satellite.
Here’s an example for Landsat 5. https://code.earthengine.google.co.in/aa3ae8052feb2335300a9b29f26b6e96
Hi Andrea, did you ever manage to tweak this? I am having a similar challenge? Please advice.
Please see this script to learn how to create landsat composites from Landsat 5,7, and 8. https://courses.spatialthoughts.com/end-to-end-gee-supplement.html#harmonized-landsat-time-series
Thanks for the code, this is very useful. I’m new in GEE. How if I want change the image with Sentinel 2 and filter with tidal height? How can I make it work? Thank you.
Using Sentinel-2 is easy. Once you learn the basics of GEE, you will be able to use the S2 collection instead of L8. Some example at https://courses.spatialthoughts.com/gee-water-resources-management.html#calculating-indices
Tidal heights is not easy and I don’t have any example code to use it.
Hi sir. First of all thank you so much for easing things out. It has been greatly useful in my project. I want to download the data for the year 2000 and 2010 as well. What changes am I supposed to do in the code? Can you guide me through it please. Thank you in anticipation for your kind gesture.
Landsat 8 (which this code uses) was launched in 2013. You need to change the collection to appropriate Landsat satellite for the time periods. See https://www.usgs.gov/landsat-missions/landsat-satellite-missions
I get an error when smoothening the output. Is there any other possibilty to straighten/smoothen the output vectors?
Dear Mr. Ujaval,
Thank you for providing the code. I attempted to modify it to use the NDWI algorithm and was successful in generating the desired output on GEE. However, I encountered an issue when attempting to download the shapefile to my drive. The error message I received was “Error: The geometry type ‘Empty’ cannot be stored in a shapefile. (Error code: 3)”. Could you please provide any guidance on how to resolve this issue?
Thanks. This is quite useful.
HI Ujval, is there anyway to use the example here- https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/main/Real_world_examples/Coastal_erosion.ipynb for tidal height correction using Landsat?
Their tidal model is from here https://github.com/digitalearthafrica/deafrica-coastlines/wiki/Setting-up-tidal-models-for-DE-Africa-Coastlines . I think it can be extracted and uploaded to GEE to be used. Worth a try.
Hello Mr. Ujaval. Thank you for sharing. This solution is impressive and a very usefull! Is there any automated process you recommend for extracting the coastline from Landsat 5? Or maybe to adapt this script in some way? I’m Interested in 1994 to 2014 period of time.
Sure. You just need to switch to appropriate Landsat collection. I have some code to get a harmonized landsat series at https://courses.spatialthoughts.com/end-to-end-gee-supplement.html#harmonized-landsat-time-series. You’ll have to adapt it and use the method shown here.
How to generate/export vector and raster data of land cover from 2012-2023 from Google Earth engine
Thanks
You can use MODIS Landcover product. Watch our video series how to work with landcover data, convert to vector and export.
(Watch this first) Landcover Analysis using ESA WorldCover: https://www.youtube.com/playlist?list=PLppGmFLhQ1HLN6ivbJJHLZbEYLdq5zgns
(Watch this next) Extracting Time Series using MODIS Landcover: https://www.youtube.com/playlist?list=PLppGmFLhQ1HLN6ivbJJHLZbEYLdq5zgns
Hello Mr Ujaval.
Thank you for your sharing.
I have a question, how if we use Landsat-7 and Landsat-5? Should use different method?
Use the appropriate collection and cloud mask. Code is available at https://courses.spatialthoughts.com/end-to-end-gee-supplement.html#harmonized-landsat-time-series
Hello, Mr Ujaval, Your blog is very helpful to me. I want to extract the river boundary by your method. What is the difference from extracting the coastline? What parts need to be modified or adjusted?
Should be mostly the same. Except the last part where you can just get the river as a polygon rather than converting it to polyline. You may also cjeck out this project that does river width extraction using GEE https://github.com/seanyx/RivWidthCloudPaper
can i have a detail of the package about OTSU thank uuuuuu!
The code is at https://code.earthengine.google.co.in/?scriptPath=users%2Fgena%2Fpackages%3Athresholding
Method is described in the paper http://www.mdpi.com/2072-4292/8/5/386