John: I found that if I just call proj with all the lats and

> lons at once (instead of once for each segment) I can speed

> it up tremendously. Here's what I tried, using the new

> LineCollection snippets you sent me, and the updated

> matplotlib snapshot:

Yep - very good idea. Your mistake with the line collection is how

you define the segments. From matplotlib.collections.LineCollection

documentation

segments is a sequence of ( line0, line1, line2), where linen =

(x0, y0), (x1, y1), ... (xm, ym). Each line can be a different

length

That is, it is a sequence of sequences of xy tuples. When you write

zip(xs.tolist(),ys.tolist()), you have a sequence of xy tuples. If

you plotted this, you would have one giant, connected line, which is

not what you want. You want a series of disconnected lines. Thus you

need to keep track of the indices where each separate segment starts,

and split out the segments, as in the code below.

For future reference, you may also want to use PolyCollections if you

want to generate a map and you have a bunch of segments defined by

sequences of xy tuples with associated face colors.

There was no significant difference between using bilinear

interpolation with antialiased drawing vs nearest neighbor

interpolation w/o aa, so I turned both back on.

And don't forget to try the new toolbar....

I noticed that the lat/lon lines don't precisely agree with the

colormap, eg around the Aleutian Islands the light blue is not

perfectly registered with the black lines. Should I be concerned that

this reflects a problem in matplotlib, or is this kind of error

standard in the data you've provided? I think this is related to the

pixel border that appears around some images, and is magnified when

interpolation is used because the top and right borders are not

defined when interpolating. I'll continue to look into this.

Would it be OK if I used this example on the web page? I like it!

Enjoy,

JDH

import cPickle

from matplotlib.matlab import *

from matplotlib.collections import LineCollection

from proj import Proj

import Numeric

# standard parallels at 50 deg N, center longitued 107 deg W.

params = {}

params['proj'] = 'lcc'

params['R'] = 63712000

params['lat_1'] = 50

params['lat_2'] = 50

params['lon_0'] = -107

proj = Proj(params)

llcornerx, llcornery = proj(-145.5,1.)

params['x_0'] = -llcornerx # add cartesian offset so lower left corner = (0,0)

params['y_0'] = -llcornery

# create a Proj instance for desired map.

proj = Proj(params)

# set the default params for imshow

rc('image', origin='lower', cmap='jet')

ax = subplot(111)

nx = 349; ny = 277

dx = 32463.41; dy = 32463.41

xmax = (nx-1)*dx; ymax = (ny-1)*dy # size of domain to plot

C = cPickle.load( file('topodata.pickle','rb') )

im = ax.imshow(C, interpolation='bilinear',

extent=(0, xmax, 0, ymax))

# ind is the index for the start of a new group

lons = []; lats = []; ind = []

i = 0 # the current ind

for line in file('wcl.txt'):

if line.startswith('# -b'):

ind.append(i)

continue

# lon/lat

lon, lat = [float(val) for val in line.split('\t')]

lons.append(lon)

lats.append(lat)

i += 1

xs, ys = proj(Numeric.array(lons),Numeric.array(lats))

#a sequence of xy tuples

segments = [zip(xs[i0:i1], ys[i0:i1]) for i0, i1 in zip(ind[:-1], ind[1:])]

collection = LineCollection(

segments,

colors = ( (0,0,0,1), ), # black

linewidths = (1.5,),

antialiaseds = (1,),) # turn off aa for speed

ax.add_collection(collection)

corners = (min(xs), min(ys)), (max(xs), max(ys))

ax.update_datalim( corners )

axis([0, xmax, 0, ymax])

ax.set_xticks([]) # no ticks

ax.set_yticks([])

title('Lambert Conformal Conic Projection')

#savefig('test')

show()