Introduction to RSGISLib Training Course

A course introducing RSGISLib was recently given in Japan by Pete Bunting. Material from the course is available to download from the following links:

• Slides (pdf)
• Worksheet (pdf)
• Package containing, slides, worksheet and data hosted on sourceforge (400 MB, tar.gz)

The course covers pre-processing of PALSAR SAR data followed by a demonstration of pixel-based and object-based classification. All the required data and scripts are included in the package.

If you have any questions on the course please use the RSGISLib support Google group: https://groups.google.com/forum/#!forum/rsgislib-support

Convert LaTeX to Word

A big problem with writing in LaTeX is collaborating with colleagues who don’t use it. One option is to generate a Word .docx version and use the comments and track changes features in Word / LibreOffice. This does require manually copying the changes back to LaTeX so isn’t quite as nice as using latexdiff (see earlier post) but is slightly easier than adding comments to a PDF.

The best program I’ve found for converting LaTeX to Word is the open source (GPL) command line tool, pandoc (http://pandoc.org/).

Basic usage is quite straightforward:

pandoc latex_document.tex -o latex_document_word_version.docx

The conversion isn’t perfect, figures and tables can get a bit mangled, but it does a good job with the text.

Pandoc can convert between many different formats, including from markdown and reStructuredText (commonly used for software documentation) so it is worth having installed.

Create a CSV with the coordinates from geotagged photos

Recently I took a lot of photos on my phone during fieldwork (survey for the NERC-ARF LiDAR and hyperspectral calibration flights) and wanted to extract the GPS coordinates from each photo so I could load them into QGIS using the add delimited text dialogue. The aim wasn’t to have precise locations for each photo (GPS positions for each of the points surveyed were recorded separately) but to give a quick idea of the locations we’d visited before the main GPS data were processed.

I while ago (2012!) I wrote a script for this task. The script (CreateJPEGKMZ) requires pillow and the imagemagick command line tools.

It code isn’t particularly tidy (despite my recent updates) but it basically pulls out the GPS location from the EXIF tags and writes this to a CSV file. To create a KMZ a thumbnail image is created and a corresponding KML file. Both the thumbnail and KML are then zipped together to generate the KMZ which can be opened in GoogleEarth.

To use the script to write a CSV file:

CreateJPEGKMZ.py --outcsv gps_points.csv input_jpeg_files

To also create KMZ files for each photo:

CreateJPEGKMZ.py --outcsv gps_points.csv \
                 --outkmz output_kmz_files 
                 input_jpeg_files

This was also a good lesson in the benefits of having scripts in version control and publicly available.

Convert EASE-2 grid cell to latitude and longitude using Python

The EASE-2 grid is used for a number of NASA datasets including SMAP. It is described in the following paper:

Brodzik, M. J., Billingsley, B., Haran, T., Raup, B., & Savoie, M. H. (2012). EASE-Grid 2.0: Incremental but Significant Improvements for Earth-Gridded Data Sets. ISPRS International Journal of Geo-Information, 1(3), 32–45. http://doi.org/10.3390/ijgi1010032

Files with the centre coordinate of each cell for EASE-2 grids at different resolutions are available from https://nsidc.org/data/ease/tools.html as well as tools for conversion.

To read these files into Python the following steps can be used:

1. Download the relevant gridloc file.

   The FTP link for the grid location files is: ftp://sidads.colorado.edu/pub/tools/easegrid2/

   For this example I’ve chosen the 36km cylindrical EASE-2 grid (gridloc.EASE2_M36km.tgz)

2. Un-tar using:

   tar -xvf gridloc.EASE2_M36km.tgz
   

3. Read the files into Python:

   import numpy
   
   # Read binary files and reshape to correct size
   # The number of rows and columns are in the file name
   lats = numpy.fromfile('EASE2_M36km.lats.964x406x1.double', 
                         dtype=numpy.float64).reshape((406,964))
   lons = numpy.fromfile('EASE2_M36km.lons.964x406x1.double', 
                         dtype=numpy.float64).reshape((406,964))
   
   # Extract latitude and longitude
   # for a given row and column 
   grid_row = 46
   grid_column = 470
   
   lat_val = lats[grid_row, grid_column]
   lon_val = lons[grid_row, grid_column]
   

Add an ENVI header to JAXA Global Mangrove Watch PALSAR tiles

This is an update to a script I wrote a while ago for working with PALSAR tiles from JAXA’s Global Mangrove Watch.

Untar files
Using GNU parallel [1] you can untar a directory of files at once:

ls *tar.gz | parallel -j 2 tar -xf 

Which will untar two files at a time (depending on how fast your disk is you may wish to set more running at once).

Add ENVI format header
Adding an ENVI header to each of the files (using information from the JAXA header file) allows the files to be easily read using GDAL. Following a similar method outlined in previous posts the script to add ENVI headers can be downloaded from add_gmw_palsar_header.py

It is run using:

python add_gmw_palsar_header.py KC_*

This will add a header for all files in each directory matching ‘KC*’

Create mosaic using RSGISLib

Using the ‘rsgislibmosaic.py’ utilities from RSGISLib, a mosaic can be created for each dataset using:

rsgislibmosaic.py -i . -s '*HH' -o palsar_hh.kea -ot UInt16
rsgislibmosaic.py -i . -s '*HV' -o palsar_hv.kea -ot UInt16
rsgislibmosaic.py -i . -s '*_linci' -o palsar_linci.kea -ot Byte
rsgislibmosaic.py -i . -s '*_mask' -o palsar_mask.kea -ot Byte
rsgislibmosaic.py -i . -s '*_date' -o palsar_date.kea -ot UInt16

A similar procedure should work for other datasets by modifying the Python script used to produce the ENVI header files.

[1] O. Tange (2011): GNU Parallel – The Command-Line Power Tool, ;login: The USENIX Magazine, February 2011:42-47.

Tracking changes in a LaTeX document

One of the problems people often have with using LaTeX for collaborative writing is that it is difficult to track changes in a document, like in Word.

As .tex files are text documents version control such as Git or Mercurial can be used to keep track of changes in the LaTeX source. However, looking at differences in .tex files is not as easy as having a formatted copy of the document with the changes marked, in particular for people not used to LaTeX.

The latexdiff command can be used to create a LaTeX document with the changes marked, from which a PDF can be created showing where the text has changed. For example:

Basic usage is:

latexdiff origionaldoc.tex changeddoc.tex > changes.tex

By default the command will print everything to the terminal so the output needs to be redirected (using >) to a file.

When writing papers we commonly have a shell script to generate a change tex file, create a PDF and remove temp files.

# Create diff file
latexdiff --exclude-textcmd "section,subsection,sub subsection" \
   201403_GOBIA_RSGISLib_RIOS_origional.tex \
   201403_GOBIA_RSGISLib_RIOS.tex > \
   201403_GOBIA_RSGISLib_RIOS_changes_temp.tex

# Create PDF
pdflatex 201403_GOBIA_RSGISLib_RIOS_changes_temp.tex
bibtex 201403_GOBIA_RSGISLib_RIOS_changes_temp
pdflatex 201403_GOBIA_RSGISLib_RIOS_changes_temp.tex
pdflatex 201403_GOBIA_RSGISLib_RIOS_changes_temp.tex

# Rename PDF and move extra files.
mv 201403_GOBIA_RSGISLib_RIOS_changes_temp.pdf 201403_GOBIA_RSGISLib_RIOS_changes.pdf
rm 201403_GOBIA_RSGISLib_RIOS_changes_temp.*

Note sometimes changes for particular tex commands can cause problems when creating a PDF, producing an error along these lines:

! Argument of \UL@word has an extra }.
    <inserted text>
                \par 
l.292 ...\DIFaddbegin \DIFadd{Methods}\DIFaddend }
                                                  \label{sec:method}

To stop this happening you can exclude particular commands (e.g., section) using the ‘–exclude-textcmd’ flag. (Thanks to this post on stack exchange for the tip).

Reading KEA files in MATLAB

The KEA file format is built on HDF5 so anything which can read HDF5 can read a KEA file, including MATLAB. I’m not quite sure why you’d want to use MATLAB instead of using Python, in particular as you can read KEA files and many other formats through GDAL. However, I’ve been claiming it’s possible to use KEA with MATLAB for a while so I though it would be a good idea to post an example!

This example uses CASI data over Injune, which is part of the test dataset for RSGISLib, you can download from here.

1. Get file information

info = h5info('injune_p142_casi_sub_ll.kea');
% Get bandnames
info.Groups.Name

2. Read in spatial information
If you don’t care about the geospatial information you can skip this step.

% Read in resolution
res = h5read('injune_p142_casi_sub_ll.kea','/HEADER/RES');
xRes = res(1);
yRes = res(2); 

% Top left coordinate
topLeft = h5read('injune_p142_casi_sub_ll.kea','/HEADER/TL');
topLeftX = topLeft(1); % Longitude
topLeftY = topLeft(2); % Latitude

