I want to segment RGB images for land cover using k means clustering in such a fashion that the different regions of the image are marked by different colors and if possible boundaries are created separating different regions. I want something like :

enter image description here

from this :

enter image description here

Is it possible to achieve this by K-means clustering? I have been searching all over internet and many tutorials do it by k means clustering but only after converting the image to grey scale. I want to do it with an RGB image only. Is there any source that could help me begin with it? Please suggest something.

  • Hi, try this link. I tried it some time ago, but only had limited success. Maybe you can get it to work a bit better. Good luck. opencv-python-tutroals.readthedocs.org/en/latest/py_tutorials/…
    – Jcstay
    Jul 1, 2015 at 7:18
  • Hi,thank You for your suggestion @Jcstay but i have already tried the link and it did not help. Thank You though.
    – rach
    Jul 1, 2015 at 7:29
  • 3
    I would point out that the K-means algorithm, like all other clustering methods, needs and optimal fit of k. Since everything in the reference data will get assigned a class, if k is not optimized, the results can be erroneous with no support for a resulting class. In these cases, a given class can represent nothing other than noise or marginal effect in the data. Commonly, margin silhouette values are used to select an optimal k. Jul 1, 2015 at 15:25

2 Answers 2


I hacked together a solution for this and wrote a blog article a while back on a very similar topic, which I will summarize here. The script is intended to extract a river from a 4-band NAIP image using an image segmentation and classification approach.

  1. Convert image to a numpy array
  2. Perform a quick shift segmentation (Image 2)
  3. Convert segments to raster format
  4. Calculate NDVI
  5. Perform mean zonal statistics using segments and NDVI to transfer NDVI values to segments (Image 3)
  6. Classify segments based on NDVI values
  7. Evaluate results (Image 4)

This example segments an image using quickshift clustering in color (x,y) space with 4-bands (red, green, blue, NIR) rather than using K-means clustering. The image segmentation was performed using the scikit-image package. More details on a variety of image segmentation algorithms in scikit-image here. For convenience sake, I used arcpy to do much of the GIS work, although this should be pretty easy to port over to GDAL.

enter image description here

from __future__ import print_function

import arcpy

import matplotlib.pyplot as plt
import numpy as np
from skimage import io
from skimage.segmentation import quickshift

# The input 4-band NAIP image
river = r'C:\path\to\naip_image.tif'

# Convert image to numpy array
img = io.imread(river)

# Run the quick shift segmentation
segments = quickshift(img, kernel_size=3, convert2lab=False, max_dist=6, ratio=0.5)
print("Quickshift number of segments: %d" % len(np.unique(segments)))

# View the segments via Python

# Get raster metrics for coordinate info
myRaster = arcpy.sa.Raster(river)

# Lower left coordinate of block (in map units)
mx = myRaster.extent.XMin
my = myRaster.extent.YMin
sr = myRaster.spatialReference

# Note the use of arcpy to convert numpy array to raster
seg = arcpy.NumPyArrayToRaster(segments, arcpy.Point(mx, my),

outRaster = r'C:\path\to\segments.tif'
seg_temp = seg.save(outRaster)
arcpy.DefineProjection_management(outRaster, sr)

# Calculate NDVI from bands 4 and 3
b4 = arcpy.sa.Raster(r'C:\path\to\naip_image.tif\Band_4')
b3 = arcpy.sa.Raster(r'C:\path\to\naip_image.tif\Band_3')
ndvi = arcpy.sa.Float(b4-b3) / arcpy.sa.Float(b4+b3)

# Extract NDVI values based on image object boundaries
zones = arcpy.sa.ZonalStatistics(outRaster, "VALUE", ndvi, "MEAN")

# Classify the segments based on NDVI values
binary = arcpy.sa.Con(zones < 20, 1, 0)
  • 2
    This is a fantastic solution and sidesteps some of the issues with k-means and finding an optimal k. Jul 1, 2015 at 15:26
  • This is very nice, great work!!
    – Jcstay
    Jul 2, 2015 at 8:41

You could look at clustering in scikit-learn. You will need to read the data into numpy arrays (I'd suggest rasterio) and from there you can manipulate the data so that each band is a variable for classification. For example, assuming you have the three bands read into python as red, green, and blue numpy arrays:

import numpy as np
import sklearn.cluster

original_shape = red.shape # so we can reshape the labels later

samples = np.column_stack([red.flatten(), green.flatten(), blue.flatten()])

clf = sklearn.cluster.KMeans(n_clusters=5)
labels = clf.fit_predict(samples).reshape(original_shape)

import matplotlib.pyplot as plt


Note that the KMeans clustering doesn't take into account connectivity within the dataset.

  • +1 Great answer. It would be especially nice to show an example of converting color images to numpy arrays using rasterio;)
    – Aaron
    Jul 1, 2015 at 14:57
  • 1
    @Aaron Thanks! I've posted a slightly longer example including reading data using rasterio.
    – om_henners
    Jul 2, 2015 at 2:17
  • @om_henners your solution is wonderful but I have a question. The segmented image returned by your program using k means clustering is 2D. Now I need to calculate dice similarity coefficient between the original image( 3D image before splitting into R,G,B bands) and the segmented image but that needs the two to have same dimensions. How do I solve this problem?
    – rach
    Jul 4, 2015 at 10:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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