Combining satellite data with tidal modelling 875a79cf7f304968a3b324696b657f29

Background

Ocean tides are the periodic rise and fall of the ocean caused by the gravitational pull of the moon and sun and the earth’s rotation. Tides in coastal areas can greatly influence how these environments appear in satellite imagery as water levels vary by up to 12 metres (e.g. in north-western Australia). To be able to study coastal environments and processes along Australia’s coastline, it is vital to obtain data on tidal conditions at the exact moment each satellite image was acquired.

Description

This notebook demonstrates how to combine remotely sensed imagery with information about ocean tides using functions from the `dea_tools.coastal module <../Tools/dea_tools/coastal.py>`__, allowing us to analyse satellite imagery by tidal stage (e.g. low, high, ebb, flow). These functions use the Finite Element Solution 2014 (FES2014) tidal model to calculate the height (relative to mean sea level, i.e. approximately equivalent to the Australian Height Datum or AHD) and stage of the tide at the exact moment each satellite image was acquired.

These tide modelling tools and tide models underpin DEA products including DEA Coastlines, DEA Intertidal Extents (ITEM), DEA High Low Tide Imagery (HLTC), and DEA Intertidal Elevation (NIDEM).

The notebook demonstrates how to:

  1. Model tide heights for specific coordinates and times using the model_tides function

  2. Model tide heights for each satellite observation using the tidal_tag function

  3. Use tide height data to produce median composites of the coast at low and high tide

  4. Swap a dataset’s dimensions to make it easier to select imagery from low to high tide.

  5. Compute ebb or flow tide phase data to determine whether water levels were rising or falling in each satellite observation

Advanced tools:

  1. Spatially model tides into each pixel of a satellite dataset using pixel_tides

  2. Evaluate potential biases in the tides observed by a satellite using tidal_stats


Getting started

To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.

Load packages

[1]:
import datacube
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt

import sys

sys.path.insert(1, "../Tools/")
from dea_tools.datahandling import load_ard
from dea_tools.plotting import rgb, display_map
from dea_tools.coastal import model_tides, tidal_tag, pixel_tides, tidal_stats

Connect to the datacube

[2]:
dc = datacube.Datacube(app='Tidal_modelling')

Modelling tide heights using model_tides

To simply model tide heights for a specific location and set of times, we can use the dea_tools.coastal.model_tides function. For example, we can model hourly tides across a one month period (September 2022):

[3]:
# Set of times to model tides for
date_list = pd.date_range(start="2022-09-01", end="2022-09-30", freq="1H")

# Run the FES2014 tidal model
tide_df = model_tides(
    x=[122.21],
    y=[-18.20],
    time=date_list,
)

# Print outputs
tide_df.head()

Modelling tides using FES2014
[3]:
tide_model tide_m
time x y
2022-09-01 00:00:00 122.21 -18.2 FES2014 -3.606788
2022-09-01 01:00:00 122.21 -18.2 FES2014 -2.470168
2022-09-01 02:00:00 122.21 -18.2 FES2014 -0.675034
2022-09-01 03:00:00 122.21 -18.2 FES2014 1.140802
2022-09-01 04:00:00 122.21 -18.2 FES2014 2.591887

Tide heights for each time and coordinate are included in the tide_m column above (representing tide height in metres relative to Mean Sea Level).

We can also plot out resulting tides to view how tides changed across this month. By looking at the y-axis, we can see that tides ranged from a minimum of ~-4 metres up to a maximum of +4 metres relative to Mean Sea Level:

[4]:
tide_df.reset_index(["x", "y"]).tide_m.plot();
../../../_images/notebooks_How_to_guides_Tidal_modelling_10_0.png

“One-to-many” and “one-to-one” modes

By default, the model_tides function operates in “one-to-many” mode, which will model tides for every requested timestep, at every requested location. For example, if we provided five x, y coordinates and five timesteps, the function would return:

5 locations * 5 timesteps = 25 modelled tides

However, often you may have a list of locations and matching timesteps. Using “one-to-one” mode, we can model tides for only these exact pairs of locations and times:

5 timesteps at 5 locations = 5 modelled tides

