Tuesday, November 26, 2013

Extracting Features from Images

Editor's Note: Republished from "Learning Geospatial Analysis with Python" with permission from Packt Publishing. This post is a follow-on to the previous post "Image Classification with Numpy and GDAL"

The ability to classify an image leads us to another remote-sensing capability. Vector GIS data such as shapefiles are typically extracted from remotely-sensed images. Extraction normally involves an analyst clicking around each object in an image and drawing the feature to save it as data. But it is also possible with good remotely-sensed data and proper pre-processing to automatically extract features from an image.

For this example we'll take a subset of our Landsat 8 thermal image to isolate a group of barrier islands as seen in the following screenshot:

You can download this image here:

Our goal with this example is to automatically extract the three islands in the image as a shapefile. But before we can do that, we need to mask out any data we aren't interested in. For example, the water has a wide range of pixel values, as do the islands themselves. If we just want to extract the islands themselves, we need to push all pixel values into just two bins to make the image black and white. This technique is called thresholding . The islands in the image have enough contrast with the water in the background such that thresholding should isolate them nicely.

In the following script we will read the image into an array and then histogram the image using only two bins. We will then use the colors black and white to color the two bins. This script is simply a modified version of our classification script with a very limited output:

import gdalnumeric

# Input file name (thermal image)
src = "islands.tif"

# Output file name
tgt = "islands_classified.tiff"

# Load the image into numpy using gdal
srcArr = gdalnumeric.LoadFile(src)

# Split the histogram into 20 bins as our classes
classes = gdalnumeric.numpy.histogram(srcArr, bins=2)[1]

lut = [[255,0,0],[0,0,0],[255,255,255]]

# Starting value for classification
start = 1

