Plotting animations of land cover 7dad63b9c97041b2b0dc4b75d4e050aa

Background

Land cover is the physical surface of the Earth, including trees, shrubs, grasses, soils, exposed rocks, water bodies, plantations, crops and built structures. Digital Earth Australia Land Cover (DEA Land Cover) is a continental dataset that maps annual land cover classifications for Australia from 1988 to the present. Detailed information about DEA Land Cover can be found in the DEA Land Cover notebook and on the DEA Land Cover product details page.

Description

This notebook introduces the lc_animation() function and demonstrates how it can be used to visualise and communicate change over time. Topics covered include:

  1. Loading a time series of DEA Land Cover data as an xarray dataset.

  2. Plotting land cover as an animation.

  3. Plotting a stacked line plot next to a land cover animation.

  4. Modifying the way land cover is plotted as an animation.

Getting started

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

Load packages

[1]:

%matplotlib inline

import math
import os
import sys

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

sys.path.insert(1, os.path.abspath("../Tools"))
from dea_tools.landcover import lc_animation
from dea_tools.plotting import display_map

Connect to the datacube

Connect to the datacube so we can access DEA data.

[2]:

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

Create query and load a time series of land cover data

In order to generate animations of DEA Land Cover, we first need to load a time series of data for an area. As an example, let’s load both the base 6-class classification (Level 3) layer and the full classification (Level 4) layer for Lake Wallawalla. This waterbody is located in the Murray River Valley in north-western Victoria, and goes through cycles of drying and filling.

[3]:

# Lake Wallawalla
point_x, point_y = (-34.1795, 141.1937)

lat = (point_x - 0.05, point_x + 0.05)
lon = (point_y - 0.05, point_y + 0.05)

# Display area on map
display_map(x=lon, y=lat)

[3]:
Make this Notebook Trusted to load map: File -> Trust Notebook
[4]:

# Build query and load data
query = {
    "y": lat,
    "x": lon,
    "time": ("1988", "2021"),
}

# Load DEA Land Cover data from the datacube
land_cover_data = dc.load(
    product="ga_ls_landcover_class_cyear_2",
    output_crs="EPSG:3577",
    measurements=["level3", "level4", "waterper_wat_cat_l4d_au"],
    resolution=(-25, 25),
    **query
)

[5]:

# See what we have loaded
land_cover_data

[5]:

<xarray.Dataset>
Dimensions:                  (time: 33, y: 471, x: 399)
Coordinates:
  * time                     (time) datetime64[ns] 1988-01-01 ... 2020-01-01
  * y                        (y) float64 -3.755e+06 -3.755e+06 ... -3.767e+06
  * x                        (x) float64 8.378e+05 8.378e+05 ... 8.478e+05
    spatial_ref              int32 3577
Data variables:
    level3                   (time, y, x) uint8 112 112 111 112 ... 111 111 216
    level4                   (time, y, x) uint8 34 34 16 34 34 ... 96 16 16 96
    waterper_wat_cat_l4d_au  (time, y, x) uint8 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
Attributes:
    crs:           EPSG:3577
    grid_mapping:  spatial_ref

Plot land cover animations

Plot default animation of Level 3 data

lc_animation() has a number of options for controlling the inclusion of colour bars and tick marks and the ability to add an animated stacked line plot beside the animated map. To begin, let’s plot the Level 3 layer with all default options. This will plot a single animated map with tick marks and no colour bar. A copy of the animation will be saved to file as default_animation.gif. This name can be changed by providing an alternative as the variable file_name as will be done in following examples.

[6]:

# Select Level 3 data
level_3 = land_cover_data["level3"]

# Generate plot
lc_animation(level_3)

[6]:

<IPython.core.display.Image object>

Plot animation with a stacked line plot

Adding the optional setting stacked_plot=True will generate a stacked line plot that is synchronised with the animated map. The stacked line plot shows what percentage of the map is taken up by each class in each year. For this example we will use the full classification, Level 4 layer.

[7]:

# Select Level 4
level_4 = land_cover_data["level4"]

# Set file name
file_name = "level_4_with_stacked_plot"

# Generate plot
lc_animation(level_4, file_name=file_name, stacked_plot=True)

[7]:

<IPython.core.display.Image object>

Modifying lc_animation()