To demonstrate “one-to-one” mode, imagine we have a pandas.Dataframe where each row contains unique site locations and times:

[5]:
sites_df = pd.DataFrame({
    "time": pd.date_range(start="2022-09-01", end="2022-09-30", periods=5),
    "x": [122.21, 122.22, 122.23, 122.24, 122.25],
    "y": [-18.20, -18.20, -18.20, -18.20, -18.20]
})

sites_df
[5]:
time x y
0 2022-09-01 00:00:00 122.21 -18.2
1 2022-09-08 06:00:00 122.22 -18.2
2 2022-09-15 12:00:00 122.23 -18.2
3 2022-09-22 18:00:00 122.24 -18.2
4 2022-09-30 00:00:00 122.25 -18.2

Using “one-to-one” mode, we can model a tide height for each row in our dataframe, and add it as a new dataframe column:

[6]:
# Model tides in "one-to-one" mode
tide_df = model_tides(
    x=sites_df.x,
    y=sites_df.y,
    time=sites_df.time,
    mode="one-to-one",
)

# Add results as a new datframe column
sites_df["tide_height"] = tide_df.tide_m.values
sites_df.style.set_properties(**{"background-color": "#FFFF8F"}, subset=["tide_height"])
Modelling tides using FES2014 in parallel
100%|██████████| 5/5 [00:06<00:00,  1.34s/it]
[6]:
  time x y tide_height
0 2022-09-01 00:00:00 122.210000 -18.200000 -3.606788
1 2022-09-08 06:00:00 122.220000 -18.200000 -1.738778
2 2022-09-15 12:00:00 122.230000 -18.200000 -3.303217
3 2022-09-22 18:00:00 122.240000 -18.200000 -0.859731
4 2022-09-30 00:00:00 122.250000 -18.200000 -3.687130

Modelling tide heights for each satellite observation using tidal_tag

However, often it is valuable to estimate the height of the tide at the exact moment satellite images were taken over a coastline. 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.

To demonstrate how this can be done, we first need to load in an example time series of satellite imagery.

Set up data query

First we set up a query to define the area, time period and other parameters required for loading data. In this example, we will load several years of Landsat 8 data for intertidal mud flats south of Broome in Western Australia. We load the 'nbart_red', 'nbart_green', 'nbart_blue' bands so that we can plot the data as true colour imagery.

The dask_chunks parameter allows us to use Dask to lazily load data rather than load data directly into memory, which can take a long time and large amounts of memory. Lazy loading can be a very useful approach for when you need to load large amounts of data without crashing your analysis. In coastal applications, it allows us to load (using either .compute() or by plotting our data) only a small subset of observations from our entire time series (e.g. only low or high tide observations) without having to load the entire dataset into memory first, which can greatly decrease processing times.

For more information about using Dask, refer to the Parallel processing with Dask notebook.

[7]:
# Set up data load query
query = {
    'x': (122.10, 122.33),
    'y': (-18.09, -18.31),
    'time': ('2019', '2021'),
    'measurements': ['nbart_red', 'nbart_green', 'nbart_blue'],
    'output_crs': 'EPSG:32755',
    'resolution': (-30, 30),
    'group_by': 'solar_day',
    'dask_chunks': {},
}

We can preview the area that we will load data for:

[8]:
display_map(x=query['x'], y=query['y'])
[8]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load satellite time-series

To obtain some satellite data to analyse, we use the load_ard function to import a time series of Landsat 8 observations as an xarray.Dataset. The input data does not need to be from Landsat: any remotely-sensed imagery with timestamps and spatial coordinates provide enough data to run the tidal model.

[9]:
# Load available data from Landsat 8
ds = load_ard(dc=dc,
              products=["ga_ls8c_ard_3"],
              min_gooddata=0.9,
              mask_filters=[("opening", 5), ("dilation", 5)],
              ls7_slc_off=False,
              **query)

# Print output data
ds
Finding datasets
    ga_ls8c_ard_3
