Combining tides with satellite data¶

This guide demonstrates how to combine tide modelling with satellite Earth observation (EO) data using the tag_tides and pixel_tides functions from eo_tides.eo.

Both these functions allow you to model the height of the tide at the exact moment of satellite image acquisition. This can then allow you to analyse satellite EO data by tidal conditions - for example, filter your data to satellite imagery collected during specific tidal stages (e.g. low or high tide).

Although both functions perform a similar function, they differ in complexity and performance. tag_tides assigns a single tide height to each timestep/satellite image, which is fast and efficient, and suitable for small-scale applications where tides are unlikely to vary across your study area. In constrast, pixel_tide models tides both through time and spatially, returning a tide height for every satellite pixel. This can be critical for producing seamless coastal EO datasets at large scale - however comes at the cost of performance.

Table 1. Comparison of tag_tides and pixel_tides

tag_tides	pixel_tides
Assigns a single tide height to each timestep/satellite image	Assigns a tide height to every individual pixel through time to capture spatial tide dynamics
🔎 Ideal for local or site-scale analysis	🌏 Ideal for regional to global-scale coastal product generation
✅ Fast, low memory use	❌ Slower, higher memory use
❌ Single tide height per image can produce artefacts in complex tidal regions	✅ Produce spatially seamless results across large extents by applying analyses at the pixel level

Getting started¶