# Set up the output image
rgb = gdalnumeric.numpy.zeros((3, srcArr.shape[0], \
# Process all classes and assign colors
for i in range(len(classes)):
    mask = gdalnumeric.numpy.logical_and(start <= srcArr, \ 
    srcArr <= classes[i])
    for j in range(len(lut[i])):
        rgb[j] = gdalnumeric.numpy.choose(mask, (rgb[j], \
    start = classes[i]+1 

# Save the image
gdalnumeric.SaveArray(rgb.astype(gdalnumeric.numpy.uint8), \
 tgt, format="GTIFF", prototype=src)
The output looks great as seen in the following screenshot:

The islands are clearly isolated so our extraction script will be able to identify them as polygons and save them to a shapefile. The GDAL library has a method called Polygonize() that does exactly that. It groups all sets of isolated pixels in an image and saves them out as a feature data set. One interesting technique we will use in this script is to use our input image as a mask. The Polygonize() method allows you to specify a mask that will use the color black as a filter that will prevent the water from being extracted as a polygon, so we'll end up with just the islands. Another area to note in the script is that we copy the georeferencing information from our source image to our shapefile to geolocate it properly:

import gdal
import ogr, osr

# Thresholded input raster name
src = "islands_classified.tiff"
# Output shapefile name
tgt = "extract.shp"
# OGR layer name
tgtLayer = "extract"
# Open the input raster
srcDS = gdal.Open(src)
# Grab the first band
band = srcDS.GetRasterBand(1)
# Force gdal to use the band as a mask
mask = band
# Set up the output shapefile
driver = ogr.GetDriverByName("ESRI Shapefile")
shp = driver.CreateDataSource(tgt)
# Copy the spatial reference
srs = osr.SpatialReference()
layer = shp.CreateLayer(tgtLayer, srs=srs)
# Set up the dbf file
fd = ogr.FieldDefn( "DN", ogr.OFTInteger )
dst_field = 0
# Automatically extract features from an image!
extract = gdal.Polygonize(band, mask, layer, dst_field, [], None)
The output shapefile is simply called extract.shp. We'll bring that script back here to look at our shapefile, but we'll add something to it. The largest island has a small lagoon which shows up as a hole in the polygon. In order to properly render it, we have to deal with parts in a shapefile record. The previous example using that script did not do that, so we'll add that piece as we loop through the shapefile features. The code comments in the following code outline the technique.  You can download the pure Python PNGCanvas module from the Tao of Mac website:

import shapefile
import pngcanvas
# Open the extracted islands
r = shapefile.Reader("extract.shp")
# Setup the world to pixels conversion
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 800
iheight = 600
xratio = iwidth/xdist
yratio = iheight/ydist
polygons = []
# Loop through all shapes
for shape in r.shapes():
  # Loop through all parts to catch
  # polygon holes!
  for i in range(len(
    pt = None
    if I < len(
      pt = shape.points[[i][i+1]]
      pt = shape.points[[i]:]
    for x,y in pt:
      px = int(iwidth - ((r.bbox[2] - x) * xratio))
      py = int((r.bbox[3] - y) * yratio)
# Set up the output canvas
c = pngcanvas.PNGCanvas(iwidth,iheight)
# Loop through the polygons and draw them
for p in polygons:
# Save the image
f = file("extract.png", "wb")
The following screenshot shows our automatically extracted island features! Commercial packages that do this kind of work can easily cost tens of thousands of dollars. While these packages are very robust, it is still fun to see how far you can get with simple Python scripts and a few open source packages. In many cases you can do everything you need to do.
 The western-most island contains the polygon hole as shown in the following screenshot, which is zoomed to that area:
Automated feature extraction is a holy grail within geospatial analysis because of the cost and tedious effort required to manually extract features. The key to feature extraction is proper image classification. Automated feature extraction works well with water bodies (and islands), roads, farm fields, buildings, and other features that tend to have high-contrast pixel values with their background.

Friday, November 22, 2013

Image Classification with Numpy and GDAL

Editor's Note:  This post is published with permission from Packt Publishing and originally published in my book. 

Automated Remote Sensing ( ARS ) is rarely ever done in the visible spectrum. The most commonly available wavelengths outside of the visible spectrum are infrared and near-infrared. The following scene is a thermal image (band 10) from a fairly recent Landsat 8 flyover of the US Gulf Coast from New Orleans, Louisiana to Mobile, Alabama. Major natural features in the image are labeled so you can orient yourself:

Because every pixel in that image has a reflectance value, it is information. Python can "see" those values and pick out features the same way we intuitively do by grouping related pixel values. We can colorize pixels based on their relation to each other to simplify the image and view related features. This technique is called classification. Classifying can range from fairly simple groupings based only on some value distribution algorithm derived from the histogram to complex methods involving training data sets and even computer learning and artificial intelligence. The simplest forms are called unsupervised classifications, whereas methods involving some sort of training data to guide the computer are called supervised. It should be noted that classification techniques are used across many fields, from medical doctors trying to spot cancerous cells in a patient's body scan, to casinos using facial-recognition software on security videos to automatically spot known con-artists at blackjack tables.
To introduce remote sensing classification we'll just use the histogram to group pixels with similar colors and intensities and see what we get. First you'll need to download the Landsat 8 scene here:
Instead of our histogram() function from previous examples, we'll use the version included with NumPy that allows you to easily specify a number of bins and returns two arrays with the frequency as well as the ranges of the bin values. We'll use the second array with the ranges as our class definitions for the image. The lut or look-up table is an arbitrary color palette used to assign colors to classes. You can use any colors you want.
Note: In the following code we use the gdalnumeric module which allows Numpy and GDAL to work together.  GDAL reads compatible image/data formats and converts them to Numpy arrays.  The Numpy module was formerly called Numeric. 
import gdalnumeric

# Input file name (thermal image)
src = "thermal.tif"

# Output file name
tgt = "classified.jpg"

# Load the image into numpy using gdal
srcArr = gdalnumeric.LoadFile(src)

# Split the histogram into 20 bins as our classes
classes = gdalnumeric.numpy.histogram(srcArr, bins=20)[1]

# Color look-up table (LUT) - must be len(classes)+1.
# Specified as R,G,B tuples 
lut = [[255,0,0],[191,48,48],[166,0,0],[255,64,64],
\ [0,133,0],[57,230,57],[103,230,103],[184,138,0]]

# Starting value for classification
start = 1

# Set up the RGB color JPEG output image
rgb = gdalnumeric.numpy.zeros((3, srcArr.shape[0],
srcArr.shape[1],), gdalnumeric.numpy.float32)
# Process all classes and assign colors
for i in range(len(classes)):
    mask = gdalnumeric.numpy.logical_and(start <= \
 srcArr, srcArr <= classes[i])
    for j in range(len(lut[i])):
      rgb[j] = gdalnumeric.numpy.choose(mask, (rgb[j], \ 
    start = classes[i]+1 

# Save the image    
gdalnumeric.SaveArray(rgb.astype(gdalnumeric.numpy.uint8), \
tgt, format="JPEG")

The following image is our classification output, which we just saved as a JPEG. We didn't specify the prototype argument when saving as an image, so it has no georeferencing information:

This result isn't bad for a very simple unsupervised classification. The islands and coastal flats show up as different shades of green. The clouds were isolated as shades of orange and dark blues. We did have some confusion inland where the land features were colored the same as the Gulf of Mexico. We could further refine this process by defining the class ranges manually instead of just using the histogram.

Thursday, November 7, 2013

Must See: Going Spatial with Python and Github

Check out this lightning talk by the talented Sara Safavi at the Austin, TX PyLadies Meetup:
Spatial Pythonista
Sara Safavi

And see the rest of her work at: 

Tuesday, November 5, 2013

Geospatial Python - the Book

I'm pleased to announce Packt Publishing officially released my new book "Learning Geospatial Analysis with Python".  Packt already has two other books related to geospatial Python.  All three books are featured on the right-side border of this page.  I went through a great deal of effort to compliment these two existing books with different but relevant information.

My goal for this book is to provide geospatial Python examples with the following priorities:
  1. Use Pure Python examples as much as possible (standard library only, or third party libraries that only use the standard library)
  2. Use PYPI-compatible libraries which can be installed using easy_install or pip
  3. If I must use Python bindings to a third-party compiled library, focus on methods which lack good examples elsewhere and avoid repeating other books and blogs
If you've ever tried to implement even a simple spatial algorithm, you'll find the vast majority of programmers rely on a relatively small set of core libraries.  These libraries are then built up into larger toolkits which tie those libraries together and add additional value.

The following diagram is a high-level summary of the geospatial software ecosystem, and thereby the Python geospatial ecosystem.  The "800-pound gorillas" including GEOS, PROJ.4 and GDAL are the foundation of nearly all other software.  In this diagram the CGAL and JTS libraries are different.  Both CGAL and JTS are the reference implementations for many of the spatial algorithms later used by GEOS:

The geospatial software industry consists of a small handful of open-source
libraries which support the vast majority of all other software, both proprietary
and open-source.

These libraries are:
  1. Actively maintained and free (open source!)
  2. Implement operations and data access complex enough that nobody wants to reinvent the wheel
I have used these libraries for years.  But I am also stubbornly curious.  I like to understand the algorithms inside the proverbial black box.  My favorite way to understand algorithms is through Python.  I do not have a computer science degree or a degree in mathematics.  I am ignorantly dangerous in almost any programming language.

But what I like about Python is the syntax. It is so clean and expressive it makes hard things easier.  I am not alone.  Scores of non-programmer biologists, meteorologists, economists, medical researchers, artists and other specialists have found Python a user-friendly way to model and solve both trivial and complicated real-world problems.  

But there is a large gap in the literature. For any given problem you can almost always find pure algorithms in research papers or implementations in C/C++ libraries.  But rarely can you find working, practical implementations in Python.  Of course many of the compiled libraries are open source but their structure is complex enough that the casual learner must invest a significant amount of time in understanding most of the code base to uncover a given algorithm.  

The most common reason cited for using compiled libraries versus pure Python is that interpreted languages like Python are far too slow to be efficient and useful for most geospatial operations.  I find this reason unsatisfying for the following reasons:
  1. It discredits Python as a learning tool.  There's no need for speed in learning.
  2. As server platforms become increasingly virtualized, multiplied, and abstracted, our view of resource constraints can continually be relaxed in exchange for faster development with more participants 
  3. There are many tools out there which make pure Python faster when you are ready to optimize
  4. As geospatial applications become more task specific, implementing a handful of algorithms instead of using a "kitchen-sink" toolkit becomes desirable
  5. As mobile devices become increasingly important, dragging around a large compiled tool chain becomes less desirable.
  6. When you bind Python to a compiled library, both of which are updated regularly you have an additional maintenance routine supporting multiple platforms.  This issue varies in difficulty but it's never easy.  This group excludes compiled libraries designed specifically for Python (Numpy, PIL, Fiona, Shapely).
  7. To me, Python's greatest value is as a meta-syntax.  The Python C interpreter alone runs on 30+ platforms.  But the Python C interpreter is often called a "reference implementation".  The Python interpreter has been implemented in dozens of other programming languages including Java (Jython) and .NET (IronPython), and Python itself (PyPy).  So why do people go through the trouble to port Python so many places?  Because the programmatic expression allowed by the syntax is just great.  I always liked the metaphor, "Python is executable pseudocode".  The excellent design of the language from a user standpoint is why nearly every geospatial platform, proprietary or open-source, has a Python API as opposed to other options.
I originally wrote the PyShp library to learn more about both Python and Shapefiles.  I kept it simple by releasing a single file using only the standard library.  Since its release PyShp has become extremely popular. It proved to me there is indeed a demand for spatial software written in pure Python without external dependencies.

"Learning Geospatial Analysis with Python" is my humble attempt at laying a plank in the bridge between abstract algorithms and as close to pure Python implementations as I could.  I took the word "learning" seriously.  I don't try to teach Python as a language because that's been done and done well already.  I do explore common geospatial concepts early on and move into algorithms later in the book which are more advanced (and useful!).

So if you're new to geospatial analysis the first few chapters try to ease you into the field from a software and data perspective before switching to Python.  If you're an experienced geospatial developer I included some interesting examples of virtually undocumented features of GDAL in Python as well as examples of using Numpy for operations that are normally handed off to GDAL or or SciPy.  This book contains more examples of Remote Sensing using Python than any book I've come across yet.

The book is available from Amazon, O'Reilly, and Barnes & Noble as both a paperback and all the major eBook formats.  I hope you find it useful and I hope it increases the use of Python for geospatial analysis and applications.