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.
http://nitace.bsd.uchicago.edu:8080/files/share/map.png
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 -
http://sourceforge.net/mailarchive/message.php?msg_id=6542965
Here's some example code:
# shapefile from
# http://edcftp.cr.usgs.gov/pub/data/nationalatlas/statesp020.tar.gz
# shapefile lib from http://thuban.intevation.org
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 = shp.info()
ax = gca()
dpi = ax.dpi
bbox = ax.bbox
transx = ax.xaxis.transData
transy = ax.yaxis.transData
left=; right=; bottom=; top=
db = dbflib.open(dbfile)
# 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
else:
seen[state] = area, verts
for state, tup in seen.items():
area, verts = tup
poly = Polygon(dpi, bbox, verts,
fill=False,
transx=transx, transy=transy)
x = [tx for tx, ty in verts]
y = [ty for tx, ty in verts]
ax.xaxis.datalim.update(x)
ax.yaxis.datalim.update(y)
ax.add_patch(poly)
set(gca(), 'xlim', ax.xaxis.datalim.bounds())
set(gca(), 'ylim', ax.yaxis.datalim.bounds())
savefig('map', dpi=150)
axis('Off')
show()