Introduction to the DEA Wetlands Insight Tool (WIT) 9422f4f51287457885febfd1ec6d545c

Note: ga_ls9c_ard_3 can be included once derivative Water Observations and Fractional Cover products have been generated

Background

Wetlands provide a wide range of ecosystem services including improving water quality, carbon sequestration, as well as providing habitat for fish, amphibians, reptiles and birds. Managing wetlands in Australia is challenging due to competing pressures for water availability and highly variable climatic settings. The DEA Wetlands Insight Tool has been developed to provide catchment managers, environmental water holders, and wetlands scientists a consistent historical baseline of wetlands dynamics from 1987 onwards.

The DEA Wetlands Insight Tool summarises how the amount of water, green vegetation, dry vegetation and bare soil varies over time within each wetland. It provides the user with the ability to compare how the wetland is behaving now with how it has behaved in the past. This allows users to identify how changes in water availability have affected the wetland. It achieves this by presenting a combined view of Water Observations from Space, Tasseled Cap Wetness and Fractional Cover measurements from the Landsat series of satellites, summarised as a stacked line plot to show how that wetland has changed over time.

Expected output

This notebook demonstrates the method for reproducing results from the Wetlands Insight Tool. At the end of the notebook you will display and save a WIT plot:

f2603f3596214b8fa21f65067479cf03

Applications

The product is designed to support wetland managers, catchment managers and environmental waterholders in understanding whether or not wetlands are changing over time. In instances where the wetlands are changing, the tool allows users to identify whether the changes are gradual, rapid, once-off or cyclical in nature. For example the response of wetlands to the following drivers can be assessed:

  • Changes in river flow volumes

  • Changes in flood frequency

  • Long term shifts in rainfall

  • Wet-season/Dry-season shifts in water availability

  • Invasive weeds

  • Environmental watering events

Care should be used when interpreting Wetlands Insight Tool results as increases/decreases in particular cover types can be associated with different processes. For example an increase in green cover could indicate canopy recovery of desirable wetland species or an increase in the amount of invasive weeds.

Publications

  • Dunn, B., Ai, E., Alger, M.J., Fanson, B., Fickas, K.C., Krause, C.E., Lymburner, L., Nanson, R., Papas, P., Ronan, M., Thomas, R.F., 2023. Wetlands Insight Tool: Characterising the Surface Water and Vegetation Cover Dynamics of Individual Wetlands Using Multidecadal Landsat Satellite Data. Wetlands 43, 37. Available: https://doi.org/10.1007/s13157-023-01682-7

  • Dunn, B., Lymburner, L., Newey, V., Hicks, A. and Carey, H., 2019. Developing a Tool for Wetland Characterization Using Fractional Cover, Tasseled Cap Wetness And Water Observations From Space. IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium, 2019, pp. 6095-6097. Available: https://doi.org/10.1109/IGARSS.2019.8897806.

  • Krause, C., Dunn, B., Bishop-Taylor, R., Adams, C., Burton, C., Alger, M., Chua, S., Phillips, C., Newey, V., Kouzoubov, K., Leith, A., Ayers, D., Hicks, A., DEA Notebooks contributors 2021. Digital Earth Australia notebooks and tools repository. Geoscience Australia, Canberra. https://doi.org/10.26186/145234

Please cite Dunn et. al 2023 and Krause et al 2021 when using the results of this notebook, and contact us as we’d love to hear about your use case! Licencing information is provided at the bottom of this notebook, and requires attribution.

WIT calculation

For each pixel, we calculate its WIT values based on the method of Dunn, et al., 2023. The calculation in this flowchart below shows the workflow implemented in this notebook.

455493d1cf584aef949f7402596e7966

Getting Started

Load packages

[1]:

import datetime
import geopandas as gpd
import itertools
import numpy as np
import pandas as pd
import xarray as xr

import datacube

import sys

sys.path.insert(1, "../Tools/")
import dea_tools.bandindices
import dea_tools.datahandling
from dea_tools.spatial import xr_rasterize
from dea_tools.dask import create_local_dask_cluster
import dea_tools.wetlands

