Scalable image segmentation using RSGISLib

To for application to very large remote sensing datasets, an approach to “Scalable image segmentation” presented in [1] using RSGISLib. In the paper a 30 m spatial resolution satellite mosaic of Australia was segmented by splitting into tiles, processing each tile on a separate node of a HPC, merging and then performing a second segmentation to remove artefacts at tile boundaries.
For a full description of the approach see section 5.3 of the open access paper (available to view online at http://www.mdpi.com/2072-4292/6/7/6111/htm).

A version of this approach, designed to be used on a single machine, is now available in RSGISLib through the performTiledSegmentation function. The function is called in a similar way to the existing runShepherdSegmentation function.

For instructions on installing the latest version of RSGISLib see this post, skipping to section 7 if you already have Linux or OS X installed.

The code below presents an example of applying either the tiled or standard segmentation through the use of a command line flag:

#!/usr/bin/env python
"""
Apply segmentation to image
"""

from __future__ import print_function
import argparse
import shutil
import sys
import tempfile
from rsgislib.segmentation import tiledsegsingle
from rsgislib.segmentation import segutils

# Set values for clustering
NUM_CLUSTERS = 60
MIN_PIXELS = 100
DIST_THRESHOLD = 100

# Set values for tiles segmentation
TILE_WIDTH = 2000
TILE_HEIGHT = 2000

parser = argparse.ArgumentParser()
parser.add_argument("inputimage", nargs=1,
                    type=str, help="Input Image")
parser.add_argument("outputclumps", nargs=1,
                    type=str, help="Output clumps")
parser.add_argument("--tiled",
                    default=False,
                    action="store_true",
                    help="Use tiled segmentation")
args = parser.parse_args()

# Make temp directory for intermediate files
temp_dir = tempfile.mkdtemp(prefix='rsgislib_seg_')

if args.tiled:
    # If requested run tiled segmentation
     tiledsegsingle.performTiledSegmentation(args.inputimage[0],
                                             args.outputclumps[0],
                                             tmpDIR=temp_dir,
                                             tileWidth=TILE_WIDTH,
                                             tileHeight=TILE_HEIGHT,
                                             validDataThreshold=0.3,
                                             numClusters=NUM_CLUSTERS,
                                             minPxls=MIN_PIXELS,
                                             distThres=DIST_THRESHOLD,
                                             sampling=100, kmMaxIter=200)
else:
    # If not run standard
     segutils.runShepherdSegmentation(args.inputimage[0],
                                      args.outputclumps[0],
                                      tmpath=temp_dir,
                                      numClusters=NUM_CLUSTERS,
                                      minPxls=MIN_PIXELS,
                                      distThres=DIST_THRESHOLD,
                                      sampling=100, kmMaxIter=200)
shutil.rmtree(temp_dir)

Note, script edited for post – full version available from GitHub.

[1] Clewley, D.; Bunting, P.; Shepherd, J.; Gillingham, S.; Flood, N.; Dymond, J.; Lucas, R.; Armston, J.; Moghaddam, M. A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables. Remote Sensing 2014, 6, 6111-6135. http://www.mdpi.com/2072-4292/6/7/6111

Copy a Shapefile to a Raster Attribute Table

The object based image analysis features in RSGISLib are based around storing object attributes as a Raster Attribute Table (RAT), as described in the following paper:

Clewley, D.; Bunting, P.; Shepherd, J.; Gillingham, S.; Flood, N.; Dymond, J.; Lucas, R.; Armston, J.; Moghaddam, M. A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables. Remote Sensing 2014, 6, 6111-6135. (open access)

However, data to be used as part of the analysis are often stored as attributes of a shapefile. A RAT can be created from a shapefile using the gdal_rasterize command and the copyShapefile2RAT function in RSGISLib using the following steps:

1. Rasterise vector
   This can done directly from the command line or by calling gdal_rasterize from a Python script using subprocess.

   import subprocess 
   
   input_vector = 'brig_re09_binary.shp'
   rasterised_vector = 'brig_re09_binary_raster.kea'
       
   rasterise_cmd = ['gdal_rasterize','-of','KEA',
                                     '-a','BGL',
                                     '-tr','30','30',
                                     input_vector,
                                     rasterised_vector]
   print('Running: ', ' '.join(rasterise_cmd))
   subprocess.call(rasterise_cmd)
   
   

   The above code sets the output raster pixel size to 30 x 30 m and uses the ‘BGL’ column from the shapefile (which is an integer) for the pixel values.

   Note this image is only used to get the extent and pixel size of the final RAT image. If you have another image you wish to match the extent and pixel size to this can be used instead of rasterising the vector, in this case skip this step and use the existing image instead of ‘rasterised_vector’.

2. Create RAT and copy shapefile attributes

   from rsgislib import vectorutils
   import rsgislib
    
   input_vector = 'brig_re09_binary.shp'
   rasterised_vector = 'brig_re09_binary_raster.kea'
   output_rat_image = 'brig_re09_attributes.kea'
       
   vectorutils.copyShapefile2RAT(input_vector, rasterised_vector,output_rat_image)
   
   

   This will copy all the attributes from the input shapefile to a RAT. Note the output format is always KEA as it has support for large attribute tables and stores them with compression.

To view the attributes open the file in TuiView. The output RAT can then be used as any other RAT, e.g., attributing objects with image statistics or using as part of a classification. For more details see other posts on object based image analysis.

Object-based classification using Random Forests

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.

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

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

gdal_rasterize -of KEA -ot Byte -a WETCLASS \
     -tr 100 100 nwi_shp.shp nwi_raster.kea

3. Stack bands

To stack the bands (assumed to have the same projection and resolution) the stack bands command, within RSGISLib, is used.

#/usr/bin/env python
import rsgislib
from rsgislib import imageutils

# Segmentation files
palsarList = ['palsar_hh.kea',
            'also_hv.kea']

# All files
datalayersList = palsarList + \
                 ['palsar_hh_tex.kea',
                  'palsar_hv_tex.kea',
                  'elevation.kea',
                  'slope.kea']

# Band names
bandNamesList = ['hh','hv',
                 'hh_tex','hv_tex'
                 'elevation', 'slope']

# Out file type and format
gdalformat = 'KEA'
dataType = rsgislib.TYPE_32FLOAT

# Create stack for segmentation
outputImage = 'palsar_hhhv.kea'
imageutils.stackImageBands(palsarList, None, \
     outputImage, None, 0, gdalformat, dataType)

# Create stack of all data layers
outputDataStack = 'classification_stack.kea'
imageutils.stackImageBands(datalayersList, bandNamesList, \
     outputDataStack, None, 0, gdalformat, dataType)

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.

4. Image segmentation

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

#/usr/bin/env python
from rsgislib.segmentation import segutils
 
inputImage = 'palsar_hhhv.kea'
clumpsFile = 'palsar_hhhv_clumps_elim_final.kea'
meanImage = 'palsar_hhhv_clumps_elim_final_mean.kea'

# Run segmentation
segutils.runShepherdSegmentation(inputImage, clumpsFile,
                    meanImage, numClusters=100, minPxls=100)

The method is similar to that described in an earlier post, with the addition of a step to include categorical data.

5. Attribute segments

#!/usr/bin/env python

from rsgislib import rastergis
from rsgislib.rastergis import ratutils
from rsgislib import imageutils
 
dataStack = 'classification_stack.kea'
classFile = 'nwi_raster.kea'
clumpsFile = 'palsar_hhhv_clumps_elim_final.kea'
 
# Attribute segments with data
ratutils.populateImageStats(dataStack, clumpsFile,
                    calcMean=True)

# Convert training data to RAT
codeStats = list()
codeStats.append(rastergis.BandAttStats(band=1, minField='wetCode'))
rastergis.populateRATWithStats(classFile, classFile, codeStats)

# Attribute segments with class
rastergis.strClassMajority(clumpsFile, classFile, \
              'wetCode', 'wetCode', False)

6. Perform classification

Much of the classification code is sorting the data from the RAT so it can be input into scikit-learn.

#!/usr/bin/env python
from rios import rat
import osgeo.gdal as gdal
import numpy as np
from sklearn.ensemble import RandomForestClassifier

# Open RAT
inRatFile = 'palsar_hhhv_clumps_elim_final.kea'
ratDataset = gdal.Open(inRatFile, gdal.GA_Update)

# Set column names
x_col_names = ['hh','hv',
               'hh_tex','hv_tex'
               'elevation', 'slope']

y_col_name = 'wetCode'

# Set up list to hold data
X = []

# Read in data from each column
for colName in x_col_names:
    X.append(rat.readColumn(ratDataset, colName))

# Read in training data
y = rat.readColumn(ratDataset, y_col_name) 
# Set NA values to 0
y = np.where(y == b'NA',0,y)
y = y.astype(np.int16)

X.append(y)

X = np.array(X)
X = X.transpose()

# Remove rows with 0 (NA) for wetCode
X_train = X[X[:,-1] != 0]

# Remove non-finite values
X_train = X_train[np.isfinite(X_train).all(axis=1)]

# Split into variables (X) and class (y)
y_train = X_train[:,-1]
X_train = X_train[:,0:-1]

# Train Random Forests
clf = RandomForestClassifier(n_estimators=500, max_features=3, \
      oob_score=True, n_jobs=6, verbose=2)

clf.fit(X_train, y_train)

# Set NaN values to 0
X = np.where(np.isfinite(X),X,0)

# Apply classification
predictClass = clf.predict(X[:,0:-1])

# Write out data to RAT
rat.writeColumn(ratDataset, 'predictClass', predictClass)
ratDataset = None

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.

References

[1] Clewley, D.; Bunting, P.; Shepherd, J.; Gillingham, S.; Flood, N.; Dymond, J.; Lucas, R.; Armston, J.; Moghaddam, M. A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables. Remote Sensing 2014, 6, 6111-6135.
[2] Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.
[3] Whitcomb, J., Moghaddam, M., McDonald, K., Podest, E., Kellndorfer, J., Wetlands Map of Alaska Using L-Band Radar Satellite Imagery, Canadian Journal of Remote Sensing, 2009, Vol. 35, pp. 54-72

A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables – Bonus Features

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.

This first post is on:

Clewley, D.; Bunting, P.; Shepherd, J.; Gillingham, S.; Flood, N.; Dymond, J.; Lucas, R.; Armston, J.; Moghaddam, M. A Python-Based Open Source System for Geographic Object-Based Image Analysis (GEOBIA) Utilizing Raster Attribute Tables. Remote Sensing 2014, 6, 6111-6135.

The article is open access and is available to download from:

http://www.mdpi.com/2072-4292/6/7/6111

1. Installing the software.

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

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

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

   RSGISLib (Python):

   from rsgislib.segmentation import segutils
    
   inputImage = '../Data/naip_newhogansouth_2012_sub.tif'
   clumpsFile = 'naip_newhogansouth_2012_clumps_elim_final.kea'
     
   # Run segmentation
   segutils.runShepherdSegmentation(inputImage, clumpsFile,
                              numClusters=60, minPxls=100, 
                              distThres=100, bands=None, 
                              sampling=100, kmMaxIter=200)
   

   OTB (Bash):

   otbcli_Segmentation -in ../Data/naip_newhogansouth_2012_sub.tif \
       -mode.vector.out naip_newhogansouth_2012_otb_seg.shp
   

   To time the segmentation the UNIX ‘time’ command was used.

   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:

   gdal_rasterize -a DN -tr 30 30 -ot UInt32 \
      -of KEA segmentation_vector.shp \
      segmentation_raster.kea
   

   To convert to a RAT the following Python function can be used:

   from rsgislib import rastergis
   clumps='segmentation_raster.kea'
   pyramids=True
   colourtable=True
   rastergis.populateStats(clumps, colourtable, pyramids)
   

   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.

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

Image segmentation & attribution utilities in RSGISLib

Included with RSGISLib are two command line tools to segment an image, and attribute each segment:

# Segmentation
rsgislibsegmentation.py --input jers1palsar_stack.kea \
--output jers1palsar_stack_clumps_elim_final.kea \
--outmeanimg jers1palsar_stack_clumps_elim_final_mean.kea \
-tmpath $PWD --numclusters 100 --minpxls 100

# Attribute segments
rsgislibattributerat.py --inimage jers1palsar_stack.kea \
--inclumps jers1palsar_stack_clumps_elim_final.kea \
--mean

To populate the image statistics, band names are used (where available), these can be set using the ‘setbandnames.py’ script from RSGIS Scripts.

These command line tools use the Python utility functions ‘runShepherdSegmentation‘ from segutils and ‘populateImageStats‘ from ratutils. These utility functions can be called directly from Python:

from rsgislib.rastergis import ratutils
from rsgislib.segmentation import segutils
from rsgislib import imageutils

inputImage = 'jers1palsar_stack.kea'
clumpsFile = 'jers1palsar_stack_clumps_elim_final.kea'
meanImage = 'jers1palsar_stack_clumps_elim_final_mean.kea'

# Set band names
bandNames = ['98_summer','98_winter','07_HH','07_HV']
imageutils.setBandNames(inputImage, bandNames)

# Run segmentation
segutils.runShepherdSegmentation(inputImage, clumpsFile,
                    meanImage, numClusters=100, minPxls=100)

# Attribute segments
ratutils.populateImageStats(inputImage, clumpsFile,
                    calcMean=True)

Note the latest version of RSGISLib (2.1.752) is required for this. Older versions don’t include the ‘setBandNames’ function and require all parameters to be set in the utility functions.