Tag Archives: RAT

Calculations on large Raster Attribute Tables using ratapplier

The RIOS Library (http://rioshome.org) provides two methods of manipulating columns within a Raster Attribute Table (RAT):

  1. The ‘rat’ module can read an entire column to memory.
    from osgeo import gdal
    from rios import rat
    
    # Open RAT dataset
    rat_dataset = gdal.Open("clumps.kea", gdal.GA_Update)
    
    # Get columns with average red and NIR for each object
    red = rat.readColumn(rat_dataset, "RedAvg")
    nir = rat.readColumn(rat_dataset, "NIR1Avg")
    
    ndvi = (nir - red) / (nir + red)
    
    # Write out column
    rat.writeColumn(rat_dataset, "NDVI", ndvi)
    
    # Close RAT dataset
    rat_dataset = None
    
  2. 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.
    from rios import ratapplier
    
    def _ratapplier_calc_ndvi(info, inputs, outputs):
        """
        Calculate NDVI from RAT.
    
        Called by ratapplier
        """
        # Get columns with average red and NIR for each object
        # within block
        red = getattr(inputs.inrat, "RedAvg")
        nir = getattr(inputs.inrat, "NIR1Avg")
    
        # Calculate NDVI
        ndvi = (nir - red) / (nir + red)
    
        # Save to 'NDVI' column (will create if doesn't exist)
        setattr(outputs.outrat,"NDVI", ndvi)
    
    if __name__ == "__main__":
    
        # Set up rat applier for input / output
        in_rats = ratapplier.RatAssociations()
        out_rats = ratapplier.RatAssociations()
    
        # Pass in clumps file
        # Same file is used for input and output to write
        # to existing RAT
        in_rats.inrat = ratapplier.RatHandle("clumps.kea")
        out_rats.outrat = ratapplier.RatHandle("clumps.kea")
    
        # Apply function to all rows in chunks
        ratapplier.apply(_ratapplier_calc_ndvi, in_rats, out_rats)
    

Although using ratapplier looks slightly more complicated at first writing scripts to use it rather than the ‘rat’ interface means they will scale much better to larger datasets.

Both these examples are available to download from https://bitbucket.org/snippets/danclewley/7E666

Advertisements

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
  2. 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
    
  3. Define training data
  4. 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
    
  5. Stack bands
  6. 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.

  7. Image segmentation
  8. 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.

  9. Attribute segments
  10. #!/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)
    
  11. Perform classification
  12. 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

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.

Convert SSURGO soil data to a Raster Attribute Table

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

To do this the following steps are required:

  1. Use gdal_rasterize to create a raster.
  2. Use RSGISLib to convert to a RAT.
  3. Add a column for each attribute using RIOS.

A Python script to perform these steps can be downloaded from https://github.com/MiXIL/SSURGO-Utilities.

An example of usage is:

python convertSSURGO2RAT.py --indir CA669 \
  --colname claytotal_ \
  --outformat KEA

A list of all available columns can be viewed using:

python convertSSURGO2RAT.py --printcols

To export the columns as a standard raster (using RSGISLib) pass in the ‘–raster’ flag.

SSURGO data can be downloaded from the USDA NRCS Geospatial Data Portal (http://datagateway.nrcs.usda.gov/)

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.