How to perform image segmentation on 4-band geotiff using Python's scikit-image?

Question

I am attempting to read a 4-band (red, green, blue, near-infrared) geotiff (example data) and perform a quickshift segmentation using the scikit-image module in Python.

I have created the following script (based on the scikit example):

from __future__ import print_function
from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np

from skimage.segmentation import felzenszwalb, slic, quickshift
from skimage.segmentation import mark_boundaries
from skimage.util import img_as_float

image = r'C:\path\to\my\geotiff.tif'
img = io.imread(image, as_grey=False, plugin="gdal")
segments_quick = quickshift(img, kernel_size=3, max_dist=6, ratio=0.5)

I get the following error:

ValueError: the input array must be have a shape == (.., ..,[ ..,] 3)), got (4, 436, 553)

I am pretty sure the numpy array needs to be reshaped somehow. How can I properly read the multiband geotiff into a numpy array and perform the image segmentation?

Answer 1

I believe your problem is that quickshift() thinks your image is rgb. I downloaded a random image from the link you provided and read it into skimage.

img = io.imread('./m_4111722_ne_11_1_20100704.tif')

I resized it to 128x128x4 (to make computation easy)

img = transform.resize(img, (128, 128, 4))

then ran quickshift()

segments = quickshift(img, kernel_size=3, max_dist=6, ratio=0.5)

and got the same error.

ValueError: the input array must be have a shape == (.., ..,[ ..,] 3)), got (128, 128, 4)

Higher up in the stack trace it says

skimage/segmentation/_quickshift.pyx inskimage.segmentation._quickshift.\
quickshift (skimage/segmentation/_quickshift.c:1710)()

/****/****/anaconda/lib/python2.7/site-packages/skimage/color/colorconv.pyc in rgb2lab(rgb)
    901     This function uses rgb2xyz and xyz2lab.
    902     """
--> 903     return xyz2lab(rgb2xyz(rgb))

So you can see _quickshift.pyx is trying to convert rgb --> xyz and then xyz --> lab. So its assuming your image is rgb. The skimage docs for quickshift() shows it has a flag convert2lab that defaults to True.

convert2lab : bool, optional (default True) Whether the input should be converted to Lab colorspace prior to segmentation. For this purpose, the input is assumed to be RGB.

If I re-run your function with that flag set to False

segments = quickshift(img, kernel_size=3, convert2lab=False, max_dist=6, ratio=0.5)

it runs.

plt.imshow(segments);

Edit:

Just as an aside, I noticed your image shape is (4, 436, 553) which is also problematic. skimage expects the color channel to be last. This can be remedied with

img = img.transpose(1, 2, 0)