As in the previous example, our first step is to tell eo-tides the location of our tide model or clipped tide model directory (if you haven't set this up, refer to the setup instructions here):

In [1]:

directory = "../../tests/data/tide_models/"

directory = "../../tests/data/tide_models/"

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 Landsat 8 and 9 satellite data from 2024 over the city of Broome, Western Australia - a macrotidal region with extensive intertidal coastal habitats. We will load this data from the Digital Earth Australia STAC catalogue.

Tip

For a more detailed guide to using STAC metadata and odc-stac to find and load satellite data, refer to the Digital Earth Australia STAC user guide.

In [2]:

import odc.stac
import pystac_client

# Connect to STAC catalog
catalog = pystac_client.Client.open("https://explorer.dea.ga.gov.au/stac")

# 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
bbox = [122.12, -18.25, 122.43, -17.93]
query = catalog.search(
    bbox=bbox,
    collections=["ga_ls8c_ard_3", "ga_ls9c_ard_3"],
    datetime="2024-01-01/2024-12-31",
    filter = "eo:cloud_cover < 5"  # Filter to images with <5% cloud
)

# Load data into xarray format
ds = odc.stac.load(
    items=list(query.items()),
    bands=["nbart_red", "nbart_green", "nbart_blue"],
    crs="utm",
    resolution=30,
    groupby="solar_day",
    bbox=bbox,
    fail_on_error=False,
    chunks={},
)

# Plot the first image
ds.isel(time=0).odc.explore(vmin=50, vmax=3000)

import odc.stac import pystac_client # Connect to STAC catalog catalog = pystac_client.Client.open("https://explorer.dea.ga.gov.au/stac") # 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 bbox = [122.12, -18.25, 122.43, -17.93] query = catalog.search( bbox=bbox, collections=["ga_ls8c_ard_3", "ga_ls9c_ard_3"], datetime="2024-01-01/2024-12-31", filter = "eo:cloud_cover < 5" # Filter to images with <5% cloud ) # Load data into xarray format ds = odc.stac.load( items=list(query.items()), bands=["nbart_red", "nbart_green", "nbart_blue"], crs="utm", resolution=30, groupby="solar_day", bbox=bbox, fail_on_error=False, chunks={}, ) # Plot the first image ds.isel(time=0).odc.explore(vmin=50, vmax=3000)

/env/lib/python3.10/site-packages/odc/geo/_rgba.py:56: RuntimeWarning: invalid value encountered in cast
  return x.astype("uint8")
/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(

Out[2]:

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

Using tag_tides¶

We can pass our satellite dataset ds to the tag_tides function to model a tide for each timestep in our dataset. This can help sort and filter images by tide height, allowing us to learn more about how coastal environments respond to the effect of changing tides.

The tag_tides function uses the time and date of acquisition and the geographic centroid of each satellite observation as inputs for the selected tide model (EOT20 by default). It returns an xarray.DataArray called tide_height, with a modelled tide for every timestep in our satellite dataset:

In [3]:

from eo_tides.eo import tag_tides

tides_da = tag_tides(
    data=ds,
    directory=directory,
)

# Print modelled tides
print(tides_da)

from eo_tides.eo import tag_tides tides_da = tag_tides( data=ds, directory=directory, ) # Print modelled tides print(tides_da)

Setting tide modelling location from dataset centroid: 122.27, -18.09
Modelling tides with EOT20
<xarray.DataArray 'tide_height' (time: 45)> Size: 180B
array([-0.25108945,  0.5275667 ,  1.5171705 ,  2.011433  ,  0.36809078,
        1.9124153 ,  1.0965405 ,  1.1305977 , -0.21305004,  1.7526643 ,
       -0.2175682 , -0.80696225, -0.9524483 ,  3.0443938 , -1.3169092 ,
        3.3964403 ,  2.625125  , -0.9654651 , -0.4368508 ,  2.581808  ,
        1.7244624 , -1.0423656 ,  0.16597499,  0.73508024, -0.33408186,
       -0.10130765,  0.7978594 , -1.8499157 ,  1.6090035 , -1.2717861 ,
       -1.8342571 ,  2.162794  ,  2.7683735 , -2.6152036 , -2.39916   ,
        2.5937896 ,  2.1230242 , -1.6377252 ,  3.2850509 , -0.3772273 ,
        0.6212255 ,  1.6580964 ,  0.71566176, -0.03352478,  1.0641807 ],
      dtype=float32)
Coordinates:
  * time        (time) datetime64[ns] 360B 2024-01-07T01:55:31.679580 ... 202...
    tide_model  <U5 20B 'EOT20'

We can easily combine these modelled tides with our original satellite data for further analysis. The code below adds our modelled tides as a new tide_height variable under Data variables.

In [4]:

ds["tide_height"] = tides_da
print(ds)

ds["tide_height"] = tides_da print(ds)

<xarray.Dataset> Size: 703MB
Dimensions:      (y: 1185, x: 1099, time: 45)
Coordinates:
  * y            (y) float64 9kB 8.017e+06 8.017e+06 ... 7.982e+06 7.982e+06
  * x            (x) float64 9kB 4.068e+05 4.068e+05 ... 4.397e+05 4.398e+05
    spatial_ref  int32 4B 32751
  * time         (time) datetime64[ns] 360B 2024-01-07T01:55:31.679580 ... 20...
    tide_model   <U5 20B 'EOT20'
Data variables:
    nbart_red    (time, y, x) float32 234MB dask.array<chunksize=(1, 1185, 1099), meta=np.ndarray>
    nbart_green  (time, y, x) float32 234MB dask.array<chunksize=(1, 1185, 1099), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 234MB dask.array<chunksize=(1, 1185, 1099), meta=np.ndarray>
    tide_height  (time) float32 180B -0.2511 0.5276 1.517 ... -0.03352 1.064

Tip

You could also model tides and insert tide heights into ds in a single step via:
ds["tide_height"] = tag_tides(ds, ...)

We can plot this new tide_height variable over time to inspect the tide heights observed by the satellites in our time series:

In [5]:

ds.tide_height.plot()

ds.tide_height.plot()

Out[5]:

[<matplotlib.lines.Line2D at 0x7f6056bfc850>]

Selecting and analysing satellite data by tide¶

Having tide_height as a variable allows us to select and analyse our satellite data using information about tides. For example, we could sort by tide_height, then plot the lowest and highest tide images in our time series:

In [6]:

# Sort by tide and plot the first and last image
ds_sorted = ds.sortby("tide_height")
ds_sorted.isel(time=[0, -1]).odc.to_rgba(vmin=50, vmax=3000).plot.imshow(col="time")

# Sort by tide and plot the first and last image ds_sorted = ds.sortby("tide_height") ds_sorted.isel(time=[0, -1]).odc.to_rgba(vmin=50, vmax=3000).plot.imshow(col="time")

Out[6]:

<xarray.plot.facetgrid.FacetGrid at 0x7f6056ae61a0>

Using pixel_tides¶

The previous examples show how to model a single tide height for each satellite image using the centroid of the image as a tide modelling location. However, in reality tides vary spatially – potentially by several metres in areas of complex tidal dynamics. This means that an individual satellite image can capture a range of tide conditions.

We can use the pixel_tides function to capture this spatial variability in tide heights. For efficient processing, this function first models tides into a low resolution grid surrounding each satellite image in our time series. This lower resolution data includes a buffer around the extent of our satellite data so that tides can be modelled seamlessly across analysis boundaries.

First, let's reload our satellite data for a fresh start:

In [7]:

# Load data into xarray format
ds = odc.stac.load(
    items=list(query.items()),
    bands=["nbart_red", "nbart_green", "nbart_blue"],
    crs="utm",
    resolution=30,
    groupby="solar_day",
    bbox=bbox,
    fail_on_error=False,
    chunks={},
)

# Load data into xarray format ds = odc.stac.load( items=list(query.items()), bands=["nbart_red", "nbart_green", "nbart_blue"], crs="utm", resolution=30, groupby="solar_day", bbox=bbox, fail_on_error=False, chunks={}, )

Now run pixel_tides, passing our satellite dataset ds as an input:

In [8]:

from eo_tides.eo import pixel_tides

# Model tides spatially
tides_lowres = pixel_tides(
    data=ds,
    resample=False,
    directory=directory,
)

# Print output
print(tides_lowres)

from eo_tides.eo import pixel_tides # Model tides spatially tides_lowres = pixel_tides( data=ds, resample=False, directory=directory, ) # Print output print(tides_lowres)

Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides with EOT20
Returning low resolution tide array
<xarray.DataArray 'tide_height' (time: 45, y: 13, x: 13)> Size: 30kB
array([[[-0.30606133, -0.30717686, -0.30773783, ..., -0.25615942,
         -0.25615942, -0.25615942],
        [-0.30401564, -0.30424863, -0.30387282, ..., -0.25615942,
         -0.25615942, -0.25615942],
        [-0.30068266, -0.29933682, -0.29704148, ..., -0.25615942,
         -0.25615942, -0.25615942],
        ...,
        [-0.27365166, -0.27505815, -0.2606054 , ..., -0.25108945,
         -0.25108945, -0.25108945],
        [-0.27030516, -0.26545134, -0.26545134, ..., -0.25108945,
         -0.25108945, -0.25108945],
        [-0.26545134, -0.26545134, -0.26545134, ..., -0.25108945,
         -0.25108945, -0.25108945]],

       [[ 0.4048146 ,  0.40729165,  0.4100286 , ...,  0.51402086,
          0.51402086,  0.51402086],
        [ 0.41507572,  0.41939604,  0.42421374, ...,  0.51402086,
          0.51402086,  0.51402086],
        [ 0.42650127,  0.43261963,  0.43943623, ...,  0.51402086,
          0.51402086,  0.51402086],
...
        [-0.05945582, -0.06040129, -0.04247594, ..., -0.03352478,
         -0.03352478, -0.03352478],
        [-0.05551685, -0.04903932, -0.04903932, ..., -0.03352478,
         -0.03352478, -0.03352478],
        [-0.04903932, -0.04903932, -0.04903932, ..., -0.03352478,
         -0.03352478, -0.03352478]],

       [[ 0.88739955,  0.90003395,  0.91273254, ...,  1.0558591 ,
          1.0558591 ,  1.0558591 ],
        [ 0.9003904 ,  0.91464955,  0.9293343 , ...,  1.0558591 ,
          1.0558591 ,  1.0558591 ],
        [ 0.9152668 ,  0.93159544,  0.9488063 , ...,  1.0558591 ,
          1.0558591 ,  1.0558591 ],
        ...,
        [ 1.0099615 ,  1.0095476 ,  1.0438287 , ...,  1.0641807 ,
          1.0641807 ,  1.0641807 ],
        [ 1.0166272 ,  1.026946  ,  1.026946  , ...,  1.0641807 ,
          1.0641807 ,  1.0641807 ],
        [ 1.026946  ,  1.026946  ,  1.026946  , ...,  1.0641807 ,
          1.0641807 ,  1.0641807 ]]], dtype=float32)
Coordinates:
  * time         (time) datetime64[ns] 360B 2024-01-07T01:55:31.679580 ... 20...
  * x            (x) float64 104B 3.925e+05 3.975e+05 ... 4.475e+05 4.525e+05
  * y            (y) float64 104B 8.028e+06 8.022e+06 ... 7.972e+06 7.968e+06
    tide_model   <U5 20B 'EOT20'
    spatial_ref  int32 4B 32751

If we plot the resulting data, we can see that we now have two-dimensional tide surfaces for each timestep in our data (instead of the single tide height per timestamp returned by the tag_tides function).

Blue values below indicate low tide pixels, while red indicates high tide pixels. If you look closely, you may see some spatial variability in tide heights within each timestep, with slight variations in tide heights along the north-west side of the study area:

In [9]:

# Plot the first four timesteps in our data
tides_lowres.isel(time=slice(0, 4)).plot.imshow(col="time", vmin=-1, vmax=1, cmap="RdBu")

# Plot the first four timesteps in our data tides_lowres.isel(time=slice(0, 4)).plot.imshow(col="time", vmin=-1, vmax=1, cmap="RdBu")

Out[9]:

<xarray.plot.facetgrid.FacetGrid at 0x7f60811c18a0>

Reprojecting into original high-resolution spatial grid¶

By setting resample=True, we can use interpolation to re-project our low resolution tide data back into the resolution of our satellite image, resulting in an individual tide height for every pixel in our dataset through time and space:

In [10]:

# Model tides spatially
tides_highres = pixel_tides(
    data=ds,
    resample=True,
    directory=directory,
)

# Plot the first four timesteps in our data
tides_highres.isel(time=slice(0, 4)).plot.imshow(col="time", vmin=-1, vmax=1, cmap="RdBu")

# Model tides spatially tides_highres = pixel_tides( data=ds, resample=True, directory=directory, ) # Plot the first four timesteps in our data tides_highres.isel(time=slice(0, 4)).plot.imshow(col="time", vmin=-1, vmax=1, cmap="RdBu")

Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides with EOT20
Reprojecting tides into original resolution

Out[10]:

<xarray.plot.facetgrid.FacetGrid at 0x7f6047d99fc0>

tides_highres will have exactly the same dimensions as ds, with a unique tide height for every satellite pixel:

In [11]:

ds.sizes

ds.sizes

Out[11]:

Frozen({'y': 1185, 'x': 1099, 'time': 45})

In [12]:

tides_highres.sizes

tides_highres.sizes

Out[12]:

Frozen({'time': 45, 'y': 1185, 'x': 1099})

Because of this, our stack of tides can be added as an additional 3D variable in our dataset:

In [13]:

ds["tide_height_pixel"] = tides_highres
print(ds)

ds["tide_height_pixel"] = tides_highres print(ds)

<xarray.Dataset> Size: 938MB
Dimensions:            (y: 1185, x: 1099, time: 45)
Coordinates:
  * y                  (y) float64 9kB 8.017e+06 8.017e+06 ... 7.982e+06
  * x                  (x) float64 9kB 4.068e+05 4.068e+05 ... 4.398e+05
    spatial_ref        int32 4B 32751
  * time               (time) datetime64[ns] 360B 2024-01-07T01:55:31.679580 ...
    tide_model         <U5 20B 'EOT20'
Data variables:
    nbart_red          (time, y, x) float32 234MB dask.array<chunksize=(1, 1185, 1099), meta=np.ndarray>
    nbart_green        (time, y, x) float32 234MB dask.array<chunksize=(1, 1185, 1099), meta=np.ndarray>
    nbart_blue         (time, y, x) float32 234MB dask.array<chunksize=(1, 1185, 1099), meta=np.ndarray>
    tide_height_pixel  (time, y, x) float32 234MB -0.3038 -0.3039 ... 1.064

Calculating tide height min/max/median/quantiles for each pixel¶

Min, max or any specific quantile of all tide heights observed over a region can be calculated for each pixel by passing in a list of quantiles/percentiles using the calculate_quantiles parameter.

This calculation is performed on the low resolution modelled tide data before reprojecting to higher resolution, so should be faster than calculating min/max/median tide at high resolution:

In [14]:

# Model tides spatially
tides_highres_quantiles = pixel_tides(
    data=ds,
    calculate_quantiles=(0, 0.5, 1),
    directory=directory,
)

# Plot quantiles
tides_highres_quantiles.plot.imshow(col="quantile")

# Model tides spatially tides_highres_quantiles = pixel_tides( data=ds, calculate_quantiles=(0, 0.5, 1), directory=directory, ) # Plot quantiles tides_highres_quantiles.plot.imshow(col="quantile")

Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides with EOT20
Computing tide quantiles
Reprojecting tides into original resolution

Out[14]:

<xarray.plot.facetgrid.FacetGrid at 0x7f6047cf2350>

Modelling custom times¶

Instead of using times contained in the time dimension of our dataset, we can also calculate pixel-based tides for a custom set of times:

In [15]:

import pandas as pd

custom_times = pd.date_range(
    start="2022-01-01", 
    end="2022-01-02", 
    freq="6H",
)

# Model tides spatially
tides_highres = pixel_tides(
    data=ds, 
    time=custom_times,
    directory=directory,
)

# Plot custom timesteps
tides_highres.plot.imshow(col="time")

import pandas as pd custom_times = pd.date_range( start="2022-01-01", end="2022-01-02", freq="6H", ) # Model tides spatially tides_highres = pixel_tides( data=ds, time=custom_times, directory=directory, ) # Plot custom timesteps tides_highres.plot.imshow(col="time")

Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides with EOT20
Reprojecting tides into original resolution

Out[15]:

<xarray.plot.facetgrid.FacetGrid at 0x7f6047be1000>

Next steps¶

Now that we have learnt to combine tide modelling with satellite data, we can learn how to calculate statistics describing local tide dynamics, as well as biases caused by interactions between tidal processes and satellite orbits.