map projections

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()

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.

I'm pretty sure it's not the data. I only see this when I use imshow, not pcolor.

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

Sure. I've reworked the example a bit - now I read in a regular lat/lon grid from a pickle and use numarray's spline interpolation function to interpolate to the native projection grid. Here's the modified example:

from matplotlib.matlab import *
from matplotlib.collections import LineCollection from proj import Proj
import numarray, cPickle
from numarray import nd_image

# example to demonstrate plotting data on a map projection.
# requires numarray, proj module (which in turn requires # proj command from http://proj.maptools.org)

# 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.)
xmax=11297266.68; ymax=8959901.16
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)

# define grid (nx x ny regularly spaced native projection grid)
nx = 349; ny = 277 dx = xmax/(nx-1); dy = ymax/(ny-1)
xgrid = dx*numarray.indices((ny,nx))[1,:,:]
ygrid = dy*numarray.indices((ny,nx))[0,:,:]
# compute lons, lats of regular projection grid.
lonout, latout = proj(xgrid, ygrid, inverse=True)
# make sure lons are between 0 and 360
lonout = numarray.where(lonout < 0, lonout+360, lonout)
latout = latout+90

# read in topo data from pickle (on a regular lat/lon grid)
topodict = cPickle.load(open('etopo20.pickle','rb')) lons = topodict['lons']
lats = topodict['lats']
topoin = topodict['topo']

# find coordinates of native projection grid.
xcoords = (len(lons)-1)*(lonout-lons[0])/(lons[-1]-lons[0])
ycoords = (len(lats)-1)*(latout-lats[0])/(lats[-1]-lats[0])
coords = [ycoords,xcoords]
# interpolate to projection grid using numarray.nd_image spline filter.
topodat = numarray.nd_image.map_coordinates(topoin,coords,mode='nearest')

ax = subplot(111) # use imshow rather than pcolor for speed # set the default params for imshow rc('image', origin='lower', cmap='jet') im = ax.imshow(topodat, interpolation='nearest',extent=(0, xmax, 0, ymax)) #pcolor(xgrid, ygrid, topodat, shading='flat')

# read in coastline data from pickle.
wcl = cPickle.load(open('wcl.pickle','rb'))
ind = wcl['segment_index']; lons = wcl['lons']; lats = wcl['lats']
# transform coastline segements to map projection coordinates. xs, ys = proj(lons,lats)
# a sequence of xy tuples segments = [zip(xs[i0:i1], ys[i0:i1]) for i0, i1 in zip(ind[:-1], ind[1:])]

# line collection
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('20-minute Topography: Lambert Conformal Conic Projection') show()

I've tarred up all the pickles and scripts to run this at http://whitaker.homeunix.org/~jeff/plotmap.tar.gz. The resulting image is at http://whitaker.homeunix.org/~jeff/plotmap.png. When I get some time I'd like to wrap some of this in some easy-to-use functions, maybe with an interface similar to the matlab mapping toolbox.

Thanks for all your help!

-Jeff

···

On Wed, 14 Jul 2004, John Hunter wrote:

  -- 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

Jeff Whitaker writes:

> When I get some time I'd like to wrap some of this in some
> easy-to-use functions, maybe with an interface similar to the
> matlab mapping toolbox.

That would be very nice -- I think the matlab mapping/colorscaling
examples at

http://www.mathworks.com/company/newsletters/digest/nov02/earth_pt4.html
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectType=file&objectId=2611

represent the vast majority of things that people would want to do
with map/image overlays.

Regards, Phil