Introduction to DEA High and Low Tide Imagery 
Sign up to the DEA Sandbox to run this notebook interactively from a browser
Compatibility: Notebook currently compatible with both the
NCIandDEA SandboxenvironmentsProducts used: high_tide_comp_20p, low_tide_comp_20p
Background
Intertidal zones are those which are exposed to the air at low tide and underwater at high tide. These include sandy beaches, tidal flats, rocky shores and reefs.
Intertidal zones form critical habitats for a wide range of organisms, but are faced with increasing threats, including coastal erosion and a rise in sea levels.
The ever-changing nature of the tides makes it hard to systematically capture consistent imagery of the intertidal zone, particularly across large regions. Pressure is mounting on this zone from sea-level rise and anthropogenic sources such as land reclamation and aquaculture.
What this product offers
The DEA High and Low Tide Imagery product is a suite of cloud-free composite imagery of the intertidal zones at high and low tide around the Australian coast. It calculates the geometric median (geomedian) of the highest and lowest 20% of the observed tidal range in Digital Earth Australia (DEA)’s archive of Landsat satellite images.
To generate these composites, the archive of Landsat images has been paired with regional tidal modelling, generated by Oregon State Tidal Prediction software. This allows the archive to be sorted by tide height rather than date, so the intertidal zone can be visualised at any stage of the tide regime.
The data has been captured on a 25m spatial scale.
Applications
Mapping cover types within the intertidal zone
Visualising the full observed extent of the tidal range around the Australian continental coastline
Publications
Sagar, S., Phillips, C., Bala, B., Roberts, D., & Lymburner, L. (2018). Generating continental scale pixel-based surface reflectance composites in coastal regions with the use of a multi-resolution tidal model. Remote Sensing, 10(3), 480.
Description
This notebook will demonstrate how to load data from the DEA High and Low Tide product suite using the Digital Earth Australia datacube. Topics covered include:
Inspecting the products and measurements available in the datacube
Loading low and high tide data for a coastal location
Plotting low and high tide data in true and false colour
Converting low and high tide data to a remote sensing water index (NDWI), and use this to map wet pixels at low and high tide
Note: Visit the DEA High and Low Tide Imagery product documentation for detailed technical information including methods, quality, and data access.
Getting started
To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.
Load packages
Import Python packages that are used for the analysis.
[1]:
import datacube
import matplotlib.pyplot as plt
import sys
sys.path.insert(1, '../Tools/')
from dea_tools.plotting import rgb
/env/lib/python3.8/site-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
warnings.warn(
Connect to the datacube
Connect to the datacube so we can access DEA data.
[2]:
dc = datacube.Datacube(app='DEA_High_and_Low_Tide_Imagery')
Available products and measurements
List products
We can use datacube’s list_products functionality to inspect the DEA High and Low Tide Imagery products that are available in the datacube. The table below shows the product names that we will use to load the data, a brief description of the data, and its coordinate reference system, spatial resolution and spatial dimensions (i.e. 25 m pixels in the Australian Albers EPSG:3577 map projection).
[3]:
# List DEA High and Low Tide Imagery products available in DEA
dc_products = dc.list_products()
dc_products.loc[['high_tide_comp_20p', 'low_tide_comp_20p']]
[3]:
| name | description | license | default_crs | default_resolution | |
|---|---|---|---|---|---|
| name | |||||
| high_tide_comp_20p | high_tide_comp_20p | High tide 20 percentage composites 25m v. 2.0.0 | None | EPSG:3577 | (-25, 25) |
| low_tide_comp_20p | low_tide_comp_20p | Low tide 20 percentage composites 25m v. 2.0.0 | None | EPSG:3577 | (-25, 25) |
List measurements
We can further inspect the data available for each DEA High and Low Tide Imagery product using datacube’s list_measurements functionality. The table below lists each of the measurements available in the data, which represent the different satellite imagery bands (e.g. red, green, blue etc). These correspond to bands in the original Landsat satellite imagery.
The table also provides information about the measurement data types, units, nodata value and other technical information about each measurement.
[4]:
dc_measurements = dc.list_measurements()
dc_measurements.loc[['high_tide_comp_20p']]
[4]:
| name | dtype | units | nodata | aliases | flags_definition | spectral_definition | ||
|---|---|---|---|---|---|---|---|---|
| product | measurement | |||||||
| high_tide_comp_20p | blue | blue | float32 | metres | -999 | NaN | NaN | NaN |
| green | green | float32 | metres | -999 | NaN | NaN | NaN | |
| red | red | float32 | metres | -999 | NaN | NaN | NaN | |
| nir | nir | float32 | metres | -999 | NaN | NaN | NaN | |
| swir1 | swir1 | float32 | metres | -999 | NaN | NaN | NaN | |
| swir2 | swir2 | float32 | metres | -999 | NaN | NaN | NaN |
Loading data
Now that we know what products and measurements are available for the products, we can load data from the datacube for an example location:
[5]:
# Set up a region to load low and high tide data
query = {
'x': (-1041555.354316434, -1025555.354316434),
'y': (-1992666.5901746228, -1976666.5901746228),
'crs': 'EPSG:3577'
}
# Load low and high tide data from the datacube
low_tide_ds = dc.load(
product='low_tide_comp_20p',
measurements=['red', 'green', 'blue', 'nir', 'swir1', 'swir2'],
**query)
high_tide_ds = dc.load(
product='high_tide_comp_20p',
measurements=['red', 'green', 'blue', 'nir', 'swir1', 'swir2'],
**query)
We can now view the data that we loaded. The satellite bands listed under Data variables should match the measurements displayed in the previous List measurements step.
[6]:
low_tide_ds
[6]:
<xarray.Dataset>
Dimensions: (time: 1, x: 641, y: 641)
Coordinates:
* time (time) datetime64[ns] 2008-06-01
* y (y) float64 -1.977e+06 -1.977e+06 ... -1.993e+06 -1.993e+06
* x (x) float64 -1.042e+06 -1.042e+06 ... -1.026e+06 -1.026e+06
spatial_ref int32 3577
Data variables:
red (time, y, x) float32 0.04964104 0.049942903 ... 0.11672047
green (time, y, x) float32 0.13809542 0.13810375 ... 0.07829063
blue (time, y, x) float32 0.142985 0.14270617 ... 0.051898982
nir (time, y, x) float32 0.019549698 0.019038271 ... 0.21109466
swir1 (time, y, x) float32 0.014030585 0.013525395 ... 0.3071515
swir2 (time, y, x) float32 0.012055442 0.011902479 ... 0.22296384
Attributes:
crs: EPSG:3577
grid_mapping: spatial_ref- time: 1
- x: 641
- y: 641
- time(time)datetime64[ns]2008-06-01
- units :
- seconds since 1970-01-01 00:00:00
array(['2008-06-01T00:00:00.000000000'], dtype='datetime64[ns]')
- y(y)float64-1.977e+06 ... -1.993e+06
- units :
- metre
- resolution :
- -25.0
- crs :
- EPSG:3577
array([-1976662.5, -1976687.5, -1976712.5, ..., -1992612.5, -1992637.5, -1992662.5]) - x(x)float64-1.042e+06 ... -1.026e+06
- units :
- metre
- resolution :
- 25.0
- crs :
- EPSG:3577
array([-1041562.5, -1041537.5, -1041512.5, ..., -1025612.5, -1025587.5, -1025562.5]) - spatial_ref()int323577
- spatial_ref :
- PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]
- grid_mapping_name :
- albers_conical_equal_area
array(3577, dtype=int32)
- red(time, y, x)float320.04964104 ... 0.11672047
- units :
- metres
- nodata :
- -999
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[0.04964104, 0.0499429 , 0.04796788, ..., 0.0698581 , 0.07028087, 0.07129946], [0.04977173, 0.05028106, 0.04845902, ..., 0.07083403, 0.07082441, 0.07164738], [0.05028646, 0.04904343, 0.04848236, ..., 0.07220325, 0.07141209, 0.07238168], ..., [0.03457597, 0.03418804, 0.03401429, ..., 0.09348667, 0.10701565, 0.12069853], [0.03427565, 0.03422121, 0.03421114, ..., 0.13420059, 0.14639324, 0.14832586], [0.0340373 , 0.03425661, 0.0341205 , ..., 0.14411034, 0.12841429, 0.11672047]]], dtype=float32) - green(time, y, x)float320.13809542 ... 0.07829063
- units :
- metres
- nodata :
- -999
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[0.13809542, 0.13810375, 0.13611284, ..., 0.16628462, 0.16715015, 0.16753429], [0.13907972, 0.13858926, 0.13676979, ..., 0.16755676, 0.16750498, 0.16788842], [0.13941373, 0.13773836, 0.13758245, ..., 0.17009951, 0.16770892, 0.16861743], ..., [0.09712204, 0.09715024, 0.09738219, ..., 0.07293889, 0.07549716, 0.07789179], [0.09719902, 0.09735054, 0.09746309, ..., 0.0832534 , 0.08605166, 0.08651371], [0.09688519, 0.09715211, 0.09706473, ..., 0.08467433, 0.0809047 , 0.07829063]]], dtype=float32) - blue(time, y, x)float320.142985 0.14270617 ... 0.051898982
- units :
- metres
- nodata :
- -999
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[0.142985 , 0.14270617, 0.14022838, ..., 0.15646422, 0.15658316, 0.15758859], [0.1431268 , 0.14306079, 0.14109887, ..., 0.15757109, 0.15733393, 0.15787108], [0.14304389, 0.14213254, 0.14189798, ..., 0.15845639, 0.15719387, 0.15804994], ..., [0.10451855, 0.10477079, 0.10463735, ..., 0.0517672 , 0.05075943, 0.05098888], [0.1042261 , 0.10487015, 0.1045967 , ..., 0.05206816, 0.05304268, 0.05325721], [0.10414909, 0.10487483, 0.1045597 , ..., 0.0537968 , 0.05303383, 0.05189898]]], dtype=float32) - nir(time, y, x)float320.019549698 ... 0.21109466
- units :
- metres
- nodata :
- -999
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[0.0195497 , 0.01903827, 0.01729411, ..., 0.01804621, 0.01742984, 0.0181765 ], [0.01860858, 0.01868952, 0.01735434, ..., 0.01726199, 0.01792473, 0.01808378], [0.01894449, 0.01818096, 0.01768152, ..., 0.0176274 , 0.01755256, 0.0177814 ], ..., [0.01788817, 0.0178327 , 0.0176415 , ..., 0.20651203, 0.2181411 , 0.23163185], [0.01711432, 0.01721582, 0.01721008, ..., 0.24189241, 0.25608 , 0.25172725], [0.01676285, 0.01663971, 0.01671357, ..., 0.2434728 , 0.22394842, 0.21109466]]], dtype=float32) - swir1(time, y, x)float320.014030585 ... 0.3071515
- units :
- metres
- nodata :
- -999
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[0.01403059, 0.01352539, 0.01283986, ..., 0.01007457, 0.00948367, 0.00991209], [0.01376353, 0.01397293, 0.01270587, ..., 0.00991901, 0.00992017, 0.00984262], [0.0141523 , 0.01298638, 0.01300007, ..., 0.01000375, 0.00972117, 0.00956158], ..., [0.01237123, 0.01203357, 0.01193276, ..., 0.26145387, 0.27997833, 0.3093873 ], [0.01203714, 0.01188824, 0.01172496, ..., 0.32827324, 0.34819764, 0.35292023], [0.01148987, 0.011447 , 0.01168679, ..., 0.34250706, 0.32114953, 0.3071515 ]]], dtype=float32) - swir2(time, y, x)float320.012055442 ... 0.22296384
- units :
- metres
- nodata :
- -999
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[0.01205544, 0.01190248, 0.01067546, ..., 0.007613 , 0.00745695, 0.00867213], [0.01175852, 0.01178023, 0.01051343, ..., 0.00758195, 0.00799356, 0.00779425], [0.01145432, 0.01081816, 0.01065391, ..., 0.00869842, 0.00774264, 0.00757237], ..., [0.01058437, 0.00979857, 0.01010403, ..., 0.16180673, 0.18422522, 0.21883045], [0.01068483, 0.00994988, 0.00974412, ..., 0.24057686, 0.2628302 , 0.2671504 ], [0.00984477, 0.00934011, 0.00959785, ..., 0.26111013, 0.23739797, 0.22296384]]], dtype=float32)
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
Plotting data
We can plot the data using the rgb function to view the loaded low and high tide imagery. By plotting the ['red', 'green', 'blue'] bands, we can view the imagery in true colour:
[7]:
# Plot low and high tide imagery side-by-side
fig, axes = plt.subplots(ncols=2, figsize=(14, 6))
rgb(low_tide_ds, bands=['red', 'green', 'blue'], ax=axes[0])
rgb(high_tide_ds, bands=['red', 'green', 'blue'], ax=axes[1])
axes[0].title.set_text('Low tide true-colour imagery')
axes[1].title.set_text('High tide true-colour imagery')
By plotting the ['swir1', 'nir', 'green'] bands, we can view the imagery in false colour. This view emphasises the changing position of the land-water boundary between low and high tide.
Note: For more information about plotting satellite imagery in true and false colour, refer to the Introduction to Plotting notebook.
[8]:
# Plot low and high tide imagery side-by-side
fig, axes = plt.subplots(ncols=2, figsize=(14, 6))
rgb(low_tide_ds, bands=['swir1', 'nir', 'green'], ax=axes[0])
rgb(high_tide_ds, bands=['swir1', 'nir', 'green'], ax=axes[1])
axes[0].title.set_text('Low tide false-colour imagery')
axes[1].title.set_text('High tide false-colour imagery')
Example application: computing remote sensing indices on low and high tide imagery
Because the DEA High and Low Tide Imagery data was produced using a geomedian approach that preserves the band relationships within the modelled spectra at each pixel, the surface reflectance values can be used for remote sensing applications (for example, image classification, remote sensing index calculation, habitat mapping) in exactly the same way as spectral bands from individual Landsat observations.
For example, we can use the data to compute simple remote sensing indices such as the Normalized Difference Water Index (NDWI). This index will have high values where a pixel is likely to be open water (e.g. NDWI > 0, or blue colours below):
[9]:
# Compute NDWI using the formula (green - nir) / (green + nir)
high_tide_ndwi = (high_tide_ds.green - high_tide_ds.nir) / \
(high_tide_ds.green + high_tide_ds.nir)
low_tide_ndwi = (low_tide_ds.green - low_tide_ds.nir) / \
(low_tide_ds.green + low_tide_ds.nir)
# Plot high and low tide NDWI side-by-side
fig, axes = plt.subplots(ncols=2, figsize=(14, 6))
low_tide_ndwi.plot(ax=axes[0], vmin=-0.5, vmax=0.5, cmap='RdYlBu')
high_tide_ndwi.plot(ax=axes[1], vmin=-0.5, vmax=0.5, cmap='RdYlBu')
axes[0].title.set_text('Low tide NDWI')
axes[1].title.set_text('High tide NDWI')
axes[1].get_yaxis().set_visible(False)
fig.tight_layout()
A possible application of this may be to map the distribution of water at low and high tide using a NDWI threshold:
[10]:
# Identify water pixels using an example threshold value of 0
high_tide_water = high_tide_ndwi > 0
low_tide_water = low_tide_ndwi > 0
# Plot distribution of water at high and low tide
fig, axes = plt.subplots(ncols=2, figsize=(14, 6))
low_tide_water.plot(ax=axes[0])
high_tide_water.plot(ax=axes[1])
axes[0].title.set_text('Low tide water')
axes[1].title.set_text('High tide water')
axes[1].get_yaxis().set_visible(False)
fig.tight_layout()
Additional information
License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.
Contact: If you need assistance, please post a question on the Open Data Cube Discord chat or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on
GitHub.
Last modified: December 2023
Compatible datacube version:
[12]:
print(datacube.__version__)
1.8.5
Tags
Tags: NCI compatible, sandbox compatible, DEA products, high_tide_comp_20p, low_tide_comp_20p, rgb, NDWI, intertidal, coastal