Re: [matplotlib-devel] image module questions

[John H. <jdh...@ac...>]

>>>>> "John" == John N S Gill <jn...@eu...> writes:

    John> Quite often I find myself with numbers for different parts
    John> of the world that I want to map, shading regions different
    John> colours according to the numbers.

    John> The problem I've tended to run into with mapping projects
    John> has been getting shape files that aren't distributed under a
    John> 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()