% Size
imageSize = h5read('injune_p142_casi_sub_ll.kea','/HEADER/SIZE');
imageSizeX = imageSize(1);
imageSizeY = imageSize(2);

% Calculate bottom right
bottomRightX = topLeftX + (double(imageSizeX) * xRes);
bottomRightY = topLeftY + (double(imageSizeY) * yRes);

lon_axis = linspace(topLeftX,bottomRightX,imageSizeX);
lat_axis = linspace(bottomRightY,topLeftY,imageSizeY);

3. Read in data
This assumes reading in bands 12, 8 and 3.

band12 = transpose(h5read('injune_p142_casi_sub_ll.kea','/BAND12/DATA'));
band8 = transpose(h5read('injune_p142_casi_sub_ll.kea','/BAND8/DATA'));
band3 = transpose(h5read('injune_p142_casi_sub_ll.kea','/BAND3/DATA'));

4. Display data
A three band composite is created, this is stretched and then displayed.

% Create composite
composite = cat(3,band12,band8,band3);

% Apply linear stretch
stretchedComposite = imadjust(composite,stretchlim(composite));

% View image
h = imagesc(lon_axis, lat_axis, stretchedComposite);

If you’ve used the image from the example the output should look something like this:

More details about the KEA format (including a description of the data structure) are available in the following paper: Peter Bunting and Sam Gillingham, The KEA image file format, Computers & Geosciences, Volume 57, August 2013, Pages 54-58, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2013.03.025.

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.

Nested proportional symbols in OpenLayers

For the visualising data from the SoilSCAPE in-situ network of soil moisture sensors I wanted to use proportional symbols, to provide a visual comparison of the soil moisture from different nodes. The data is received in real time, so I wanted the maps to automatically update and show the current soil moisture. Additionally, as there were three sensors attached to each node I wanted to use three circles, nested within each other, with a different colour for each circle, to show the current moisture for all three sensors, which made it slightly more tricky.

This is the final map for Tonzi Ranch.

To produce the maps three main components were used: OpenLayers, imagemagick, and a Python script running as a cron job.

Maps – OpenLayers

All the maps on the SoilSCAPE site use OpenLayers, which allows you to display a text file using:

    var gsat = new OpenLayers.Layer.Google(
        "Google Satellite",
        {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22});

   // Read soil moisture sensor information (location, icon, description etc., from text file)
   var layer = new OpenLayers.Layer.Vector("Soil moisture sensors", {
        strategies: [new OpenLayers.Strategy.BBOX({resFactor: 1.1})],
        protocol: new OpenLayers.Protocol.HTTP({
    url: './all/all-points_sm.txt',
    format: new OpenLayers.Format.Text()
        })
    });

    // Add controls to map
    map = new OpenLayers.Map( 'map');

    // Add layers to map
    map.addLayers([gsat, layer]);
    map.addControl(new OpenLayers.Control.Zoom());
    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.zoomToExtent(mapBounds.transform(proj, map.getProjectionObject()));

The format for the text file is:

