Monitoring coastal erosion along Africa’s coastline

Keywords: data used; landsat 8, water; tide modelling, water; waterline extraction, band index; MNDWI; data methods; compositing

Background

Over 40% of the world’s population lives within 100 km of the coastline. However, coastal environments are constantly changing, with erosion and coastal change presenting a major challenge to valuable coastal infrastructure and important ecological habitats. Up-to-date data on coastal change and erosion is essential for coastal managers to be able to identify and minimise the impacts of coastal change and erosion.

Monitoring coastlines and rivers using field surveys can be challenging and hazardous, particularly at regional or national scale. Aerial photography and LiDAR can be used to monitor coastal change, but this is often expensive and requires many repeated flights over the same areas of coastline to build up an accurate history of how the coastline has changed across time.

Digital Earth Africa use case

Imagery from satellites such as the NASA/USGS Landsat program is available for free for the entire planet, making satellite imagery a powerful and cost-effective tool for monitoring coastlines and rivers at regional or national scale.

By identifying and extracting the precise boundary between water and land based on satellite data, it is possible to extract accurate shorelines that can be compared across time to reveal hotspots of erosion and coastal change.

The usefulness of satellite imagery in the coastal zone can be affected by the presence of clouds, sun-glint over water, poor water quality (e.g. sediment) and the influence of tides. The effect of these factors can be reduced by combining individual noisy images into cleaner “summary” or composite layers, and filtering the data to focus only on images taken at certain tidal conditions (e.g. high tide).

Description

In this example, we combine data from the Landsat 5, 7 and 8 satellites with image compositing and tide filtering techniques to accurately map shorelines across time, and identify areas that have changed significantly between 1987 and 2019. The worked example demonstrates how to:

  1. Load in a cloud-free Landsat time series

  2. Compute a water index (MNDWI)

  3. Filter images by tide height

  4. Create “summary” or composite images for given time periods

  5. Extract and visualise shorelines across time

Getting started

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

After finishing the analysis, return to the “Analysis parameters” cell, modify some values (e.g. choose a different location or time period to analyse) and re-run the analysis. There are additional instructions on modifying the notebook at the end.

Load packages

Load key Python packages and supporting functions for the analysis.

[1]:
%matplotlib inline

import datacube
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os

