Plotting animations from time series imagery

What does this notebook do?

This notebook demonstrates how to import a time series of DEA cloud-free Landsat imagery from multiple sensors (i.e. Landsat 5, 7 and 8) as an xarray dataset, and then plot the data as a time series animation with one or two panels that can be either one or three-band images (e.g. true and false colour), or with a line graph (e.g. a pixel drill) in the right column. Animations can be produced as either GIFs, MP4s or WMVs for any area in Australia using a standard datacube query.

Requirements:

You need to run the following commands from the command line prior to launching jupyter notebooks from the same terminal so that the required libraries and paths are set:

module use /g/data/v10/public/modules/modulefiles

module load dea

module load ffmpeg

This notebook also uses three external functions called load_clearlandsat, animated_timeseries, animated_doubletimeseries and animated_timeseriesline. These functions are available in the 10_Scripts folder of the dea-notebooks Github repository. Note that these functions have been developed by DEA users, not the DEA development team, and so are provided without warranty. If you find an error or bug in the functions, please either create an ‘Issue’ in the Github repository, or fix it yourself and create a ‘Pull’ request to contribute the updated function back into the repository (See the repository README for instructions on creating a Pull request).

Date: October 2018

Author: Robbi Bishop-Taylor

Tags: plot, animation, time_series, gif, mp4, wmv, movie, Landsat5, Landsat7, Landsat8, load_clearlandsat, animated_timeseries, animated_doubletimeseries’, :index:`DEAPlotting, DEADataHandling, Scripts, pixeldrill, Outputting data, animated_timeseriesline

[1]:
# Import modules
import datacube
import sys
import os
import pandas as pd
from IPython.display import Image
import matplotlib.pyplot as plt

# Import external dea-notebooks functions using relative link to Scripts directory
sys.path.append('../10_Scripts')
import DEADataHandling
import DEAPlotting

# Set up datacube instance
dc = datacube.Datacube(app='Time series animation')

%load_ext autoreload
%autoreload 2

../10_Scripts/DEAPlotting.py:23: UserWarning: matplotlib.pyplot as already been imported, this call will have no effect.
  matplotlib.use('agg')

Import satellite data to animate

Define the query bounds for datacube extraction using a dict. This should include x and y limits, potentially a list of measurements (i.e. the bands you want to extract like ‘red’, ‘green’, ‘blue’; this significantly speeds up the import) and a time extent. If no time is given, the function defaults to all timesteps available to all sensors (e.g. 1987-2018).

[2]:
# Set up spatial and temporal query.
query = {'x': (969476, 988476),
         'y': (-3568950, -3551951),
         'time': ('2015-01-01', '2017-01-01'),
         'crs': 'EPSG:3577'}

Extract cloud-free clear Landsat observations from all sensors

Use the load_clearlandsat function to load Landsat observations and PQ data for multiple sensors (i.e. ls5, ls7, ls8), and return a single xarray dataset containing only observations that contain greater than a specified proportion of clear pixels. This uses dask to only load in the filtered observations, and results in a visually appealing time series of observations that are not affected by cloud!

[3]:
# Set the minimum proportion of clear pixels (pixels with no clouds or other nodata)
masked_prop=0.99

# Load in only clear Landsat observations with < 1% unclear values
ds = DEADataHandling.load_clearlandsat(dc=dc, query=query,
                                       bands_of_interest=['red', 'green', 'blue', 'nir', 'swir1'],
                                       masked_prop=masked_prop)
print(ds)

# Plot first five time series images:
ds[['red', 'green', 'blue']].isel(time=[1, 2, 3, 4, 5]).to_array().plot.imshow(col='time', robust=True)
Loading ls5 pixel quality
    Skipping ls5; no valid data for query
Ignoring SLC-off observations for ls7
Loading ls7 pixel quality
    Loading 0 filtered ls7 timesteps
Loading ls8 pixel quality
    Loading 17 filtered ls8 timesteps
Combining and sorting ls7, ls8 data
    Replacing invalid -999 values with NaN (data will be coerced to float64)
<xarray.Dataset>
Dimensions:    (time: 17, x: 760, y: 680)
Coordinates:
  * y          (y) float64 -3.552e+06 -3.552e+06 ... -3.569e+06 -3.569e+06
  * x          (x) float64 9.695e+05 9.695e+05 9.695e+05 ... 9.884e+05 9.885e+05
  * time       (time) datetime64[ns] 2015-01-22T00:20:21 ... 2016-12-10T00:20:41.500000
Data variables:
    red        (time, y, x) float64 1.157e+03 1.155e+03 ... 2.034e+03 2.574e+03
    green      (time, y, x) float64 818.0 798.0 830.0 ... 1.681e+03 2.039e+03
    blue       (time, y, x) float64 563.0 553.0 557.0 ... 1.244e+03 1.499e+03
    nir        (time, y, x) float64 1.939e+03 1.945e+03 ... 3.107e+03 3.5e+03
    swir1      (time, y, x) float64 3.164e+03 3.186e+03 ... 4.049e+03 4.208e+03
    data_perc  (time) float64 1.0 0.9999 0.9996 0.9928 ... 0.9996 0.9999 1.0
Attributes:
    crs:      EPSG:3577
[3]:
<xarray.plot.facetgrid.FacetGrid at 0x7fb1febffef0>
../../_images/notebooks_08_Outputting_data_PlottingAnimatedGifs_7_2.png

Plot entire time series as a one-panel animated GIF

The animated_timeseries function takes an xarray time series and exports a one band or three band (e.g. true or false colour) GIF, MP4 or WMV animation showing changes in the landscape across time. Here, we plot the xarray as an animated GIF that includes a date annotation for each frame. We set the interval between the animation frames to 80 milliseconds. For three-band RGB animations like this, the function will automatically select an appropriate colour stretch by clipping the data to remove values smaller or greater than the 2 and 98th percentiles. This can be controlled further with the percentile_stretch parameter.

[4]:
# Produce time series animation of red, green and blue bands
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                interval=80, width_pixels=350)

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[4]:
<IPython.core.display.Image object>

We can also use different band combinations (e.g. false colour), add a title, and change the font size using annotation_kwargs, which passes a dictionary of values to the matplotlib plt.annotate function (see https://matplotlib.org/api/_as_gen/matplotlib.pyplot.annotate.html for options):

[5]:
# Produce time series animation of red, green and blue bands
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                width_pixels=350,
                                bands=['swir1', 'nir', 'green'],
                                title='Time-series animation',
                                interval=80, percentile_stretch=[0.01, 0.99],
                                annotation_kwargs={'fontsize': 33})

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[5]:
<IPython.core.display.Image object>

It is also possible to plot a single band image instead of a three band. For example, we could plot an index like the Normalized Difference Water Index (NDWI), which has high values where a pixel is likely to be open water (e.g. NDWI > 0). By default the colour bar limits are set based on percentile_stretch; set percentile_stretch=(0.0, 1.00) to show the full data range, or reduce it to e.g. percentile_stretch=(0.01, 0.99) to optimise colours by discarding outliers/extreme values.

[6]:
# Compute NDWI using the formula (green - nir) / (green + nir). This will calculate
# NDWI for every time-step in the dataset:
ds_ndwi = (ds.green - ds.nir) / (ds.green + ds.nir)

# We can now add this back into our dataset as a new data variable:
ds['NDWI'] = ds_ndwi

# Produce time series animation of NDWI:
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                width_pixels=350, bands=['NDWI'],
                                title='NDWI', interval=80,
                                percentile_stretch=(0.0, 1.00))

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')

Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[6]:
<IPython.core.display.Image object>

We can customise animations based on a single band like NDWI by specifying parameters using the onebandplot_kwargs which will be sent to the matplotlib plt.imshow function (see https://matplotlib.org/api/_as_gen/matplotlib.pyplot.imshow.html for options). For example, we can use a more appropriate blue colour scheme with 'cmap':'Blues', and set 'vmin':0.1, 'vmax':0.8 to overrule the default colour bar limits with manually specified values:

[7]:
# Produce time series animation using a custom colour scheme and limits:
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                width_pixels=350, bands=['NDWI'],
                                title='NDWI', interval=80,
                                onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8})

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[7]:
<IPython.core.display.Image object>

Two special kwargs (tick_fontsize, tick_colour) can be used to control the tick labels on the colourbar. This can be useful for example when the tick labels are difficult to see against a dark background:

[8]:
# Produce time series animation using a custom colour scheme and limits:
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                width_pixels=350, bands=['NDWI'],
                                title='NDWI', interval=80,
                                onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8,
                                                    'tick_fontsize': 25, 'tick_colour': 'grey'})

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[8]:
<IPython.core.display.Image object>

One band animations show a colour bar by default, but this can be disabled:

[9]:
# Produce time series animation without a colour bar:
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                width_pixels=350, bands=['NDWI'],
                                title='NDWI', interval=80,
                                onebandplot_cbar=False,
                                onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8})

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[9]:
<IPython.core.display.Image object>

Plot a two panel side-by-side animation

We can use the animated_doubletimeseries function to plot an animation that compares two different xarray datasets or band combinations at the same time:

[10]:
# Animate datasets
DEAPlotting.animated_doubletimeseries(ds1=ds, ds2=ds,
                                      output_path='animated_doubletimeseries.gif',
                                      width_pixels=500,
                                      bands1=['red', 'green', 'blue'],
                                      bands2=['swir1', 'nir', 'green'],
                                      percentile_stretch2=[0.01, 0.99],
                                      title1='True colour', title2='False colour')

# Plot animated gif
plt.close()
Image(filename='animated_doubletimeseries.gif')

Generating 17 frame animation (i.e. timesteps in shortest dataset)
    Exporting animation to animated_doubletimeseries.gif
[10]:
<IPython.core.display.Image object>

We can also plot a three-band image on one side of the plot, and a single band image on the other side:

[11]:
# Animate datasets
DEAPlotting.animated_doubletimeseries(ds1=ds, ds2=ds,
                                      output_path='animated_doubletimeseries.gif',
                                      width_pixels=500,
                                      bands1=['red', 'green', 'blue'], bands2=['NDWI'],
                                      title1='True colour', title2='NDWI',
                                      onebandplot_kwargs2={'cmap':'Blues'})

# Plot animated gif
plt.close()
Image(filename='animated_doubletimeseries.gif')
Generating 17 frame animation (i.e. timesteps in shortest dataset)
    Exporting animation to animated_doubletimeseries.gif
[11]:
<IPython.core.display.Image object>

Annotations can be customised for each panel individually:

[12]:
# Animate datasets
DEAPlotting.animated_doubletimeseries(ds1=ds, ds2=ds,
                                      output_path='animated_doubletimeseries.gif',
                                      width_pixels=500,
                                      bands1=['red', 'green', 'blue'], bands2=['NDWI'],
                                      title1='True colour', title2='NDWI',
                                      onebandplot_kwargs2={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8, 'tick_fontsize': 15},
                                      annotation_kwargs1={'fontsize': 25},
                                      annotation_kwargs2={'fontsize': 18})

# Plot animated gif
plt.close()
Image(filename='animated_doubletimeseries.gif')

Generating 17 frame animation (i.e. timesteps in shortest dataset)
    Exporting animation to animated_doubletimeseries.gif
[12]:
<IPython.core.display.Image object>

Plotting a side-by-side time series and line graph

The final animated_timeseriesline function allows you to plot an xarray dataset in the left panel of the animation, and a set of line graphs from a pandas dataframe in the right hand panel. This may be useful for comparing satellite imagery to time series data like river guages, water quality indices, or pixel drills across time.

First, you need time series data in pandas dataframe format. The dataframe must contain time steps in a DatetimeIndex column, and one or more numeric data columns to plot as lines in the right panel. Column names are used to label the lines on the plot, so assign them informative names. Lines are plotted by showing all parts of the line with dates on or before the current timestep (i.e. for a 2006 time step, only the portion of the lines with dates on or before 2006 will be plotted for that frame. Here, we will extract a pixel drill of surface water data from the centre of the dataset:

[13]:
# Extract time series of NDWI values at location (340, 340)
line_df = ds.NDWI.isel(x=340, y=340).to_dataframe()

# Drop all unnecessary columns, and rename to give informative name on plot
line_df.drop(columns=['y', 'x'], inplace=True)
line_df.rename(columns={'NDWI':'Normalized Difference Water Index (NDWI)'}, inplace = True)

# Print first 10 rows
line_df.head(10)

[13]:
Normalized Difference Water Index (NDWI)
time
2015-01-22 00:20:21.000 -0.158416
2015-02-07 00:20:18.000 -0.237256
2015-03-11 00:20:02.000 -0.298246
2015-03-27 00:19:53.500 -0.317982
2015-07-01 00:19:47.000 -0.105911
2015-09-19 00:20:21.000 -0.165733
2015-11-22 00:20:32.500 -0.198198
2015-12-24 00:20:32.000 -0.197972
2016-01-09 00:20:27.500 -0.198205
2016-02-10 00:20:23.000 -0.243345
[14]:
# Animate line plot in right panel of animation
DEAPlotting.animated_timeseriesline(ds=ds, df=line_df,
                                    output_path='animated_timeseriesline.gif',
                                    width_pixels=500,
                                    bands=['NDWI'],
                                    percentile_stretch=(0.00, 1.00))

# Plot animated gif
plt.close()
Image(filename='animated_timeseriesline.gif')

/g/data/v10/public/modules/dea-env/20181015/lib/python3.6/site-packages/matplotlib/figure.py:2362: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  warnings.warn("This figure includes Axes that are not compatible "
    Exporting animation to animated_timeseriesline.gif
[14]:
<IPython.core.display.Image object>

To demonstrate how multiple lines can be added to the same plot:

[15]:
# Add new column to dataframe
line_df['Inverted NDWI'] = line_df['Normalized Difference Water Index (NDWI)'] * -1.0

# Animate line plot in right panel of animation
DEAPlotting.animated_timeseriesline(ds=ds, df=line_df,
                                    output_path='animated_timeseriesline.gif',
                                    width_pixels=500,
                                    bands=['NDWI'],
                                    percentile_stretch=(0.00, 1.00))

# Plot animated gif
plt.close()
Image(filename='animated_timeseriesline.gif')

/g/data/v10/public/modules/dea-env/20181015/lib/python3.6/site-packages/matplotlib/figure.py:2362: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  warnings.warn("This figure includes Axes that are not compatible "
    Exporting animation to animated_timeseriesline.gif
[15]:
<IPython.core.display.Image object>

Lines do not need to have the same timestamps. To combine two different time series, you can simply concatenate the two dataframes and sort by date, and the function will take care of nodata values.

[16]:
# Create example dataframe
secondline_df = pd.DataFrame(data={'Test data': [0, 0.2, 0.4]},
                             index=[pd.Timestamp('2015-03-01'),
                                    pd.Timestamp('2016-01-01'),
                                    pd.Timestamp('2016-06-01')])

secondline_df
[16]:
Test data
2015-03-01 0.0
2016-01-01 0.2
2016-06-01 0.4
[17]:
# Concatenate into one dataframe and sort by date
combinedline_df = pd.concat([line_df, secondline_df])
combinedline_df = combinedline_df.sort_index()

# Print the first 10 rows
combinedline_df.head(10)
/g/data/v10/public/modules/dea-env/20181015/lib/python3.6/site-packages/ipykernel/__main__.py:2: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  from ipykernel import kernelapp as app
[17]:
Inverted NDWI Normalized Difference Water Index (NDWI) Test data
2015-01-22 00:20:21.000 0.158416 -0.158416 NaN
2015-02-07 00:20:18.000 0.237256 -0.237256 NaN
2015-03-01 00:00:00.000 NaN NaN 0.0
2015-03-11 00:20:02.000 0.298246 -0.298246 NaN
2015-03-27 00:19:53.500 0.317982 -0.317982 NaN
2015-07-01 00:19:47.000 0.105911 -0.105911 NaN
2015-09-19 00:20:21.000 0.165733 -0.165733 NaN
2015-11-22 00:20:32.500 0.198198 -0.198198 NaN
2015-12-24 00:20:32.000 0.197972 -0.197972 NaN
2016-01-01 00:00:00.000 NaN NaN 0.2
[18]:
# Animate line plot in right panel of animation, using dots on the lines to highlight each observation point
DEAPlotting.animated_timeseriesline(ds=ds, df=combinedline_df,
                                    output_path='animated_timeseriesline.gif',
                                    width_pixels=500,
                                    bands=['red', 'green', 'blue'],
                                    pandasplot_kwargs={'style':'.-', 'ms': 10, 'ylim': (-1.0, 1.0)})

# Plot animated gif
plt.close()
Image(filename='animated_timeseriesline.gif')

    Exporting animation to animated_timeseriesline.gif
[18]:
<IPython.core.display.Image object>

Shapefile overlay

All three animation functions support plotting shapefiles over satellite imagery. To do this, simply pass in a shapefile path with shapefile_path. The shapefile must be in the same projection system as the satellite data (e.g. Australian Albers EPSG:3577 below):

[19]:
# Produce time series animation of red, green and blue bands
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                interval=150, width_pixels=350,
                                bands=['NDWI'], onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8},
                                shapefile_path='../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_poly.shp')

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[19]:
<IPython.core.display.Image object>

You can customise styling for the shapefile overlays using the shapefile_kwargs argument:

[20]:
# Produce time series animation of red, green and blue bands
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                interval=150, width_pixels=350,
                                bands=['NDWI'], onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8},
                                shapefile_path='../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_poly.shp',
                                shapefile_kwargs = {'linewidth': 2, 'edgecolor': 'blue', 'facecolor': 'red'})

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[20]:
<IPython.core.display.Image object>

Multiple shapefiles can be overlaid by providing a list of shapefile paths:

[21]:
# Produce time series animation of red, green and blue bands
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                interval=150, width_pixels=350,
                                bands=['NDWI'], onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8},
                                shapefile_path=['../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_poly.shp',
                                                '../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_point.shp',
                                                '../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_line.shp'])

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[21]:
<IPython.core.display.Image object>

To customise styling for each shapefile individually, you can pass in a list of dictionaries to shapefile_kwargs. This list should be equal in length to the number of shapefile paths you pass to shapefile_path (e.g. one dictionary corresponding to each path). Note that at the current time this functionality is only available for the animated_timeseries function (not animated_timeseriesline or animated_doubletimeseries):

[22]:
# Produce time series animation of red, green and blue bands
DEAPlotting.animated_timeseries(ds=ds, output_path='animated_timeseries.gif',
                                interval=150, width_pixels=350,
                                bands=['NDWI'], onebandplot_kwargs={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8},
                                shapefile_path=['../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_poly.shp',
                                                '../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_point.shp',
                                                '../Supplementary_data/Files/PlottingAnimatedGifs/PlottingAnimatedGifs_line.shp'],
                                shapefile_kwargs = [{'linewidth': 2, 'edgecolor': 'blue', 'facecolor': 'red'},
                                                    {'markersize': 100, 'edgecolor': 'black', 'facecolor': 'black'},
                                                    {'linewidth': 4, 'edgecolor': 'orange'},])

# Plot animated gif
plt.close()
Image(filename='animated_timeseries.gif')
Generating 17 frame animation
    Exporting animation to animated_timeseries.gif
[22]:
<IPython.core.display.Image object>

Available output formats

The above examples have focused on exporting animated GIFs, but MP4 and WMV files can also be generated: * .mp4: fast to generate, smallest file sizes and often highest quality; suitable for Twitter/social media and recent versions of Powerpoint * .wmv: slow to generate, large file sizes; suitable for Powerpoint 2010 * .gif: slow to generate, large file sizes; suitable for all versions of Powerpoint and Twitter/social media

[23]:
# Animate datasets as a WMV file
DEAPlotting.animated_doubletimeseries(ds1=ds, ds2=ds,
                                      output_path='animated_doubletimeseries.wmv',
                                      bands1=['red', 'green', 'blue'], bands2=['NDWI'],
                                      title1='True colour', title2='NDWI',
                                      onebandplot_kwargs2={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8},
                                      annotation_kwargs1={'fontsize': 25},
                                      annotation_kwargs2={'fontsize': 13})

# Animate datasets as a MP4 file
DEAPlotting.animated_doubletimeseries(ds1=ds, ds2=ds,
                                      output_path='animated_doubletimeseries.mp4',
                                      bands1=['red', 'green', 'blue'], bands2=['NDWI'],
                                      title1='True colour', title2='NDWI',
                                      onebandplot_kwargs2={'cmap':'Blues', 'vmin':0.1, 'vmax':0.8},
                                      annotation_kwargs1={'fontsize': 25},
                                      annotation_kwargs2={'fontsize': 13})

Generating 17 frame animation (i.e. timesteps in shortest dataset)
    Exporting animation to animated_doubletimeseries.wmv
Generating 17 frame animation (i.e. timesteps in shortest dataset)
    Exporting animation to animated_doubletimeseries.mp4
../../_images/notebooks_08_Outputting_data_PlottingAnimatedGifs_44_1.png
../../_images/notebooks_08_Outputting_data_PlottingAnimatedGifs_44_2.png
[ ]: