Re: [Matplotlib-users] map projections

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

>>>>> "Jeff" == Jeff Whitaker <js...@fa...> writes:

    Jeff> John: I found that if I just call proj with all the lats and
    Jeff> lons at once (instead of once for each segment) I can speed
    Jeff> it up tremendously. Here's what I tried, using the new
    Jeff> LineCollection snippets you sent me, and the updated
    Jeff> 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()