lc_animation() contains many optional parameters allowing customisation of the appearance of the land cover animations. Below are examples of how some of them operate. For the full details of the parameters please refer to the landcover.py file.

Plot animation with colour bar

A colour bar can be added to the stand alone animation by including the optional setting colour_bar=True.
> Note: The colour bar does not work with the stacked line plot option.

Let’s generate an example using the water persistence layer.

[8]:

# Select water persistance layer
water_per = land_cover_data["waterper_wat_cat_l4d_au"]

# Set file name
file_name = "water_persistence_animation"

# Generate plot
lc_animation(water_per, file_name=file_name, colour_bar=True)

[8]:

<IPython.core.display.Image object>

Plot with no axes ticks

If you want to make an animation for communication purposes, perhaps you don’t want to have any tick marks on the axes. In this case you include label_ax=False.

[9]:

# Select Level 4
level_4 = land_cover_data["level4"]

# Set file name
file_name = "level_4_no_ticks_example"

# Generate plot
lc_animation(level_4, file_name=file_name, label_ax=False)

[9]:

<IPython.core.display.Image object>

Why is my animation producing an error or plotting no data?

lc_animation() uses the name of the xarray.DataArray to identify what colour map should be used to create the animation. If, in your analysis, you have created a new xarray.DataArray object you will need to ensure that you either: 1. Name the xarray.DataArray after the layer being plotted. 2. Use the measurement setting for lc_animation() to specify the layer.

Plotting a layer with the incorrect colour map will result in an error or the data being plotted as nodata.

In the next cell we do some analysis on the Level 3 classification to isolate one class. When looking at this new xarray.DataArray, we can see that the name is now different ("Bare_only" compared to "level3").

[10]:

# Select Level 3
level_3 = land_cover_data["level3"]

# Do some data manipulation then create a new Xarray.DataArray
# with a new layer name "Bare_only"
data = np.where(level_3 == 216, level_3, 0)
xr_bare = xr.DataArray(
    data=data,
    coords=level_3.coords,
    dims=level_3.dims,
    name="Bare_only",
    attrs=None,
)

xr_bare

[10]:

<xarray.DataArray 'Bare_only' (time: 33, y: 471, x: 399)>
array([[[  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        ...,
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0]],

       [[  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        ...,
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0]],

       [[  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0],
        ...,
...
        ...,
        [  0,   0,   0, ...,   0,   0, 216],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0,   0]],

       [[  0, 216, 216, ..., 216, 216, 216],
        [216, 216, 216, ...,   0,   0,   0],
        [216, 216, 216, ...,   0,   0,   0],
        ...,
        [  0,   0,   0, ...,   0, 216, 216],
        [  0,   0,   0, ...,   0,   0, 216],
        [  0,   0,   0, ...,   0,   0, 216]],

       [[  0, 216, 216, ..., 216, 216, 216],
        [216, 216, 216, ..., 216,   0, 216],
        [216, 216, 216, ...,   0,   0,   0],
        ...,
        [  0,   0,   0, ...,   0, 216, 216],
        [  0,   0,   0, ...,   0,   0,   0],
        [  0,   0,   0, ...,   0,   0, 216]]], dtype=uint8)
Coordinates:
  * time         (time) datetime64[ns] 1988-01-01 1989-01-01 ... 2020-01-01
  * y            (y) float64 -3.755e+06 -3.755e+06 ... -3.767e+06 -3.767e+06
  * x            (x) float64 8.378e+05 8.378e+05 ... 8.477e+05 8.478e+05
    spatial_ref  int32 3577

If we try to generate an animation for this without specifying which DEA Land Cover classification should be used to colour the plot, we will get an error. In order to avoid getting this error, we provide the name of the classification being used to the ‘measurement’ variable. In this case, the measurement being used is "level3".

[11]:

# Set file name
file_name = "bare_surface_example"

# Generate plot
lc_animation(xr_bare, measurement='level3', file_name=file_name, colour_bar=True, width_pixels=15)

[11]:

<IPython.core.display.Image object>

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

Compatible datacube version:

[12]:

print(datacube.__version__)

1.8.11

Tags

**Tags**: :index:`sandbox compatible`, :index:`landsat 5`, :index:`landsat 7`, :index:`landsat 8`, :index:`DEA Land Cover`, :index:`time series`, :index: `LCCS`, :index:`colour maps`, :index:`data visualisation`