Combining satellite data with tidal modelling using OTPS

• Products used: ls5_sr, ls7_sr, ls8_sr

Keywords data used; landsat 5, data used; landsat 7, data used; landsat 8, tidal; modelling, intertidal, dask

Contexte

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. To be able to study environmental processes along Africa’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 tidally tag remotely sensed imagery using functions from the deafrica_tools.coastal package so that images can be extracted or analysed by tidal stage (e.g. low, high, ebb, flow). These functions use the OTPS TPXO8 tidal model to calculate the height (relative to mean sea level) and stage of the tide at the exact moment each satellite image was acquired.

The notebook demonstrates how to:

1. Load an example time series of satellite data

2. Use the tidal_tag function from deafrica_tools.coastal to model tide heights for each satellite observation

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

4. Swap a dataset’s dimensions to compute a rolling median along the tide_height dimension

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

6. Use the tidal_stats function to evaluate any biases in the tidal conditions observed by a satellite

---

Getting started

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

Load packages

%matplotlib inline

import datacube
import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt
from datacube.utils.geometry import Geometry

from deafrica_tools.plotting import rgb, display_map
from deafrica_tools.datahandling import load_ard, mostcommon_crs
from deafrica_tools.coastal import tidal_tag, tidal_stats
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.areaofinterest import define_area

Set up a Dask cluster

Dask can be used to better manage memory use down and conduct the analysis in parallel. For an introduction to using Dask with Digital Earth Africa, see the Dask notebook.

Note: We recommend opening the Dask processing window to view the different computations that are being executed; to do this, see the Dask dashboard in DE Africa section of the Dask notebook.

To use Dask, set up the local computing cluster using the cell below.

create_local_dask_cluster()

Connect to the datacube

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

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 30 years of Landsat 5, 7 and 8 data for the Geba River Estuary in southern Guinea-Bissau. We load the 'red', 'green', 'blue' bands so that we can plot the data as true colour imagery.

Note: The dask_chunks parameter allows us to use Dask to lazily load data rather than load data directly into memory, as the standard loading method 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.

Define location

To define the area of interest, there are two methods available:

1. By specifying the latitude, longitude, and buffer. This method requires you to input the central latitude, central longitude, and the buffer value in square degrees around the center point you want to analyze. For example, lat = 10.338, lon = -1.055, and buffer = 0.1 will select an area with a radius of 0.1 square degrees around the point with coordinates (10.338, -1.055).

2. By uploading a polygon as a GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.

To use one of these methods, you can uncomment the relevant line of code and comment out the other one. To comment out a line, add the "#" symbol before the code you want to comment out. By default, the first option which defines the location using latitude, longitude, and buffer is being used.

# Define the location
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=11.689, lon=-15.674, buffer=0.15)

# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')

#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)

# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])

# Create a reusable query
query = {
    'x': lon_range,
    'y': lat_range,
    'time': ('1988-01-01', '2018-12-31'),
    'measurements': ['red', 'green', 'blue'],
    'resolution': (-30, 30),
    'dask_chunks': {}
}

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

display_map(x=query['x'], y=query['y'])

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 5, 7 and 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.

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

# Load available data from all three Landsat satellites
ds = load_ard(dc=dc,
              products=['ls5_sr', 'ls7_sr', 'ls8_sr'],
              output_crs=output_crs,
              align=(15, 15),
              ls7_slc_off=False,
              group_by='solar_day',
              **query)

# Print output data
print(ds)

Using pixel quality parameters for USGS Collection 2
Finding datasets
    ls5_sr
    ls7_sr
    Ignoring SLC-off observations for ls7
    ls8_sr
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Returning 227 time steps as a dask array
<xarray.Dataset>
Dimensions:      (time: 227, y: 1109, x: 1093)
Coordinates:
  * time         (time) datetime64[ns] 1988-07-17T10:52:47.658038 ... 2018-12...
  * y            (y) float64 1.309e+06 1.309e+06 ... 1.276e+06 1.276e+06
  * x            (x) float64 4.102e+05 4.102e+05 ... 4.429e+05 4.429e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

