Basic Image Change Detection in QGIS and Python

Change detection allows you to automatically highlight the differences between two images in the same area if they are properly orthorectified. We'll do a simple difference change detection on two images, which are several years apart, to see the differences in urban development and the natural environment.  Change detection algorithms can become quite sophisticated.  However at the most basic level, you can do a simple, literal mathematical difference between the pixel values in the two images.  That basic image difference is exactly what we'll do in this example.   This recipe is from my book, the"QGIS Python Programming Cookbook".   I show how do the same operation using the gdalnumeric module outside of QGIS in my other book "Learning Geospatial Analysis with Python".

We're going to start with two aerial images taken several years apart.  We'll subtract the before image from the after image and and visualize the result using a color ramp.

The before and after aerial photos we'll use for this recipe.
Click for  a larger version.

Getting ready

You will need a data directory to store the source images.  I usually create a directory named qgis_data for these types of experiments.  You can download the two images used in this example from https://github.com/GeospatialPython/qgis/blob/gh-pages/change-detection.zip?raw=true and put them in a directory named change-detection in the rasters directory of your qgis_data directory. Note that the file is 55 megabytes, so it may take several minutes to download.  Also, we'll step through this script a few lines at a time.  It's meant to be run in the QGIS Python Console.  You can download the entire script here: https://github.com/GeospatialPython/qgis/raw/master/4985_08_15-change.py

How to do it...

We'll use the QGIS raster calculator to subtract the images in order to get the difference, which will highlight significant changes. We'll also add a color ramp shader to the output in order to visualize the changes. To do this, we need to perform the following steps:

First, we need to import the libraries that we need in to the QGIS console:

from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.analysis import *

Now, we'll set up the path names and raster names for our images:

before = "/Users/joellawhead/qgis_data/rasters/change-detection/before.tif"
|after = "/Users/joellawhead/qgis_data/rasters/change-detection/after.tif"
beforeName = "Before"
afterName = "After"

Next, we'll establish our images as raster layers:

beforeRaster = QgsRasterLayer(before, beforeName)
afterRaster = QgsRasterLayer(after, afterName)

Then, we can build the calculator entries:

beforeEntry = QgsRasterCalculatorEntry()
afterEntry = QgsRasterCalculatorEntry()
beforeEntry.raster = beforeRaster
afterEntry.raster = afterRaster
beforeEntry.bandNumber = 1
afterEntry.bandNumber = 2
beforeEntry.ref = beforeName + "@1"
afterEntry.ref = afterName + "@2"
entries = [afterEntry, beforeEntry]

Now, we'll set up the simple expression that does the math for remote sensing:

exp = "%s - %s" % (afterEntry.ref, beforeEntry.ref)

Then, we can set up the output file path, the raster extent, and pixel width and height:

output = "/Users/joellawhead/qgis_data/rasters/change-detection/change.tif"
e = beforeRaster.extent()
w = beforeRaster.width()
h = beforeRaster.height()

Now, we perform the calculation:

change = QgsRasterCalculator(exp, output, "GTiff", e, w, h, entries)
change.processCalculation()

Finally, we'll load the output as a layer, create the color ramp shader, apply it to the layer, and add it to the map, as shown here:

lyr = QgsRasterLayer(output, "Change")
algorithm = QgsContrastEnhancement.StretchToMinimumMaximum
limits = QgsRaster.ContrastEnhancementMinMax
lyr.setContrastEnhancement(algorithm, limits)
s = QgsRasterShader()
c = QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
i = []
qri = QgsColorRampShader.ColorRampItem
i.append(qri(0, QColor(0,0,0,0), 'NODATA'))
i.append(qri(-101, QColor(123,50,148,255), 'Significant Itensity Decrease'))
i.append(qri(-42.2395, QColor(194,165,207,255), 'Minor Itensity Decrease'))
i.append(qri(16.649, QColor(247,247,247,0), 'No Change'))
i.append(qri(75.5375, QColor(166,219,160,255), 'Minor Itensity Increase'))
i.append(qri(135, QColor(0,136,55,255), 'Significant Itensity Increase'))
c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
lyr.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(lyr)

