Pages

Thursday, December 2, 2010

Dot Density Maps with Python and OGR

If you use Python for GIS sooner or later you'll use GDAL for manipulating raster data and its vector cousin OGR for working with vector data. OGR has a Python API for most of the methods in the C++ library and even provides some basic geometry analysis. And most importantly it can read/write and therefore convert data in a variety of vector file and database formats.

OGR provides a fast way to create dot density maps. A dot density maprepresents statistical information about an area as mathematically distributed points. Areas with higher values have a higher concentration of points. This is one of my favorite types of maps because it is a great example of GIS - visualizing geographic data in a way that is instantly comprehensible.

I'm using OGR in this example because it can read and write shapefiles. But unlike the Python Shapefile Library it can also perform basic geometry operations needed for this sample. Most GIS programs would display the population information on some type of memory layer instead of actually outputting a shapefile for the density layer as demonstrated here. But we're going to keep things simple for this example and just create a shapefile.

Assuming you have Python installed, here are some basic gdal/ogr installation instructions.
1. Go to http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries and download the gdal binary for your platform
2. Extract the directory to your hard drive
3. Add the "bin" directory within the gdal folder to your system shell path
4. Set the path to the "data" directory in the gdal folder to an environment variable called "GDAL_DATA"
5. Install the appropriate python module for your Python version and platform from here: http://pypi.python.org/pypi/GDAL/1.6.0#downloads

There is no error handling in this sample so if you run it multiple times you need delete the output dot density shapefile.

Note that this type of rendering only works when you have one polygon representing each data value. For example you couldn't do this operation with a world country boundary shapefile because islands like Hawaii associated with a country would force an inaccurate representation. For that type of map you need to use a choropleth map.

Also note that when you use OGR for shapefile editing you must specify a "layer" after opening a file. This extra step is necessary because OGR handles dozens of formats, some of which are layered vector formats such as DWG using the same API. Also because OGR is a wrapped C library you have to adjust to explicitly destroying objects and extreme camel casing on method calls usually not found in Python.

OGR and the raster equivalent GDAL are two very powerful libraries which dominate the open source geospatial world. They are also included in several well-known commercial packages thanks to the commercial-friendly MIT license.

3 comments:

Hopefully you will gimme some tips and/or insights. I have hundreds of shapefiles, all of them contain only polygons. I need to extract the attributes at fixed distances, say every 1000 meters. It is a bit like rasterize: add a layer of points on top of the polygons layer and get the attributes at each point. Unless I am missing something, neither your nice library or the OGR one let you read the attributes at specific locations, I mean is there a way to know the attributes of this or that shapefile at this lat and lon?

And when you're not online...

Search GeospatialPython.com

@SpatialPython on Twitter

About Me

My name is Joel Lawhead and I'm the CIO and lead Python geek at NVisionSolutions.com. NVision is a geospatial technology and engineering company on the Mississippi Gulf Coast at the NASA Stennis Space Center.
I have been using Python for 13 years including authoring "Learning Geospatial Analysis with Python" from Packt Publishing and contributions to the O'Reilly "Python Cookbook" and have worked in the geospatial industry for 13 years. I actually started out in the newspaper industry with a degree in Journalism.