point	title	description	icon	iconSize	iconOffset
36.0021,-98.63108333	Node#101	Descriptions	icon.png	27	-2.5,-2.5

Real time updates – Python

The measurements from the sensors are sent to our lab via SMS and added to a MySQL database. The open layers text file is generated at regular intervals using a Python script to pull data from the database running as a cronjob. This approach is used to create maps of sensor status, with the link to the appropriate colour symbol saved in the open layers file based on the last received measurement. A description is also saved in the file with a link to a plot of measurements (also generated using a Python script running as a cronjob).

Symbols – imagemagick

To create the nested symbols imagemagick is used to create a different symbol for each node, by adjusting the size of three symbols and stacking them together. A Python script is used for this which can be downloaded from SoilSCAPE Utilities on Bitbucket. The script is run using:

python createPropSymbol.py outDIR outSymbolName.png Val1 Val2 Val3

To embed the maps within our Drupal site, an iframe is used.

More information about the SoilSCAPE project is available on the project website. Data can also be downloaded using the links at the bottom of each map. For more details on the wireless sensor aspect of the project see some of Agnelo Silva’s publications. The project is funded by NASA ESTO.

Version control with Mercurial

I’ve been quite enthusiastically advocating version control software (specifically mercurial) recently amongst some colleagues for managing code. Whilst the principles of version control are pretty familiar to most people, there are better ways than creating multiple copies with different names such as _final_final_verson2!

Version control software provides a better way of not only retaining different versions of files but allowing you to compare changes between them so you can easily find the lines that have been changed, particularly useful if you’ve just broken something! Mercurial and Git are two popular examples of version controls software. What’s nice about these (compared to other version control software such as subversion) is that all the changes are kept locally so you just need the software installed on your computer to have a folder under version control. For small projects, with only a couple of files, I just have local mercurial repositories set up.

When multiple people are working on the same code, these software work best when combined with hosting such as that provided by Bitbucket or github. My personal preference is Bitbucket as it offers mercurial (whereas github uses only git), allows private repositories for free (with unlimited users if you’re in academia) and I think it has a nicer interface. A lot of the software we’re developing / using is on Bitbucket (see bitbucket.org/petebunting/ and bitbucket.org/chchrsc/).

The command line interface for mercurial can be downloaded for Windows and OS X from http://mercurial.selenic.com, or can be installed from the package manager if you’re using linux:

sudo apt-get install mercurial

If you’d rather use a GUI than the command line TortoiseHG is a good option for Windows or sourcetree for Windows or OS X.

There are loads of good resources for learning mercurial on the internet, hginit.com is a particularly well written one and a good place to start. I’ve also put together a quick ‘cheat sheet’ with the commands I use most frequently.

Mercurial Cheat Sheet

 hg clone REMOTE LOCAL # Create local copy.
 hg pull -u # Pull changes from remote server (bitbucket) and update local copy.
 hg merge # Merge changes from remote server with local changes (sometimes required after pulling).
 hg commit # Commit changes locally (see setting up hgrc).
 hg push # Push changes to remote server (bitbucket).
 hg add # Add new file to version control (need to commit after).
 hg remove # Remove a file from filesystem and version control (need to commit after).
 hg mv # Move a file in the filesystem and version control (need to commit after).
 hg tag # Create name for revision (e.g., version number).
 hg update -r REV/TAG # Revert to revision or tag.
 hg revert FILE # Revert file to last commit (can also revert to revision)
 hg log # List all commits (good idea to pipe to less under UNIX).
 hg status # Show status of all files, can use -m to only show modified files.
 hg init # Initialise mercurial repository in directory (for new repositories only).

Setting up hgrc
The first time you commit you will likely an error ‘no username supplied’ to fix this you need to add the following lines to .hg/hgrc

 [ui]
 username = Firstname Lastname <email>
 verbose = True
 [paths]
 default = PATH TO REPOSITORY

The tag command is particularly useful for creating version numbers (e.g., version1.1) or noting the version of code used for publications etc. It does take some getting into but you soon get used to running ‘hg commit’ whenever you make changes.