Model tide heights for each observation

We use the tidal_tag function from deafrica_coastaltools to associate each satellite observation in our timeseries with a tide height relative to mean sea level. This function uses the time and date of acquisition and the geographic location of each satellite observation as inputs to the OSU Tidal Prediction Software (OTPS) tidal model. From Sagar et al. 2015:

Note: The OTPS TPX08 tidal model consists of a multi-resolution bathymetric grid solution, with a 1/6° solution in the global open ocean, and a 1/30° local resolution solution to improve modelling in complex shallow water environments. The OTPS model is based on a system of linear partial differential equations, called Laplace’s tidal equations, parametrised with nine harmonic tidal constituents. The model is fitted to track-averaged TOPEX/Poseidon altimeter data collected from 1992 to 2016 and Jason-1 (Poseidon 2) altimeter data from 2002 to 2013, enabling estimation of the tidal height and harmonic constituents at discrete temporal epochs and spatial locations.

# Model tide heights
ds_tidal = tidal_tag(ds)

# Print output data
print(ds_tidal)

Setting tide modelling location from dataset centroid: -15.67, 11.69
<xarray.Dataset>
Dimensions:      (time: 227, y: 1109, x: 1093)
Coordinates:
  * time         (time) datetime64[ns] 1988-07-17T10:52:47.658038 ... 2018-12...
  * y            (y) float64 1.309e+06 1.309e+06 ... 1.276e+06 1.276e+06
  * x            (x) float64 4.102e+05 4.102e+05 ... 4.429e+05 4.429e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    tide_height  (time) float64 -0.194 -0.605 0.168 ... 0.535 -0.217 -0.821
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

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_height variable in the xarray.Dataset (the variable should appear under Data variables above).

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

ds_tidal.tide_height.plot(linewidth=0.5)

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

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 5% and highest 5% percent of tide heights, and use these to filter our observations. We can then combine and plot these filtered observations to visualise what the landscape looks like at low and high tide:

Note: We recommend opening the Dask processing window to view the different computations that are being executed; to do this, see the Dask dashboard in DE Africa section of the Dask notebook.

Note: An alternative approach to combining observations into a composite can be achieved using a geomedian, which ensures band relationships are kept consistent. More information is provided in the Generating geomedian composites notebook.

# Calculate the lowest and highest 5% of tides
lowest_5, highest_5 = ds_tidal.tide_height.quantile([0.05, 0.95]).values

