Mapping the intertidal zone¶
The intertidal zone (i.e. the land along the coast that is periodically inundated by the tide) 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 mapping data describing the spatial extents of the intertidal zone are essential for managing these environments, and predicting when and where these threats will have the greatest impact. However, the intertidal zone is challenging and expensive to map at large scale using intensive manual survey methods - particularly across large coastal regions.
Satellite Earth observation (EO) data is freely available for the entire planet, making satellite imagery a powerful and cost-effective tool for mapping the intertidal zone at regional, national scale or global scale.
This case study will demonstrate a simple intertidal mapping workflow that combines free and open Landsat satellite data with tide modelling from eo-tides.
The workflow includes:
- Loading a time-series of cloud-free satellite data from the cloud using
odc-stac - Converting our satellite data to a remote sensing water index (NDWI)
- Modelling tides for each satellite image and inspecting how these observed tides match up to the full local astronomical tide range
- Filtering our satellite imagery to low and high tide observations
- Combining noisy individual images into clean low and high tide median NDWI composites
- Using these composites to extract the extent of the intertidal zone
More information
For more information about the workflows described below, refer to Sagar et al. 2017, Sagar et al. 2018, and Bishop-Taylor et al. 2019.
Getting started¶
Import any Python packages we need for the analysis:
import odc.stac
import pystac_client
import planetary_computer
import matplotlib.pyplot as plt
from eo_tides.eo import tag_tides
from eo_tides.stats import tide_stats
from eo_tides.utils import list_models
Tide model directory¶
We need to tell eo-tides the location of our tide model directory (if you haven't set this up, refer to the setup instructions here):
directory = "../../tests/data/tide_models"
# Confirm we have model data
list_models(directory=directory, show_supported=False);
──────────────────────────────────────────────────────────────────────────────── 🌊 | Model | Expected path ──────────────────────────────────────────────────────────────────────────────── ✅ │ EOT20 │ ../../tests/data/tide_models/EOT20/ocean_tides ✅ │ GOT5.5 │ ../../tests/data/tide_models/GOT5.5/ocean_tides ✅ │ HAMTIDE11 │ ../../tests/data/tide_models/hamtide ──────────────────────────────────────────────────────────────────────────────── Summary: Available models: 3/50
Analysis parameters¶
To make our analysis more re-usable, we can define some important parameters up-front. The default will load Landsat 8 and 9 satellite data from 2022-23 over the city of Broome, Western Australia - a macrotidal region with extensive intertidal coastal habitats.
Tip
Leave the defaults below unchanged the first time you run through this notebook.
# Set the study area (xmin, ymin, xmax, ymax)
bbox = [122.12, -18.25, 122.43, -17.93]
# Set the time period
start_date = "2022-01-01"
end_date = "2023-12-31"
# Satellite products and bands to load
satellite_sensors = ["landsat-c2-l2"]
bands = ["green", "nir08"]
# Tide model to use
tide_model = "EOT20"
Load satellite data using odc-stac¶
Now we can load a time-series of satellite data over our area of interest using the Open Data Cube's odc-stac package.
This powerful package allows us to load open satellite data (e.g ESA Sentinel-2 or NASA/USGS Landsat) for any time period and location on the planet, and load our data into a multi-dimensional xarray.Dataset format dataset.
In this example, we will load our data from the Microsoft Planetary Computer STAC catalogue, and restrict our data to cloud-free images with less than 10% cloud (eo:cloud_cover": {"lt": 10}"), and load our data at low resolution (resolution=50) to improve load times.
Tip
For a more detailed guide to using STAC metadata and odc-stac to find and load satellite data, refer to this STAC user guide.
# Connect to STAC catalog
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
# Set cloud access defaults
odc.stac.configure_rio(
cloud_defaults=True,
aws={"aws_unsigned": True},
)
# Build a query and search the STAC catalog for all matching items
query = catalog.search(
bbox=bbox,
collections=satellite_sensors,
datetime=f"{start_date}/{end_date}",
query={
"eo:cloud_cover": {"lt": 10}, # Filter to images with <5% cloud
"platform": {"in": ["landsat-8", "landsat-9"]}, # No Landsat 7
},
)
# Load data into xarray format
ds = odc.stac.load(
items=query.item_collection(),
bands=bands,
crs="utm",
resolution=50,
groupby="solar_day",
bbox=bbox,
fail_on_error=False,
chunks={},
)
# Apply USGS Landsat Collection 2 scaling factors to convert
# surface reflectance to between 0.0 and 1.0. See:
# https://www.usgs.gov/faqs/how-do-i-use-a-scale-factor-landsat-level-2-science-products
ds = (ds.where(ds != 0) * 0.0000275 + -0.2).clip(0, 1)
print(ds)
<xarray.Dataset> Size: 454MB
Dimensions: (time: 121, y: 712, x: 659)
Coordinates:
* y (y) float64 6kB 8.017e+06 8.017e+06 ... 7.982e+06 7.982e+06
* x (x) float64 5kB 4.068e+05 4.069e+05 ... 4.397e+05 4.397e+05
spatial_ref int32 4B 32751
* time (time) datetime64[ns] 968B 2022-01-01T01:55:34.757654 ... 20...
Data variables:
green (time, y, x) float32 227MB dask.array<chunksize=(1, 712, 659), meta=np.ndarray>
nir08 (time, y, x) float32 227MB dask.array<chunksize=(1, 712, 659), meta=np.ndarray>
Converting satellite data to NDWI water index¶
To help us map the distribution of tidal inundation across our study area, we can convert our spectral bands to a remote sensing water index that enhances contrast between water and dry land. We will use the Normalised Difference Water Index (NDWI), which is particularly useful for intertidal mapping as it is less vulnerable to misclassifications caused by wet sand and mud after high tide:
$$ \text{NDWI} = \frac{\text{Green} - \text{NIR}}{\text{Green} + \text{NIR}} $$
In the image below, red pixels represent dry land pixels, while blue pixels represent water:
# Calculate NDWI
ds[["ndwi"]] = (ds.green - ds.nir08) / (ds.green + ds.nir08)
# Plot a single timestep
ds.ndwi.isel(time=1).odc.explore(vmin=-0.5, vmax=0.5, cmap="RdBu")
/env/lib/python3.10/site-packages/dask/core.py:133: RuntimeWarning: invalid value encountered in divide return func(*(_execute_task(a, cache) for a in args)) /env/lib/python3.10/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. dest = _reproject(
Modelling tide heights for each satellite image¶
We can now use the tag_tides function to model tides for each image in our satellite timeseries:
ds["tide_height"] = tag_tides(
data=ds,
model=tide_model,
directory=directory,
)
Setting tide modelling location from dataset centroid: 122.27, -18.09 Modelling tides using EOT20
We can plot a histogram of the tide heights of our satellite images. This shows that we have decent coverage of the tide range between approximately -2.8 to 3.5 m Above Mean Sea Level (AMSL):
ds["tide_height"].plot.hist(bins=10, figsize=(5, 3))
plt.title("Histogram of tide heights in satellite timeseries");
Exploring tide biases¶
It is also important to compare our satellite tide observations against the full range of astronomical tides at our location. This lets us become aware of any major biases in our data. For example, here we can see that our data is biased away from low tide observations. This means that our final intertidal maps will underestimate the lower extent of the intertidal zone.
tide_stats(
data=ds,
model=tide_model,
directory=directory,
);
Using tide modelling location: 122.27, -18.09 Modelling tides using EOT20 Modelling tides using EOT20 🌊 Modelled astronomical tide range: 9.46 metres. 🛰️ Observed tide range: 6.37 metres. 🔴 67% of the modelled astronomical tide range was observed at this location. 🟡 The highest 12% (1.09 metres) of the tide range was never observed. 🔴 The lowest 21% (2.00 metres) of the tide range was never observed. 🌊 Mean modelled astronomical tide height: -0.00 metres. 🛰️ Mean observed tide height: 0.65 metres. ⬆️ The mean observed tide height was 0.65 metres higher than the mean modelled astronomical tide height.
Extracting low and high tide images¶
Now that we have satellite images and associated tide heights, we can filter our data to extract only images from specific tidal stages.
In this example, we will focus on extracting low tide and high tide images. We can identify low and high tide images by calculating quantiles of tide heights across time, selecting images from the lowest and highest 20% of tide heights.
Tip
Try modifying the lowtide_cutoff and hightide_cutoff tide cutoffs below to customise the tidal stages we will analyse.
# Calculate low and high tide height thresholds
lowtide_cutoff = 0.2
hightide_cutoff = 0.8
lowtide_thresh, hightide_thresh = ds.tide_height.quantile([lowtide_cutoff, hightide_cutoff])
print(f"Low tide threshold: {lowtide_thresh:.2f} metres AMSL")
print(f"High tide threshold: {hightide_thresh:.2f} metres AMSL")
Low tide threshold: -1.00 metres AMSL High tide threshold: 2.51 metres AMSL
We can now use these thresholds to select just the subset of our data with tide heights lower and higher than our thresholds:
# Extract subset of low and high tide images
ds_lowtide = ds.sel(time=ds.tide_height <= lowtide_thresh)
ds_hightide = ds.sel(time=ds.tide_height >= hightide_thresh)
# Plot extracted images over all images
ds.tide_height.plot(marker="o", linewidth=0, label="Other satellite images")
ds_hightide.tide_height.plot(marker="o", linewidth=0, label="High tide images")
ds_lowtide.tide_height.plot(marker="o", linewidth=0, label="Low tide images")
plt.axhline(lowtide_thresh, color="black", linestyle="dashed")
plt.axhline(hightide_thresh, color="black", linestyle="dashed")
plt.legend(bbox_to_anchor=(1.01, 0.6))
plt.title("Low and high tide satellite images");
Creating low and high tide NDWI composites¶
Individual satellite images can be affected by many sources of noise, including clouds, saltspray, whitewater and sunglint. These sources of noise can make it difficult to map coastal environments consistently across time and along the coastline.
To produce consistent analysis outputs, it can be useful to combine multiple satellite images into a single, cleaner "composite" image. This can improve the signal-to-noise ratio of our data, and allow us to map coastal features and change more accurately (at the cost of lower temporal resolution).
In this example, we will combine our low and high tide NDWI data into two clean, relatively noise-free low tide and high tide median composites. We use a median calculation for this analysis as medians are extremely robust to noisy data.
Important
Note the use of .load() below. Up to this point, our entire analysis has been "lazy", which means we haven't loaded the majority of our satellite data - we have simply "queued" up our analysis to run in a single step. This makes it quick and easy to write code without having to wait for every step of our workflow to run every time, or ever worrying about running out of memory. Running .load() triggers our entire analysis to run, and then load our final outputs into memory for further use. For more information about lazy loading and processing, see Parallel processing with Dask.
Important
Be patient; this step may take several minutes to complete.
# Combine NDWI into single median composites
ndwi_lowtide_median = ds_lowtide.ndwi.median(dim="time")
ndwi_hightide_median = ds_hightide.ndwi.median(dim="time")
# Process our data and load into memory
ndwi_lowtide_median.load()
ndwi_hightide_median.load();
Now that we have processed our median NDWI composites, we can plot them on a map:
ndwi_lowtide_median.odc.explore(vmin=-0.5, vmax=0.5, cmap="RdBu")
ndwi_hightide_median.odc.explore(vmin=-0.5, vmax=0.5, cmap="RdBu")
Or plot them together for a direct comparison. We can see that the low tide NDWI composite contains additional exposed intertidal terrain along the coastline that becomes inundated at high tide:
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ndwi_lowtide_median.plot.imshow(ax=axes[0], vmin=-0.5, vmax=0.5, cmap="RdBu")
ndwi_hightide_median.plot.imshow(ax=axes[1], vmin=-0.5, vmax=0.5, cmap="RdBu")
axes[0].set_title("Low tide NDWI composite")
axes[1].set_title("High tide NDWI composite");
Mapping intertidal extent¶
We can now perform a simple comparison of our two outputs to identify intertidal pixels as any pixels that are dry in our low tide composite, but wet in our high tide composite.
First, let's classify our NDWI outputs into binary dry (True or 1.0) vs. wet (False or 0.0) using a standard 0.0 NDWI threshold. In the images below, yellow pixels represent water, and purple pixels represent land:
# Classify our NDWI layers
ndwi_lowtide_bool = ndwi_lowtide_median > 0.0
ndwi_hightide_bool = ndwi_hightide_median > 0.0
# Plot outputs
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ndwi_lowtide_bool.plot.imshow(ax=axes[0])
ndwi_hightide_bool.plot.imshow(ax=axes[1])
axes[0].set_title("Low tide classified NDWI")
axes[1].set_title("High tide classified NDWI");
Now, we can identify pixels that are True (dry) in our low tide data, and False (wet) in our high tide data.
Yellow pixels in our resulting map represent likely intertidal pixels that are exposed at low tide and inundated at high tide!
intertidal = ~ndwi_lowtide_bool & ndwi_hightide_bool
intertidal.odc.explore()
Exporting data¶
As a final step, we can export our classified intertidal map as a GeoTIFF raster that we can load into a GIS like QGIS or ArcGIS Pro:
intertidal.astype("int16").odc.write_cog("intertidal_map.tif", overwrite=True);
Next steps¶
Now that you have completed running this example, here's some possible next steps:
- Download the exported
intertidal_map.tifand load it into a GIS software (QGIS, ArcGIS Pro) to inspect the output classification. - Return to Analysis parameters, and re-run the analysis for a different location (
bbox) or time period (start_date,end_date). - Try modifying Analysis parameters to select a different tide model (e.g.
tide_model="GOT5.5"ortide_model="HAMTIDE11"). - Return to Extracting low and high tide images, and change the percentage thresholds (
lowtide_cutoff,hightide_cutoff) used to identify low and high tide images. - Advanced: In this simple analysis, we filtered to cloud-free images by discarding entire images with more than 10% cloud. Rather than filtering cloud-free images, consider masking clouds at the pixel-level using cloud masking bands that are packaged with the data (e.g.
qa_pixel). - Advanced: NDWI is just one of many possible remote sensing water indices that can be used for coastal analysis. Update the workflow to use a different water index (e.g. MNDWI), ensuring that you load any new bands required for the index calculation.
In addition, consider the following questions:
- What implications do the tide biases calculated in Exploring tide biases have for the outputs of this analysis?
- Are our outputs likely to fully map the entire intertidal zone? What areas of the intertidal zone are likely to be better or poorly mapped?
- If you experimented with running the analysis using a different tide model, how does this influence our results?