image module questions

Quite often I find myself with numbers for different parts

    > of the world that I want to map, shading regions different
    > colours according to the numbers.

    > The problem I've tended to run into with mapping projects
    > has been getting shape files that aren't distributed under a
    > restrictive licence.

I've experimented with this a bit using shapefiles from the national
atlas (I think they are distributed under a permissive license).
thuban has a nice python interface for reading shape files and their
associated db files. I've held off on pursuing this functionality
until we get an efficient path/polygon extension, which has been
discussed a number of times but is still on the TODO list.

The example below renders the US map from 2000 polygons, and is slow
for interactive use with that many polygons.

For nice map navigation, you would probably want a better navigation
toolbar (hand pan, zoom to region) which was discussed many moons ago
but has also languished due to lack of time -

Here's some example code:

# shapefile from

# shapefile lib from

import shapelib, dbflib, shptree
from matplotlib.patches import Polygon
from matplotlib.matlab import *
filename = 'statesp020.shp'
dbfile = 'statesp020.dbf'
shp = shapelib.ShapeFile(filename)

numShapes, type, smin, smax =
ax = gca()
dpi = ax.dpi
bbox = ax.bbox
transx = ax.xaxis.transData
transy = ax.yaxis.transData

left=[]; right=[]; bottom=[]; top=[]

db =

# just get the main polys for each state
seen = {}
for i in range(numShapes):
    rec = db.read_record(i)
    state = rec['STATE']
    area = rec['AREA']
    obj = shp.read_object(i)
    verts = obj.vertices()[0]
    if seen.has_key(state):
        have, tmp = seen[state]
        if area>have:
            seen[state] = area, verts
        seen[state] = area, verts
for state, tup in seen.items():
    area, verts = tup
    poly = Polygon(dpi, bbox, verts,
                   transx=transx, transy=transy)

    x = [tx for tx, ty in verts]
    y = [ty for tx, ty in verts]

set(gca(), 'xlim', ax.xaxis.datalim.bounds())
set(gca(), 'ylim', ax.yaxis.datalim.bounds())
savefig('map', dpi=150)