Getting started with Python-based open source geoprocessing

Well, I've had this Blog configured for some time and never found the inspiration that I needed to get started until today. In my profession, I seem to be constantly evaluating new geospatial technologies to find the right tool for the right job and as it would happen, many of these tools lack documentation so I am left to my own vices. Often I can find a good tutorial or some example code but sometimes I can't. My new-found inspiration for this blog is to share my misery with others and hopefully help others find the information they need to make progress.

This week's exploration is all about the SWIG-based Python bindings for GDAL or the Geospatial Data Abstraction Library (http://gdal.org/). GDAL is a miracle tool for handling raster-based spatial data. You can translate data between various formats such as GeoTIFF, GRID, etc..., and you can also reproject the rasters to a new coordinate system. Combining GDAL and Python is a natural marriage because it allows you to script these processes to simplify repetitive tasks. The main problem is that there is little documentation for beginners such as myself. The following is a simple script that I have been working through slowly to learn the ins and outs of the Python GDAL interface. I know it must be painfully simple for some but this is how I work things out. As I learn more, I'll keep posting things here so that others can learn along with me.

# Import the gdal and osr modules
from osgeo import gdal,osr

# GDAL works on the concept of the "Driver" which automagically recognizes the type of file you are trying to open.
# The simplest method is to register all the drivers so that it has the maximum amount of options for files that you are opening.

gdal.AllRegister()

# Open a new GDAL dataset
dataset = gdal.Open("c://tmp//gdal//usgsicurrent.bil")
if dataset is None:
print "Could not open dataset"

# Print some information about the file
print 'Driver: ', dataset.GetDriver().ShortName,'/',dataset.GetDriver().LongName
print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize,'x',dataset.RasterCount
print 'Projection is ',dataset.GetProjection()

# Get some information about the image like the corner coordinates and cell sizes
geotransform = dataset.GetGeoTransform()
if not geotransform is None:
print 'Origin = (',geotransform[0], ',',geotransform[3],')'
print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'

xOrigin = geotransform[0]
yOrigin = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

# Create new Spatial References using OSR
outproj = osr.SpatialReference()
ingeog = osr.SpatialReference()

# Set the output projection to the US National Atlas Lambert Equal Area Azimuthal Projection
outproj.ImportFromEPSG(2163)

# Set the input projection to WGS84
ingeog.SetWellKnownGeogCS("WGS84")

# Use the ExportToWkt() to set the input raster's coordinate system using the spatial reference created earlier.
dataset.SetProjection(outproj.ExportToWkt())

# Create a new Coordinate Transformation object
ct = osr.CoordinateTransformation(ingeog, outproj)

# Transform a point from Geographic to Projected Coordinates
(MyX,MyY,MyZ) = ct.TransformPoint(-100.0,45.00)

# Compute the offset
xOffset = int((MyX - xOrigin) / pixelWidth)
yOffset = int((MyY - yOrigin) / pixelHeight)

# Read the raster value for the tranformed coordinate
band = dataset.GetRasterBand(1) # 1-based index
data = band.ReadAsArray(xOffset, yOffset, 1, 1)
value = data[0,0]

# Print the value
print value

I know, it doesn't seem like much but this little program has many uses and the individual components are are useful are the whole. As I learn more and have more examples, I'll continue to share them.

This site helped with the coordinate system tranformation: http://help.maptiler.org/coordinates/batch-transformation

Here is a GREAT tutorial on Python-based geoprocessing using GDAL/OGR: http://www.gis.usu.edu/~chrisg/python/

Here is the official GDAL API tutorial with some information about the Python implementation: http://www.gdal.org/gdal_tutorial.html