# Create local dask cluster to improve data load time
client = create_local_dask_cluster(return_client=True)

Client

Client-7e5ff440-2d10-11ee-842b-42d3dbc549f9

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /user/bex.dunn@ga.gov.au/proxy/8787/status

Cluster Info

Connect to the datacube

Connect to the datacube so we can access DEA data. The app parameter is a unique name for the analysis which is based on the notebook file name.

[2]:

dc = datacube.Datacube(app="DEA_Wetlands_Insight_Tool")

Analysis parameters

  • poly: The polygon of your wetland of interest. This can be a shapefile or a GeoJSON. User will need to drag and drop the file into the sandbox and update the file path. Your polygon should be in EPSG:3577 GDA94 Australian Albers.

Load a shapefile that describes your wetland study area

[3]:

# This cell loads and plots the wetlands polygon. If no output, check your polygon
poly = gpd.read_file("../Supplementary_data/DEA_Wetlands_Insight_Tool/hird.shp")
poly.geometry[0]

[3]:
../../../_images/notebooks_DEA_products_DEA_Wetlands_Insight_Tool_12_0.svg
[4]:

# Specifying coordinate reference system of the polygon.
gpgon = datacube.utils.geometry.Geometry(poly.geometry[0], crs=poly.crs)

[5]:

# Select time period. Consistent data is available from 01-09-1987.
time = ("2014-01-01", "2016-01-01")
# time = ('1985-01-01', '2024-01-01')

Load Landsat, Water Observations, and Fractional Cover:

Load and mask Landsat data

Load Landsat 5, 7 and 8 data. Not including Landsat 7 SLC off period (31-05-2003 to 06-04-2022). If the following cell successfully loads data, it will generate output text, similar to:

Finding datasets
    ga_ls8c_ard_3
    ga_ls7e_ard_3 (ignoring SLC-off observations)
    ga_ls5t_ard_3
Applying fmask pixel quality/cloud mask
Loading 8 time steps

For further information on the functions used here, go to the function documentation.

[6]:

# Define which spectral bands are being used in the analysis
bands = [
    f"nbart_{band}" for band in ("blue", "green", "red", "nir", "swir_1", "swir_2")
]

# Load Landsat 5, 7 and 8 data. Not including Landsat 7 SLC off period (31-05-2003 to 06-04-2022).
ds = dea_tools.datahandling.load_ard(
    dc,
    products=["ga_ls8c_ard_3", "ga_ls7e_ard_3", "ga_ls5t_ard_3"],
    ls7_slc_off=False,
    measurements=bands,
    geopolygon=gpgon,
    output_crs="EPSG:3577",
    resolution=(-30, 30),
    resampling={"fmask": "nearest", "*": "bilinear"},
    time=time,
    group_by="solar_day",
    dask_chunks={},
)

# Load into memory using Dask
ds.load()

Finding datasets
    ga_ls8c_ard_3
    ga_ls7e_ard_3 (ignoring SLC-off observations)
    ga_ls5t_ard_3
Applying fmask pixel quality/cloud mask
Returning 89 time steps as a dask array

[6]:

<xarray.Dataset>
Dimensions:       (time: 89, y: 127, x: 61)
Coordinates:
  * time          (time) datetime64[ns] 2014-01-05T00:10:35.094128 ... 2015-1...
  * y             (y) float64 -3.968e+06 -3.968e+06 ... -3.972e+06 -3.972e+06
  * x             (x) float64 1.09e+06 1.09e+06 1.09e+06 ... 1.091e+06 1.091e+06
    spatial_ref   int32 3577
Data variables:
    nbart_blue    (time, y, x) float32 791.0 728.0 707.0 754.0 ... nan nan nan
    nbart_green   (time, y, x) float32 1.097e+03 1.025e+03 990.0 ... nan nan nan
    nbart_red     (time, y, x) float32 1.422e+03 1.339e+03 1.293e+03 ... nan nan
    nbart_nir     (time, y, x) float32 2.151e+03 2.009e+03 1.957e+03 ... nan nan
    nbart_swir_1  (time, y, x) float32 3.506e+03 3.378e+03 3.314e+03 ... nan nan
    nbart_swir_2  (time, y, x) float32 2.41e+03 2.291e+03 2.245e+03 ... nan nan
