Raster Images: Extracting Features of Interest Based on Color

While I am always interested in remote sensing, my focus has been on GIS/GPS data and applications, largely due to the nature of my day-to-day routines at work which include GIS analysis and map production using almost exclusively vector data.

Working with Google Earth the other day, I wondered how easy it would be to select an area of a raster image based on certain color, i.e. to select or identify features characterized by certain color ranges of each raster bands. A simple quality control measure would be a comparison of the raster selection with a shapefile with point data for a few sample features.

Since I don’t have a Spatial Analyst license, I tried to attack this with only ArcView…

Step 1
Find a raster file and add individual bands as separate layers to your map. In my case, this was a (more or less random) 6BM (approx. 3000x2000px) GeoTIFF with 3 bands.

Step 2
Find the color range in each band that is typical for the type of feature. – By default, one-band rasters will be displayed as grayscale images, using a stretched color ramp. A neat little tool I wasn’t familiar with before is the Pixel Inspector. By moving your cursor over the raster, you get a look at the 8-bit color value of each pixel and thus a sense of the range values for your feature.

With a bit of trial and error, you can narrow down the range of pixel values that best outlines the features you’re interested in, and by using the classified color ramp on the Symbology tab, you’re able to check what areas will be included when using as certain range of values.

Step 3
Now it’s time to create a new raster file that contains only those colors making up the features of interest. Of course, there will be some other noise of similarly colored features. But I am hoping that upon merging the three bands, only those areas that are part of all three rasters will be left over.

Turns out that to turn my original one-band raster into one containing only the pixels of interest, I need the Raster Calculator, which, in understand, is a feature of Spatial Analyst.

Time to look for Open Source help. And I quickly found OpenEv, which is part of FWTools, used GDAL, has Python bindings, and provides all kinds of raster manipulation capabilities. All I had to do was open my one-band raster and use some Python to run my own – VERY CRUDE – raster calculation. I don’t know much about numpy (Python’s module for all kinds of numerical/mathematical tricks) but this was pretty straight forward and goes something like this:

#Band1.img is the single band raster saved in ArcGIS
#columns and rows are the number of columns/rows respectively
# 100 and 180 are the min/max pixel values of interest

array1 = LoadFile(`Band1.img`)

for c in range(columns):
    for r in range(rows):
        pixel = array1[c:c+1,r:r+1]
        if pixel < 100 or pixel > 180:
             array1[c:c+1,r:r+1]=0

display(array1)

Now, I can export array1 as another GeoTIFF and open it in ArcGIS. Basically, I get a raster where the pixels of interest are unchanged, and everything else has been set to zero or BLACK. That’s as far as I’ve gotten. Now, I just need to the same for all three bands and work on recombing them.

PS: I stumbled across another open source project that might help with raster file – TerraLib/TerraView. But I haven’t had a chance to try it out.

Advertisements

3 Comments

Filed under Uncategorized

3 responses to “Raster Images: Extracting Features of Interest Based on Color

  1. Dear Arne,

    I found the very good post! Work for some time with the TerraLib which has a fantastic potential, but has not yet TerraView tools for classification, I believe that soon we are implementing a tool for digital image processing, so you can test done with other projects as TerraLib the InterIMAGE, which supports object-oriented classification. It also has TerraPixel, very interesting project! It is worth studying TerraLib!

    On the topic of the post in my blog I have some tutorials and articles, whether your interest …

    http://geotecnologias.wordpress.com/

  2. gui

    I haven’t used OpenEV, but I use numpy quite a bit (along with other vectorial langage as R/matlab).

    According to the OpenEV documentation, the LoadFile function returns a numpy array. Your two imbricated loop can then be written :

    import numpy as np
    array[np.logical_or(array100180)] = 0

    Nice blog anyway and keep posting about your many experiments!

  3. gui

    I mixed up the writing… It should be :

    array[np.logical_or(array1180)] = 0

    this page is great to quickly have a grasp with numpy :

    http://www.scipy.org/Cookbook

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s