Modelling intertidal elevation using tidal data
Products used: s2_l2a
Prerequisites: For more information about the tidal modelling step of this analysis, refer to the Tidal modelling notebook
Keywords: data used; Sentinel-2, water; tide modelling, water; waterline extraction, band index; NDWI, data methods; compositing
Background
Intertidal environments support important ecological habitats (e.g. sandy beaches and shores, tidal flats and rocky shores and reefs), and provide many valuable benefits such as storm surge protection, carbon storage and natural resources for recreational and commercial use. However, intertidal zones are faced with increasing threats from coastal erosion, land reclamation (e.g. port construction), and sea level rise. Accurate elevation data describing the height and shape of the coastline is needed to help predict when and where these threats will have the greatest impact. However, this data is expensive and challenging to map across the entire intertidal zone of a continent the size of Africa.
Digital Earth Africa use case
The rise and fall of the tide can be used to reveal the three-dimensional shape of the coastline by mapping the boundary betweeen water and land across a range of known tides (e.g. from low tide to high tide). Assuming that the land-water boundary is a line of constant height relative to mean sea level (MSL), elevations can be modelled for the area of coastline located between the lowest and highest observed tide.
Imagery from satellites such as the Copernicus Sentinel-2 program is available for free for the entire planet, making satellite imagery a powerful and cost-effective tool for modelling the 3D shape and structure of the intertidal zone at regional or national scale.
Reference: Bishop-Taylor et al. 2019, the first 3D model of Australia’s entire coastline: the National Intertidal Digital Elevation Model.
Description
In this example, we demonstrate a simplified version of the National Intertidal Digital Elevation Model (NIDEM) method that combines data from Sentinel-2 satellite with tidal modelling, image compositing and spatial interpolation techniques. We first map the boundary between land and water from low to high tide, and use this information to generate smooth, continuous 3D elevation maps of the intertidal zone. The resulting data may assist in mapping the habitats of threatened coastal species, identifying areas of coastal erosion, planning for extreme events such as storm surges and flooding, and improving models of how sea level rise will affect the coastline. This worked example takes users through the code required to:
Load in a cloud-free Sentinel-2 time series.
Compute a water index (NDWI).
Tag and sort satellite images by tide height.
Create “summary” or composite images that show the distribution of land and water at discrete intervals of the tidal range (e.g. at low tide, high tide).
Extract and visualise the topography of the intertidal zone as depth contours.
Interpolate depth contours into a smooth, continuous Digital Elevation Model (DEM) of the intertidal zone.
Getting started
To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.
After finishing the analysis, return to the “Analysis parameters” cell, modify some values (e.g. choose a different location or time period to analyse) and re-run the analysis. There are additional instructions on modifying the notebook at the end.
Load packages
Load key Python packages and supporting functions for the analysis.
[1]:
%matplotlib inline
import os
import datacube
import xarray as xr
import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from odc.geo.cog import write_cog
from odc.geo.geom import Geometry
from deafrica_tools.datahandling import load_ard, mostcommon_crs
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.coastal import tidal_tag
from deafrica_tools.spatial import subpixel_contours, interpolate_2d, contours_to_arrays
from deafrica_tools.plotting import display_map, map_shapefile, rgb
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.areaofinterest import define_area
Set up a Dask cluster
Dask can be used to better manage memory use down and conduct the analysis in parallel. For an introduction to using Dask with Digital Earth Africa, see the Dask notebook.
Note: We recommend opening the Dask processing window to view the different computations that are being executed; to do this, see the Dask dashboard in DE Africa section of the Dask notebook.
To use Dask, set up the local computing cluster using the cell below.
[2]:
create_local_dask_cluster()
Client
Client-546f4068-d3dd-11ef-9a93-1e01ca941b00
Connection method: Cluster object | Cluster type: distributed.LocalCluster |
Dashboard: /user/victoria@kartoza.com/proxy/8787/status |
Cluster Info
LocalCluster
771e0b39
Dashboard: /user/victoria@kartoza.com/proxy/8787/status | Workers: 1 |
Total threads: 7 | Total memory: 59.21 GiB |
Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-ab91ef3e-9e49-4210-90ef-2ca341f41cc9
Comm: tcp://127.0.0.1:40365 | Workers: 1 |
Dashboard: /user/victoria@kartoza.com/proxy/8787/status | Total threads: 7 |
Started: Just now | Total memory: 59.21 GiB |
Workers
Worker: 0
Comm: tcp://127.0.0.1:37259 | Total threads: 7 |
Dashboard: /user/victoria@kartoza.com/proxy/33033/status | Memory: 59.21 GiB |
Nanny: tcp://127.0.0.1:44349 | |
Local directory: /tmp/dask-scratch-space/worker-a5m2u_2z |
Connect to the datacube
Activate the datacube database, which provides functionality for loading and displaying stored Earth observation data.
[3]:
dc = datacube.Datacube(app='Intertidal_elevation_Sentinel-2')
Analysis parameters
The following cell sets the parameters, which define the area of interest and the length of time to conduct the analysis over. The parameters are
lat
: The central latitude to analyse (e.g.-11.384
).lon
: The central longitude to analyse (e.g.43.287
).buffer
: The number of square degrees to load around the central latitude and longitude. For reasonable loading times, set this as0.1
or lower.time_range
: The date range to analyse (e.g.('2017', '2022')
)
Select location
To define the area of interest, there are two methods available:
By specifying the latitude, longitude, and buffer. This method requires you to input the central latitude, central longitude, and the buffer value in square degrees around the center point you want to analyze. For example,
lat = 10.338
,lon = -1.055
, andbuffer = 0.1
will select an area with a radius of 0.1 square degrees around the point with coordinates (10.338, -1.055).Alternatively, you can provide separate buffer values for latitude and longitude for a rectangular area. For example,
lat = 10.338
,lon = -1.055
, andlat_buffer = 0.1
andlon_buffer = 0.08
will select a rectangular area extending 0.1 degrees north and south, and 0.08 degrees east and west from the point(10.338, -1.055)
.For reasonable loading times, set the buffer as
0.1
or lower.By uploading a polygon as a
GeoJSON or Esri Shapefile
. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files buttonin the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files
(.cpg, .dbf, .shp, .shx)
. Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.
To use one of these methods, you can uncomment the relevant line of code and comment out the other one. To comment out a line, add the "#"
symbol before the code you want to comment out. By default, the first option which defines the location using latitude, longitude, and buffer is being used.
If running the notebook for the first time, keep the default settings below. This will demonstrate how the analysis works and provide meaningful results. The example generates an intertidal elevation model for an area off the coast of Guinea Bissau.
To ensure that the tidal modelling part of this analysis works correctly, please make sure the centre of the study area is located over water when setting lat
and lon
.
[4]:
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-19.25612, lon=44.32470, buffer=0.04)
# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
#aoi = define_area(vector_path='aoi.shp')
#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)
# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])
# Set the range of dates for the analysis
time_range = ('2017', '2022')
View the selected location
The next cell will display the selected area on an interactive map. Feel free to zoom in and out to get a better understanding of the area you’ll be analysing. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point.
[5]:
display_map(x=lon_range, y=lat_range)
[5]:
Load cloud-masked Sentinel-2 data
The first step in this analysis is to load in Sentinel-2 data for the lat_range
, lon_range
and time_range
we provided above. The code below uses the load_ard
function to load in Sentinel-2 analysis-ready data for the area and time specified. For more infmation, see the Using load_ard notebook. The function will also automatically mask out clouds from the dataset, allowing us to focus on pixels that contain useful data:
[6]:
# Create the 'query' dictionary object, which contains the longitudes,
# latitudes and time provided above
query = {
'x': lon_range,
'y': lat_range,
'time': time_range,
'measurements': ['red', 'green', 'blue', 'nir'],
'resolution': (-10, 10),
}
# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='s2_l2a', query=query)
# Load available data
s2_ds = load_ard(dc=dc,
products=['s2_l2a'],
output_crs=output_crs,
min_gooddata=0.5,
group_by='solar_day',
dask_chunks={},
**query)
/home/jovyan/dev/deafrica-sandbox-notebooks/Tools/deafrica_tools/datahandling.py:244: UserWarning: Setting 'min_gooddata' percentage to > 0.0 will cause dask arrays to compute when loading pixel-quality data to calculate 'good pixel' percentage. This can slow the return of your dataset.
warnings.warn(
Using pixel quality parameters for Sentinel 2
Finding datasets
s2_l2a
Counting good quality pixels for each time step
Filtering to 642 out of 821 time steps with at least 50.0% good quality pixels
Applying pixel quality/cloud mask
Returning 642 time steps as a dask array
Once the load is complete, examine the data by printing it in the next cell. The Dimensions
argument revels the number of time steps in the data set, as well as the number of pixels in the x
(longitude) and y
(latitude) dimensions.
[7]:
print(s2_ds)
<xarray.Dataset> Size: 8GB
Dimensions: (time: 642, y: 889, x: 845)
Coordinates:
* time (time) datetime64[ns] 5kB 2017-01-07T07:23:27 ... 2022-12-27...
* y (y) float64 7kB 7.875e+06 7.875e+06 ... 7.866e+06 7.866e+06
* x (x) float64 7kB 4.248e+05 4.248e+05 ... 4.332e+05 4.333e+05
spatial_ref int32 4B 32738
Data variables:
red (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
green (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
blue (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
nir (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
Attributes:
crs: epsg:32738
grid_mapping: spatial_ref
Plot example timestep in true colour
To visualise the data, use the pre-loaded rgb
utility function to plot a true colour image for a given time-step. White areas indicate where clouds or other invalid pixels in the image have been masked.
Change the value for timestep
and re-run the cell to plot a different timestep (timesteps are numbered from 0
to n_time - 1
where n_time
is the total number of timesteps; see the time
listing under the Dimensions
category in the dataset print-out above).
[8]:
# Set the timesteps to visualise
timestep = 12
# Generate RGB plots at each timestep
rgb(s2_ds, index=timestep)

Compute Normalised Difference Water Index
To extract intertidal depth contours, we need to be able to seperate water from land in our study area. To do this, we can use our Sentinel-2 data to calculate a water index called the Normalised Difference Water Index
, or NDWI. This index uses the ratio of green and near-infrared radiation to identify the presence of water. The formula is:
where Green
is the green band and NIR
is the near-infrared band.
When it comes to interpreting the index, High values (greater than 0, blue colours) typically represent water pixels, while low values (less than 0, red colours) represent land. You can use the cell below to calculate and plot one of the images after calculating the index.
[9]:
# Calculate the water index
s2_ds = calculate_indices(s2_ds, index='NDWI',
satellite_mission ='s2')
# Plot the resulting image for the same timestep selected above
s2_ds.NDWI.isel(time=timestep).plot(cmap='RdBu',
size=6,
vmin=-0.8,
vmax=0.8)
plt.show()

How does the plot of the index compare to the optical image from earlier? Was there water or land anywhere you weren’t expecting?
Model tide heights
The location of the shoreline can vary greatly from low to high tide. In the code below, we aim to calculate the height of the tide at the exact moment each Sentinel-2 image was acquired. This will allow us to built a sorted time series of images taken at low tide to high tide, which we will use to generate the intertidal elevation model.
The tidal_tag
function below uses the OTPS TPXO8 tidal model to calculate the height of the tide at the exact moment each satellite image in our dataset was taken, and adds this as a new tide_m
attribute in our dataset (for more information about this function, refer to the Tidal modelling notebook).
Note: this function can only model tides correctly if the centre of your study area is located over water. If this isn’t the case, you can specify a custom tide modelling location by passing a coordinate to
tidepost_lat
andtidepost_lon
(e.g.tidepost_lat=11.228, tidepost_lon=-15.860
).
[10]:
# Calculate tides for each timestep in the satellite dataset
s2_ds = tidal_tag(ds=s2_ds, tidepost_lat=None, tidepost_lon=None)
# Print the output dataset with new `tide_m` variable
print(s2_ds)
Setting tide modelling location from dataset centroid: 44.32, -19.26
Modelling tides using FES2014 tidal model
<xarray.Dataset> Size: 10GB
Dimensions: (time: 642, y: 889, x: 845)
Coordinates:
* time (time) datetime64[ns] 5kB 2017-01-07T07:23:27 ... 2022-12-27...
* y (y) float64 7kB 7.875e+06 7.875e+06 ... 7.866e+06 7.866e+06
* x (x) float64 7kB 4.248e+05 4.248e+05 ... 4.332e+05 4.333e+05
spatial_ref int32 4B 32738
Data variables:
red (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
green (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
blue (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
nir (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
NDWI (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
tide_m (time) float64 5kB 0.6238 -0.9891 -0.7224 ... -1.41 0.1797
Attributes:
crs: epsg:32738
grid_mapping: spatial_ref
[11]:
# Plot the resulting tide heights for each Sentinel-2 image:
s2_ds.tide_m.plot()
plt.show()

Create water index summary images from low to high tide
Using these tide heights, we can sort our Sentinel-2 dataset by tide height to reveal which parts of the landscape are inundated or exposed from low to high tide.
Individual remote sensing images can be affected by noise, including clouds, sunglint and poor water quality conditions (e.g. sediment). To produce cleaner images that can be compared more easily between tidal stages, we can create ‘summary’ images or composites that combine multiple images into one image to reveal the ‘typical’ or median appearance of the landscape at different tidal stages. In this case, we use the median as the summary statistic because it prevents strong outliers (like stray clouds) from skewing the data, which would not be the case if we were to use the mean.
In the code below, we take the time series of images, sort by tide and categorise each image into 9 discrete tidal intervals, ranging from the lowest (tidal interval 1) to the highest tides observed by Sentinel-2 (tidal interval 9). For more information on this method, refer to Sagar et al. 2018.
[12]:
# Sort every image by tide height
s2_ds = s2_ds.sortby('tide_m')
# Bin tide heights into 9 tidal intervals from low (1) to high tide (9)
binInterval = np.linspace(s2_ds.tide_m.min(),
s2_ds.tide_m.max(),
num=10)
tide_intervals = pd.cut(s2_ds.tide_m,
bins=binInterval,
labels=range(1, 10),
include_lowest=True)
# Add interval to dataset
s2_ds['tide_interval'] = xr.DataArray(tide_intervals.astype(int),
[('time', s2_ds.time.data)])
print(s2_ds)
<xarray.Dataset> Size: 10GB
Dimensions: (time: 642, y: 889, x: 845)
Coordinates:
* time (time) datetime64[ns] 5kB 2019-09-29T07:25:07 ... 2020-07-...
* y (y) float64 7kB 7.875e+06 7.875e+06 ... 7.866e+06 7.866e+06
* x (x) float64 7kB 4.248e+05 4.248e+05 ... 4.332e+05 4.333e+05
spatial_ref int32 4B 32738
Data variables:
red (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
green (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
blue (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
nir (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
NDWI (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
tide_m (time) float64 5kB -1.995 -1.911 -1.901 ... 1.001 1.013
tide_interval (time) int64 5kB 1 1 1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9 9 9
Attributes:
crs: epsg:32738
grid_mapping: spatial_ref
[13]:
s2_ds.sortby('time').tide_m.plot()
for i in binInterval: plt.axhline(i, c='black', alpha=0.5)

Now that we have a dataset where each image is classified into a discrete range of the tide, we can combine our images into a set of nine individual images that show where land and water is located from low to high tide. This step can take several minutes to process.
Note: We recommend opening the Dask processing window to view the different computations that are being executed; to do this, see the Dask dashboard in DE Africa section of the Dask notebook.
[14]:
# For each interval, compute the median water index and tide height value
s2_intervals = (s2_ds[['tide_interval', 'NDWI', 'tide_m']]
.compute()
.groupby('tide_interval')
.median(dim='time'))
# Plot the resulting set of tidal intervals
s2_intervals.NDWI.plot(col='tide_interval', col_wrap=5, cmap='RdBu')
plt.show()
/opt/venv/lib/python3.12/site-packages/xarray/core/concat.py:540: UserWarning: No index created for dimension tide_interval because variable tide_interval is not a coordinate. To create an index for tide_interval, please first call `.set_coords('tide_interval')` on this object.
ds.expand_dims(dim_name, create_index_for_new_dim=create_index_for_new_dim)

The plot above should make it clear how the shape and structure of the coastline changes significantly from low to high tide as low-lying tidal flats are quickly inundated by increasing water levels.
Extract depth contours from imagery
We now want to extract an accurate boundary between land and water for each of the tidal intervals above. The code below identifies the depth contours based on the boundary between land and water by tracing a line along pixels with a water index value of 0
(the boundary between land and water water index values). It returns a geopandas.GeoDataFrame
with one depth contour for each tidal interval that is labelled with tide heights in metres relative to Mean Sea Level.
[15]:
# Set up attributes to assign to each waterline
attribute_df = pd.DataFrame({'tide_m': s2_intervals.tide_m.values})
# Extract waterlines
contours_gdf = subpixel_contours(da=s2_intervals.NDWI,
z_values=0,
crs=s2_ds.geobox.crs,
attribute_df=attribute_df,
min_vertices=20,
dim='tide_interval')
# Plot output shapefile over the top of the first tidal interval water index
fig, ax = plt.subplots(1, 1, figsize=(15, 10))
s2_intervals.NDWI.sel(tide_interval=1).plot(ax=ax,
cmap='Greys',
add_colorbar=False)
contours_gdf.plot(ax=ax, column='tide_m', cmap='YlOrRd', legend=True)
plt.show()
Operating in single z-value, multiple arrays mode

The above plot is a basic visualisation of the depth contours returned by the subpixel_contours
function. Deeper contours (in m relative to Mean Sea Level) are coloured in yellow; more shallow contours are coloured in red. Now have the shapefile, we can use a more complex function to make an interactive plot for viewing the topography of the intertidal zone.
Plot interactive map of depth contours coloured by time
The next cell provides an interactive map with an overlay of the depth contours identified in the previous cell. Run it to view the map.
Zoom in to the map below to explore the resulting set of depth contours. Using this data, we can easily identify areas of the coastline which are only exposed in the lowest of tides, or other areas that are only covered by water during high tides.
[16]:
map_shapefile(gdf=contours_gdf, attribute='tide_m', continuous=True)
While the contours above provide valuable information about the topography of the intertidal zone, we can extract additional information about the 3D structure of the coastline by converting them into an elevation raster (i.e. a Digital Elevation Model or DEM).
In the cell below, we convert the shapefile above into an array of points with X, Y and Z coordinates, where the Z coordinate is the point’s elevation relative to Mean Sea Level. We then use these XYZ points to interpolate smooth, continuous elevations across the intertidal zone using linear interpolation.
[17]:
# First convert our contours shapefile into an array of XYZ points
xyz_array = contours_to_arrays(contours_gdf, 'tide_m')
# Interpolate these XYZ points over the spatial extent of the Sentinel-2 dataset
intertidal_dem = interpolate_2d(ds=s2_intervals,
x_coords=xyz_array[:, 0],
y_coords=xyz_array[:, 1],
z_coords=xyz_array[:, 2])
# Plot the output
intertidal_dem.plot(cmap='viridis', size=8, robust=True)
plt.show()

You can see in the output above that our interpolation results are very messy. This is because the interpolation extends across areas of our study area that are not affected by tides (e.g. areas of water located beyond the lowest observed tide, and on land). To clean up the data, we can restrict the DEM to only the area between the lowest and highest observed tides:
[18]:
# Identify areas that are always wet (e.g. below low tide), or always dry
above_lowest = s2_intervals.isel(tide_interval=0).NDWI < 0
below_highest = s2_intervals.isel(tide_interval=-1).NDWI > 0
# Keep only pixels between high and low tide
intertidal_dem_clean = intertidal_dem.where(above_lowest & below_highest)
# Plot the cleaned dataset
intertidal_dem_clean.plot(cmap='viridis', size=8, robust=True)
plt.show()

Export intertidal DEM as a GeoTIFF
As a final step, we can take the intertidal DEM we created and export it as a GeoTIFF that can be loaded in GIS software like QGIS or ArcMap (to download the dataset from the Sandbox, locate it in the file browser to the left, right click on the file, and select “Download”).
[19]:
# Export as a GeoTIFF
write_cog(intertidal_dem_clean,
fname='intertidal_dem_s2.tif',
overwrite=True)
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option WIDTH
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option HEIGHT
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option COUNT
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option DTYPE
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option CRS
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option TRANSFORM
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option NODATA
[19]:
PosixPath('intertidal_dem_s2.tif')
Additional information
License The code in this notebook is licensed under the Apache License, Version 2.0.
Digital Earth Africa data is licensed under the Creative Commons by Attribution 4.0 license.
Contact If you need assistance, please post a question on the DE Africa Slack channel or on the GIS Stack Exchange using the open-data-cube
tag (you can view previously asked questions here).
If you would like to report an issue with this notebook, you can file one on Github.
Compatible datacube version
[20]:
print(datacube.__version__)
1.8.20
[21]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[21]:
'2025-01-16'