# Filter our data to low and high tide observations
filtered_low = ds_tidal.where(ds_tidal.tide_height <= lowest_5, drop=True)
filtered_high = ds_tidal.where(ds_tidal.tide_height >= highest_5, 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_height')
ds_highlow['tide_height'] = ['Low tide', 'High tide']

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

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_height variable.

# Model tide heights
ds_tidal = tidal_tag(ds, swap_dims=True)

# Print output data
print(ds_tidal)

Setting tide modelling location from dataset centroid: -15.67, 11.69
<xarray.Dataset>
Dimensions:      (tide_height: 227, y: 1109, x: 1093)
Coordinates:
  * y            (y) float64 1.309e+06 1.309e+06 ... 1.276e+06 1.276e+06
  * x            (x) float64 4.102e+05 4.102e+05 ... 4.429e+05 4.429e+05
    spatial_ref  int32 32628
  * tide_height  (tide_height) float64 -0.978 -0.978 -0.952 ... 1.284 1.325
Data variables:
    red          (tide_height, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    green        (tide_height, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    blue         (tide_height, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

The dataset now contains three dimensions: tide_height, x and y.

This can make it easier to analyse the data with respect to tide, e.g. computing a rolling median by tide height (e.g. along the tide_height dimension), as is done in the following cell.

Note: We recommend opening the Dask processing window to view the different computations that are being executed; to do this, see the Dask dashboard in DE Africa section of the Dask notebook.

# First we need to update the chunks used by Dask to allow us to do a
# rolling median without having to load all our data into memory first
ds_rechunked = ds_tidal.chunk(chunks={'tide_height': 13})

# Compute a rolling median that will go through every satellite
# observation, and take the median of that timestep and its 15 neighbours
ds_rolling = ds_rechunked.rolling(tide_height=13,
                                  center=True,
                                  min_periods=1).median()

# Plot the lowest, median and highest tide rolling median image
rgb(ds_rolling, index_dim='tide_height', index=[0, len(ds_tidal.tide_height)//2, -1])

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

# Model tide heights
ds_tidal = tidal_tag(ds, ebb_flow=True)

# Print output data
print(ds_tidal)

Setting tide modelling location from dataset centroid: -15.67, 11.69
Modelling tidal phase (e.g. ebb or flow)
<xarray.Dataset>
Dimensions:      (time: 227, y: 1109, x: 1093)
Coordinates:
  * time         (time) datetime64[ns] 1988-07-17T10:52:47.658038 ... 2018-12...
  * y            (y) float64 1.309e+06 1.309e+06 ... 1.276e+06 1.276e+06
  * x            (x) float64 4.102e+05 4.102e+05 ... 4.429e+05 4.429e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    tide_height  (time) float64 -0.194 -0.605 0.168 ... 0.535 -0.217 -0.821
    ebb_flow     (time) <U4 'Flow' 'Flow' 'Ebb' 'Ebb' ... 'Flow' 'Flow' 'Flow'
Attributes:
    crs:           epsg:32628
    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:

ds_tidal[['time', 'tide_height', 'ebb_flow']].to_dataframe()

	tide_height	ebb_flow	spatial_ref
time
1988-07-17 10:52:47.658038	-0.194	Flow	32628
1988-08-18 10:52:49.179013	-0.605	Flow	32628
1988-10-05 10:52:43.680006	0.168	Ebb	32628
1988-10-21 10:52:37.662000	0.334	Ebb	32628
1988-12-08 10:52:31.815006	0.805	Flow	32628
...	...	...	...
2018-10-24 11:21:58.388427	1.094	Flow	32628
2018-11-09 11:22:01.151598	0.730	Flow	32628
2018-11-25 11:22:00.936674	0.535	Flow	32628
2018-12-11 11:21:57.773631	-0.217	Flow	32628
2018-12-27 11:21:57.806437	-0.821	Flow	32628

227 rows × 3 columns

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

ds_tidal.where(ds_tidal.ebb_flow == 'Ebb', drop=True)

<xarray.Dataset>
Dimensions:      (time: 106, y: 1109, x: 1093)
Coordinates:
  * time         (time) datetime64[ns] 1988-10-05T10:52:43.680006 ... 2018-10...
  * y            (y) float64 1.309e+06 1.309e+06 ... 1.276e+06 1.276e+06
  * x            (x) float64 4.102e+05 4.102e+05 ... 4.429e+05 4.429e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>
    tide_height  (time) float64 0.168 0.334 -0.648 -0.395 ... 0.065 0.724 1.284
    ebb_flow     (time) object 'Ebb' 'Ebb' 'Ebb' 'Ebb' ... 'Ebb' 'Ebb' 'Ebb'
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

xarray.Dataset

• Dimensions:
  • time: 106
  • y: 1109
  • x: 1093
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    1988-10-05T10:52:43.680006 ... 2...

    array(['1988-10-05T10:52:43.680006000', '1988-10-21T10:52:37.662000000',
           '1989-09-22T10:47:48.881000000', '1989-10-08T10:47:35.757019000',
           '1989-11-25T10:46:03.037063000', '1989-12-11T10:45:39.717013000',
           '1990-10-11T10:41:41.717013000', '1990-10-27T10:41:35.028050000',
           '1990-11-12T10:41:35.413013000', '1990-11-28T10:41:21.799050000',
           '1990-12-14T10:41:27.273031000', '1990-12-30T10:41:52.516081000',
           '1991-11-15T10:46:17.278000000', '1994-07-18T10:40:08.911025000',
           '1995-06-19T10:28:05.602013000', '1996-03-01T10:28:32.486056000',
           '1998-02-19T10:58:17.962063000', '1998-03-07T10:58:33.143063000',
           '1998-03-23T10:58:47.923088000', '1999-11-29T11:14:51.752254000',
           '2000-02-01T11:14:52.634172000', '2000-06-24T11:13:50.212411000',
           '2000-07-10T11:13:35.865891000', '2000-08-27T11:13:01.880397000',
           '2001-01-02T11:12:27.903011000', '2001-02-03T11:12:32.034612000',
           '2001-02-19T11:12:29.745162000', '2001-07-29T11:11:16.260082000',
           '2001-08-14T11:10:55.610240000', '2001-09-15T11:10:36.015143000',
           '2002-01-05T11:10:46.454847000', '2002-01-21T11:10:56.736123000',
           '2002-02-22T11:10:57.914884000', '2002-08-01T11:10:18.347957000',
           '2002-09-18T11:09:51.064379000', '2002-10-04T11:09:48.127048000',
           '2003-02-09T11:10:44.661508000', '2003-02-25T11:10:44.083847000',
           '2003-03-13T11:10:49.369281000', '2003-04-14T11:10:50.641171000',
           '2007-01-27T11:17:01.835063000', '2007-02-12T11:17:04.228056000',
           '2007-02-28T11:17:04.037038000', '2007-03-16T11:17:01.474038000',
           '2007-07-22T11:15:41.294019000', '2007-08-07T11:15:32.020056000',
           '2007-08-23T11:15:21.541031000', '2007-09-24T11:15:06.340069000',
           '2009-09-29T11:12:01.100056000', '2009-10-31T11:12:19.951069000',
           '2010-11-03T11:11:54.285038000', '2010-11-19T11:11:53.312019000',
           '2011-10-21T11:10:04.063031000', '2013-05-03T11:23:53.310888000',
           '2013-05-19T11:24:06.077361000', '2013-07-06T11:24:04.405354000',
           '2013-11-11T11:23:49.100948000', '2013-11-27T11:23:41.849117000',
           '2013-12-13T11:23:38.930442000', '2013-12-29T11:23:30.292655000',
           '2014-01-14T11:23:15.376568000', '2014-05-22T11:21:34.999331000',
           '2014-06-07T11:21:45.055068000', '2014-06-23T11:21:47.968851000',
           '2014-07-09T11:21:56.081333000', '2014-07-25T11:21:58.141203000',
           '2014-11-14T11:22:18.466748000', '2014-11-30T11:22:16.880842000',
           '2014-12-16T11:22:12.037765000', '2015-01-01T11:22:06.615041000',
           '2015-01-17T11:22:05.793171000', '2015-02-02T11:22:01.893389000',
           '2015-05-25T11:21:09.083851000', '2015-06-10T11:21:19.878233000',
           '2015-06-26T11:21:25.078335000', '2015-07-12T11:21:36.708205000',
           '2015-08-13T11:21:46.824490000', '2015-12-03T11:22:13.740165000',
           '2015-12-19T11:22:13.664612000', '2016-01-04T11:22:10.202099000',
           '2016-01-20T11:22:10.392217000', '2016-02-05T11:22:06.795677000',
           '2016-06-12T11:21:46.786090000', '2016-06-28T11:21:54.490874000',
           '2016-07-14T11:22:02.039857000', '2016-08-15T11:22:09.140616000',
           '2016-12-21T11:22:19.326793000', '2017-01-06T11:22:16.948666000',
           '2017-01-22T11:22:12.191980000', '2017-02-07T11:22:04.129007000',
           '2017-02-23T11:21:58.867480000', '2017-07-01T11:21:51.930455000',
           '2017-07-17T11:21:55.828621000', '2017-08-02T11:22:03.551977000',
           '2017-09-03T11:22:11.680478000', '2017-09-19T11:22:14.104889000',
           '2018-01-09T11:22:11.054208000', '2018-01-25T11:22:02.430782000',
           '2018-02-10T11:21:55.226966000', '2018-02-26T11:21:48.758176000',
           '2018-03-14T11:21:39.531096000', '2018-07-20T11:21:17.343248000',
           '2018-08-05T11:21:25.659501000', '2018-08-21T11:21:34.094337000',
           '2018-09-22T11:21:44.685387000', '2018-10-08T11:21:52.881616000'],
          dtype='datetime64[ns]')

  • y
    (y)
    float64
    1.309e+06 1.309e+06 ... 1.276e+06

    units :

    metre

    resolution :

    -30.0

    crs :

    epsg:32628

    array([1308870., 1308840., 1308810., ..., 1275690., 1275660., 1275630.])

  • x
    (x)
    float64
    4.102e+05 4.102e+05 ... 4.429e+05

    units :

    metre

    resolution :

    30.0

    crs :

    epsg:32628

    array([410160., 410190., 410220., ..., 442860., 442890., 442920.])

  • spatial_ref
    ()
    int32
    32628

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]

    grid_mapping_name :

    transverse_mercator

    array(32628, dtype=int32)

• Data variables: (5)
  • red
    (time, y, x)
    float32
    dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>

    crs :

    epsg:32628

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 490.14 MiB 4.62 MiB Shape (106, 1109, 1093) (1, 1109, 1093) Dask graph 106 chunks in 13 graph layers Data type float32 numpy.ndarray

  • green
    (time, y, x)
    float32
    dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>

    crs :

    epsg:32628

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 490.14 MiB 4.62 MiB Shape (106, 1109, 1093) (1, 1109, 1093) Dask graph 106 chunks in 13 graph layers Data type float32 numpy.ndarray

  • blue
    (time, y, x)
    float32
    dask.array<chunksize=(1, 1109, 1093), meta=np.ndarray>

    crs :

    epsg:32628

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 490.14 MiB 4.62 MiB Shape (106, 1109, 1093) (1, 1109, 1093) Dask graph 106 chunks in 13 graph layers Data type float32 numpy.ndarray

  • tide_height
    (time)
    float64
    0.168 0.334 -0.648 ... 0.724 1.284

    array([ 0.168,  0.334, -0.648, -0.395,  0.66 ,  0.995, -0.633, -0.324,
            0.28 ,  0.202,  0.618,  0.918, -0.326,  0.038, -0.816,  0.441,
           -0.676, -0.175, -0.128, -0.896,  0.249, -0.648, -0.235,  0.876,
           -0.678, -0.35 ,  0.316, -0.218, -0.195,  0.908, -0.978, -0.675,
           -0.291, -0.563,  0.66 ,  0.929, -0.655, -0.382, -0.184,  0.713,
           -0.56 , -0.315,  0.321,  0.49 , -0.568, -0.455, -0.093,  0.753,
            0.193,  0.735,  0.564,  0.689, -0.36 , -0.596, -0.394,  0.793,
           -0.586, -0.37 ,  0.223,  0.232,  0.684, -0.63 , -0.316,  0.303,
            0.343,  0.873, -0.622, -0.658, -0.322,  0.212,  0.263,  0.752,
           -0.647, -0.66 , -0.234,  0.304,  0.97 , -0.651, -0.713, -0.264,
            0.178,  0.325, -0.652, -0.666, -0.144,  0.512, -0.682, -0.733,
           -0.195,  0.134,  0.411, -0.661, -0.638, -0.042,  0.62 ,  1.193,
           -0.722, -0.717, -0.118,  0.1  ,  0.507, -0.675, -0.577,  0.065,
            0.724,  1.284])

  • ebb_flow
    (time)
    object
    'Ebb' 'Ebb' 'Ebb' ... 'Ebb' 'Ebb'

    array(['Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb',
           'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb', 'Ebb'], dtype=object)

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1988-10-05 10:52:43.680006', '1988-10-21 10:52:37.662000',
                   '1989-09-22 10:47:48.881000', '1989-10-08 10:47:35.757019',
                   '1989-11-25 10:46:03.037063', '1989-12-11 10:45:39.717013',
                   '1990-10-11 10:41:41.717013', '1990-10-27 10:41:35.028050',
                   '1990-11-12 10:41:35.413013', '1990-11-28 10:41:21.799050',
                   ...
                   '2018-01-09 11:22:11.054208', '2018-01-25 11:22:02.430782',
                   '2018-02-10 11:21:55.226966', '2018-02-26 11:21:48.758176',
                   '2018-03-14 11:21:39.531096', '2018-07-20 11:21:17.343248',
                   '2018-08-05 11:21:25.659501', '2018-08-21 11:21:34.094337',
                   '2018-09-22 11:21:44.685387', '2018-10-08 11:21:52.881616'],
                  dtype='datetime64[ns]', name='time', length=106, freq=None))

  • y
    PandasIndex

    PandasIndex(Index([1308870.0, 1308840.0, 1308810.0, 1308780.0, 1308750.0, 1308720.0,
           1308690.0, 1308660.0, 1308630.0, 1308600.0,
           ...
           1275900.0, 1275870.0, 1275840.0, 1275810.0, 1275780.0, 1275750.0,
           1275720.0, 1275690.0, 1275660.0, 1275630.0],
          dtype='float64', name='y', length=1109))

  • x
    PandasIndex

    PandasIndex(Index([410160.0, 410190.0, 410220.0, 410250.0, 410280.0, 410310.0, 410340.0,
           410370.0, 410400.0, 410430.0,
           ...
           442650.0, 442680.0, 442710.0, 442740.0, 442770.0, 442800.0, 442830.0,
           442860.0, 442890.0, 442920.0],
          dtype='float64', name='x', length=1093))

• Attributes: (2)

  crs :

  epsg:32628

  grid_mapping :

  spatial_ref

Evaluating observed vs. all modelled tide heights

The complex 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 the OTPS tidal model to model tide heights at a regular interval (every two hours by default) across the entire time period covered by the input satelliter timeseries 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.

Note: 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.

out_stats = tidal_stats(ds)

Setting tide modelling location from dataset centroid: -15.67, 11.69

79% of the full 2.91 m modelled tidal range is observed at this location.
The lowest 14% and highest 7% of tides are never observed.

Observed tides do not increase or decrease significantly over the ~31 year period.
All tides do not increase or decrease significantly over the ~31 year period.

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)

• observed_slope: slope of any relationship between observed tide heights and time

• all_slope: slope of any relationship between all modelled tide heights and time

• observed_pval: significance/p-value of any relationship between observed tide heights and time

• all_pval: significance/p-value of any relationship between all modelled tide heights and time

print(out_stats)

tidepost_lat        11.689
tidepost_lon       -15.674
observed_mean_m      0.086
all_mean_m           0.000
observed_min_m      -0.978
all_min_m           -1.393
observed_max_m       1.325
all_max_m            1.518
observed_range_m     2.303
all_range_m          2.911
spread               0.791
low_tide_offset      0.143
high_tide_offset     0.066
observed_slope       0.005
all_slope           -0.000
observed_pval        0.287
all_pval             0.370
dtype: float64

---

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Africa 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.

Compatible datacube version:

print(datacube.__version__)

1.8.15

Last Tested:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

'2023-08-14'