from deafrica_tools.datahandling import load_ard, mostcommon_crs
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.coastal import tidal_tag
from deafrica_tools.spatial import subpixel_contours
from deafrica_tools.plotting import display_map, rgb, map_shapefile
from deafrica_tools.dask import create_local_dask_cluster
/env/lib/python3.6/site-packages/geopandas/_compat.py:88: UserWarning: The Shapely GEOS version (3.7.2-CAPI-1.11.0 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow.
  shapely_geos_version, geos_capi_version_string

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.

[2]:
create_local_dask_cluster()

Client

Cluster

  • Workers: 1
  • Cores: 2
  • Memory: 13.11 GB

Connect to the datacube

Activate the datacube database, which provides functionality for loading and displaying stored Earth observation data.

[3]:
dc = datacube.Datacube(app="Coastal_erosion")

Analysis parameters

The following cell set important parameters for the analysis:

  • lat: The central latitude to analyse (e.g. 14.283).

  • lon: The central longitude to analyse (e.g. -16.921).

  • buffer: The number of square degrees to load around the central latitude and longitude. For reasonable loading times, set this as 0.1 or lower.

  • time_range: The date range to analyse (e.g. ('2013', '2020'))

  • time_step: This parameter allows us to choose the length of the time periods we want to compare: e.g. shorelines for each year, or shorelines for each six months etc. 1Y will generate one coastline for every year in the dataset; 6M will produce a coastline for every six months, etc.

  • tide_range: The minimum and maximum proportion of the tidal range to include in the analysis. For example, tide_range = (0.50, 1.00) will select all satellite images taken when the tide was greater than the median (i.e. 50th percentile) of all tide heights and less than the maximum (i.e. 100th percentile) of all tide heights. This allows you to seperate the effect of erosion from the influence of tides by producing shorelines for specific tidal conditions (e.g. low tide, average tide, high tide shorelines etc).

If running the notebook for the first time, keep the default settings below. This will demonstrate how the analysis works and provide meaningful results. The example explores coastal change in Ponto, Senegal.

To run the notebook for a different area, make sure Landsat 5, 7 and 8 data is available for the new location, which you can check at the DE Africa Explorer (use the drop-down menu to view all Landsat products).

To ensure that the tidal modelling part of this analysis works correctly, please make sure the centre of the study area is located over water when setting lat_range and lon_range.

[4]:
# Define the area of interest
lat = 14.283
lon = -16.921
buffer = 0.015

# Combine central lat,lon with buffer to get area of interest
lat_range = (lat-buffer, lat+buffer)
lon_range = (lon-buffer, lon+buffer)

# Set the range of dates for the analysis, time step and tide range
time_range = ('2013', '2020')
time_step = '1Y'
tide_range = (0.50, 1.00)

View the selected location

The next cell will display the selected area on an interactive map. Feel free to zoom in and out to get a better understanding of the area you’ll be analysing. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point.

[5]:
display_map(x=lon_range, y=lat_range)
[5]:

Load cloud-masked Landsat data

The first step in this analysis is to load in Landsat data for the lat_range, lon_range and time_range we provided above. The code below uses the load_ard function to load in data from the Landsat 5, 7 and 8 satellites for the area and time specified. For more infmation, see the Using load_ard notebook. The function will also automatically mask out clouds from the dataset, allowing us to focus on pixels that contain useful data:

[6]:
# Create the 'query' dictionary object, which contains the longitudes,
# latitudes and time provided above
query = {
    'x': lon_range,
    'y': lat_range,
    'time': time_range,
    'measurements': ['red', 'green', 'blue', 'swir_1'],
    'resolution': (-30, 30),
}

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

# Load available data Landsat 8
landsat_ds = load_ard(dc=dc,
                      products=['ls8_sr'],
                      output_crs=output_crs,
                      align=(15, 15),
                      dask_chunks={'time': 1},
                      group_by='solar_day',
                      **query)

Using pixel quality parameters for USGS Collection 2
Finding datasets
    ls8_sr
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Returning 169 time steps as a dask array

Once the load is complete, examine the data by printing it in the next cell. The Dimensions argument revels the number of time steps in the data set, as well as the number of pixels in the x (longitude) and y (latitude) dimensions.

[7]:
print(landsat_ds)
<xarray.Dataset>
Dimensions:      (time: 169, x: 109, y: 113)
Coordinates:
  * time         (time) datetime64[ns] 2013-03-26T11:26:35.081094 ... 2020-12...
  * y            (y) float64 1.582e+06 1.582e+06 ... 1.578e+06 1.578e+06
  * x            (x) float64 2.912e+05 2.912e+05 ... 2.944e+05 2.944e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    swir_1       (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

Plot example timestep in true colour

To visualise the data, use the pre-loaded rgb utility function to plot a true colour image for a given time-step. White areas indicate where clouds or other invalid pixels in the image have been masked.

Change the value for timestep and re-run the cell to plot a different timestep (timesteps are numbered from 0 to n_time - 1 where n_time is the total number of timesteps; see the time listing under the Dimensions category in the dataset print-out above).

[8]:
# Set the timesteps to visualise
timestep = 40

# Generate RGB plots at each timestep
rgb(landsat_ds,
    index=timestep,
    percentile_stretch=(0.05, 0.95))

../../../_images/sandbox_notebooks_Real_world_examples_Coastal_erosion_19_0.png

Compute Modified Normalised Difference Water Index

To extract shoreline locations, we need to be able to seperate water from land in our study area. To do this, we can use our Landsat data to calculate a water index called the Modified Normalised Difference Water Index, or MNDWI. This index uses the ratio of green and mid-infrared radiation to identify the presence of water (Xu 2006). The formula is:

\[\begin{aligned} \text{MNDWI} &= \frac{(\text{Green} - \text{MIR})}{(\text{Green} + \text{MIR})} \end{aligned}\]

where Green is the green band and MIR is the mid-infrared band. For Landsat, we can use the Short-wave Infrared (SWIR) 1 band as our measure for MIR.

When it comes to interpreting the index, High values (greater than 0, blue colours) typically represent water pixels, while low values (less than 0, red colours) represent land. You can use the cell below to calculate and plot one of the images after calculating the index.

[9]:
# Calculate the water index
landsat_ds = calculate_indices(landsat_ds, index='MNDWI',
                               collection='c1')

# Plot the resulting image for the same timestep selected above
landsat_ds.MNDWI.isel(time=timestep).plot(cmap='RdBu',
                                          size=6,
                                          vmin=-0.5,
                                          vmax=0.5)
plt.show()

../../../_images/sandbox_notebooks_Real_world_examples_Coastal_erosion_21_0.png

How does the plot of the index compare to the optical image from earlier? Was there water or land anywhere you weren’t expecting?

Model tide heights

The location of the shoreline can vary greatly from low to high tide. In the code below, we aim to reduce the effect of tides by modelling tide height data, and keeping only the satellite images that were taken at specific tidal conditions. For example, if tide_range = (0.50, 1.00), we are telling the analysis to focus only on satellite images taken when the tide was between median (50th percentile) and maximum conditions (100th percentile).

The tidal_tag function below uses the OTPS TPXO8 tidal model to calculate the height of the tide at the exact moment each satellite image in our dataset was taken, and adds this as a new tide_height attribute in our dataset.

Important note: this function can only model tides correctly if the centre of your study area is located over water. If this isn’t the case, you can specify a custom tide modelling location by passing a coordinate to tidepost_lat and tidepost_lon (e.g. tidepost_lat=14.283, tidepost_lon=-16.921).

[10]:
# Calculate tides for each timestep in the satellite dataset
landsat_ds = tidal_tag(ds=landsat_ds, tidepost_lat=None, tidepost_lon=None)

# Print the output dataset with new `tide_height` variable
print(landsat_ds)

Setting tide modelling location from dataset centroid: -16.92, 14.28
<xarray.Dataset>
Dimensions:      (time: 169, x: 109, y: 113)
Coordinates:
  * time         (time) datetime64[ns] 2013-03-26T11:26:35.081094 ... 2020-12...
  * y            (y) float64 1.582e+06 1.582e+06 ... 1.578e+06 1.578e+06
  * x            (x) float64 2.912e+05 2.912e+05 ... 2.944e+05 2.944e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    swir_1       (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    MNDWI        (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    tide_height  (time) float64 0.008 -0.157 0.538 0.58 ... 0.048 0.057 -0.27
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

Now that we have modelled tide heights, we can plot them to visualise the range of tide that was captured by Landsat across our time series. In the plot below, red dashed lines also show the subset of the tidal range we selected using the tide_range parameter. The plot should make it clear that limiting the range of the tides for the analysis should give you more consistent results. A large variance in the tide height could obscure your results, so consistency is critical as you want to compare the change in the shoreline from year to year.

[11]:
# Calculate the min and max tide heights to include based on the % range
min_tide, max_tide = landsat_ds.tide_height.quantile(tide_range)

# Plot the resulting tide heights for each Landsat image:
landsat_ds.tide_height.plot()
plt.axhline(min_tide, c='red', linestyle='--')
plt.axhline(max_tide, c='red', linestyle='--')
plt.show()

../../../_images/sandbox_notebooks_Real_world_examples_Coastal_erosion_26_0.png

Filter Landsat images by tide height

Here we take the Landsat dataset and only keep the images with tide heights we want to analyse (i.e. tides within the heights given by tide_range). This will result in a smaller number of images (e.g. ~70 images compared to ~140):

[12]:
# Keep timesteps larger than the min tide, and smaller than the max tide
landsat_filtered = landsat_ds.sel(time=(landsat_ds.tide_height > min_tide) &
                                       (landsat_ds.tide_height <= max_tide))
print(landsat_filtered)
<xarray.Dataset>
Dimensions:      (time: 84, x: 109, y: 113)
Coordinates:
  * time         (time) datetime64[ns] 2013-05-26T11:29:31.343448 ... 2020-10...
  * y            (y) float64 1.582e+06 1.582e+06 ... 1.578e+06 1.578e+06
  * x            (x) float64 2.912e+05 2.912e+05 ... 2.944e+05 2.944e+05
    spatial_ref  int32 32628
Data variables:
    red          (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    green        (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    swir_1       (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    MNDWI        (time, y, x) float32 dask.array<chunksize=(1, 113, 109), meta=np.ndarray>
    tide_height  (time) float64 0.538 0.58 0.805 0.578 ... 0.35 0.537 0.659
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

Combine observations into noise-free summary images

Individual remote sensing images can be affected by noisy data, including clouds, sunglint and poor water quality conditions (e.g. sediment). To produce cleaner images that can be compared more easily across time, we can create ‘summary’ images or composites that combine multiple images into one image to reveal the median or ‘typical’ appearance of the landscape for a certain time period. In this case, we use the median as the summary statistic because it prevents strong outliers (like masked cloud values) from skewing the data, which would not be the case if we were to use the mean.

In the code below, we take the time series of images and combine them into single images for each time_step. For example, if time_step = '2Y', the code will produce one new image for each two-year period in the dataset. This step can take several minutes to load if the study area is large.

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.

[13]:
# Combine into summary images by `time_step`
landsat_summaries = (landsat_filtered.MNDWI
                     .compute()
                     .resample(time=time_step, closed='left')
                     .median('time'))

# Plot the output summary images
landsat_summaries.plot(col='time',
                       cmap='RdBu',
                       col_wrap=4,
                       vmin=-0.5,
                       vmax=0.5)
plt.show()

../../../_images/sandbox_notebooks_Real_world_examples_Coastal_erosion_30_0.png

Extract shorelines from imagery

We now want to extract an accurate shoreline for each each of the summary images above (e.g. 2014, 2015 etc. summaries). The code below identifies the boundary between land and water by tracing a line along pixels with a water index value of 0 (halfway between land and water water index values). It returns a shapefile with one line for each time step:

[14]:
# Set up attributes to assign to each waterline
attribute_data = {'time': [str(i)[0:10] for i in landsat_summaries.time.values]}
attribute_dtypes = {'time': 'str'}

# Extract waterline contours for the '0' water index threshold:
contour_gdf = subpixel_contours(da=landsat_summaries,
                                 z_values=0,
                                 crs=landsat_ds.geobox.crs,
                                 affine=landsat_ds.geobox.transform,
                                 output_path=f'output_waterlines.geojson',
                                 min_vertices=50)

# Plot output shapefile over the first MNDWI layer in the time series
landsat_summaries.isel(time=0).plot(size=12,
                                    cmap='Greys',
                                    add_colorbar=False)

# add the contours to the plot
contour_gdf.plot(ax=plt.gca(),
                 column='time',
                 cmap='YlOrRd',
                 legend=True,
                 legend_kwds={'loc': 'lower right'})
plt.show()

Operating in single z-value, multiple arrays mode
Writing contours to output_waterlines.geojson
/env/lib/python3.6/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))
../../../_images/sandbox_notebooks_Real_world_examples_Coastal_erosion_32_2.png

The above plot is a basic visualisation of the contours returned by the contour_extract function. Given we now have the shapefile, we can use a more complex function to make an interactive plot for viewing the change in shoreline over time below.

Plot interactive map of output shorelines coloured by time

The next cell provides an interactive map with an overlay of the shorelines identified in the previous cell. Run it to view the map (this step can take several minutes to load if the study area is large).

Zoom in to the map below to explore the resulting set of shorelines. Older shorelines are coloured in yellow, and more recent shorelines in red. Hover over the lines to see the time period for each shoreline printed above the map. Using this data, we can easily identify areas of coastline or rivers that have changed significantly over time, or areas that have remained stable over the entire time period.

[15]:
map_shapefile(gdf=contour_gdf, attribute='time', hover_col='time')

Drawing conclusions

Here are some questions to think about: * What can you conclude about the change in the shoreline? * Which sections of the shoreline have seen the most change? * Is the change consistent with erosion? * What other information might you need to draw additional conclusions about the cause of the change?

Next steps

When you are done, return to the “Set up analysis” cell, modify some values (e.g. time_range, tide_range, time_step or lat/lon) and rerun the analysis.

If you’re going to change the location, you’ll need to make sure Landsat data is available for the new location, which you can check at the DE Africa Explorer (use the drop-down menu to view all Landsat products).


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:

[16]:
print(datacube.__version__)
1.8.4.dev52+g07bc51a5

Last Tested:

[17]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[17]:
'2021-05-19'