Tag Archives: KEA

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

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:

casi_data_in_matlab

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.