6

Is there an equivalent in GDAL to the Arcpy ISO data unsupervised classification tool, or a series of methods using GDAL/python that can accomplish this?

  • 1
    try this cs.umd.edu/~mount/Projects/ISODATA or the kMean from OTB (but the latter is not exactly the same as ISODATA). – radouxju Nov 10 '14 at 14:34
  • The Python package pyradar includes an ISOdata classifier in Python. GDAL, as with most Python scripts, is used to import and export your image data to/from numpy arrays. – Kersten Nov 10 '14 at 15:17
  • Thanks for the feedback, both of those look like good possibilities. – Brian Nov 11 '14 at 14:02
5

Scikit-learn has some excellent unsupervised classification/clustering algorithms. The batched K-means algorithm works quickly with large datasets.

Here is an example using the KEA file format. You will have to modify this slightly to work with whatever raster format you use.

import rsgislib
from rsgislib import imageutils
from osgeo import gdal
import numpy
import glob
import sys
from sklearn.cluster import MiniBatchKMeans

# Search current directory for KEA rasters
InputImages = glob.glob('*.kea')

# Define the number of spectral classes
SpectralClasses = 6

# Define format and datatype of output raster
gdalformat = 'KEA'
gdaldatatype = gdal.GDT_Byte

# Define the classifier
clf = MiniBatchKMeans(n_clusters=SpectralClasses, init='k-means++', max_iter=10, batch_size=10000, verbose=0, compute_labels=True, random_state=None, tol=0.0, max_no_improvement=100, init_size=2000, n_init=10, reassignment_ratio=0.05)

# Terminate the script if the returned list is empty
if not InputImages:
  sys.exit("Error: No input images provided.")

for Raster in sorted(InputImages):
    OutImage = Raster.replace('.kea','_Classified.kea')

    print("Reading " + Raster)
    Image = gdal.Open(Raster, gdal.GA_ReadOnly)

    # Set up empty list to hold data
    TestData = []

    # Read data from each band
    for band in range(Image.RasterCount):
        band += 1

        B = Image.GetRasterBand(band)

        Array = B.ReadAsArray()

        # Get shape of array
        Shape = numpy.ma.shape(Array)

        # Flatten to 1D array
        Array = Array.flatten()

        TestData.append(Array)

        del Array

    TestData = numpy.array(TestData, dtype=numpy.dtype('float32')) # Convert to float to prevent sklearn error/warning message 
    TestData = numpy.transpose(TestData)

    print("Performing K-means classification...")
    clf.fit(TestData, y=None)
    predictedClass = clf.predict(TestData)

    del TestData

    predictedClass = predictedClass + 1 #Add 1 to exclude zeros in output raster
    predictedClass = numpy.reshape(predictedClass, Shape) # Reshape the numpy array to match the original image

    # Create an output raster the same size as the input image
    driver = gdal.GetDriverByName(gdalformat)
    metadata = driver.GetMetadata()
    output = driver.Create(OutImage, Image.RasterXSize, Image.RasterYSize, 1, gdaldatatype)

    # Create projection info for the output raster
    output.SetProjection(Image.GetProjectionRef())
    output.SetGeoTransform(Image.GetGeoTransform()) 

    # Write classification to band 1
    output_band = output.GetRasterBand(1)
    output_band.WriteArray(predictedClass)

    # Close datasets
    output_band = None
    output = None
    Image = None
    del predictedClass

    # Build image overviews
    imageutils.popImageStats(OutImage, True, 0, True)
    print("Done." + '\n')

print("All images processed.")

You will also want to create a binary mask image so you can clip/remove nodata regions from the output. I have skipped this step for brevity.

If you don't want to code all of this yourself through Python, you can undertake the k-means and ISOdata classifications using RSGISLib. Type "rsgiscmd -h classification" in the terminal for a list of options.

2

Brian, you could ask this question at http://lists.osgeo.org/mailman/listinfo/gdal-dev/ and perhaps get a more specific target audience. Also, while not providing a GDAL way to accomplishment this page (http://www.yale.edu/ceo/Projects/swap/landcover/Unsupervised_classification.htm) gives a nice summary of the differences between KMEANS and ISO DATA perhaps providing some insight into how you might write a custom function similar to http://geoexamples.blogspot.com/2012/02/raster-classification-with-gdal-python.html

Best, Derek

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.