Modelling intertidal elevation using tidal data a92a1db208ba45bfbedef29287932c45

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 Australia.

Digital Earth Australia 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 NASA/USGS Landsat 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. Recently, Geoscience Australia combined 30 years of Landsat data from the Digital Earth Australia archive with tidal modelling to produce the first 3D model of Australia’s entire coastline: the DEA Intertidal Elevation dataset (also known as NIDEM; for more information, see Bishop-Taylor et al. 2019 or access the dataset here).

Description

In this example, we demonstrate a simplified version of the NIDEM method that combines data from the Landsat 5, 7 and 8 satellites 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 Australian coastline. This worked example takes users through the code required to:

  1. Load in a cloud-free Landsat time series

  2. Compute a water index (NDWI)

  3. Tag and sort satellite images by tide height

  4. 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)

  5. Extract and visualise the topography of the intertidal zone as depth contours

  6. 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 datacube
import sys
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datacube.utils.cog import write_cog

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.datahandling import load_ard, mostcommon_crs
from dea_tools.plotting import rgb, display_map
from dea_tools.bandindices import calculate_indices
from dea_tools.spatial import subpixel_contours, interpolate_2d, contours_to_arrays
from dea_tools.coastal import tidal_tag
from dea_tools.dask import create_local_dask_cluster

# Create local dask cluster to improve data load time
client = create_local_dask_cluster(return_client=True)

