Create a GIF from a series of images

Creating a video is a cool way to visualise a time series of remote sensing data / outputs. Rather than creating a full video you can also use a GIF. One advantage of GIFs is you can embed them in Powerpoint presentations and they should will play in presentation mode (I’ve found them more reliable than videos in other formats).

You can use the ImageMagick convert command to create a GIF:

convert -delay 20 -loop 0 in_files/*.png out_animation.gif

However, depending on how your files are named it might not sort the numbers as expected leading to jumps in the video. This problem was addressed in the following post on the ‘remotelysensible’ blog:

Alphanumeric string sorting in Python

This can be used to sort the files in Python then call the convert command using subprocess:

#!/usr/bin/env python

import argparse
import re
import subprocess

re_natural = re.compile('[0-9]+|[^0-9]+')

def natural_key(s):
   """
   Natural sorting function from http://stackoverflow.com/questions/12184015/in-python-how-can-i-naturally-sort-a-list-of-alphanumeric-strings-such-that-alp
   """
   return [(1, int(c)) if c.isdigit() else (0, c.lower()) for c in re_natural.findall(s)] + [s]

if __name__ == '__main__':
   parser = argparse.ArgumentParser(description='Create an animated gif from a series of images')
   parser.add_argument("images", nargs='+',type=str, help="Input images")
   parser.add_argument('-o', '--outgif',
                       metavar ='Output gif',
                       help ='Output name for gif',
                       required=True)
   args=parser.parse_args()

   sorted_png_files = sorted(args.images, key=natural_key)

   convert_cmd = ['convert','-delay','20','-loop','0']
   convert_cmd.extend(sorted_png_files)
   convert_cmd.extend([args.outgif])

   subprocess.call(convert_cmd)

The script is run using:

python create_gif.py -o out_animation.gif in_files/*.png 

As well as creating GIFS of a time series you can use SPD 3D Points Viewer to export images of a point cloud as it rotates and create a gif from these similar to this video of airborne LiDAR data over Lake Vyrnwy on YouTube: https://www.youtube.com/watch?v=P_9998KJib4

The script can be downloaded from GitHub. Note the delay is hardcoded at 20/100th of a second but this can be changed if needed.

Apply the same stretch to multiple files in TuiView

There are some cases when it is necessary to play around with the stretch used to display data in TuiView (edit -> stretch). For example, if the scene contains very bright or dark areas (e.g., clouds/shadows) this can skew the statistics and make the images display too light/dark. One option is to use a local stretch (round house icon), which will only use values from the current display to calculate statistics. You can also manually set the minimum and maximum values.

Once a good stretch has been obtained this can be saved to the image for formats which support it (e.g., KEA) using the ‘Save Stretch and Lookup Table’ button (disk icon). If multiple files need to be opened with the same stretch applied you can pass in a file with a saved stretch when opening TuiView e.g.,

tuiview --stretchfromgdal image_with_saved_stretch.kea *kea

This will open all files matching the pattern ‘*bil’, applying the stretch saved in ‘image_with_saved_stretch.kea’

You can also save a stretch to a text file using the ‘Export Stretch and Lookup Table to Text File’ button (disk with ABC icon), this can be passed in when opening TuiView using

tuiview --stretchfromtext casi.stretch *bil

This will open all files matching the pattern ‘*bil’, applying the stretch saved in the text file ‘casi.stretch’

Viewing hyperspectral data in TuiView

As TuiView is based on GDAL, it has always been able to open the ENVI files commonly used for hyperspectral data. However, there have been some recent enhancements to TuiView which have improved handling of hyperspectral data.

TuiView can be installed for Linux, OS X and Windows through conda using:

conda install -c rios tuiview

The source is available to download from https://bitbucket.org/chchrsc/tuiview/.

As hyperspectral data contain many bands, specifying the ones to use for display when you load the file in makes things easier. You can do this from the command line using the ‘–bands’ argument. For example:

tuiview --rgb --stddev --bands 280,460,160 f159223b_mapped_osng.bil

For the Specim Fenix dataset used as an example the command will display a colour composite of NIR, SWIR and Red at approximately the centre wavelengths of the Landsat 8 bands with a standard deviation stretch applied to the pixel values. Opening the image may take a while as the statistics need to be calculated for the the three bands selected for display in order to apply a standard deviation stretch.

Once the file has opened you can use the ‘Query tool’ to display spectral profiles for each pixel.

If you get an error that TuiView is unable to open the file it could be that ‘data ignore value = 0’ is set in the header and the GDAL driver is unable to handle it. You can export the environmental variable:

export GDAL_PAM_ENABLED=ON

on Linux/OS X, or:

set GDAL_PAM_ENABLED=ON

on Windows. This will create a ‘*.aux.xml’ metadata file with additional features the driver doesn’t support when TuiView opens the file. More information is available in the GDALPamDataset documentation.

If you have hyperspectral data which isn’t mapped (e.g., level1b / level2 data in the NASA Data processing levels) TuiView will not open these by default and will give an error about only allowing north-up images. If you know the data don’t have geospatial information you can force TuiView to open them by setting:

export TUIVIEW_ALLOW_NOGEO=YES

on Linux/OS X, or:

set TUIVIEW_ALLOW_NOGEO=YES

on Windows.

More details about available environmental variables used by TuiView are available on the TuiView Wiki.

Hyperspectral data from the NERC ARSF used for this post were provided courtesy of NERC via the NERC Earth Observation Data Centre (NEODC).

Add a colour table in RSGISLib

For thematic rasters (e.g., classification) it is useful to save a colour scheme with the data for visualisation. This can be accomplished using a colour table saved as fields in an attribute table. There is a function in RSGISLib which will add colours to an existing RAT. To use the function on rasters which don’t already have an attribute table (i.e., the class is stored as the pixel value) one must be created and the pixel values copied to a column.

To accomplish this in RSGISLib the following is used:

import collections
from rsgislib import rastergis

classification='brigalow_regrowth_classification.kea'

# Add histogram (automatically adds attribute table)
rastergis.populateStats(classification, addclrtab=False, \
        calcpyramids=False, ignorezero=False)

# Add pixel values to attribute table
bandStats = []
bandStats.append(rastergis.BandAttStats(band=1, maxField='Class'))

rastergis.populateRATWithStats(classification, \
                                classification, bandStats)

# Add colour table
classcolours = {}
colourCat = collections.namedtuple('ColourCat', \
                        ['red', 'green', 'blue', 'alpha'])
classcolours[0] = colourCat(red=0, green=0, blue=0, alpha=0)
classcolours[1] = colourCat(red=255, green=0, blue=0, alpha=255)
classcolours[2] = colourCat(red=0, green=0, blue=255, alpha=255)
classcolours[3] = colourCat(red=0, green=200, blue=0, alpha=255)
classcolours[4] = colourCat(red=0, green=100, blue=0, alpha=255)
rastergis.colourClasses(classification, 'Class', classcolours)

# Add pyramids (for fast display)
rastergis.populateStats(classification, addclrtab=False, \
          calcpyramids=True, ignorezero=False)

If you don’t have an existing colour scheme you can add a random colour to each class by running ‘rastergis.populateStats’ and setting ‘addclrtab’ to True.

Note: if you add a colour table to a KEA file it will be recognised in TuiView and ArcMap (using the KEA driver) but not QGIS. Colour tables in ERDAS Imagine format work fine in QGIS.

Plot 2D Histogram of Pixel Values from Two Images

Recently I wanted to plot the pixel values of two images against each other. I though it would be good to combine extracting the pixel values (using RIOS) and plotting the data (using matplotlib) in a single script.

The script I wrote (two_band_scatter_plot.py) is available to download from the RSGIS Scripts repository

There are three main steps:

1. RIOS is used to iterate through the image in small blocks, for each block only pixels that have data in each band are retained.
2. The numpy.histogram2d is used to produce a 2D histogram.
3. The histogram is plotted using matplotlib.

The 2D histogram was based on code from http://oceanpython.org/2013/02/25/2d-histogram/. As there are a lot of points, displaying as a 2D histogram makes the information easier to see.

The simplest usage of the script is:

python two_band_scatter_plot.py \
--imageA imageA.kea \
--imageB imageB.kea \
--plot outplot.pdf 

Which will plot the first band of each image. To change the bands used, labels for the plot and scaling additional flags can be passed in.

python two_band_scatter_plot.py \
--imageA imageA.kea \
--imageB imageB.kea \
--bandA 1 --bandB 1 \
--plot outplot.pdf \
--plotMin 0 --plotMax 0.5 \ 
--labelA "13157 - 004 (HH)" \
--labelB "13157 - 006 (HH)" \

The output looks something like this:

The script assumes both images are the same projection and resolution, and RIOS takes care of the spatial information. It would be possible to adapt so RIOS resamples the data, see the RIOS Wiki for more details.

Create an animated GIF from a series of images

ImageMagick is a collection of command line tools which allows you to do lots of cool things with images. I use it a lot to create animations, in particular to add multiple images and text to a single slide. If you want to make something simpler than a movie, it is possible to generate a gif animation using ImageMagick. The simplest command is:

convert -delay 100 *files.png animation.gif

Which will create an animation from all files matching ‘*files.png’ with a 100 x 1/100 s (1 s) delay between images.

If you want to add a title for each image you can also do this with ImageMagick, using something similar to:

convert image.png -gravity south \
 -strokewidth 10 -pointsize 100 \
 -fill white -annotate 0 "Title" \
 image_with_title.png 

An example script, using RSGISLib to colour images and imagemagick to create a GIF, is available on github.

Export quicklooks with vector overlay using TuiView

Often when you have a large stack of images (e.g., a time series) it’s useful to create quick looks so you can flick through them. Normally I use gdal_translate with the ‘-scale’ and ‘-outsize’ flags, as described in a previous post.

Recently I wanted to quickly export all bands of an image for a presentation, with some points overlain. The TuiView manual describes how it’s possible to automate the ‘Save Current Display’ menu option using Python, which sounded ideal for this task.

Following the instructions to export a single image in the manual, I added a loop to export each band and name the output file with the band name, which produced exactly the result I needed. As this is something I’m likely to need to do again, I tidied the script up, made it a bit more general and added some more of the options available. The script can be downloaded from RSGIS Scripts on Bitbucket.

To export all bands in an image (original task):

python export_quicklook_tuiview.py --inimage tonzi_data_layers.kea \
       --outimage tonzi_data_layers_ql.png \
       --invector sensor_locations.shp \
       --allbands

I also added the options to export an image using the default stretch:

python export_quicklook_tuiview.py --inimage tonzi_data_layers.kea \
       --invector sensor_locations.shp \
       --outimage tonzi_data_layers_ql.png

Or if the data has already been stretched and RGB mapped to bands 1,2,3

python export_quicklook_tuiview.py --inimage tonzi_data_2sd_rgb.kea \
       --invector sensor_locations.shp \
       --outimage tonzi_data_rgb_2sd_ql.png --rgb --nostretch

To add text to the quick looks (e.g., the filename) imagemagick can be used:

convert in_image.png -gravity north -strokewidth 10 \
-pointsize 100 -annotate 0 "Title" out_image_with_title.png

To run for all images in a directory GNU Parallel can be used:

ls *kea | parallel \
python export_quicklook_tuiview.py --inimage {} \
       --invector sensor_locations.shp \
       --outimage {.}_ql.png