Attributes:
    crs:           EPSG:3577
    grid_mapping:  spatial_ref

Load Water Observations and Fractional Cover datasets

Load Fractional Cover and Water Observations data into the same spatial grid and resolution as the loaded Landsat dataset.

[7]:

ds_wo = dc.load(
    "ga_ls_wo_3", resampling="nearest", group_by="solar_day", like=ds, dask_chunks={}
)
ds_fc = dc.load(
    "ga_ls_fc_3", resampling="nearest", group_by="solar_day", like=ds, dask_chunks={}
)

# Load data into memory
ds_wo.load()
ds_fc.load()

[7]:

<xarray.Dataset>
Dimensions:      (time: 163, y: 127, x: 61)
Coordinates:
  * time         (time) datetime64[ns] 2014-01-05T00:10:35.094128 ... 2015-12...
  * y            (y) float64 -3.968e+06 -3.968e+06 ... -3.972e+06 -3.972e+06
  * x            (x) float64 1.09e+06 1.09e+06 1.09e+06 ... 1.091e+06 1.091e+06
    spatial_ref  int32 3577
Data variables:
    bs           (time, y, x) uint8 12 9 11 15 19 16 16 14 ... 19 7 2 2 2 2 2 2
    pv           (time, y, x) uint8 7 5 7 13 10 6 7 8 2 ... 14 14 11 9 6 9 6 8 8
    npv          (time, y, x) uint8 80 85 81 71 70 77 76 ... 89 92 89 92 90 90
    ue           (time, y, x) uint8 9 10 10 9 8 8 9 10 10 ... 12 8 7 6 6 5 5 5