Counting good quality pixels for each time step using fmask
/env/lib/python3.8/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(
Filtering to 36 out of 136 time steps with at least 90.0% good quality pixels
Applying morphological filters to pixel quality mask: [('opening', 5), ('dilation', 5)]
Applying fmask pixel quality/cloud mask
Returning 36 time steps as a dask array
[9]:
<xarray.Dataset>
Dimensions:      (time: 36, y: 1003, x: 1002)
Coordinates:
  * time         (time) datetime64[ns] 2019-02-02T01:55:03.684910 ... 2021-12...
  * y            (y) float64 7.812e+06 7.812e+06 ... 7.782e+06 7.782e+06
  * x            (x) float64 -2.204e+06 -2.204e+06 ... -2.174e+06 -2.174e+06
    spatial_ref  int32 32755
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
Attributes:
    crs:           EPSG:32755
    grid_mapping:  spatial_ref

“Tagging” each satellite image with tide height metadata

We can now use the tidal_tag function from dea_tools.coastal to “tag” or associate each satellite observation in our time series with a tide height relative to mean sea level (i.e. approximately equivalent to the Australian Height Datum or AHD). This function uses the time and date of acquisition and the geographic location of each satellite observation as inputs to the Finite Element Solution 2014 (FES2014) tidal model:

FES2014 is the last version of the FES (Finite Element Solution) tide model developed in 2014-2016. It is an improved version of the FES2012 model. This new FES2014 model has been developed, implemented and validated by the LEGOS, NOVELTIS and CLS, within a CNES funded project. FES2014 takes advantage of longer altimeter time series and better altimeter standards, improved modelling and data assimilation techniques, a more accurate ocean bathymetry and a refined mesh in most of the shallow water regions. Special efforts have been dedicated to address the major non-linear tides issue and to the determination of accurate tidal currents.

The function will automatically select a tide modelling location based on the dataset centroid. It will then output modelled tide heights as a new tide_m variable in the xarray.Dataset (the variable should appear under Data variables below):

[10]:
# Model tide heights
ds_tidal = tidal_tag(ds)

# Print output data
ds_tidal

Setting tide modelling location from dataset centroid: 122.22, -18.20
Modelling tides using FES2014
[10]:
<xarray.Dataset>
Dimensions:      (time: 36, y: 1003, x: 1002)
Coordinates:
  * time         (time) datetime64[ns] 2019-02-02T01:55:03.684910 ... 2021-12...
  * y            (y) float64 7.812e+06 7.812e+06 ... 7.782e+06 7.782e+06
  * x            (x) float64 -2.204e+06 -2.204e+06 ... -2.174e+06 -2.174e+06
    spatial_ref  int32 32755
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    tide_m       (time) float32 1.491 2.212 2.516 ... 3.096 1.183 -0.5053
Attributes:
    crs:           EPSG:32755
    grid_mapping:  spatial_ref

We can easily plot this new variable to inspect the range of tide heights observed by the satellites in our time series. In this example, our observed tide heights range from approximately -2.0 to 4.0 m relative to Mean Sea Level:

[11]:
ds_tidal.tide_m.plot();
../../../_images/notebooks_How_to_guides_Tidal_modelling_24_0.png

Example tide height analysis

To demonstrate how tidally tagged images can be used to produce composites of high and low tide imagery, we can compute the lowest 20% and highest 20% of tide heights, and use these to filter our observations. We can then combine and plot these filtered observations to visualise how the landscape appears at low and high tide:

[12]:
# Calculate the lowest and highest 10% of tides
lowest_20, highest_20 = ds_tidal.tide_m.quantile([0.2, 0.8]).values

# Filter our data to low and high tide observations
filtered_low = ds_tidal.where(ds_tidal.tide_m <= lowest_20, drop=True)
filtered_high = ds_tidal.where(ds_tidal.tide_m >= highest_20, drop=True)

# Take the simple median of each set of low and high tide observations to
# produce a composite (alternatively, observations could be combined
# using a geomedian to keep band relationships consistent)
median_low = filtered_low.median(dim='time', keep_attrs=True)
median_high = filtered_high.median(dim='time', keep_attrs=True)

# Combine low and high tide medians into a single dataset and give
# each layer a meaningful name
ds_highlow = xr.concat([median_low, median_high], dim='tide_m')
ds_highlow['tide_m'] = ['Low tide', 'High tide']

# Plot low and high tide medians side-by-side
rgb(ds_highlow, col='tide_m')

../../../_images/notebooks_How_to_guides_Tidal_modelling_26_0.png

Swapping dimensions

The tidal_tag function allows you to use tide heights as the primary dimension in the dataset, rather than time. Setting swap_dims=True will swap the time dimension in the original xarray.Dataset to the new tide_m variable.

[13]:
# Model tide heights
ds_tidal = tidal_tag(ds, swap_dims=True)

# Print output data
ds_tidal
Setting tide modelling location from dataset centroid: 122.22, -18.20
Modelling tides using FES2014
[13]:
<xarray.Dataset>
Dimensions:      (tide_m: 36, y: 1003, x: 1002)
Coordinates:
  * y            (y) float64 7.812e+06 7.812e+06 ... 7.782e+06 7.782e+06
  * x            (x) float64 -2.204e+06 -2.204e+06 ... -2.174e+06 -2.174e+06
    spatial_ref  int32 32755
  * tide_m       (tide_m) float32 -2.372 -1.68 -0.7717 ... 3.243 3.272 3.562
Data variables:
    nbart_red    (tide_m, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_green  (tide_m, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_blue   (tide_m, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
Attributes:
    crs:           EPSG:32755
    grid_mapping:  spatial_ref

The dataset now contains three dimensions: tide_m, x and y. This can make it easier to analyse the data with respect to tide, e.g. plotting the lowest and highest tide image in our dataset:

[14]:
# Plot first (lowest tide) and last (highest tide) images
rgb(ds_tidal, index_dim='tide_m', index=[0, -1])
../../../_images/notebooks_How_to_guides_Tidal_modelling_30_0.png

Note: The white dots in the images above are false positives in the automated cloud mask caused by bright sandy shorelines being mistaken for clouds.

Modelling ebb and flow tidal phases

The tidal_tag function also allows us to determine whether each satellite observation was taken while the tide was rising/incoming (flow tide) or falling/outgoing (ebb tide) by setting ebb_flow=True. This is achieved by comparing tide heights 15 minutes before and after the observed satellite observation.

Ebb and flow data can provide valuable contextual information for interpreting satellite imagery, particularly in tidal flat or mangrove forest environments where water may remain in the landscape for considerable time after the tidal peak.

[15]:
# Model tide heights
ds_tidal = tidal_tag(ds, ebb_flow=True)

# Print output data
ds_tidal
Setting tide modelling location from dataset centroid: 122.22, -18.20
Modelling tides using FES2014
Modelling tidal phase (e.g. ebb or flow)
Modelling tides using FES2014
[15]:
<xarray.Dataset>
Dimensions:      (time: 36, y: 1003, x: 1002)
Coordinates:
  * time         (time) datetime64[ns] 2019-02-02T01:55:03.684910 ... 2021-12...
  * y            (y) float64 7.812e+06 7.812e+06 ... 7.782e+06 7.782e+06
  * x            (x) float64 -2.204e+06 -2.204e+06 ... -2.174e+06 -2.174e+06
    spatial_ref  int32 32755
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    tide_m       (time) float32 1.491 2.212 2.516 ... 3.096 1.183 -0.5053
    ebb_flow     (time) <U4 'Flow' 'Flow' 'Flow' 'Ebb' ... 'Flow' 'Flow' 'Flow'
Attributes:
    crs:           EPSG:32755
    grid_mapping:  spatial_ref

We now have data giving us the both the tide height and tidal phase (‘ebb’ or ‘flow’) for every satellite image:

[16]:
ds_tidal[['time', 'tide_m', 'ebb_flow']].to_dataframe().drop('spatial_ref',
                                                             axis=1)
[16]:
tide_m ebb_flow
time
2019-02-02 01:55:03.684910 1.490903 Flow
2019-02-18 01:55:01.585984 2.211541 Flow
2019-03-06 01:54:56.923564 2.516171 Flow
2019-07-12 01:55:13.930435 -0.324644 Ebb
2019-07-28 01:55:19.106108 -0.132350 Ebb
2019-08-13 01:55:25.096137 1.924859 Ebb
2019-08-29 01:55:29.115074 2.744836 Ebb
2019-09-14 01:55:33.783014 2.975977 Flow
2019-09-30 01:55:38.714933 3.113354 Flow
2019-10-16 01:55:41.223276 1.808664 Flow
2019-11-01 01:55:41.070156 -0.087908 Flow
2019-11-17 01:55:38.012099 -0.449818 Flow
2019-12-03 01:55:36.758377 -1.679660 Flow
2019-12-19 01:55:34.623293 -2.372203 Flow
2020-03-08 01:55:12.515274 2.302160 Flow
2020-03-24 01:55:04.353944 2.814702 Flow
2020-04-09 01:54:56.077227 3.272317 Flow
2020-04-25 01:54:48.941936 2.009573 Flow
2020-08-15 01:55:18.231542 0.121439 Ebb
2020-08-31 01:55:26.558107 2.191674 Ebb
2020-09-16 01:55:32.698914 2.813782 Flow
2020-10-02 01:55:36.635460 3.098040 Flow
2020-10-18 01:55:37.892585 3.140818 Flow
2020-12-05 01:55:38.895041 -0.627862 Flow
2021-03-27 01:55:05.454732 2.624109 Flow
2021-04-12 01:55:00.410953 3.218125 Flow
2021-04-28 01:54:52.230407 3.561958 Flow
2021-05-14 01:54:52.033925 1.869562 Flow
2021-06-15 01:55:07.799884 -0.100711 Flow
2021-08-18 01:55:27.286586 -0.771717 Ebb
2021-09-03 01:55:31.946655 0.578454 Ebb
2021-10-05 01:55:40.578676 3.021137 Flow
2021-10-21 01:55:43.617074 3.243399 Flow
2021-11-06 01:55:42.131253 3.096294 Flow
2021-11-22 01:55:38.241313 1.183015 Flow
2021-12-08 01:55:39.297251 -0.505341 Flow

We could for example use this data to filter our observations to keep ebbing phase observations only:

[17]:
ds_tidal.where(ds_tidal.ebb_flow == 'Ebb', drop=True)
[17]:
<xarray.Dataset>
Dimensions:      (time: 8, y: 1003, x: 1002)
Coordinates:
  * time         (time) datetime64[ns] 2019-07-12T01:55:13.930435 ... 2021-09...
  * y            (y) float64 7.812e+06 7.812e+06 ... 7.782e+06 7.782e+06
  * x            (x) float64 -2.204e+06 -2.204e+06 ... -2.174e+06 -2.174e+06
    spatial_ref  int32 32755
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    tide_m       (time) float32 -0.3246 -0.1324 1.925 ... 2.192 -0.7717 0.5785
    ebb_flow     (time) object 'Ebb' 'Ebb' 'Ebb' 'Ebb' 'Ebb' 'Ebb' 'Ebb' 'Ebb'
Attributes:
    crs:           EPSG:32755
    grid_mapping:  spatial_ref

Advanced

Modelling tides for each pixel in a satellite time series 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 up to several metres across a distance of kilometres in areas of complex tidal dynamics. This means that an individual satellite image can contain a range of tide height conditions.

To capture this spatial variability in tide heights, we can use the pixel_tides function from dea_tools.coastal. For efficient processing, this function first models tides into a low resolution 5000 m grid surrounding each satellite image in our time series. This lower resolution data will also include a buffer around the extent of our data tides can be modelled seamlessly across analysis boundaries:

[18]:
# Extract a subset of our satellite data above to make plotting easier
ds_subset = ds.isel(time=slice(9, 13))

# Model tides spatially using `pixel_tides`
tides_lowres = pixel_tides(ds_subset, resample=False)
Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides using FES2014 in parallel
100%|██████████| 5/5 [00:06<00:00,  1.21s/it]
Returning low resolution tide array

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

Blue values below indicate low tide pixels, while red indicates high tide pixels. If you look closely at the second and third timestep below, you may also be able to see some spatial variability in tide heights within each timestep, with slightly higher tides in the north-west of the images:

[19]:
tides_lowres.plot.imshow(col='time', vmin=-1, vmax=1, cmap='RdBu')
[19]:
<xarray.plot.facetgrid.FacetGrid at 0x7fad57c011f0>
../../../_images/notebooks_How_to_guides_Tidal_modelling_41_1.png

Reprojecting tide heights back 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 value for every single pixel in our dataset through time and space:

[20]:
# Model tides spatially using `pixel_tides`
tides_highres, tides_lowres = pixel_tides(ds_subset, resample=True)

tides_highres.plot.imshow(col='time', vmin=-1, vmax=1, cmap='RdBu')
Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides using FES2014 in parallel
100%|██████████| 5/5 [00:05<00:00,  1.18s/it]
Reprojecting tides into original array
[20]:
<xarray.plot.facetgrid.FacetGrid at 0x7fad57c11d60>
../../../_images/notebooks_How_to_guides_Tidal_modelling_43_4.png

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

[21]:
ds_subset.sizes
[21]:
Frozen({'time': 4, 'y': 1003, 'x': 1002})
[22]:
tides_highres.sizes
[22]:
Frozen({'time': 4, 'y': 1003, 'x': 1002})

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

[23]:
ds_subset['tide_m'] = tides_highres
ds_subset
[23]:
<xarray.Dataset>
Dimensions:      (time: 4, y: 1003, x: 1002)
Coordinates:
  * time         (time) datetime64[ns] 2019-10-16T01:55:41.223276 ... 2019-12...
  * y            (y) float64 7.812e+06 7.812e+06 ... 7.782e+06 7.782e+06
  * x            (x) float64 -2.204e+06 -2.204e+06 ... -2.174e+06 -2.174e+06
    spatial_ref  int32 32755
    tide_model   <U7 'FES2014'
Data variables:
    nbart_red    (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_green  (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    nbart_blue   (time, y, x) float32 dask.array<chunksize=(1, 1003, 1002), meta=np.ndarray>
    tide_m       (time, y, x) float32 1.847 1.847 1.847 ... -1.68 -1.68 -1.68
    ebb_flow     (time) <U4 'Flow' 'Flow' 'Flow' 'Flow'
Attributes:
    crs:           EPSG:32755
    grid_mapping:  spatial_ref

Calculating min/max/median/quantiles of tide heights 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.

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:

[24]:
tides_highres_quantiles, tides_lowres_quantiles = pixel_tides(
    ds_subset, resample=True, calculate_quantiles=(0, 0.5, 1))
Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides using FES2014 in parallel
100%|██████████| 5/5 [00:05<00:00,  1.19s/it]
Computing tide quantiles
Reprojecting tides into original array
[25]:
tides_highres_quantiles.plot.imshow(col='quantile')
[25]:
<xarray.plot.facetgrid.FacetGrid at 0x7fadfc4324c0>
../../../_images/notebooks_How_to_guides_Tidal_modelling_51_1.png

Pixel-based tides for 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:

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

tides_highres, tides_lowres = pixel_tides(ds_subset,
                                          resample=True,
                                          times=custom_times)
Creating reduced resolution 5000 x 5000 metre tide modelling array
Modelling tides using FES2014 in parallel
100%|██████████| 5/5 [00:06<00:00,  1.22s/it]
Reprojecting tides into original array
[27]:
tides_highres.plot.imshow(col='time')
[27]:
<xarray.plot.facetgrid.FacetGrid at 0x7fad602cdca0>
../../../_images/notebooks_How_to_guides_Tidal_modelling_54_1.png

Evaluating tidal biases using tidal_stats

The complex temporal behaviour of tides mean that a sun synchronous sensor like Landsat does not observe the full range of the tidal cycle at all locations. Biases in the proportion of the tidal range observed by satellites can prevent us from obtaining data on areas of the coastline exposed or inundated at the extremes of the tidal range. This can risk gaining misleading insights into the true extent of the area of the coastline affected by tides, and make it difficult to compare high or low tide images fairly in different locations.

The tidal_stats function can assist in evaluating how the range of tides observed by satellites compare to the full tidal range. It works by using our tidal model to to model all available tide heights at a regular interval (every two hours by default) across the entire time period covered by the input satellite time series dataset. This is then compared against the tide heights in observed by the satellite and used to calculate a range of statistics and a plot that summarises potential biases in the data.

For a more detailed discussion of the issue of tidal bias in sun-synchronous satellite observations of the coastline, refer to the ‘Limitations and future work’ section in Bishop-Taylor et al. 2018.

[28]:
out_stats = tidal_stats(ds)
Setting tide modelling location from dataset centroid: 122.22, -18.20
Modelling tides using FES2014
Modelling tides using FES2014
                                                 tide_model    tide_m
time                       x          y
2019-02-02 01:55:03.684910 122.215013 -18.200098    FES2014  1.490903
2019-02-02 03:55:03.684910 122.215013 -18.200098    FES2014  0.854936
2019-02-02 05:55:03.684910 122.215013 -18.200098    FES2014 -0.925995
2019-02-02 07:55:03.684910 122.215013 -18.200098    FES2014 -1.479792
2019-02-02 09:55:03.684910 122.215013 -18.200098    FES2014 -0.741872
...                                                     ...       ...
2021-12-07 17:55:03.684910 122.215013 -18.200098    FES2014  2.950938
2021-12-07 19:55:03.684910 122.215013 -18.200098    FES2014 -0.416146
2021-12-07 21:55:03.684910 122.215013 -18.200098    FES2014 -3.455013
2021-12-07 23:55:03.684910 122.215013 -18.200098    FES2014 -3.563971
2021-12-08 01:55:03.684910 122.215013 -18.200098    FES2014 -0.522441

[12481 rows x 2 columns]

59% of the 10.04 m modelled astronomical tidal range is observed at this location.
The lowest 26% and highest 15% of astronomical tides are never observed.

../../../_images/notebooks_How_to_guides_Tidal_modelling_56_1.png

The function also outputs a pandas.Series object containing a set of statistics that compare the observed vs. full modelled tidal ranges. These statistics include:

  • tidepost_lat: latitude used for modelling tide heights

  • tidepost_lon: longitude used for modelling tide heights

  • observed_min_m: minimum tide height observed by the satellite (in metre units)

  • all_min_m: minimum tide height from full modelled tidal range (in metre units)

  • observed_max_m: maximum tide height observed by the satellite (in metre units)

  • all_max_m: maximum tide height from full modelled tidal range (in metre units)

  • observed_range_m: tidal range observed by the satellite (in metre units)

  • all_range_m: full modelled tidal range (in metre units)

  • spread_m: proportion of the full modelled tidal range observed by the satellite (see Bishop-Taylor et al. 2018)

  • low_tide_offset: proportion of the lowest tides never observed by the satellite (see Bishop-Taylor et al. 2018)

  • high_tide_offset: proportion of the highest tides never observed by the satellite (see Bishop-Taylor et al. 2018)

[29]:
out_stats
[29]:
tidepost_lat        -18.200
tidepost_lon        122.215
observed_mean_m       1.553
all_mean_m           -0.000
observed_min_m       -2.372
all_min_m            -5.019
observed_max_m        3.562
all_max_m             5.025
observed_range_m      5.934
all_range_m          10.044
spread                0.591
low_tide_offset       0.263
high_tide_offset      0.146
dtype: float64

Additional information

FES2014: Tidal modelling is provided by the FES2014 global tidal model, implemented using functions from the pyTMD Python package. FES2014 was produced by NOVELTIS, LEGOS, CLS Space Oceanography Division and CNES. It is distributed by AVISO, with support from CNES (https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes/description-fes2014.html).

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:

[30]:
print(datacube.__version__)
1.8.13

Tags

Tags: sandbox compatible, NCI compatible, landsat 8, display_map, load_ard, rgb, tidal_tag, tidal_stats, model_tides, pixel_tides, tide modelling, intertidal, Dask, lazy loading, rolling window