The newer ‘ratapplier’ module is modelled after the ‘applier’ module for images and allows a function to be applied to chunks of rows, making it particularly useful for a RAT which is too large to load to memory.

In our recent paper on an open source system for object based image classification [1] we mentioned linking with scikit-learn [2] to apply different classification algorithms. This post presents an example using Random Forests to give an idea of all the steps required. Random Forests is an ensamble learning algorithm which utilises many decision trees, each of which ‘vote’ for the final class. To reduce correlation between trees only a subset of predictors (data layers) and training samples are used for each class. Random forests is increasing in popularity within remote sensing, an example of usage is the pixel-based classification of Alaska by Whitcombe et al 2009 [3]

Taking [3] as an example this post demonstrates how the classification can be applied at the object level.

Install software

As detailed in previous posts the software required can be installed under Linux/OS X through conda using

conda install -c osgeo \
rsgislib rios tuiview scikit-learn

Define training data

Random Forests works best with a large number of training samples. For this example we’re using the National Wetlands Inventory (NWI) which are available as an ESRI Shapefile. For each class an integer code (WETCLASS) has been assigned and stored within the attribute table. The polygons are converted to a raster, where the pixel ID is the class using:

Two stacks are created, one for the classification, which contains all data layers, and a second for segmentation which contains only SAR data. As the training data raster is categorical it is kept as a separate layer, because a separate function is required to attribute the segments.

Image segmentation

The following code is used for segmentation (for more detail see earlier post)

There are other algorithms in scikit-learn which can also be applied instead of Random Forests once the data is in the correct format. The big advantage of this system is the entire process can be applied within a single Python script so multiple algorithms / parameters can be easily tested and the performance evaluated.

There is an option to pass in a scale to convert the horizontal spacing from degrees to metres. However, if the DEM covers a range of latitudes it is desirable to convert the horizontal spacing on a per-pixel basis. To accomplish this I wrote a Python script, with RIOS used to read and write out data in blocks, making it efficient to process large datasets.

Implementation As the slope calculation requires looping through a block of data, a task which is slow in pure Python, two strategies are used to improve the speed:

Numba is used to compile the Python code – which is a lot faster than pure Python

Fortran code is used to perform the actual slope calculation, compiled as a Python module using f2py

To provide a comparison of the time required for each implementation slope was calculated from a 1800 x 1800 pixel DEM using each method:

So the Fortran version is slightly faster than the Numba version but only slightly and required a lot more effort to implement. Both are significantly faster than the pure Python version. The code provides a nice example of applying more complicated algorithms to images with RIOS used to handle the I/O.

To account for noise in the DEM there is also a version of the slope calculation which will use least squares fitting to fit a plane over a window of pixels and calculate the slope from this. As the Python version requires calls to the NumPy linear fitting code there is no improvement using Numba. The Fortran version uses Accelerate under OS X or ATLAS under Linux for the plane fitting.

This is the first post, of what I hope will be a regular feature, with some ‘Bonus Features’ of recently published papers I have been involved with. The idea is to provide some of the details considered too technical (or trivial) for an academic publication, code snippets, a bit of backstory and other things I think might be of interest.

Binaries for OS X and Linux are made available through conda. See the software page for more instructions.

RFC40 and RSGISLib 2.1 / 2.2 differences

Whilst we were writing the paper there were some major changes in RSGISLib due to the release of GDAL 1.11 which included the RFC40 changes proposed and implemented by Sam and Pete. These changes meant that rather than loading the entire RAT to memory (which was fast but the size of the RAT which could be processed was limited by the amount of available RAM), rows were accessed from disk (which was slower but removed the RAM limitation). For RSGISLib Pete and myself removed all the RAT functions after the 2.1 release and started adding them back to take advantage of the new RAT interface in GDAL 1.11. For some functions this required reworking the algorithm. After a pretty intense couple of days coding we managed to port most of the main functions across, and we’re gradually adding the rest. Therefore, some of the functions listed in the paper as available in 2.1 aren’t available in the latest release yet.

In RIOS, Neil and Sam added a ‘ratapplier’ interface to process large RATs chunks at a time, similar to the ‘applier’ interface for processing images. More information is in the RIOS documentation. The ratapplier interface is backwards compatible with pre-RFC40 versions of GDAL (although all rows are loaded at once).

The RFC40 changes in TuiView allowed massive RATs (e.g., segmentations for Australia and Alaska) to be easily visualised on moderate specification machines. Being able to open and navigate a RAT with 10s of millions of rows on a laptop really is an impressive feat!

Software Comparison

To run the segmentation in RSGISLib and OTB the following scripts were used (note OTB isn’t currently available in conda, for the paper we installed through the ubuntugis-unstable package repository).

There are lots of algorithms in OTB for segmentation so if you find the algorithm in RSGSILib doesn’t quite fit your requirements I’d highly recommend trying OTB. Although OTB can produce a raster output the vector output algorithm is able to process larger datasets and utilise multiple cores. To convert to a raster gdal_rasterize can be used. For example:

This can then be used exactly the same as the segmentation produced in RSGISLib. We think the ability to use different segmentation algorithms, from different packages, is a real benefit of the system.

Examples

The ‘Change in Mangroves Extent’ example was part of a course Pete and Myself gave at JAXA’s 20th Kyoto and Carbon Meeting (agenda and presentations available here, full workshop available to download from SourceForge). For segmenting the image and attributing objects there are two utility functions available in RSGISLib, described in an earlier post. For the course and paper we used the old RIOS RAT interface. However, it is recommended to use the new ratapplier interface. As an example the classification of water would be:

from rios import ratapplier
import numpy
def classifyWater(info, inputs, outputs):
# Read the 1996 JERS-1 Mean dB values for the clumps
# The column is represented as a numpy array of size
# block length x 1.
HH96MeandB = inputs.inrat.HH96MeandB
# Create a new numpy array with the same dimensions (i.e., length)
# as the 'HH96MeandB' array. The data type has been defined as
# an 8 bit integer (i.e., values from -128 to 128).
# All pixel values will be initialised to zero
Water96 = numpy.zeros_like(HH96MeandB, dtype=numpy.int8)
# Similar to an SQL where selection the where numpy where function
# allows a selection to be made. In this case all array elements
# with a 1996 HH value less than -12 dB are being selected and
# the corresponding elements in the Water96 array will be set to 1.
Water96 = numpy.where((HH96MeandB &amp;amp;lt; -12), 1, Water96)
# Save out to column 'Water96'
outputs.outrat.Water96 = Water96
# Set up inputs and outputs for ratapplier
inRats = ratapplier.RatAssociations()
outRats = ratapplier.RatAssociations()
inRats.inrat = ratapplier.RatHandle('N06W053_96-10_segs.kea')
outRats.outrat = ratapplier.RatHandle('N06W053_96-10_segs.kea')
print('Classifying water')
ratapplier.apply(classifyWater, inRats, outRats)

If you have any questions or comments about the system described in the paper email the RSGISLib support group (rsgislib-support@googlegroups.com).

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 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.

SSURGO (Soil Survey Geographic database) provides soil information across the United States. The data is provide as Shapefiles with the mapping units. The attributes for each polygon are stored as a text files, which need to be imported into an Access database and linked with the shapefile.

An alternative for working with SSURGO data is to convert the shapefile to a raster, parse the text files and store the attributes for each mapping unit as a Raster Attribute Table (RAT).