/env/lib/python3.8/site-packages/distributed/node.py:180: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 45189 instead
  warnings.warn(

Client

Client-1d05aeee-63c4-11ed-844a-62d0465acdaa

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /user/robbi.bishoptaylor@ga.gov.au/proxy/45189/status

Cluster Info

Connect to the datacube

Activate the datacube database, which provides functionality for loading and displaying stored Earth observation data.

[2]:

dc = datacube.Datacube(app='Intertidal_elevation')

Analysis parameters

The following cell set important parameters for the analysis:

  • lat_range: The latitude range to analyse (e.g. (-27.60, -27.665)). For reasonable load times, keep this to a range of ~0.1 degrees or less.

  • lon_range: The longitude range to analyse (e.g. (153.33, 153.425)). For reasonable load times, keep this to a range of ~0.1 degrees or less.

  • time_range: The date range to analyse (e.g. ('2000', '2019'))

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 intertidal flats in Moreton Bay south of Brisbane, Australia.

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_range and lon_range.

[3]:

lat_range = (-18.14, -18.26)
lon_range = (122.15, 122.28)
time_range = ('2010', '2020')

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.

[4]:

display_map(x=lon_range, y=lat_range)

[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load cloud-masked Landsat data

The first step in this analysis is to load in Landsat data for the lat_range, lon_range and time_range we provided above. The code below first connects to the datacube database, and then uses the load_ard function to load in data from the Landsat 5, 7 and 8 satellites for the area and time included in lat_range, lon_range and time_range. We set min_gooddata=0.9 to load only satellite images containing less than 10% cloud, allowing us to focus on satellite images that contain useful data.

[5]:

# Create the 'query' dictionary object, which contains the longitudes,
# latitudes and time provided above
query = {
    'y': lat_range,
    'x': lon_range,
    'time': time_range,
    'measurements': ['nbart_red', 'nbart_green', 'nbart_blue', 'nbart_nir'],
    'resolution': (-30, 30),
}

# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='ga_ls8c_ard_3', query=query)

# Load available data from all three Landsat satellites
landsat_ds = load_ard(dc=dc,
                      products=['ga_ls5t_ard_3',
                                'ga_ls7e_ard_3',
                                'ga_ls8c_ard_3'],
                      output_crs=output_crs,
                      align=(15, 15),
                      group_by='solar_day',
                      min_gooddata=0.9,
                      dask_chunks={'time': 1, 'x': 2000, 'y': 2000},
                      **query)

/home/jovyan/Robbi/dea-notebooks/Real_world_examples/../Tools/dea_tools/datahandling.py:229: 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("Setting 'min_gooddata' percentage to > 0.0 "

Finding datasets
    ga_ls5t_ard_3
    ga_ls7e_ard_3
    ga_ls8c_ard_3
Counting good quality pixels for each time step

CPLReleaseMutex: Error = 1 (Operation not permitted)

Filtering to 120 out of 829 time steps with at least 90.0% good quality pixels
Applying pixel quality/cloud mask
Returning 120 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.

[6]:

landsat_ds

[6]:

<xarray.Dataset>
Dimensions:      (time: 120, y: 445, x: 461)
Coordinates:
  * time         (time) datetime64[ns] 2010-01-17T01:40:18.801940 ... 2020-12...
  * y            (y) float64 -2.006e+06 -2.006e+06 ... -2.019e+06 -2.019e+06
  * x            (x) float64 4.101e+05 4.101e+05 ... 4.239e+05 4.239e+05
    spatial_ref  int32 32651
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_nir    (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
Attributes:
    crs:           epsg:32651
    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).

[7]:

# Set the timesteps to visualise
timestep = 1

# Generate RGB plots at each timestep
rgb(landsat_ds, index=timestep)
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_16_0.png

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 Landsat 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:

\[\begin{aligned} \text{NDWI} &= \frac{(\text{Green} - \text{NIR})}{(\text{Green} + \text{NIR})} \end{aligned}\]

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.

[8]:

# Calculate the water index
landsat_ds = calculate_indices(landsat_ds,
                               index='NDWI',
                               collection='ga_ls_3')

# Plot the resulting image for the same timestep selected above
landsat_ds.NDWI.isel(time=timestep).plot(cmap='RdBu',
                                         size=6,
                                         vmin=-0.8,
                                         vmax=0.8)

[8]:

<matplotlib.collections.QuadMesh at 0x7f8742b98a00>
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_18_1.png

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 Landsat 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 FES2014 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 and tidepost_lon (e.g. tidepost_lat=-27.73, tidepost_lon=153.46).

[9]:

# Calculate tides for each timestep in the satellite dataset
landsat_ds = tidal_tag(ds=landsat_ds, tidepost_lat=None, tidepost_lon=None)

# Print the output dataset with new `tide_m` variable
landsat_ds

Setting tide modelling location from dataset centroid: 122.22, -18.20
Modelling tides using FES2014 tidal model

[9]:

<xarray.Dataset>
Dimensions:      (time: 120, y: 445, x: 461)
Coordinates:
  * time         (time) datetime64[ns] 2010-01-17T01:40:18.801940 ... 2020-12...
  * y            (y) float64 -2.006e+06 -2.006e+06 ... -2.019e+06 -2.019e+06
  * x            (x) float64 4.101e+05 4.101e+05 ... 4.239e+05 4.239e+05
    spatial_ref  int32 32651
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_nir    (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    NDWI         (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    tide_m       (time) float64 0.5997 -0.4498 2.84 ... 3.098 3.141 -0.6279
Attributes:
    crs:           epsg:32651
    grid_mapping:  spatial_ref

Now that we have modelled tide heights, we can plot them to visualise the range of tide that was captured by Landsat across our time series:

[10]:

# Plot the resulting tide heights for each Landsat image:
landsat_ds.sortby('time').tide_m.plot()
plt.show()
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_23_0.png

Create water index summary images from low to high tide

Using these tide heights, we can sort our Landsat 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 Landsat (tidal interval 9). For more information on this method, refer to Sagar et al. 2018.

[11]:

# Sort every image by tide height
landsat_ds = landsat_ds.sortby('tide_m')

# Bin tide heights into 9 tidal intervals from low (1) to high tide (9)
binInterval = np.linspace(landsat_ds.tide_m.min(),
                          landsat_ds.tide_m.max(),
                          num=10)
tide_intervals = pd.cut(landsat_ds.tide_m,
                        bins=binInterval,
                        labels=range(1, 10),
                        include_lowest=True)

# Add interval to dataset
landsat_ds['tide_interval'] = xr.DataArray(tide_intervals,
                                           coords=[landsat_ds.time])

landsat_ds

[11]:

<xarray.Dataset>
Dimensions:        (time: 120, y: 445, x: 461)
Coordinates:
  * time           (time) datetime64[ns] 2016-10-23T01:55:38.539139 ... 2016-...
  * y              (y) float64 -2.006e+06 -2.006e+06 ... -2.019e+06 -2.019e+06
  * x              (x) float64 4.101e+05 4.101e+05 ... 4.239e+05 4.239e+05
    spatial_ref    int32 32651
Data variables:
    nbart_red      (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_green    (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_blue     (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    nbart_nir      (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    NDWI           (time, y, x) float32 dask.array<chunksize=(1, 445, 461), meta=np.ndarray>
    tide_m         (time) float64 -2.187 -2.118 -2.094 ... 3.741 3.754 3.76
    tide_interval  (time) int64 1 1 1 1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9 9 9 9
Attributes:
    crs:           epsg:32651
    grid_mapping:  spatial_ref

We can plot the boundaries between the nine tidal intervals on the same plot we generated earlier:

[12]:

landsat_ds.sortby('time').tide_m.plot()
for i in binInterval: plt.axhline(i, c='black', alpha=0.5)
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_27_0.png

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.

[13]:

# For each interval, compute the median water index and tide height value
landsat_intervals = (landsat_ds[['tide_interval', 'NDWI', 'tide_m']]
                     .compute()
                     .groupby('tide_interval')
                     .median(dim='time'))

# Shut down Dask client now that we have processed the data we need
client.close()

# Plot the resulting set of tidal intervals
landsat_intervals.NDWI.plot(col='tide_interval', col_wrap=5, cmap='RdBu')

[13]:

<xarray.plot.facetgrid.FacetGrid at 0x7f87430a9160>
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_29_1.png

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.

[14]:

# Set up attributes to assign to each waterline
attribute_df = pd.DataFrame({'tide_m': landsat_intervals.tide_m.values})

# Extract waterlines
contours_gdf = subpixel_contours(da=landsat_intervals.NDWI,
                                 z_values=0,
                                 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))
landsat_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)


Operating in single z-value, multiple arrays mode

[14]:

<AxesSubplot:title={'center':'tide_interval = 1, spatial_ref = 32651'}, xlabel='x [metre]', ylabel='y [metre]'>
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_32_2.png

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.

[15]:

# Plot on interactive map
contours_gdf.explore(column='tide_m', cmap='YlOrRd')

[15]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Interpolate contours into a Digital Elevation Model (DEM)

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.

[16]:

# 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 Landsat dataset
intertidal_dem = interpolate_2d(ds=landsat_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)
plt.show()
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_38_0.png

You can see in the output above that our interpolation includes areas of ocean and land that are outside of the area affected by tides. To clean up the data, we can restrict the DEM to only the area between the lowest and highest observed tides:

[17]:

# Identify areas that are always wet (e.g. below low tide), or always dry
above_lowest = landsat_intervals.isel(tide_interval=0).NDWI < 0
below_highest = landsat_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)
plt.show()
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_40_0.png

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 DEA Sandbox, locate it in the file browser to the left, right click on the file, and select “Download”).

[18]:

# Export as a GeoTIFF
write_cog(geo_im=intertidal_dem_clean,
          fname='intertidal_dem.tif',
          overwrite=True)

[18]:

PosixPath('intertidal_dem.tif')

Next steps

When you are done, return to the “Set up analysis” cell, modify some values (e.g. time_range, lat_range, lon_range) and rerun the analysis.

If you’re going to change the location, you’ll need to make sure Landsat 5, 7 and 8 data is available for the new location, which you can check at the DEA Explorer (use the drop-down menu to view all Landsat products).

DEA Intertidal Elevation (NIDEM)

For more information about the science behind this notebook, please refer to the scientific article outlining the application of this approach to the entire Australian coastline: Bishop-Taylor et al. 2019 Between the tides: Modelling the elevation of Australia’s exposed intertidal zone at continental scale.

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube 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.

Last modified: December 2023

Compatible datacube version:

[19]:

print(datacube.__version__)

1.8.6

Tags

Tags: NCI compatible, sandbox compatible, landsat 5, landsat 7, landsat 8, load_ard, mostcommon_crs, calculate_indices, tidal_tag, subpixel_contours, display_map, rgb, interpolate_2d, contours_to_arrays, NDWI, tide modelling, image compositing, waterline extraction, spatial interpolation, intertidal, real world, GeoTIFF, exporting data