Attributes:
    crs:           PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["G...
    grid_mapping:  spatial_ref

Locate and remove any observations which aren’t in all three datasets

[8]:

missing = set()
for t1, t2 in itertools.product(
    [ds_fc.time.values, ds_wo.time.values, ds.time.values], repeat=2
):
    missing_ = set(t1) - set(t2)
    missing |= missing_

ds_fc = ds_fc.sel(time=[t for t in ds_fc.time.values if t not in missing])
ds = ds.sel(time=[t for t in ds.time.values if t not in missing])
ds_wo = ds_wo.sel(time=[t for t in ds_wo.time.values if t not in missing])

Calculate Tasseled Cap Wetness from the Landsat data

We use Tasseled Cap Wetness (TCW) to identify areas in wetlands that are wet, but not identified by the WOfS algorithm as open water. This includes areas of mixed vegetation and water like in palustrine wetlands. TCW is the third component of a Tasseled Cap Index analysis and can be thresholded to classify wet pixels (Fisher et al., 2016). Tasseled Cap Index analysis is a linear principal component analysis of Landsat imagery with a Procrustes’ Rotation (Kauth and Thomas, 1976), producing three components roughly corresponding to brightness (TCB), greenness (TCG) and wetness (TCW; Roberts et al., 2018). The coefficients of Crist (1985) were used to calculate TCW.

TCW has values between −12,915 and 7032 when applied to DEA Landsat data. As TCW increases, the pixel becomes wetter, with open water having TCW values near and above zero (Pasquarella et al., 2016). We threshold TCW at −350, with values above this threshold used to characterise ‘wet’ pixels. We developed the TCW threshold to identify ‘wet’ pixels by comparing TCW to known inundation extents in the Macquarie Marshes, a large floodplain wetland in the Murray-Darling Basin of south-eastern Australia. We chose the threshold for the WIT workflow to minimise the absolute difference in extent between the observed inundated area and the area of thresholded TCW as described further here: https://link.springer.com/article/10.1007/s13157-023-01682-7#appendices

[9]:

tcw = dea_tools.bandindices.calculate_indices(
    ds, index="TCW", collection="ga_ls_3", normalise=False, drop=True, inplace=False
)

Dropping bands ['nbart_blue', 'nbart_green', 'nbart_red', 'nbart_nir', 'nbart_swir_1', 'nbart_swir_2']

Divide Fractional Cover by 100

Divide the FC values by 100 to keep them in \([0, 1]\) (ignoring for now that values can exceed 100 due to how FC is calculated) Keeps data types the same in the output raster.

[10]:

bs = ds_fc.bs / 100
pv = ds_fc.pv / 100
npv = ds_fc.npv / 100

Generate the WIT raster bands:

Create an empty dataset called output_rast and populate with values from input datasets.

[11]:

rast_names = ["pv", "npv", "bs", "wet", "water"]
output_rast = {n: xr.zeros_like(bs) for n in rast_names}

output_rast["bs"].values[:] = bs
output_rast["pv"].values[:] = pv
output_rast["npv"].values[:] = npv

Masking

  1. Mask pixels not in the polygon

  2. Use the Water Observations bit flags to mask noncontiguous data, low solar incidence angle, cloud, and water out of the wet category. We distinguish open water from wet as mutually exclusive categories.

Thresholding

We threshold TCW at −350, with values above this threshold used to characterise ‘wet’ pixels (Dunn, et al., 2023).

[12]:

# Rasterise the shapefile where poly is the vector data and pv is the xarray template
poly_raster = xr_rasterize(poly, pv) > 0

# Mask includes No data, Non contiguous data, Cloud shadow, Cloud, and water.
# See /notebooks/DEA_products/DEA_Water_Observations.html#Understanding-WOs-bit-flags for more detail.
mask = (ds_wo.water & 0b01100011) == 0
mask &= poly_raster

# Set open water to water present and classified as water as per Water Observations and bit flags
open_water = ds_wo.water & (1 << 7) > 0

# Set wet pixels where not masked and above threshold of -350
wet = tcw.where(mask).TCW > -350

[13]:

# Adding wet and water values to output raster

# TCW
output_rast["wet"].values[:] = wet.values.astype(float)
for name in rast_names[:3]:
    output_rast[name].values[wet.values] = 0

# WO
output_rast["water"].values[:] = open_water.values.astype(float)
for name in rast_names[:4]:
    output_rast[name].values[open_water.values] = 0

[14]:

# Masking again
ds_wit = xr.Dataset(output_rast).where(mask)

Mask entire observations where the polygon is more than 10% masked:

[15]:

# Calculate percentage missing
pc_missing = (~mask).where(poly_raster).mean(dim=["x", "y"])

ds_wit = ds_wit.where(pc_missing < 0.1)

The WIT results are now computed. All that’s left is to normalise them.

Normalise Fractional Cover Values in WIT result

Users are provided the option to normalise the WIT result. Normalisation only affects vegetation Fractional Cover values. We rescale Fractional Cover values by the area of the vegetation so that the total percentage adds to 1. This makes the graph more interpretable and is done as the Fractional Cover algorithm sometimes returns values which add to greater than 1.

We suggest normalisation when displaying the WIT plots but perhaps not if calculating population statistics.

[16]:

# Convert ds_wit: XArray.Dataset to polygon_base_df: pandas.DataFrame

polygon_base_df = pd.DataFrame()
polygon_base_df["date"] = ds_wit.time.values

for band in rast_names:
    polygon_base_df[band] = ds_wit[band].mean(dim=["x", "y"])

[17]:

polygon_base_df = dea_tools.wetlands.normalise_wit(polygon_base_df)

Create WIT Plot

Edit the polygon name in the cell below, if required. A png file will be created in this directory with the same name as the polygon, unless the png_name variable below is edited.

[18]:

polygon_name = "Hird Swamp"
png_name = polygon_name  # file will be png_name.png

[19]:

dea_tools.wetlands.display_wit_stack_with_df(polygon_base_df, polygon_name, png_name)
../../../_images/notebooks_DEA_products_DEA_Wetlands_Insight_Tool_39_0.png

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

[20]:

print(datacube.__version__)

1.8.13

Tags

**Tags**: :index:`calculate_indices`, :index:`sandbox compatible`, :index:`landsat 5`, :index:`landsat 7`, :index:`landsat 8`, :index:`time series`, :index:`Wetlands_Insight_Tool_(WIT)`, :index:`water observations`, :index:`fractional cover`, :index:`tasseled cap`