You should end up with the following image.  As a general rule of thumb, green shows area of change where something was built or trees were removed, and purple areas show where a road or structure was likely moved.  There is lots of noise and false positives.  But the algorithm does very well in areas of significant change.

BAM! Change detection! Click for larger version

How it works...

If a building is added in the new image, it will be brighter than its surroundings. If a building is removed, the new image will be darker in that area. The same holds true for vegetation, to some extent.

Here is a close up of an area of intense change that demonstrates what we're talking about.  In the following image, a single softball field was expanded into an entire softball field complex.  In the upper left of the image you can see where the trees were cleared and bright-green grass and dirt infields were added.  In the purple area in the lower right of the image you can see where the original field was rotated.  Further in that direction, you can see some of the noise where the intensity changed between the images but it doesn't have much meaning. However if you zoomed out, the areas of significant change stand out very clearly which is usually the point of this kind of analysis:

Looking at these softball fields close up really
shows what's happening. Click for larger version

Summary

The concept is simple. We subtract the older image data from the new image data. Concentrating on urban areas tends to be highly-reflective and results in higher image pixel values.

Clipping a Shapefile in Pure Python

In a previous post, I showed how to use pure Python to subset a shapefile by attributes.  That method is fast if the attributes define the data you want to clip.  But sometimes you need to clip data spatially.

In this post, we'll clip a road shapefile using bounding boxes.  This is the simplest type of clip.  In a future post, I'll demonstrate how to use irregular polygons which is more complex.  However the ingredients I'll use for that recipe are already on this blog (this post + PyShp + the ray-casting point-in-polygon algorithm).

Here's our goal for this example. We have a shapefile with the major roads of the United States which you can download a sample dataset here. We want to use a simple bounding box to subset the roads for Puerto Rico. Here's a view of part of that shapefile including Florida and Puerto Rico in the lower right of the image.  It gives you a good idea of the density and number of records in this data.

Each road in the shapefile has its own shape record with the bounding box of the road in the record header.  We'll just do a simple bounding box comparison to isolate the roads within our desired bounding box and then write those results to a new shapefile.  You'll also need the latest version of PyShp which provides the newer iterShapeRecords() method to provide and efficient iterator for both the shapes and the dbf attributes in a single object:

import shapefile

# Create a reader instance for our US Roads shapefile
r = shapefile.Reader("roadtrl020") 

# Create a writer instance copying the reader's shapefile type
w = shapefile.Writer(r.shapeType)

# Copy the database fields to the writer
w.fields = list(r.fields)

# Our selection box that contains Puerto Rico
xmin = -67.5
xmax = -65.0
ymin = 17.8
ymax = 18.6

# Iterate through the shapes and attributes at the same time
for road in r.iterShapeRecords():
    # Shape geometry
    geom = road.shape
    # Database attributes 
    rec = road.record
    # Get the bounding box of the shape (a single road)
    sxmin, symin, sxmax, symax = geom.bbox
    # Compare it to our Puerto Rico bounding box.
    # go to the next road as soon as a coordinate is outside the box
    if sxmin <  xmin: continue
    elif sxmax > xmax: continue
    elif symin < ymin: continue
    elif symax > ymax: continue
    # Road is inside our selection box.
    # Add it to the new shapefile
    w._shapes.append(geom)
    w.records.append(rec)
    
# Save the new shapefile! (.shp, .shx, .dbf)
w.save("Puerto_Rico_Roads")

And here's our result!

Yes, there are dozens of tools and libraries that will do the same thing.  But this is the only solution in the entire (known) universe that is implemented in pure Python.  You don't need any compiled C code.  If you have the ability to run Python 2 or 3, you can run this without any additional system permissions.  It will run where most other libraries won't and when that's what your requirements are, nothing else will do!

Geolocating Photos in QGIS

Photos taken with GPS-enabled cameras, including smartphones, store location information in the header of the file, in a format called EXIF tags. These tags are largely based on the same header tags used by the TIFF image standard. If you have a directory full of photos, you can easily sort them by date or by name.  But if you want to sort them by location, only a map will do.  In this recipe, we'll use these tags to create locations on a map in QGIS for some photos and provide links to open them. 

