Re: [Matplotlib-users] map projections

[Jeff W. <js...@fa...>]

On Tue, 13 Jul 2004, John Hunter wrote:

>
> OK, made some headway here.
>
....
>
> For the line part of the map, I extended the
> matplotlib.collections.LineCollection class to handle a sequence of
> lines, where each line is defined by a list of tuples (x0,y0), (x1,
> y1), ...  Thus all of your lines are handled by a single object,
> rather than having 1800+ separate line objects created in plot.
> Again, no python loops required.
>
> In the current form, the code takes about 1.15s to run on my machine
> and is about 30x faster than the original code you posted which
> includes the data loading part.  Nonetheless, the matplotlib part is
> much faster too, as you'll see when you interact with the data.
>

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:


from matplotlib.matlab import *
from matplotlib.collections import LineCollection 
from proj import Proj
import Numeric, time

# open file with coastline data (from world coastline database).
wcl = open('wcl.txt')

# set up map projection parameters (lambert conformal conic,
# 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

t1 = time.clock()
segnum = 0
lons = []
lats = []
for line in wcl:
     splitline = line.split()
     if splitline[0] != '#':
         lon = float(splitline[0])
 	lat = float(splitline[1])
 	lons.append(lon)
 	lats.append(lat)
xs, ys = proj(Numeric.array(lons),Numeric.array(lats))
minx, maxx = min(xs), max(xs) 
miny, maxy = min(ys), max(ys) 
xypts = zip(xs.tolist(),ys.tolist())

wcl.close() # close coastline file
# all args are sequences, length 1 in case of linewidths and 
# antialiased 
collection = LineCollection(segments     = xypts,
                             colors       = ( (0,0,0,1), ), # black
                             linewidths   = (1.5,),
                             antialiaseds = (0,),)  # turn off aa for speed 
ax.add_collection(collection) 
ax.update_datalim( ((minx, miny), (maxx, maxy)))

axis([0, xmax, 0, ymax]) 
ax.set_xticks([]) # no ticks 
ax.set_yticks([]) 
title('Lambert Conformal Conic Projection') 
t2 = time.clock()
print t2-t1
show()

But when I run it, I get the following error:

Traceback (most recent call last):
   File 
"/sw/lib/python2.3/site-packages/matplotlib/backends/backend_gtkagg.py", 
line 75, in callback
     self.draw()
   File 
"/sw/lib/python2.3/site-packages/matplotlib/backends/backend_gtkagg.py", 
line 42, in draw
     agg.draw()
   File 
"/sw/lib/python2.3/site-packages/matplotlib/backends/backend_agg.py", line 
291, in draw
     self.figure.draw(self.renderer)
   File "/sw/lib/python2.3/site-packages/matplotlib/figure.py", line 236, 
in draw
     for a in self.axes: a.draw(renderer)
   File "/sw/lib/python2.3/site-packages/matplotlib/axes.py", line 668, in 
draw
     c.draw(renderer)
   File "/sw/lib/python2.3/site-packages/matplotlib/collections.py", line 
286, in draw
     self._transOffset)
TypeError: CXX: type error.


-Jeff

  -- 
Jeffrey S. Whitaker         Phone : (303)497-6313
NOAA/OAR/CDC  R/CDC1        FAX   : (303)497-6449
325 Broadway                Web   : http://www.cdc.noaa.gov/~jsw
Boulder, CO, USA 80305-3328 Office: Skaggs Research Cntr 1D-124