Modelling intertidal elevation using tidal data 7b904a90cfb94561afcab45d0e9ea078

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

Description

In this example, we demonstrate an intertidal elevation extraction method that combines data from the Landsat 8 and 9 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

Note: The intertidal elevation method presented here has recently been superseded by an improved method published in the DEA Intertidal product. Visit the DEA Intertidal product documentation for detailed technical information about this new product, or view the updated code in the DEA Intertidal Github repository.

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

import datacube
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

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, extract_vertices, xr_interpolate
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)

Client

Client-16f41e64-39c3-11ef-8d01-229fa11ce3cd

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /user/robbi.bishoptaylor@ga.gov.au/proxy/8787/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. For reasonable load times, keep this to a range of ~0.1 degrees or less.

  • lon_range: The longitude range to analyse. For reasonable load times, keep this to a range of ~0.1 degrees or less.

  • time_range: The date range to analyse.

By default, this example will generate an intertidal elevation model for tidal flats near Roebuck Bay near Broome, Western Australia.

If running the notebook for the first time, keep the default settings below. This will demonstrate how the analysis works and provide meaningful results.

[3]:

lat_range = (-18.15, -18.25)
lon_range = (122.16, 122.27)
time_range = ("2020", "2023")

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 8 and 9 satellites for the area and time included in lat_range, lon_range and time_range.

We set cloud_cover=(0, 5) to load only satellite scenes that contain less than 5% 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 Landsat
landsat_ds = load_ard(
    dc=dc,
    products=["ga_ls8c_ard_3", "ga_ls9c_ard_3"],
    output_crs=output_crs,
    align=(15, 15),
    group_by="solar_day",
    dask_chunks={"x": 2048, "y": 2048},
    cloud_cover=(0, 5),
    **query
)

Finding datasets
    ga_ls8c_ard_3
    ga_ls9c_ard_3
Applying fmask pixel quality/cloud mask
Returning 167 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> Size: 387MB
Dimensions:      (time: 167, y: 371, x: 390)
Coordinates:
  * time         (time) datetime64[ns] 1kB 2020-01-29T01:49:37.148233 ... 202...
  * y            (y) float64 3kB -2.007e+06 -2.007e+06 ... -2.018e+06 -2.018e+06
  * x            (x) float64 3kB 4.112e+05 4.112e+05 ... 4.228e+05 4.228e+05
    spatial_ref  int32 4B 32651
Data variables:
    nbart_red    (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_green  (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_nir    (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), 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 (e.g. timestep=10) and re-run the cell to plot a different timestep.

[7]:

# Set the timesteps to visualise
timestep = 0

# 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 0x7fe35f9c1360>
../../../_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.21, -18.20
Modelling tides using FES2014

[9]:

<xarray.Dataset> Size: 483MB
Dimensions:      (time: 167, y: 371, x: 390)
Coordinates:
  * time         (time) datetime64[ns] 1kB 2020-01-29T01:49:37.148233 ... 202...
  * y            (y) float64 3kB -2.007e+06 -2.007e+06 ... -2.018e+06 -2.018e+06
  * x            (x) float64 3kB 4.112e+05 4.112e+05 ... 4.228e+05 4.228e+05
    spatial_ref  int32 4B 32651
Data variables:
    nbart_red    (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_green  (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_nir    (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    NDWI         (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    tide_m       (time) float32 668B -0.9305 1.639 2.303 ... -0.8275 0.267 0.175
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.tide_m.sortby("time").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 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 discrete tidal intervals, ranging from the lowest to the highest tides observed by Landsat. For more information on this method, refer to Sagar et al. 2018.

[11]:

# Set the number of tide intervals. More intervals can produce more
# detail, but add noise
tide_intervals = 5

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

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

landsat_ds

[11]:

<xarray.Dataset> Size: 483MB
Dimensions:        (time: 167, y: 371, x: 390)
Coordinates:
  * time           (time) datetime64[ns] 1kB 2020-01-29T01:49:37.148233 ... 2...
  * y              (y) float64 3kB -2.007e+06 -2.007e+06 ... -2.018e+06
  * x              (x) float64 3kB 4.112e+05 4.112e+05 ... 4.228e+05 4.228e+05
    spatial_ref    int32 4B 32651
Data variables:
    nbart_red      (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_green    (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_blue     (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    nbart_nir      (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    NDWI           (time, y, x) float32 97MB dask.array<chunksize=(1, 371, 390), meta=np.ndarray>
    tide_m         (time) float32 668B -0.9305 1.639 2.303 ... 0.267 0.175
    tide_interval  (time) category 207B <class 'xarray.core.extension_array.P...
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.tide_m.sortby("time").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"]]
    .groupby("tide_interval")
    .median(dim="time")
    .compute()
)

# 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")

/env/lib/python3.10/site-packages/dask/utils.py:78: RuntimeWarning: All-NaN slice encountered
  return func(*args, **kwargs)

[13]:

<xarray.plot.facetgrid.FacetGrid at 0x7fe35c39f8b0>
../../../_images/notebooks_Real_world_examples_Intertidal_elevation_29_2.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

/home/jovyan/Robbi/dea-notebooks/Real_world_examples/../Tools/dea_tools/spatial.py:536: UserWarning: The `squeeze` kwarg to GroupBy is being removed.Pass .groupby(..., squeeze=False) to disable squeezing, which is the new default, and to silence this warning.
  contour_arrays = {
../../../_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 first extract vertices from the contours above to get a dataset of points. We then use the elevation of these points to interpolate smooth, continuous elevations across the intertidal zone using linear interpolation.

[16]:

# Convert contour lines into point vertices for interpolation
gdf_points = extract_vertices(contours_gdf)

# Interpolate points over the spatial extent of the Landsat dataset
intertidal_dem = xr_interpolate(
    ds=landsat_intervals,
    gdf=gdf_points,
).tide_m

# 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 on an interactive map
intertidal_dem_clean.odc.explore()

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

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

intertidal_dem_clean.odc.write_cog(fname='intertidal_dem.tif', overwrite=True);

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

A more accurate and detailed intertidal elevation modelling method has recently been published in the DEA Intertidal product. Visit the DEA Intertidal product documentation for detailed technical information about this new product, or view the updated code in the DEA Intertidal Github repository.

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 Discord chat 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: July 2024

Compatible datacube version:

[19]:

print(datacube.__version__)

1.8.19

Tags

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