Getting ready

You need to make sure you have QGIS installed as well as the Python Imaging Library (PIL) for the Python distribution that QGIS uses.  QGIS should already have PIL included. You can download some sample geotagged images here or you can use any photo taken outdoors with your smartphone.  Unzip these photos and place them in a directory accessible to QGIS.

How to do it...

PIL can parse EXIF tags. We will gather the file names of the photos, parse the location information, convert it to decimal degrees, create the point vector layer, add the photo locations, and add a QGIS action link to the attributes. To do this, we need to perform the following steps:

In the QGIS Python Console, import the libraries that we'll need, including the Image module (PIL) and its ExifTags module for parsing image data and the glob module for doing wildcard file searches:

import glob
import Image
from ExifTags import TAGS

Next, we'll create a function that can parse the header data:

def exif(img):
   exif_data = {}
   try:  
       i = Image.open(img)
       tags = i._getexif()
       for tag, value in tags.items():
         decoded = TAGS.get(tag, tag)
           exif_data[decoded] = value
   except:
       pass
return exif_data

Now, we'll create a function that can convert degrees-minute-seconds to decimal degrees, which is how coordinates are stored in JPEG images:

def dms2dd(d, m, s, i):
   sec = float((m * 60) + s)
   dec = float(sec / 3600)
   deg = float(d + dec)
   if i.upper() == 'W':
       deg = deg * -1
   elif i.upper() == 'S':
       deg = deg * -1
   return float(deg)

Next, we'll define a function to parse the location data from the header data:

def gps(exif):
   lat = None
   lon = None
   if exif['GPSInfo']:      
       # Lat
       coords = exif['GPSInfo']
       i = coords[1]
       d = coords[2][0][0]
       m = coords[2][1][0]
       s = coords[2][2][0]
       lat = dms2dd(d, m ,s, i)
       # Lon
       i = coords[3]
       d = coords[4][0][0]
       m = coords[4][1][0]
       s = coords[4][2][0]
       lon = dms2dd(d, m ,s, i)
return lat, lon

Next, we'll loop through the photos directory, get the file names, parse the location information, and build a simple dictionary to store the information, as follows:

photos = {}
photo_dir = "/Users/joellawhead/qgis_data/photos/"
files = glob.glob(photo_dir + "*.jpg")
for f in files:
   e = exif(f)
   lat, lon = gps(e)
 photos[f] = [lon, lat]

Now, we'll set up the vector layer for editing:

lyr_info = "Point?crs=epsg:4326&field=photo:string(75)"
vectorLyr = QgsVectorLayer(lyr_info, \"Geotagged Photos" , "memory")
vpr = vectorLyr.dataProvider()

We'll add the photo details to the vector layer:

features = []
for pth, p in photos.items():
   lon, lat = p
   pnt = QgsGeometry.fromPoint(QgsPoint(lon,lat))
   f = QgsFeature()
   f.setGeometry(pnt)
   f.setAttributes([pth])
   features.append(f)
vpr.addFeatures(features)
vectorLyr.updateExtents()

Now, we can add the layer to the map and make the active layer:

QgsMapLayerRegistry.instance().addMapLayer(vectorLyr)
iface.setActiveLayer(vectorLyr)
activeLyr = iface.activeLayer()

Finally, we'll add an action that allows you to click on it and open the photo.  Actions are a simple yet powerful feature of QGIS that let you define the result of certain user actions like clicking on features:

actions = activeLyr.actions()
actions.addAction(QgsAction.OpenUrl, "Photos", \'[% "photo" %]')

How it works...

Using the included PIL EXIF parser, getting location information and adding it to a vector layer is relatively straight forward. This action is a default option for opening a URL. However, you can also use Python expressions as actions to perform a variety of tasks. The following screenshot shows an example of the data visualization and photo popup:

There's more...

You can download the complete code sample here. Another plugin called Photo2Shape is available that performs a similar function, but it requires you to install an external EXIF tag parser while this approach uses the parser in PIL included with QGIS.