from matplotlib import rcParams
from matplotlib import __version__ as matplotlib_version
# check to make sure matplotlib is not too old.
mpl_required_version = '0.87.3'
if matplotlib_version < mpl_required_version:
    raise ImportError('your matplotlib is too old - basemap '
                      'requires at least %s'%mpl_required_version)
from matplotlib.collections import LineCollection
from matplotlib.patches import Polygon
from matplotlib.lines import Line2D
import pyproj, sys, os, math, dbflib
from proj import Proj
from greatcircle import GreatCircle, vinc_dist, vinc_pt
import matplotlib.numerix as NX
from matplotlib.numerix import ma, isnan
from matplotlib.mlab import linspace
from matplotlib.numerix.mlab import squeeze
from matplotlib.cbook import popd, is_scalar
from shapelib import ShapeFile

# look in sys.prefix for directory containing basemap files if
# BASEMAP_DATA_PATH env var not set.
_datadir = os.environ.get('BASEMAP_DATA_PATH')
if not _datadir:
   _datadir = os.path.join(sys.prefix,'share','py-basemap-data')

__version__ = '0.9.2'
__revision__ = '20060831'

# test to see if array indexing is supported
# (it is not for Numeric, but is for numarray and numpy)
try:
    NX.ones(10)[[1,2]]
    _has_arrindexing = True
except:
    _has_arrindexing = False

class Basemap(object):

    """
 Set up a basemap with one of 18 supported map projections
 (cylindrical equidistant, mercator, polyconic, oblique mercator,
 transverse mercator, miller cylindrical, lambert conformal conic,
 azimuthal equidistant, equidistant conic, lambert azimuthal equal area,
 albers equal area conic, gnomonic, orthographic, sinusoidal, mollweide, 
 robinson, cassini-soldner or stereographic).
 Doesn't actually draw anything, but sets up the map projection class and
 creates the coastline, lake river and political boundary data
 structures in native map projection coordinates.
 Uses a pyrex interface to C-code from proj.4 (http://proj.maptools.org).

 Useful instance variables:

 projection - map projection ('cyl','merc','mill','lcc','eqdc','aea',
  'laea', 'nplaea', 'splaea', 'tmerc', 'omerc', 'cass', 'gnom', 'poly', 
  'sinu', 'moll', 'ortho', 'robin', 'aeqd', 'npaeqd', 'spaeqd', 'stere',
  'npstere' or 'spstere')
 (projections prefixed with 'np' or 'sp' are special case polar-centric
  versions of the parent projection)
 aspect - map aspect ratio (size of y dimension / size of x dimension).
 llcrnrlon - longitude of lower left hand corner of the desired map domain.
 llcrnrlon - latitude of lower left hand corner of the desired map domain.
 urcrnrlon - longitude of upper right hand corner of the desired map domain.
 urcrnrlon - latitude of upper right hand corner of the desired map domain.
 llcrnrx,llcrnry,urcrnrx,urcrnry - corners of map domain in projection coordinates.
 rmajor,rminor - equatorial and polar radii of ellipsoid used (in meters).
 resolution - resolution of boundary dataset being used ('c' for crude,
   'l' for low, etc.). If None, no boundary dataset is associated with the
   Basemap instance.  

 Example Usage:

>>> from matplotlib.toolkits.basemap import Basemap
>>> from pylab import *
>>> # read in topo data (on a regular lat/lon grid)
>>> etopo = load('etopo20data.gz')
>>> lons  = load('etopo20lons.gz')
>>> lats  = load('etopo20lats.gz')
>>> # create Basemap instance for Robinson projection.
>>> m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
>>> # compute native map projection coordinates for lat/lon grid.
>>> x, y = m(*meshgrid(lons,lats))
>>> # make filled contour plot.
>>> cs = m.contourf(x,y,etopo,30,cmap=cm.jet)
>>> m.drawcoastlines() # draw coastlines
>>> m.drawmapboundary() # draw a line around the map region
>>> m.drawparallels(arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels
>>> m.drawmeridians(arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
>>> title('Robinson Projection') # add a title
>>> show()

 [this example (simpletest.py) plus many others can be found in the
  examples directory of source distribution.  The "OO" version of this
  example (which does not use pylab) is called "simpletest_oo.py".]

 Contact: Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
    """

    def __init__(self,llcrnrlon=None,llcrnrlat=None,
       urcrnrlon=None,urcrnrlat=None,\
       width=None,height=None,\
       projection='cyl',resolution='c',area_thresh=None,rsphere=6370997.0,\
       lat_ts=None,lat_1=None,lat_2=None,lat_0=None,lon_0=None,\
       lon_1=None,lon_2=None,suppress_ticks=True,\
       boundinglat=None,anchor='C',ax=None):
        """
 create a Basemap instance.

 arguments:

 projection - map projection.  'cyl' - cylindrical equidistant, 'merc' -
  mercator, 'lcc' - lambert conformal conic, 'stere' - stereographic,
  'npstere' - stereographic, special case centered on north pole.
  'spstere' - stereographic, special case centered on south pole,
  'aea' - albers equal area conic, 'tmerc' - transverse mercator,
  'aeqd' - azimuthal equidistant, 'mill' - miller cylindrical,
  'npaeqd' - azimuthal equidistant, special case centered on north pole,
  'spaeqd' - azimuthal equidistant, special case centered on south pole,
  'eqdc' - equidistant conic, 'laea' - lambert azimuthal equal area,
  'nplaea' - lambert azimuthal, special case centered on north pole,
  'splaea' - lambert azimuthal, special case centered on south pole,
  'cass' - cassini-soldner (transverse cylindrical equidistant),
  'poly' - polyconic, 'omerc' - oblique mercator, 'ortho' - orthographic,
  'sinu' - sinusoidal, 'moll' - mollweide, 'robin' - robinson,
  and 'gnom' - gnomonic are currently available.  Default 'cyl'.

 The map projection region can either be specified by setting these keywords:

 llcrnrlon - longitude of lower left hand corner of the desired map domain (degrees).
 llcrnrlat - latitude of lower left hand corner of the desired map domain (degrees).
 urcrnrlon - longitude of upper right hand corner of the desired map domain (degrees).
 urcrnrlat - latitude of upper right hand corner of the desired map domain (degrees).

 or these keywords:

 width  - width of desired map domain in projection coordinates (meters).
 height - height of desired map domain in projection coordinates (meters).
 lon_0  - center of desired map domain (in degrees).
 lat_0  - center of desired map domain (in degrees).

 For 'ortho', 'sinu', 'moll', 'npstere', 'spstere', 'nplaea', 'splaea', 'nplaea', 
 'splaea', 'npaeqd', 'spaeqd' or 'robin', the values of 
 llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,width and height are ignored (because 
 either they are computed internally, or entire globe is always plotted). For the 
 cylindrical projections ('cyl','merc' and 'mill'), the default is to use 
 llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other 
 projections, either the lat/lon values of the corners or width and height must be 
 specified by the user.  

 resolution - resolution of boundary database to use. Can be 'c' (crude),
  'l' (low), 'i' (intermediate), 'h' (high), or None. Default is 'c'.
  If None, no boundary data will be read in (and class methods
  such as drawcoastlines will raise an exception if invoked).
  Resolution drops off by roughly 80%
  between datasets.  Higher res datasets are much slower to draw.
  Default 'c'. Coastline data is from the GSHHS
  (http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html).
  State, country and river datasets from the Generic Mapping
  Tools (http://gmt.soest.hawaii.edu).

 area_thresh - coastline or lake with an area smaller than area_thresh
  in km^2 will not be plotted.  Default 10000,1000,100,10 for resolution
  'c','l','i','h'.

 rsphere - radius of the sphere used to define map projection (default
  6370997 meters, close to the arithmetic mean radius of the earth). If
  given as a sequence, the first two elements are interpreted as
  the the radii of the major and minor axes of an ellipsoid. Note: sometimes
  an ellipsoid is specified by the major axis and an 'inverse flattening
  parameter' (if).  The minor axis (b) can be computed from the major axis (a)
  and the inverse flattening parameter using the formula if = a/(a-b).

 suppress_ticks - suppress automatic drawing of axis ticks and labels
 in map projection coordinates.  Default False, so parallels and meridians
 can be labelled instead. If parallel or meridian labelling is requested
 (using drawparallels and drawmeridians methods), automatic tick labelling
 will be supressed even is suppress_ticks=False.  Typically, you will
 only want to override the default if you want to label the axes in meters
 using native map projection coordinates.

 anchor - determines how map is placed in axes rectangle (passed to
 axes.set_aspect). Default is 'C', which means map is centered.
 Allowed values are ['C', 'SW', 'S', 'SE', 'E', 'NE', 'N', 'NW', 'W'].

 ax - set default axes instance (default None - pylab.gca() may be used
 to get the current axes instance).  If you don't want pylab to be imported,
 you can either set this to a pre-defined axes instance, or use the 'ax'
 keyword in each Basemap method call that does drawing. In the first case,
 all Basemap method calls will draw to the same axes instance.  In the
 second case, you can draw to different axes with the same Basemap instance.
 You can also use the 'ax' keyword in individual method calls to
 selectively override the default axes instance.

 The following parameters are map projection parameters which all default to
 None.  Not all parameters are used by all projections, some are ignored.

 lat_ts - latitude of natural origin (used for mercator, and 
  optionally for stereographic projection).
 lat_1 - first standard parallel for lambert conformal, albers
  equal area projection and equidistant conic projections. Latitude of one
  of the two points on the projection centerline for oblique mercator.
  If lat_1 is not given, but lat_0 is, lat_1 is set to lat_0 for
  lambert conformal, albers equal area and equidistant conic.
 lat_2 - second standard parallel for lambert conformal, albers
  equal area projection and equidistant conic projections. Latitude of one
  of the two points on the projection centerline for oblique mercator.
  If lat_2 is not given, it is set to lat_1 for 
  lambert conformal, albers equal area and equidistant conic.
 lon_1 - longitude of one of the two points on the projection centerline
  for oblique mercator.
 lon_2 - longitude of one of the two points on the projection centerline
  for oblique mercator.
 lat_0 - central latitude (y-axis origin) - used by all projections, 
 lon_0 - central meridian (x-axis origin) - used by all projections,
 boundinglat - bounding latitude for pole-centered projections (npstere,spstere,
  nplaea,splaea,npaeqd,spaeqd).  These projections are square regions centered
  on the north or south pole.  The longitude lon_0 is at 6-o'clock, and the
  latitude circle boundinglat is tangent to the edge of the map at lon_0.

        """

        # where to put plot in figure (default is 'C' or center)
        self.anchor = anchor
        # map projection.
        self.projection = projection

        # set up projection parameter dict.
        projparams = {}
        projparams['proj'] = projection
        try:
            if rsphere[0] > rsphere[1]:
                projparams['a'] = rsphere[0]
                projparams['b'] = rsphere[1]
            else:
                projparams['a'] = rsphere[1]
                projparams['b'] = rsphere[0]
        except:
            projparams['R'] = rsphere
            # this is a workaround for bug in proj 4.4.9 tmerc where a=b.
            #projparams['a'] = rsphere
            #projparams['b'] = rsphere - 0.01
        # set units to meters.
        projparams['units']='m'
        # check for sane values of lon_0, lat_0, lat_ts, lat_1, lat_2
        if lat_0 is not None:
            if lat_0 > 90. or lat_0 < -90.:
                raise ValueError, 'lat_0 must be between -90 and +90 degrees'
            else:
                projparams['lat_0'] = lat_0
        if lon_0 is not None:
            if lon_0 < -720. or lon_0 > 720.:
                raise ValueError, 'lon_0 must be between -720 and +720 degrees'
            else:
                projparams['lon_0'] = lon_0
        if lon_1 is not None:
            if lon_1 < -720. or lon_1 > 720.:
                raise ValueError, 'lon_1 must be between -720 and +720 degrees'
            else:
                projparams['lon_1'] = lon_1
        if lon_2 is not None:
            if lon_2 < -720. or lon_2 > 720.:
                raise ValueError, 'lon_2 must be between -720 and +720 degrees'
            else:
                projparams['lon_2'] = lon_2
        if lat_1 is not None:
            if lat_1 > 90. or lat_1 < -90.:
                raise ValueError, 'lat_1 must be between -90 and +90 degrees'
            else:
                projparams['lat_1'] = lat_1
        if lat_2 is not None:
            if lat_2 > 90. or lat_2 < -90.:
                raise ValueError, 'lat_2 must be between -90 and +90 degrees'
            else:
                projparams['lat_2'] = lat_2
        if lat_ts is not None:
            if lat_ts > 90. or lat_ts < -90.:
                raise ValueError, 'lat_ts must be between -90 and +90 degrees'
            else:
                projparams['lat_ts'] = lat_ts

        if None not in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
            # make sure lat/lon limits are converted to floats.
            self.llcrnrlon = float(llcrnrlon)
            self.llcrnrlat = float(llcrnrlat)
            self.urcrnrlon = float(urcrnrlon)
            self.urcrnrlat = float(urcrnrlat)
            # check values of urcrnrlon,urcrnrlat and llcrnrlon,llcrnrlat
            if self.urcrnrlat > 90.0 or self.llcrnrlat > 90.0:
                raise ValueError, 'urcrnrlat and llcrnrlat must be less than 90'
            if self.urcrnrlat < -90.0 or self.llcrnrlat < -90.0:
                raise ValueError, 'urcrnrlat and llcrnrlat must be greater than -90'
            if self.urcrnrlon > 720. or self.llcrnrlon > 720.:
                raise ValueError, 'urcrnrlon and llcrnrlon must be less than 720'
            if self.urcrnrlon < -360. or self.llcrnrlon < -360.:
                raise ValueError, 'urcrnrlon and llcrnrlon must be greater than -360'
        # for each of the supported projections, compute lat/lon of domain corners 
        # and set values in projparams dict as needed.
        if projection == 'lcc':
            # if lat_0 is given, but not lat_1,
            # set lat_1=lat_0
            if lat_1 is None and lat_0 is not None:
                lat_1 = lat_0
                projparams['lat_1'] = lat_1
            if lat_1 is None or lon_0 is None:
                raise ValueError, 'must specify lat_1 or lat_0 and lon_0 for Lambert Conformal basemap (lat_2 is optional)'
            if lat_2 is None:
               projparams['lat_2'] = lat_1
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
                    
        elif projection == 'eqdc':
            # if lat_0 is given, but not lat_1,
            # set lat_1=lat_0
            if lat_1 is None and lat_0 is not None:
                lat_1 = lat_0
                projparams['lat_1'] = lat_1
            if lat_1 is None or lon_0 is None:
                raise ValueError, 'must specify lat_1 or lat_0 and lon_0 for Equidistant Conic basemap (lat_2 is optional)'
            if lat_2 is None:
               projparams['lat_2'] = lat_1
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'aea':
            # if lat_0 is given, but not lat_1,
            # set lat_1=lat_0
            if lat_1 is None and lat_0 is not None:
                lat_1 = lat_0
                projparams['lat_1'] = lat_1
            if lat_1 is None or lon_0 is None:
                raise ValueError, 'must specify lat_1 or lat_0 and lon_0 for Albers Equal Area basemap (lat_2 is optional)'
            if lat_2 is None:
               projparams['lat_2'] = lat_1
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'stere':
            if lat_0 is None or lon_0 is None:
                raise ValueError, 'must specify lat_0 and lon_0 for Stereographic basemap (lat_ts is optional)'
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'spstere':
            if boundinglat is None or lon_0 is None:
                raise ValueError, 'must specify boundinglat and lon_0 for South-Polar Stereographic basemap'
            projparams['lat_ts'] = -90.
            projparams['lat_0'] = -90.
            projparams['proj'] = 'stere'
            self.llcrnrlon = lon_0+45.
            self.urcrnrlon = lon_0-135.
            proj = pyproj.Proj(projparams)
            x,y = proj(lon_0,boundinglat)
            lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
            self.urcrnrlat = self.llcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'npstere':
            if boundinglat is None or lon_0 is None:
                raise ValueError, 'must specify boundinglat and lon_0 for North-Polar Stereographic basemap'
            projparams['lat_ts'] = 90.
            projparams['lat_0'] = 90.
            projparams['proj'] = 'stere'
            self.llcrnrlon = lon_0-45.
            self.urcrnrlon = lon_0+135.
            proj = pyproj.Proj(projparams)
            x,y = proj(lon_0,boundinglat)
            lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
            self.urcrnrlat = self.llcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'splaea':
            if boundinglat is None or lon_0 is None:
                raise ValueError, 'must specify boundinglat and lon_0 for South-Polar Lambert Azimuthal basemap'
            projparams['lat_0'] = -90.
            projparams['proj'] = 'laea'
            self.llcrnrlon = lon_0+45.
            self.urcrnrlon = lon_0-135.
            proj = pyproj.Proj(projparams)
            x,y = proj(lon_0,boundinglat)
            lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
            self.urcrnrlat = self.llcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'nplaea':
            if boundinglat is None or lon_0 is None:
                raise ValueError, 'must specify boundinglat and lon_0 for North-Polar Lambert Azimuthal basemap'
            projparams['lat_0'] = 90.
            projparams['proj'] = 'laea'
            self.llcrnrlon = lon_0-45.
            self.urcrnrlon = lon_0+135.
            proj = pyproj.Proj(projparams)
            x,y = proj(lon_0,boundinglat)
            lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
            self.urcrnrlat = self.llcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'spaeqd':
            if boundinglat is None or lon_0 is None:
                raise ValueError, 'must specify boundinglat and lon_0 for South-Polar Azimuthal Equidistant basemap'
            projparams['lat_0'] = -90.
            projparams['proj'] = 'aeqd'
            self.llcrnrlon = lon_0+45.
            self.urcrnrlon = lon_0-135.
            proj = pyproj.Proj(projparams)
            x,y = proj(lon_0,boundinglat)
            lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
            self.urcrnrlat = self.llcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'npaeqd':
            if boundinglat is None or lon_0 is None:
                raise ValueError, 'must specify boundinglat and lon_0 for North-Polar Azimuthal Equidistant basemap'
            projparams['lat_0'] = 90.
            projparams['proj'] = 'aeqd'
            self.llcrnrlon = lon_0-45.
            self.urcrnrlon = lon_0+135.
            proj = pyproj.Proj(projparams)
            x,y = proj(lon_0,boundinglat)
            lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
            self.urcrnrlat = self.llcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'laea':
            if lat_0 is None or lon_0 is None:
                raise ValueError, 'must specify lat_0 and lon_0 for Lambert Azimuthal basemap'
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'merc':
            if lat_ts is None:
                raise ValueError, 'must specify lat_ts for Mercator basemap'
            # clip plot region to be within -89.99S to 89.99N
            # (mercator is singular at poles)
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
            if llcrnrlat < -89.99: llcrnrlat = -89.99
            if llcrnrlat > 89.99: llcrnrlat = 89.99
            if urcrnrlat < -89.99: urcrnrlat = -89.99
            if urcrnrlat > 89.99: urcrnrlat = 89.99
            self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
            self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection in ['tmerc','gnom','cass','poly'] :
            if lat_0 is None or lon_0 is None:
                raise ValueError, 'must specify lat_0 and lon_0 for Transverse Mercator, Gnomonic, Cassini-Soldner, Orthographic or Polyconic basemap'
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
                    
        elif projection == 'ortho':
            if lat_0 is None or lon_0 is None:
                raise ValueError, 'must specify lat_0 and lon_0 for Transverse Mercator, Gnomonic, Cassini-Soldner, Orthographic or Polyconic basemap'
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection in ['moll','robin','sinu']:
            if lon_0 is None:
                raise ValueError, 'must specify lon_0 for Robinson, Mollweide, or Sinusoidal basemap'
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'omerc':
            if lat_1 is None or lon_1 is None or lat_2 is None or lon_2 is None:
                raise ValueError, 'must specify lat_1,lon_1 and lat_2,lon_2 for Oblique Mercator basemap'
            projparams['lat_1'] = lat_1
            projparams['lon_1'] = lon_1
            projparams['lat_2'] = lat_2
            projparams['lon_2'] = lon_2
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'aeqd':
            if lat_0 is None or lon_0 is None:
                raise ValueError, 'must specify lat_0 and lon_0 for Azimuthal Equidistant basemap'
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                if width is None or height is None:
                    raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
                else:
                    if lon_0 is None or lat_0 is None:
                        raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'mill':
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'cyl':
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
        elif projection == 'sinu':
            if lon_0 is None:
                raise ValueError, 'must specify lon_0 for Robinson, Mollweide, or Sinusoidal basemap'
            if width is not None or height is not None:
                print 'warning: width and height keywords ignored for %s projection' % self.projection
            if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        else:
            raise ValueError, """
  unsupported projection, use 'cyl' - cylindrical equidistant, 'merc' -
  mercator, 'lcc' - lambert conformal conic, 'stere' - stereographic,
  'npstere' - stereographic, special case centered on north pole.
  'spstere' - stereographic, special case centered on south pole,
  'aea' - albers equal area conic, 'tmerc' - transverse mercator,
  'aeqd' - azimuthal equidistant, 'mill' - miller cylindrical,
  'npaeqd' - azimuthal equidistant, special case centered on north pole,
  'spaeqd' - azimuthal equidistant, special case centered on south pole,
  'eqdc' - equidistant conic, 'laea' - lambert azimuthal equal area,
  'nplaea' - lambert azimuthal, special case centered on north pole,
  'splaea' - lambert azimuthal, special case centered on south pole,
  'cass' - cassini-soldner (transverse cylindrical equidistant),
  'poly' - polyconic, 'omerc' - oblique mercator, 'ortho' - orthographic,
  'sinu' - sinusoidal, 'moll' - mollweide, 'robin' - robinson,
  or 'gnom' - gnomonic.  You tried '%s'""" % projection

        # initialize proj4
        proj = Proj(projparams,self.llcrnrlon,self.llcrnrlat,self.urcrnrlon,self.urcrnrlat)
        
        # make sure axis ticks are suppressed.
        self.noticks = suppress_ticks

        # make Proj instance a Basemap instance variable.
        self.projtran = proj
        # copy some Proj attributes.
        atts = ['rmajor','rminor','esq','flattening','ellipsoid','projparams']
        for att in atts:
            self.__dict__[att] = proj.__dict__[att]
        # set instance variables defining map region.
        self.xmin = proj.xmin
        self.xmax = proj.xmax
        self.ymin = proj.ymin
        self.ymax = proj.ymax
        if projection == 'cyl':
            self.aspect = (self.urcrnrlat-self.llcrnrlat)/(self.urcrnrlon-self.llcrnrlon)
        else:
            self.aspect = (proj.ymax-proj.ymin)/(proj.xmax-proj.xmin)
        self.llcrnrx = proj.llcrnrx
        self.llcrnry = proj.llcrnry
        self.urcrnrx = proj.urcrnrx
        self.urcrnry = proj.urcrnry

        # set min/max lats for projection domain.
        if projection in ['mill','cyl','merc']:
            self.latmin = self.llcrnrlat
            self.latmax = self.urcrnrlat
        elif projection in ['ortho','moll','robin','sinu']:
            self.latmin = -90.
            self.latmax = 90.
        else:
            lons, lats = self.makegrid(101,101)
            self.latmin = min(NX.ravel(lats))
            self.latmax = max(NX.ravel(lats))

        # if ax == None, pylab.gca may be used.
        self.ax = ax
        self.lsmask = None

        # set defaults for area_thresh.
        self.resolution = resolution
        # if no boundary data needed, we are done.
        if self.resolution is None:
            return
        if area_thresh is None:
            if resolution == 'c':
                area_thresh = 10000.
            elif resolution == 'l':
                area_thresh = 1000.
            elif resolution == 'i':
                area_thresh = 100.
            elif resolution == 'h':
                area_thresh = 10.
            else:
                raise ValueError, "boundary resolution must be one of 'c','l','i' or 'h'"
        # read in coastline data (only those polygons whose area > area_thresh).
        coastlons = []; coastlats = []; coastsegind = []; coastsegtype = []
        msg = """
Unable to open boundary dataset file. Only the 'crude' and 'low'
resolution datasets are installed by default. If you
are requesting 'intermediate' or 'high' resolution datasets, you
need to install the basemap_data package.  If the boundary
datasets are not installed in sys.prefix/share, the BASEMAP_DATA_PATH
environment variable must be set."""
        try:
            bdatfile = open(os.path.join(_datadir,'gshhs_'+resolution+'.txt'))
        except:
            raise IOError, msg
        for line in bdatfile:
            linesplit = line.split()
            if line.startswith('P'):
                area = float(linesplit[5])
                west,east,south,north = float(linesplit[6]),float(linesplit[7]),float(linesplit[8]),float(linesplit[9])
                typ = int(linesplit[3])
                useit = self.latmax>=south and self.latmin<=north and area>area_thresh
                if useit:
                    coastsegind.append(len(coastlons))
                    coastsegtype.append(typ)
                continue
            # lon/lat
            if useit:
                lon, lat = [float(val) for val in linesplit]
                coastlons.append(lon)
                coastlats.append(lat)
        coastsegtype.append(typ)
        coastsegind.append(len(coastlons))

        # read in country boundary data.
        cntrylons = []; cntrylats = []; cntrysegind = []
        try:
            bdatfile = open(os.path.join(_datadir,'countries_'+resolution+'.txt'))
        except:
            raise IOError, msg
        for line in bdatfile:
            linesplit = line.split()
            if line.startswith('>'):
                west,east,south,north = float(linesplit[7]),float(linesplit[8]),float(linesplit[9]),float(linesplit[10])
                useit = self.latmax>=south and self.latmin<=north
                if useit: cntrysegind.append(len(cntrylons))
                continue
            # lon/lat
            if useit:
                lon, lat = [float(val) for val in linesplit]
                cntrylons.append(lon)
                cntrylats.append(lat)
        cntrysegind.append(len(cntrylons))

        # read in state boundaries (Americas only).
        statelons = []; statelats = []; statesegind = []
        try:
            bdatfile = open(os.path.join(_datadir,'states_'+resolution+'.txt'))
        except:
            raise IOError, msg
        for line in bdatfile:
            linesplit = line.split()
            if line.startswith('>'):
                west,east,south,north = float(linesplit[7]),float(linesplit[8]),float(linesplit[9]),float(linesplit[10])
                useit = self.latmax>=south and self.latmin<=north
                if useit: statesegind.append(len(statelons))
                continue
            # lon/lat
            if useit:
                lon, lat = [float(val) for val in linesplit]
                statelons.append(lon)
                statelats.append(lat)
        statesegind.append(len(statelons))

        # read in major rivers.
        riverlons = []; riverlats = []; riversegind = []
        try:
            bdatfile = open(os.path.join(_datadir,'rivers_'+resolution+'.txt'))
        except:
            raise IOError, msg
        for line in bdatfile:
            linesplit = line.split()
            if line.startswith('>'):
                west,east,south,north = float(linesplit[7]),float(linesplit[8]),float(linesplit[9]),float(linesplit[10])
                useit = self.latmax>=south and self.latmin<=north
                if useit: riversegind.append(len(riverlons))
                continue
            # lon/lat
            if useit:
                lon, lat = [float(val) for val in linesplit]
                riverlons.append(lon)
                riverlats.append(lat)
        riversegind.append(len(riverlons))

        # extend longitudes around the earth a second time
        # so valid longitudes can range from -360 to 720.
        # This means a lot of redundant processing is done when
        # creating the class instance, but it a lot easier to figure
        # out what to do when the projection domain straddles the 
        # Greenwich meridian.
        coastlons2 = [lon+360. for lon in coastlons]
        cntrylons2 = [lon+360. for lon in cntrylons]
        statelons2 = [lon+360. for lon in statelons]
        riverlons2 = [lon+360. for lon in riverlons]
        coastlons3 = [lon-360. for lon in coastlons]
        cntrylons3 = [lon-360. for lon in cntrylons]
        statelons3 = [lon-360. for lon in statelons]
        riverlons3 = [lon-360. for lon in riverlons]

        # transform coastline polygons to native map coordinates.
        xc,yc = proj(NX.array(coastlons),NX.array(coastlats))
        xc = xc.tolist(); yc = yc.tolist()
        xc2,yc2 = proj(NX.array(coastlons2),NX.array(coastlats))
        xc3,yc3 = proj(NX.array(coastlons3),NX.array(coastlats))
        xc2 = xc2.tolist(); yc2 = yc2.tolist()
        xc3 = xc3.tolist(); yc3 = yc3.tolist()

        # set up segments in form needed for LineCollection,
        # ignoring 'inf' values that are off the map.
        segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:])]
        segmentsll = [zip(coastlons[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:])]
        segtypes = [i for i in coastsegtype[:-1]]
        segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
        segmentsll2 = [zip(coastlons2[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
        segtypes2 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
        segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
        segmentsll3 = [zip(coastlons3[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
        segtypes3 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
        self.coastsegs = segments+segments2+segments3
        self.coastsegsll = segmentsll+segmentsll2+segmentsll3
        self.coastsegtypes = segtypes+segtypes2+segtypes3

        # same as above for country segments.
        xc,yc = proj(NX.array(cntrylons),NX.array(cntrylats))
        xc = xc.tolist(); yc = yc.tolist()
        xc2,yc2 = proj(NX.array(cntrylons2),NX.array(cntrylats))
        xc3,yc3 = proj(NX.array(cntrylons3),NX.array(cntrylats))
        xc2 = xc2.tolist(); yc2 = yc2.tolist()
        xc3 = xc3.tolist(); yc3 = yc3.tolist()
        segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:])]
        segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
        segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
        self.cntrysegs = segments+segments2+segments3

        # same as above for state segments.
        xc,yc = proj(NX.array(statelons),NX.array(statelats))
        xc = xc.tolist(); yc = yc.tolist()
        xc2,yc2 = proj(NX.array(statelons2),NX.array(statelats))
        xc3,yc3 = proj(NX.array(statelons3),NX.array(statelats))
        xc2 = xc2.tolist(); yc2 = yc2.tolist()
        xc3 = xc3.tolist(); yc3 = yc3.tolist()
        segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:])]
        segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
        segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
        self.statesegs = segments+segments2+segments3

        # same as above for river segments.
        xc,yc = proj(NX.array(riverlons),NX.array(riverlats))
        xc = xc.tolist(); yc = yc.tolist()
        xc2,yc2 = proj(NX.array(riverlons2),NX.array(riverlats))
        xc3,yc3 = proj(NX.array(riverlons3),NX.array(riverlats))
        xc2 = xc2.tolist(); yc2 = yc2.tolist()
        xc3 = xc3.tolist(); yc3 = yc3.tolist()
        segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:])]
        segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
        segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
        self.riversegs = segments+segments2+segments3

        # store coast polygons for filling.
        self.coastpolygons = []
        coastpolygonsll = []
        self.coastpolygontypes = []
        if projection in ['merc','mill']:
            xsp,ysp = proj(0.,-89.9) # s. pole coordinates.
            xa,ya = proj(0.,-68.) # edge of antarctica.
            x0,y0 = proj(0.,0.)
            xm360,ym360 = proj(-360.,0.)
            x360,y360 = proj(360.,0.)
            x720,y720 = proj(720.,0.)
        for seg,segtype,segll in zip(self.coastsegs,self.coastsegtypes,self.coastsegsll):
            x = [lon for lon,lat in seg]
            y = [lat for lon,lat in seg]
            lons = [lon for lon,lat in segll]
            lats = [lat for lon,lat in segll]
            # the antarctic polygon is a nuisance, since it
            # spans all longitudes, it's not closed and will not be filled
            # without some projection dependant tweaking.
            if projection == 'cyl':
                if x[-1] == 0.000 and y[-1] < -68.: # close antarctica
                    x.append(0.)
                    y.append(-90.0000)
                    x.insert(0,360.)
                    y.insert(0,-90)
                    lons.append(0.)
                    lats.append(-90.)
                    lons.insert(0,360.)
                    lats.insert(0,-90.)
                if x[-1] == 360.000 and y[-1] < -68.:
                    x.append(360.)
                    y.append(-90)
                    x.insert(0,720.)
                    y.insert(0,-90)
                    lons.append(360.)
                    lats.append(-90.)
                    lons.insert(0,720.)
                    lats.insert(0,-90.)
                if x[-1] == -360.000 and y[-1] < -68.:
                    x.append(-360.)
                    y.append(-90)
                    x.insert(0,0.)
                    y.insert(0,-90)
                    lons.append(-360.)
                    lats.append(-90.)
                    lons.insert(0,0.)
                    lats.insert(0,-90.)
            elif projection in ['merc','mill']:
                if math.fabs(x[-1]-x0) < 1. and y[-1] < ya: # close antarctica
                    x.append(x0)
                    y.append(ysp)
                    x.insert(0,x360)
                    y.insert(0,ysp)
                    lons.append(0.)
                    lats.append(-90.)
                    lons.insert(0,360.)
                    lats.insert(0,-90.)
                if math.fabs(x[-1]-x360) < 1. and y[-1] < ya:
                    x.append(x360)
                    y.append(ysp)
                    x.insert(0,x720)
                    y.insert(0,ysp)
                    lons.append(360.)
                    lats.append(-90.)
                    lons.insert(0,720.)
                    lats.insert(0,-90.)
                if math.fabs(x[-1]-xm360) < 1. and y[-1] < ya:
                    x.append(xm360)
                    y.append(ysp)
                    x.insert(0,x0)
                    y.insert(0,ysp)
                    lons.append(-360.)
                    lats.append(-90.)
                    lons.insert(0,0.)
                    lats.insert(0,-90.)
            self.coastpolygons.append((x,y))
            coastpolygonsll.append((lons,lats))
            self.coastpolygontypes.append(segtype)

        # remove those segments/polygons that don't intersect map region.
        coastsegs = []
        coastsegtypes = []
        for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
            if self._insidemap_seg(seg):
                coastsegs.append(seg)
                coastsegtypes.append(segtype)
        self.coastsegs = coastsegs
        self.coastsegtypes = coastsegtypes
        polygons = []
        polygonsll = []
        polygontypes = []
        for poly,polytype,polyll in zip(self.coastpolygons,self.coastpolygontypes,coastpolygonsll):
            if self._insidemap_poly(poly,polyll):
                polygons.append(poly)
                polygontypes.append(polytype)
                polygonsll.append(polyll)
        self.coastpolygons = polygons
        coastpolygonsll = polygonsll
        self.coastpolygontypes = polygontypes
        states = []; rivers = []; countries = []
        for seg in self.cntrysegs:
            if self._insidemap_seg(seg):
                countries.append(seg)
        for seg in self.statesegs:
            if self._insidemap_seg(seg):
                states.append(seg)
        for seg in self.riversegs:
            if self._insidemap_seg(seg):
                rivers.append(seg)
        self.statesegs = states
        self.riversegs = rivers
        self.cntryegs = countries

        # split up segments that go outside projection limb
        coastsegs = []
        coastsegtypes = []
        for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
            xx = NX.array([x for x,y in seg],NX.Float32)
            yy = NX.array([y for x,y in seg],NX.Float32)
            i1,i2 = self._splitseg(xx,yy)
            if i1 and i2:
                for i,j in zip(i1,i2):
                    segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
                    coastsegs.append(segment)
                    coastsegtypes.append(segtype)
            else:
                coastsegs.append(seg)
                coastsegs.append(segtype)
        self.coastsegs = coastsegs
        self.coastsegtypes = coastsegtypes
        states = []
        for seg in self.statesegs:
            xx = NX.array([x for x,y in seg],NX.Float32)
            yy = NX.array([y for x,y in seg],NX.Float32)
            i1,i2 = self._splitseg(xx,yy)
            if i1 and i2:
                for i,j in zip(i1,i2):
                    segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
                    states.append(segment)
            else:
                states.append(seg)
        self.statesegs = states
        countries = []
        for seg in self.cntrysegs:
            xx = NX.array([x for x,y in seg],NX.Float32)
            yy = NX.array([y for x,y in seg],NX.Float32)
            i1,i2 = self._splitseg(xx,yy)
            if i1 and i2:
                for i,j in zip(i1,i2):
                    segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
                    countries.append(segment)
            else:
                countries.append(seg)
        self.cntrysegs = countries
        rivers = []
        for seg in self.riversegs:
            xx = NX.array([x for x,y in seg],NX.Float32)
            yy = NX.array([y for x,y in seg],NX.Float32)
            i1,i2 = self._splitseg(xx,yy)
            if i1 and i2:
                for i,j in zip(i1,i2):
                    segment = zip(xx[i:j].tolist(),yy[i:j].tolist())
                    rivers.append(segment)
            else:
                rivers.append(seg)
        self.riversegs = rivers

        # split coastline segments that jump across entire plot.
        coastsegs = []
        coastsegtypes = []
        for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
            xx = NX.array([x for x,y in seg],NX.Float32)
            yy = NX.array([y for x,y in seg],NX.Float32)
            xd = (xx[1:]-xx[0:-1])**2
            yd = (yy[1:]-yy[0:-1])**2
            dist = NX.sqrt(xd+yd)
            split = dist > 5000000.
            if NX.sum(split) and self.projection not in ['merc','cyl','mill']:
               ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
               iprev = 0
               ind.append(len(xd))
               for i in ind:
                   coastsegs.append(zip(xx[iprev:i],yy[iprev:i]))
                   coastsegtypes.append(segtype)
                   iprev = i
            else:
                coastsegs.append(seg)
                coastsegtypes.append(segtype)
        self.coastsegs = coastsegs
        self.coastsegtypes = coastsegtypes

        # special treatment of coastline polygons for
        # orthographic, sinusoidal, mollweide and robinson.
        # (polygon clipping along projection limb)
        if self.projection == 'ortho':
            lat_0 = math.radians(self.projparams['lat_0'])
            lon_0 = math.radians(self.projparams['lon_0'])
            rad = (2.*self.rmajor + self.rminor)/3.
            dtheta = 0.01
            coastpolygons = []
            coastpolygontypes = []
            for poly,polytype,polyll in zip(self.coastpolygons,self.coastpolygontypes,coastpolygonsll):
                x = poly[0]
                y = poly[1]
                lons = polyll[0]
                lats = polyll[1]
                mask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
                # replace values in polygons that are over the horizon.
                xsave = False
                ysave = False
                if NX.sum(mask):
                    i1,i2 = self._splitseg(x,y,mask=mask)
                    for i,j in zip(i1,i2):
                        if i and j != len(x):
                            dist,az1,alpha21=vinc_dist(self.flattening,rad,lat_0,lon_0,math.radians(lats[i]),math.radians(lons[i]))
                            lat1,lon1,az=vinc_pt(self.flattening,rad,lat_0,lon_0,az1,0.5*math.pi*rad)
                            dist,az2,alpha21=vinc_dist(self.flattening,rad,lat_0,lon_0,math.radians(lats[j]),math.radians(lons[j]))
                            lat2,lon2,az=vinc_pt(self.flattening,rad,lat_0,lon_0,az2,0.5*math.pi*rad)
                            gc = GreatCircle(self.rmajor,self.rminor,math.degrees(lon2),math.degrees(lat2),math.degrees(lon1),math.degrees(lat1))
                            npoints = int(gc.gcarclen/dtheta)+1
                            if npoints < 2: npoints=2
                            lonstmp, latstmp = gc.points(npoints)
                            xx, yy = self(lonstmp, latstmp)
                            xnew = x[i:j]
                            ynew = y[i:j]
                            xnew = x[i:j] + xx
                            ynew = y[i:j] + yy
                            coastpolygons.append((xnew,ynew))
                            coastpolygontypes.append(polytype)
                        elif i == 0:
                            xsave = x[0:j]
                            ysave = y[0:j]
                            lats_save = lats[0:j]
                            lons_save = lons[0:j]
                        elif j == len(x):
                            xnew = x[i:j] + xsave
                            ynew = y[i:j] + ysave
                            lonsnew = lons[i:j] + lons_save
                            latsnew = lats[i:j] + lats_save
                            dist,az1,alpha21=vinc_dist(self.flattening,rad,lat_0,lon_0,math.radians(latsnew[0]),math.radians(lonsnew[0]))
                            lat1,lon1,az=vinc_pt(self.flattening,rad,lat_0,lon_0,az1,0.5*math.pi*rad)
                            dist,az2,alpha21=vinc_dist(self.flattening,rad,lat_0,lon_0,math.radians(latsnew[-1]),math.radians(lonsnew[-1]))
                            lat2,lon2,az=vinc_pt(self.flattening,rad,lat_0,lon_0,az2,0.5*math.pi*rad)
                            gc = GreatCircle(self.rmajor,self.rminor,math.degrees(lon2),math.degrees(lat2),math.degrees(lon1),math.degrees(lat1))
                            npoints = int(gc.gcarclen/dtheta)+1
                            if npoints < 2: npoints=2
                            lonstmp, latstmp = gc.points(npoints)
                            xx, yy = self(lonstmp, latstmp)
                            xnew = xnew + xx
                            ynew = ynew + yy
                            coastpolygons.append((xnew,ynew))
                            coastpolygontypes.append(polytype)
                else:
                    coastpolygons.append(poly)
                    coastpolygontypes.append(polytype)
            self.coastpolygons = coastpolygons
            self.coastpolygontypes = coastpolygontypes
        elif self.projection in ['moll','robin','sinu']:
            lon_0 = self.projparams['lon_0']
            coastpolygons=[]
            for poly,polytype,polyll in zip(self.coastpolygons,self.coastpolygontypes,coastpolygonsll):
                x = poly[0]
                y = poly[1]
                lons = polyll[0]
                lats = polyll[1]
                xn=[]
                yn=[]
                # antarctic segment goes from 360 back to 0
                # reorder to go from lon_0-180 to lon_0+180.
                if lats[-1] < -68.0:
                    lons.reverse()
                    lats.reverse()
                    xx,yy = self(lons,lats)
                    xx = NX.array(xx); yy = NX.array(yy)
                    xdist = NX.fabs(xx[1:]-xx[0:-1])
                    if max(xdist) > 1000000:
                        nmin = NX.argmax(xdist)+1
                        xnew = NX.zeros(len(xx),NX.Float64)
                        ynew = NX.zeros(len(xx),NX.Float64)
                        lonsnew = len(xx)*[0]
                        latsnew = len(xx)*[0]
                        xnew[0:len(xx)-nmin] = xx[nmin:]
                        ynew[0:len(xx)-nmin] = yy[nmin:]
                        xnew[len(xx)-nmin:] = xx[0:nmin]
                        ynew[len(xx)-nmin:] = yy[0:nmin]
                        lonsnew[0:len(xx)-nmin] = lons[nmin:]
                        latsnew[0:len(xx)-nmin] = lats[nmin:]
                        lonsnew[len(xx)-nmin:] = lons[0:nmin]
                        latsnew[len(xx)-nmin:] = lats[0:nmin]
                        x = xnew.tolist(); y = ynew.tolist()
                        lons = lonsnew; lats = latsnew
                    else:
                        x.reverse()
                        y.reverse()
                    # close polygon (add lines along left and right edges down to S pole)
                    for phi in NX.arange(-89.999,lats[0],0.1):
                        xx,yy = self(lon_0-179.99,phi)
                        xn.append(xx); yn.append(yy)
                    xn = xn+x
                    yn = yn+y
                    for phi in NX.arange(lats[-1],-89.999,-0.1):
                        xx,yy = self(lon_0+179.99,phi)
                        xn.append(xx); yn.append(yy)
                # move points outside map to edge of map
                # along constant latitude.
                else:
                    for x,y,lon,lat in zip(x,y,lons,lats):
                        if lon > lon_0+180 or lon < lon_0-180:
                            if lon >= lon_0+180: lon=lon_0+180.
                            if lon <= lon_0-180: lon=lon_0-180.
                            xx,yy = self(lon,lat)
                            xn.append(xx); yn.append(yy)
                        else:
                            xn.append(x); yn.append(y)
                coastpolygons.append((xn,yn))
            self.coastpolygons = coastpolygons

    def _splitseg(self,xx,yy,mask=None):
        """split segment up around missing values (outside projection limb)"""
        if mask is None:
            mask = (NX.logical_or(NX.greater_equal(xx,1.e20),NX.greater_equal(yy,1.e20))).tolist()
        i1=[]; i2=[]
        mprev = 1
        for i,m in enumerate(mask):
            if not m and mprev:
                i1.append(i)
            if m and not mprev:
                i2.append(i)
            mprev = m
        if not mprev: i2.append(len(mask))
        if len(i1) != len(i2):
            raise ValueError,'error in splitting boundary segments'
        return i1,i2

    def _insidemap_seg(self,seg):
        """returns True if any point in segment is inside map region"""
        xx = [x for x,y in seg]
        yy = [y for x,y in seg]
        isin = False
        for x,y in zip(xx,yy):
            if x >= self.xmin and x <= self.xmax and y >= self.ymin and y <= self.ymax:
                isin = True
                break
        return isin

    def _insidemap_poly(self,poly,polyll):
        """returns True if any point in polygon is inside map region"""
        isin = False
        xx = poly[0]; yy = poly[1]
        if self.projection in ['moll','robin','sinu']:
            lon_0 = self.projparams['lon_0']
            lons = polyll[0]
            for lon in lons:
                if lon < lon_0+180 and lon > lon_0-180:
                    isin = True
                    break
        else:
            for x,y in zip(xx,yy):
                if x >= self.xmin and x <= self.xmax and y >= self.ymin and y <= self.ymax:
                    isin = True
                    break
        return isin

    def __call__(self,x,y,inverse=False):
        """
 Calling a Basemap class instance with the arguments lon, lat will
 convert lon/lat (in degrees) to x/y native map projection
 coordinates (in meters).  If optional keyword 'inverse' is
 True (default is False), the inverse transformation from x/y
 to lon/lat is performed.

 For cylindrical equidistant projection ('cyl'), this
 does nothing (i.e. x,y == lon,lat).

 For non-cylindrical projections, the inverse transformation
 always returns longitudes between -180 and 180 degrees. For
 cylindrical projections (self.projection == 'cyl','mill' or 'merc')
 the inverse transformation will return longitudes between
 self.llcrnrlon and self.llcrnrlat.

 input arguments lon, lat can be either scalar floats or N arrays.
        """
        return self.projtran(x,y,inverse=inverse)

    def makegrid(self,nx,ny,returnxy=False):
        """
 return arrays of shape (ny,nx) containing lon,lat coordinates of
 an equally spaced native projection grid.
 if returnxy = True, the x,y values of the grid are returned also.
        """
        return self.projtran.makegrid(nx,ny,returnxy=returnxy)

    def drawmapboundary(self,color='k',linewidth=1.0,ax=None):
        """
 draw boundary around map projection region. If ax=None (default),
 default axis instance is used, otherwise specified axis instance is used.
        """
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        x = []
        y = []
        dtheta = 0.01
        if self.projection == 'ortho': # circular region.
            r = (2.*self.rmajor+self.rminor)/3.
            for az in NX.arange(0.,2.*math.pi+dtheta,dtheta):
                x.append(r*math.cos(az)+0.5*self.xmax)
                y.append(r*math.sin(az)+0.5*self.ymax)
            xy = zip(x,y)
            bound = Polygon(xy,edgecolor=color,linewidth=linewidth)
            ax.add_collection(bound)
            bound.set_fill(False)
            bound.set_clip_on(False)
        elif self.projection in ['moll','robin','sinu']:  # elliptical region.
            # left side
            lats = NX.arange(-89.9,89.9+dtheta,dtheta).tolist()
            lons = len(lats)*[self.projparams['lon_0']-179.9]
            x,y = self(lons,lats)
            # top.
            lons = NX.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
            lats = len(lons)*[89.9]
            xx,yy = self(lons,lats)
            x = x+xx; y = y+yy
            # right side
            lats = NX.arange(89.9,-89.9-dtheta,-dtheta).tolist()
            lons = len(lats)*[self.projparams['lon_0']+179.9]
            xx,yy = self(lons,lats)
            x = x+xx; y = y+yy
            # bottom.
            lons = NX.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
            lats = len(lons)*[-89.9]
            xx,yy = self(lons,lats)
            x = x+xx; y = y+yy
            xy = zip(x,y)
            poly = Polygon(xy,edgecolor=color,linewidth=linewidth)
            ax.add_patch(poly)
            poly.set_fill(False)
            poly.set_clip_on(False)
        else: # all other projections are rectangular.
            ax.axesPatch.set_linewidth(linewidth)
            ax.axesPatch.set_facecolor(ax.get_axis_bgcolor())
            ax.axesPatch.set_edgecolor(color)
            ax.set_frame_on(True)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def fillcontinents(self,color='0.8',ax=None):
        """
 Fill continents.

 color - color to fill continents (default gray).
 ax - axes instance (overrides default axes instance)

 After filling continents, lakes are re-filled with axis background color.
        """
        if self.resolution is None:
            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance' 
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        # get axis background color.
        axisbgc = ax.get_axis_bgcolor()
        np = 0
        for x,y in self.coastpolygons:
            xa = NX.array(x,NX.Float32)
            ya = NX.array(y,NX.Float32)
        # check to see if all four corners of domain in polygon (if so,
        # don't draw since it will just fill in the whole map).
            delx = 10; dely = 10
            if self.projection in ['cyl']:
                delx = 0.1
                dely = 0.1
            test1 = NX.fabs(xa-self.urcrnrx) < delx
            test2 = NX.fabs(xa-self.llcrnrx) < delx
            test3 = NX.fabs(ya-self.urcrnry) < dely
            test4 = NX.fabs(ya-self.llcrnry) < dely
            hasp1 = sum(test1*test3)
            hasp2 = sum(test2*test3)
            hasp4 = sum(test2*test4)
            hasp3 = sum(test1*test4)
            if not hasp1 or not hasp2 or not hasp3 or not hasp4:
                xy = zip(xa.tolist(),ya.tolist())
                if self.coastpolygontypes[np] != 2:
                    poly = Polygon(xy,facecolor=color,edgecolor=color,linewidth=0)
                else: # lakes filled with background color.
                    poly = Polygon(xy,facecolor=axisbgc,edgecolor=axisbgc,linewidth=0)
                ax.add_patch(poly)
            np = np + 1
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def drawcoastlines(self,linewidth=1.,color='k',antialiased=1,ax=None):
        """
 Draw coastlines.

 linewidth - coastline width (default 1.)
 color - coastline color (default black)
 antialiased - antialiasing switch for coastlines (default True).
 ax - axes instance (overrides default axes instance)
        """
        if self.resolution is None:
            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance' 
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        coastlines = LineCollection(self.coastsegs,antialiaseds=(antialiased,))
        coastlines.set_color(color)
        coastlines.set_linewidth(linewidth)
        ax.add_collection(coastlines)
        # make sure axis ticks are turned off
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def drawcountries(self,linewidth=0.5,color='k',antialiased=1,ax=None):
        """
 Draw country boundaries.

 linewidth - country boundary line width (default 0.5)
 color - country boundary line color (default black)
 antialiased - antialiasing switch for country boundaries (default True).
 ax - axes instance (overrides default axes instance)
        """
        if self.resolution is None:
            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance' 
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        coastlines = LineCollection(self.cntrysegs,antialiaseds=(antialiased,))
        coastlines.set_color(color)
        coastlines.set_linewidth(linewidth)
        ax.add_collection(coastlines)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def drawstates(self,linewidth=0.5,color='k',antialiased=1,ax=None):
        """
 Draw state boundaries in Americas.

 linewidth - state boundary line width (default 0.5)
 color - state boundary line color (default black)
 antialiased - antialiasing switch for state boundaries (default True).
 ax - axes instance (overrides default axes instance)
        """
        if self.resolution is None:
            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance' 
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        coastlines = LineCollection(self.statesegs,antialiaseds=(antialiased,))
        coastlines.set_color(color)
        coastlines.set_linewidth(linewidth)
        ax.add_collection(coastlines)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def drawrivers(self,linewidth=0.5,color='k',antialiased=1,ax=None):
        """
 Draw major rivers.

 linewidth - river boundary line width (default 0.5)
 color - river boundary line color (default black)
 antialiased - antialiasing switch for river boundaries (default True).
 ax - axes instance (overrides default axes instance)
        """
        if self.resolution is None:
            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance' 
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        coastlines = LineCollection(self.riversegs,antialiaseds=(antialiased,))
        coastlines.set_color(color)
        coastlines.set_linewidth(linewidth)
        ax.add_collection(coastlines)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])

    def readshapefile(self,shapefile,name,drawbounds=True,
                      linewidth=0.5,color='k',antialiased=1,ax=None):
        """
 read in shape file, draw boundaries on map.

 Restrictions:
  - Assumes shapes are 2D
  - vertices must be in geographic (lat/lon) coordinates.

 shapefile - path to shapefile components.  Example:
  shapefile='/home/jeff/esri/world_borders' assumes that
  world_borders.shp, world_borders.shx and world_borders.dbf
  live in /home/jeff/esri.
 name - name for Basemap attribute to hold the shapefile
  vertices in native map projection coordinates.
  Class attribute name+'_info' is a list of dictionaries, one
  for each shape, containing attributes of each shape from dbf file.
  For example, if name='counties', self.counties
  will be a list of vertices for each shape in map projection
  coordinates and self.counties_info will be a list of dictionaries
  with shape attributes. Rings in individual shapes are split out
  into separate polygons.  Additional keys
  'RINGNUM' and 'SHAPENUM' are added to shape attribute dictionary.
 drawbounds - draw boundaries of shapes (default True)
 linewidth - shape boundary line width (default 0.5)
 color - shape boundary line color (default black)
 antialiased - antialiasing switch for shape boundaries (default True).
 ax - axes instance (overrides default axes instance)

 returns a tuple (num_shapes, type, min, max) containing shape file info.
 num_shapes is the number of shapes, type is the type code (one of
 the SHPT* constants defined in the shapelib module, see
 http://shapelib.maptools.org/shp_api.html) and min and
 max are 4-element lists with the minimum and maximum values of the
 vertices.
        """
        # open shapefile, read vertices for each object, convert
        # to map projection coordinates (only works for 2D shape types).
        try:
            shp = ShapeFile(shapefile)
        except:
            raise IOError, 'error reading shapefile %s.shp' % shapefile
        try:
            dbf = dbflib.open(shapefile)
        except:
            raise IOError, 'error reading dbffile %s.dbf' % shapefile
        info = shp.info()
        if info[1] not in [1,3,5,8]:
            raise ValueError, 'readshapefile can only handle 2D shape types'
        shpsegs = []
        shpinfo = []
        for npoly in range(shp.info()[0]):
            shp_object = shp.read_object(npoly)
            verts = shp_object.vertices()
            rings = len(verts)
            for ring in range(rings):
                lons, lats = zip(*verts[ring])
                if max(lons) > 721. or min(lons) < -721. or max(lats) > 91. or min(lats) < -91:
                    msg="""
shapefile must have lat/lon vertices  - it looks like this one has vertices
in map projection coordinates. You can convert the shapefile to geographic
coordinates using the shpproj utility from the shapelib tools
(http://shapelib.maptools.org/shapelib-tools.html)"""
                    raise ValueError,msg
                x, y = self(lons, lats)
                shpsegs.append(zip(x,y))
                if ring == 0:
                    shapedict = dbf.read_record(npoly)
                # add information about ring number to dictionary.
                shapedict['RINGNUM'] = ring+1
                shapedict['SHAPENUM'] = npoly+1
                shpinfo.append(shapedict)
        # draw shape boundaries using LineCollection.
        if drawbounds:
            # get current axes instance (if none specified).
            if ax is None and self.ax is None:
                try:
                    ax = pylab.gca()
                except:
                    import pylab
                    ax = pylab.gca()
            elif ax is None and self.ax is not None:
                ax = self.ax
            # make LineCollections for each polygon.
            lines = LineCollection(shpsegs,antialiaseds=(1,))
            lines.set_color(color)
            lines.set_linewidth(linewidth)
            ax.add_collection(lines)
            # make sure axis ticks are turned off
            if self.noticks == True:
                ax.set_xticks([])
                ax.set_yticks([])
            # set axes limits to fit map region.
            self.set_axes_limits(ax=ax)
        # save segments/polygons and shape attribute dicts as class attributes.
        self.__dict__[name]=shpsegs
        self.__dict__[name+'_info']=shpinfo
        shp.close()
        dbf.close()
        return info

    def drawparallels(self,circles,color='k',linewidth=1., \
                      linestyle='--',dashes=[1,1],labels=[0,0,0,0], \
                      fmt='%g',xoffset=None,yoffset=None,ax=None,**kwargs):
        """
 draw parallels (latitude lines).

 circles - list containing latitude values to draw (in degrees).
 color - color to draw parallels (default black).
 linewidth - line width for parallels (default 1.)
 linestyle - line style for parallels (default '--', i.e. dashed).
 dashes - dash pattern for parallels (default [1,1], i.e. 1 pixel on,
  1 pixel off).
 labels - list of 4 values (default [0,0,0,0]) that control whether
  parallels are labelled where they intersect the left, right, top or
  bottom of the plot. For example labels=[1,0,0,1] will cause parallels
  to be labelled where they intersect the left and bottom of the plot,
  but not the right and top.
 fmt is a format string to format the parallel labels (default '%g').
 xoffset - label offset from edge of map in x-direction
  (default is 0.01 times width of map in map projection coordinates).
 yoffset - label offset from edge of map in y-direction
  (default is 0.01 times height of map in map projection coordinates).
 ax - axes instance (overrides default axes instance)

 additional keyword arguments control text properties for labels (see
  pylab.text documentation)
        """
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        # don't draw meridians past latmax, always draw parallel at latmax.
        latmax = 80.
        # offset for labels.
        if yoffset is None:
            yoffset = (self.urcrnry-self.llcrnry)/100.
            if self.aspect > 1:
                yoffset = self.aspect*yoffset
            else:
                yoffset = yoffset/self.aspect
        if xoffset is None:
            xoffset = (self.urcrnrx-self.llcrnrx)/100.

        if self.projection in ['merc','cyl','mill','moll','robin','sinu']:
            lons = NX.arange(self.llcrnrlon,self.urcrnrlon+0.1,0.1).astype(NX.Float32)
        else:
            lons = NX.arange(0,360.1,0.1).astype(NX.Float32)
        # make sure latmax degree parallel is drawn if projection not merc or cyl or miller
        try:
            circlesl = circles.tolist()
        except:
            circlesl = circles
        if self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
            if max(circlesl) > 0 and latmax not in circlesl:
                circlesl.append(latmax)
            if min(circlesl) < 0 and -latmax not in circlesl:
                circlesl.append(-latmax)
        xdelta = 0.01*(self.xmax-self.xmin)
        ydelta = 0.01*(self.ymax-self.ymin)
        for circ in circlesl:
            lats = circ*NX.ones(len(lons),NX.Float32)
            x,y = self(lons,lats)
            # remove points outside domain.
            testx = NX.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta)
            x = NX.compress(testx, x)
            y = NX.compress(testx, y)
            testy = NX.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta)
            x = NX.compress(testy, x)
            y = NX.compress(testy, y)
            if len(x) > 1 and len(y) > 1:
                # split into separate line segments if necessary.
                # (not necessary for mercator or cylindrical or miller).
                xd = (x[1:]-x[0:-1])**2
                yd = (y[1:]-y[0:-1])**2
                dist = NX.sqrt(xd+yd)
                split = dist > 500000.
                if NX.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
                   ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
                   xl = []
                   yl = []
                   iprev = 0
                   ind.append(len(xd))
                   for i in ind:
                       xl.append(x[iprev:i])
                       yl.append(y[iprev:i])
                       iprev = i
                else:
                    xl = [x]
                    yl = [y]
                # draw each line segment.
                for x,y in zip(xl,yl):
                    # skip if only a point.
                    if len(x) > 1 and len(y) > 1:
                        l = Line2D(x,y,linewidth=linewidth,linestyle=linestyle)
                        l.set_color(color)
                        l.set_dashes(dashes)
                        ax.add_line(l)
        # draw labels for parallels
        # parallels not labelled for orthographic, robinson, 
        # sinusoidal or mollweide.
        if self.projection in ['ortho'] and max(labels):
            print 'Warning: Cannot label parallels on Orthographic basemap'
            labels = [0,0,0,0]
        # search along edges of map to see if parallels intersect.
        # if so, find x,y location of intersection and draw a label there.
        if self.projection == 'cyl':
            dx = 0.01; dy = 0.01
        else:
            dx = 1000; dy = 1000
        if self.projection in ['moll','robin','sinu']:
            lon_0 = self.projparams['lon_0']
        for dolab,side in zip(labels,['l','r','t','b']):
            if not dolab: continue
            # for cylindrical projections, don't draw parallels on top or bottom.
            if self.projection in ['cyl','merc','mill','moll','robin','sinu'] and side in ['t','b']: continue
            if side in ['l','r']:
                nmax = int((self.ymax-self.ymin)/dy+1)
                yy = linspace(self.llcrnry,self.urcrnry,nmax)
                # mollweide inverse transform undefined at South Pole
                if self.projection == 'moll' and yy[0] < 1.e-4:
                    yy[0] = 1.e-4
                if side == 'l':
                    lons,lats = self(self.llcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(self.urcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                if max(lons) > 1.e20 or max(lats) > 1.e20:
                    raise ValueError,'inverse transformation undefined - please adjust the map projection region'
                # adjust so 0 <= lons < 360
                lons = [(lon+360) % 360 for lon in lons]
            else:
                nmax = int((self.xmax-self.xmin)/dx+1)
                xx = linspace(self.llcrnrx,self.urcrnrx,nmax)
                if side == 'b':
                    lons,lats = self(xx,self.llcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(xx,self.urcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                if max(lons) > 1.e20 or max(lats) > 1.e20:
                    raise ValueError,'inverse transformation undefined - please adjust the map projection region'
                # adjust so 0 <= lons < 360
                lons = [(lon+360) % 360 for lon in lons]
            for lat in circles:
                # find index of parallel (there may be two, so
                # search from left and right).
                nl = _searchlist(lats,lat)
                nr = _searchlist(lats[::-1],lat)
                if nr != -1: nr = len(lons)-nr-1
                if lat<0:
                    if rcParams['text.usetex']:
                        latlabstr = r'${%s\/^{\circ}\/S}$'%fmt
                    else:
                        latlabstr = u'%s\N{DEGREE SIGN}S'%fmt
                    latlab = latlabstr%NX.fabs(lat)
                elif lat>0:
                    if rcParams['text.usetex']:
                        latlabstr = r'${%s\/^{\circ}\/N}$'%fmt
                    else:
                        latlabstr = u'%s\N{DEGREE SIGN}N'%fmt
                    latlab = latlabstr%lat
                else:
                    if rcParams['text.usetex']:
                        latlabstr = r'${%s\/^{\circ}}$'%fmt
                    else:
                        latlabstr = u'%s\N{DEGREE SIGN}'%fmt
                    latlab = latlabstr%lat
                # parallels can intersect each map edge twice.
                for i,n in enumerate([nl,nr]):
                    # don't bother if close to the first label.
                    if i and abs(nr-nl) < 100: continue
                    if n >= 0:
                        if side == 'l':
                            if self.projection in ['moll','robin','sinu']:
                                xlab,ylab = self(lon_0-179.9,lat)
                            else:
                                xlab = self.llcrnrx
                            xlab = xlab-xoffset
                            ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='center',**kwargs)
                        elif side == 'r':
                            if self.projection in ['moll','robin','sinu']:
                                xlab,ylab = self(lon_0+179.9,lat)
                            else:
                                xlab = self.urcrnrx
                            xlab = xlab+xoffset
                            ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='center',**kwargs)
                        elif side == 'b':
                            ax.text(xx[n],self.llcrnry-yoffset,latlab,horizontalalignment='center',verticalalignment='top',**kwargs)
                        else:
                            ax.text(xx[n],self.urcrnry+yoffset,latlab,horizontalalignment='center',verticalalignment='bottom',**kwargs)

        # make sure axis ticks are turned off is parallels labelled.
        if self.noticks or max(labels):
            ax.set_xticks([])
            ax.set_yticks([])
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def drawmeridians(self,meridians,color='k',linewidth=1., \
                      linestyle='--',dashes=[1,1],labels=[0,0,0,0],\
                      fmt='%g',xoffset=None,yoffset=None,ax=None,**kwargs):
        """
 draw meridians (longitude lines).

 meridians - list containing longitude values to draw (in degrees).
 color - color to draw meridians (default black).
 linewidth - line width for meridians (default 1.)
 linestyle - line style for meridians (default '--', i.e. dashed).
 dashes - dash pattern for meridians (default [1,1], i.e. 1 pixel on,
  1 pixel off).
 labels - list of 4 values (default [0,0,0,0]) that control whether
  meridians are labelled where they intersect the left, right, top or
  bottom of the plot. For example labels=[1,0,0,1] will cause meridians
  to be labelled where they intersect the left and bottom of the plot,
  but not the right and top.
 fmt is a format string to format the meridian labels (default '%g').
 xoffset - label offset from edge of map in x-direction
  (default is 0.01 times width of map in map projection coordinates).
 yoffset - label offset from edge of map in y-direction
  (default is 0.01 times height of map in map projection coordinates).
 ax - axes instance (overrides default axes instance)

 additional keyword arguments control text properties for labels (see
  pylab.text documentation)
        """
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        # don't draw meridians past latmax, always draw parallel at latmax.
        latmax = 80. # not used for cyl, merc or miller projections.
        # offset for labels.
        if yoffset is None:
            yoffset = (self.urcrnry-self.llcrnry)/100.
            if self.aspect > 1:
                yoffset = self.aspect*yoffset
            else:
                yoffset = yoffset/self.aspect
        if xoffset is None:
            xoffset = (self.urcrnrx-self.llcrnrx)/100.

        if self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
            lats = NX.arange(-latmax,latmax+0.1,0.1).astype(NX.Float32)
        else:
            lats = NX.arange(-90,90.1,0.1).astype(NX.Float32)
        xdelta = 0.1*(self.xmax-self.xmin)
        ydelta = 0.1*(self.ymax-self.ymin)
        for merid in meridians:
            lons = merid*NX.ones(len(lats),NX.Float32)
            x,y = self(lons,lats)
            # remove points outside domain.
            testx = NX.logical_and(x>=self.xmin-xdelta,x<=self.xmax+xdelta)
            x = NX.compress(testx, x)
            y = NX.compress(testx, y)
            testy = NX.logical_and(y>=self.ymin-ydelta,y<=self.ymax+ydelta)
            x = NX.compress(testy, x)
            y = NX.compress(testy, y)
            if len(x) > 1 and len(y) > 1:
                # split into separate line segments if necessary.
                # (not necessary for mercator or cylindrical or miller).
                xd = (x[1:]-x[0:-1])**2
                yd = (y[1:]-y[0:-1])**2
                dist = NX.sqrt(xd+yd)
                split = dist > 500000.
                if NX.sum(split) and self.projection not in ['merc','cyl','mill','moll','robin','sinu']:
                   ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
                   xl = []
                   yl = []
                   iprev = 0
                   ind.append(len(xd))
                   for i in ind:
                       xl.append(x[iprev:i])
                       yl.append(y[iprev:i])
                       iprev = i
                else:
                    xl = [x]
                    yl = [y]
                # draw each line segment.
                for x,y in zip(xl,yl):
                    # skip if only a point.
                    if len(x) > 1 and len(y) > 1:
                        l = Line2D(x,y,linewidth=linewidth,linestyle=linestyle)
                        l.set_color(color)
                        l.set_dashes(dashes)
                        ax.add_line(l)
        # draw labels for meridians.
        # meridians not labelled for orthographic, sinusoidal,
        # robinson or mollweide
        if self.projection in ['sinu','ortho','moll'] and max(labels):
            print 'Warning: Cannot label meridians on Sinusoidal, Mollweide or Orthographic basemap'
            labels = [0,0,0,0]
        # search along edges of map to see if parallels intersect.
        # if so, find x,y location of intersection and draw a label there.
        if self.projection == 'cyl':
            dx = 0.01; dy = 0.01
        else:
            dx = 1000; dy = 1000
        if self.projection in ['moll','sinu','robin']:
            lon_0 = self.projparams['lon_0']
            xmin,ymin = self(lon_0-179.9,-90)
            xmax,ymax = self(lon_0+179.9,90)
        for dolab,side in zip(labels,['l','r','t','b']):
            if not dolab: continue
            # for cylindrical projections, don't draw meridians on left or right.
            if self.projection in ['cyl','merc','mill','sinu','robin','moll'] and side in ['l','r']: continue
            if side in ['l','r']:
                nmax = int((self.ymax-self.ymin)/dy+1)
                yy = linspace(self.llcrnry,self.urcrnry,nmax)
                if side == 'l':
                    lons,lats = self(self.llcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(self.urcrnrx*NX.ones(yy.shape,NX.Float32),yy,inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                if max(lons) > 1.e20 or max(lats) > 1.e20:
                    raise ValueError,'inverse transformation undefined - please adjust the map projection region'
                # adjust so 0 <= lons < 360
                lons = [(lon+360) % 360 for lon in lons]
            else:
                nmax = int((self.xmax-self.xmin)/dx+1)
                xx = linspace(self.llcrnrx,self.urcrnrx,nmax)
                if side == 'b':
                    lons,lats = self(xx,self.llcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(xx,self.urcrnry*NX.ones(xx.shape,NX.Float32),inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                if max(lons) > 1.e20 or max(lats) > 1.e20:
                    raise ValueError,'inverse transformation undefined - please adjust the map projection region'
                # adjust so 0 <= lons < 360
                lons = [(lon+360) % 360 for lon in lons]
            for lon in meridians:
                # adjust so 0 <= lon < 360
                lon = (lon+360) % 360
                # find index of meridian (there may be two, so
                # search from left and right).
                nl = _searchlist(lons,lon)
                nr = _searchlist(lons[::-1],lon)
                if nr != -1: nr = len(lons)-nr-1
                if lon>180:
                    if rcParams['text.usetex']:
                        lonlabstr = r'${%s\/^{\circ}\/W}$'%fmt
                    else:
                        lonlabstr = u'%s\N{DEGREE SIGN}W'%fmt
                    lonlab = lonlabstr%NX.fabs(lon-360)
                elif lon<180 and lon != 0:
                    if rcParams['text.usetex']:
                        lonlabstr = r'${%s\/^{\circ}\/E}$'%fmt
                    else:
                        lonlabstr = u'%s\N{DEGREE SIGN}E'%fmt
                    lonlab = lonlabstr%lon
                else:
                    if rcParams['text.usetex']:
                        lonlabstr = r'${%s\/^{\circ}}$'%fmt
                    else:
                        lonlabstr = u'%s\N{DEGREE SIGN}'%fmt
                    lonlab = lonlabstr%lon
                # meridians can intersect each map edge twice.
                for i,n in enumerate([nl,nr]):
                    lat = lats[n]/100.
                    # no meridians > latmax for projections other than merc,cyl,miller.
                    if self.projection not in ['merc','cyl','mill'] and lat > latmax: continue
                    # don't bother if close to the first label.
                    if i and abs(nr-nl) < 100: continue
                    if n >= 0:
                        if side == 'l':
                            ax.text(self.llcrnrx-xoffset,yy[n],lonlab,horizontalalignment='right',verticalalignment='center',**kwargs)
                        elif side == 'r':
                            ax.text(self.urcrnrx+xoffset,yy[n],lonlab,horizontalalignment='left',verticalalignment='center',**kwargs)
                        elif side == 'b':
                            if self.projection != 'robin' or (xx[n] > xmin and xx[n] < xmax):
                                ax.text(xx[n],self.llcrnry-yoffset,lonlab,horizontalalignment='center',verticalalignment='top',**kwargs)
                        else:
                            if self.projection != 'robin' or (xx[n] > xmin and xx[n] < xmax):
                                ax.text(xx[n],self.urcrnry+yoffset,lonlab,horizontalalignment='center',verticalalignment='bottom',**kwargs)

        # make sure axis ticks are turned off if meridians labelled.
        if self.noticks or max(labels):
            ax.set_xticks([])
            ax.set_yticks([])
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)

    def gcpoints(self,lon1,lat1,lon2,lat2,npoints):
        """
 compute npoints points along a great circle with endpoints
 (lon1,lat1) and (lon2,lat2).  Returns arrays x,y
 with map projection coordinates.
        """
        gc = GreatCircle(self.rmajor,self.rminor,lon1,lat1,lon2,lat2)
        lons, lats = gc.points(npoints)
        x, y = self(lons, lats)
        return x,y

    def drawgreatcircle(self,lon1,lat1,lon2,lat2,dtheta=0.02,**kwargs):
        """
 draw a great circle on the map.

 lon1,lat1 - longitude,latitude of one endpoint of the great circle.
 lon2,lat2 - longitude,latitude of the other endpoint of the great circle.
 dtheta - points on great circle computed every dtheta radians (default 0.02).

 Other keyword arguments (**kwargs) control plotting of great circle line,
 see pylab.plot documentation for details.

 Note:  cannot handle situations in which the great circle intersects
 the edge of the map projection domain, and then re-enters the domain.
        """
        # use great circle formula for a perfect sphere.
        gc = GreatCircle(self.rmajor,self.rminor,lon1,lat1,lon2,lat2)
        if gc.antipodal:
            raise ValueError,'cannot draw great circle whose endpoints are antipodal'
        # points have spacing of dtheta radians.
        npoints = int(gc.gcarclen/dtheta)+1
        lons, lats = gc.points(npoints)
        x, y = self(lons, lats)
        self.plot(x,y,**kwargs)

    def transform_scalar(self,datin,lons,lats,nx,ny,returnxy=False,checkbounds=False,order=1,masked=False):
        """
 interpolate a scalar field (datin) from a lat/lon grid with longitudes =
 lons and latitudes = lats to a (ny,nx) native map projection grid.
 Typically used to transform data to map projection coordinates
 so it can be plotted on the map with imshow.

 lons, lats must be rank-1 arrays containing longitudes and latitudes
 (in degrees) of datin grid in increasing order
 (i.e. from dateline eastward, and South Pole northward).
 For non-cylindrical projections (those other than
 cylindrical equidistant, mercator and miller) 
 lons must fit within range -180 to 180.

 if returnxy=True, the x and y values of the native map projection grid
 are also returned.

 If checkbounds=True, values of lons and lats are checked to see that
 they lie within the map projection region.  Default is False.
 If checkbounds=False, points outside map projection region will
 be clipped to values on the boundary if masked=False. If masked=True,
 the return value will be a masked array with those points masked.

 The order keyword can be 0 for nearest-neighbor interpolation,
 or 1 for bilinear interpolation (default 1).
        """
        # check that lons, lats increasing
        delon = lons[1:]-lons[0:-1]
        delat = lats[1:]-lats[0:-1]
        if min(delon) < 0. or min(delat) < 0.:
            raise ValueError, 'lons and lats must be increasing!'
        # check that lons in -180,180 for non-cylindrical projections.
        if self.projection not in ['cyl','merc','mill']:
            lonsa = NX.array(lons)
            count = NX.sum(lonsa < -180.00001) + NX.sum(lonsa > 180.00001)
            if count > 1:
                raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)'
            # allow for wraparound point to be outside.
            elif count == 1 and math.fabs(lons[-1]-lons[0]-360.) > 1.e-4:
                raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)'
        if returnxy:
            lonsout, latsout, x, y = self.makegrid(nx,ny,returnxy=True)
        else:
            lonsout, latsout = self.makegrid(nx,ny)
        if _has_arrindexing:
            datout = interp(datin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
        else:
            datout = interp_numeric(datin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
        if returnxy:
            return datout, x, y
        else:
            return datout

    def transform_vector(self,uin,vin,lons,lats,nx,ny,returnxy=False,checkbounds=False,order=1,masked=False):
        """
 rotate and interpolate a vector field (uin,vin) from a lat/lon grid
 with longitudes = lons and latitudes = lats to a
 (ny,nx) native map projection grid.

 lons, lats must be rank-1 arrays containing longitudes and latitudes
 (in degrees) of datin grid in increasing order
 (i.e. from dateline eastward, and South Pole northward).
 For non-cylindrical projections (those other than
 cylindrical equidistant, mercator and miller) 
 lons must fit within range -180 to 180.

 The input vector field is defined in spherical coordinates (it
 has eastward and northward components) while the output
 vector field is rotated to map projection coordinates (relative
 to x and y). The magnitude of the vector is preserved.

 if returnxy=True, the x and y values of the native map projection grid
 are also returned (default False).

 If checkbounds=True, values of lons and lats are checked to see that
 they lie within the map projection region.  Default is False.
 If checkbounds=False, points outside map projection region will
 be clipped to values on the boundary if masked=False. If masked=True,
 the return value will be a masked array with those points masked.

 The order keyword can be 0 for nearest-neighbor interpolation,
 or 1 for bilinear interpolation (default 1).
        """
        # check that lons, lats increasing
        delon = lons[1:]-lons[0:-1]
        delat = lats[1:]-lats[0:-1]
        if min(delon) < 0. or min(delat) < 0.:
            raise ValueError, 'lons and lats must be increasing!'
        # check that lons in -180,180 for non-cylindrical projections.
        if self.projection not in ['cyl','merc','mill']:
            lonsa = NX.array(lons)
            count = NX.sum(lonsa < -180.00001) + NX.sum(lonsa > 180.00001)
            if count > 1:
                raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)'
            # allow for wraparound point to be outside.
            elif count == 1 and math.fabs(lons[-1]-lons[0]-360.) > 1.e-4:
                raise ValueError,'grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)'
        lonsout, latsout, x, y = self.makegrid(nx,ny,returnxy=True)
        # interpolate to map projection coordinates.
        if _has_arrindexing:
            uin = interp(uin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
            vin = interp(vin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
        else:
            uin = interp_numeric(uin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
            vin = interp_numeric(vin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
        # rotate from geographic to map coordinates.
        delta = 0.1 # incement in latitude used to estimate derivatives.
        xn,yn = self(lonsout,NX.where(latsout+delta<90.,latsout+delta,latsout-delta))
        dxdlat = NX.where(latsout+delta<90.,(xn-x)/(latsout+delta),(x-xn)/(latsout+delta))
        dydlat = NX.where(latsout+delta<90.,(yn-y)/(latsout+delta),(y-yn)/(latsout+delta))
        # northangle is the angle between true north and the y axis.
        northangle = NX.arctan2(dxdlat,dydlat)
        uout = uin*NX.cos(northangle) + vin*NX.sin(northangle)
        vout = vin*NX.cos(northangle) - uin*NX.sin(northangle)
        if returnxy:
            return uout,vout,x,y
        else:
            return uout,vout

    def rotate_vector(self,uin,vin,lons,lats,returnxy=False):
        """
 rotate a vector field (uin,vin) on a rectilinear lat/lon grid
 with longitudes = lons and latitudes = lats from geographical into map
 projection coordinates.

 Differs from transform_vector in that no interpolation is done,
 the vector is returned on the same lat/lon grid, but rotated into
 x,y coordinates.

 lons, lats must be rank-2 arrays containing longitudes and latitudes
 (in degrees) of grid.

 if returnxy=True, the x and y values of the lat/lon grid
 are also returned (default False).

 The input vector field is defined in spherical coordinates (it
 has eastward and northward components) while the output
 vector field is rotated to map projection coordinates (relative
 to x and y). The magnitude of the vector is preserved.
        """
        x, y = self(lons, lats)
        # rotate from geographic to map coordinates.
        delta = 0.1 # incement in latitude used to estimate derivatives.
        xn,yn = self(lons,NX.where(lats+delta<90.,lats+delta,lats-delta))
        dxdlat = NX.where(lats+delta<90.,(xn-x)/(lats+delta),(x-xn)/(lats+delta))
        dydlat = NX.where(lats+delta<90.,(yn-y)/(lats+delta),(y-yn)/(lats+delta))
        # northangle is the angle between true north and the y axis.
        northangle = NX.arctan2(dxdlat,dydlat)
        uout = uin*NX.cos(northangle) + vin*NX.sin(northangle)
        vout = vin*NX.cos(northangle) - uin*NX.sin(northangle)
        if returnxy:
            return uout,vout,x,y
        else:
            return uout,vout

    def set_axes_limits(self,ax=None):
        """
 Set axis limits, fix aspect ratio for map domain using current
 or specified axes instance.
        """
        # get current axes instance (if none specified).
        if ax is None and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif ax is None and self.ax is not None:
            ax = self.ax
        # update data limits for map domain.
        corners = ((self.llcrnrx,self.llcrnry), (self.urcrnrx,self.urcrnry))
        ax.update_datalim( corners )
        ax.set_xlim((self.llcrnrx, self.urcrnrx))
        ax.set_ylim((self.llcrnry, self.urcrnry))
        # turn off axes frame for non-rectangular projections.
        if self.projection in ['ortho','moll','robin','sinu']:
            ax.set_frame_on(False)
        # make sure aspect ratio of map preserved.
        # plot is re-centered in bounding rectangle.
        # (anchor instance var determines where plot is placed)
        ax.set_aspect('equal',adjustable='box',anchor=self.anchor)
        ax.apply_aspect()

    def scatter(self, *args, **kwargs):
        """
 Plot points with markers on the map (see pylab.scatter documentation).
 extra keyword 'ax' can be used to override the default axes instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            ret =  ax.scatter(*args, **kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        return ret

    def plot(self, *args, **kwargs):
        """
 Draw lines and/or markers on the map (see pylab.plot documentation).
 extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            ret =  ax.plot(*args, **kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        return ret

    def imshow(self, *args, **kwargs):
        """
 Display an image over the map (see pylab.imshow documentation).
 extent and origin keywords set automatically so image will be drawn
 over map region.
 extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        kwargs['extent']=(self.llcrnrx,self.urcrnrx,self.llcrnry,self.urcrnry)
        # use origin='lower', unless overridden.
        if not kwargs.has_key('origin'):
            kwargs['origin']='lower'
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            ret =  ax.imshow(*args, **kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # reset current active image (only if pylab is imported).
        try:
            pylab.gci._current = ret
        except:
            pass
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        return ret

    def pcolor(self,x,y,data,**kwargs):
        """
 Make a pseudo-color plot over the map.
 see pylab.pcolor documentation for definition of **kwargs
 If x or y are outside projection limb (i.e. they have values > 1.e20)
 they will be convert to masked arrays with those values masked.
 As a result, those values will not be plotted.
 extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        kwargs['extent']=(self.llcrnrx,self.urcrnrx,self.llcrnry,self.urcrnry)
        # make x,y masked arrays
        # (masked where data is outside of projection limb)
        x = ma.masked_values(NX.where(x > 1.e20,1.e20,x), 1.e20)
        y = ma.masked_values(NX.where(y > 1.e20,1.e20,y), 1.e20)
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            ret =  ax.pcolor(x,y,data,**kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # reset current active image (only if pylab is imported).
        try:
            pylab.gci._current = ret
        except:
            pass
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        return ret

    def pcolormesh(self,x,y,data,**kwargs):
        """
 Make a pseudo-color plot over the map.
 see pylab.pcolormesh documentation for definition of **kwargs
 Unlike pcolor, pcolormesh cannot handle masked arrays, and so
 cannot be used to plot data when the grid lies partially outside
 the projection limb (use pcolor or contourf instead).
 extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        kwargs['extent']=(self.llcrnrx,self.urcrnrx,self.llcrnry,self.urcrnry)
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            ret =  ax.pcolormesh(x,y,data,**kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # reset current active image (only if pylab is imported).
        try:
            pylab.gci._current = ret
        except:
            pass
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        return ret

    def contour(self,x,y,data,*args,**kwargs):
        """
 Make a contour plot over the map (see pylab.contour documentation).
 extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        # mask for points outside projection limb.
        xymask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
        data = ma.asarray(data)
        # combine with data mask.
        mask = NX.logical_or(ma.getmaskarray(data),xymask)
        data = ma.masked_array(data,mask=mask)
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            CS = ax.contour(x,y,data,*args,**kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        # reset current active image (only if pylab is imported).
        try:
            try: # new contour.
                if CS._A is not None: pylab.gci._current = CS  
            except: # old contour.
                if CS[1].mappable is not None: pylab.gci._current = CS[1].mappable
        except:
            pass
        return CS

    def contourf(self,x,y,data,*args,**kwargs):
        """
 Make a filled contour plot over the map (see pylab.contourf documentation).
 If x or y are outside projection limb (i.e. they have values > 1.e20),
 the corresponing data elements will be masked.
 Extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        # mask for points outside projection limb.
        xymask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
        data = ma.asarray(data)
        # combine with data mask.
        mask = NX.logical_or(ma.getmaskarray(data),xymask)
        data = ma.masked_array(data,mask=mask)
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            CS = ax.contourf(x,y,data,*args,**kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        # reset current active image (only if pylab is imported).
        try:
            try: # new contour.
                if CS._A is not None: pylab.gci._current = CS
            except: # old contour.
                if CS[1].mappable is not None: pylab.gci._current = CS[1].mappable
        except:
            pass
        return CS

    def quiver(self, x, y, u, v, *args, **kwargs):
        """
 Make a vector plot (u, v) with arrows on the map.

 Extra arguments (*args and **kwargs) passed to quiver Axes method (see
 pylab.quiver documentation for details).
 extra keyword 'ax' can be used to override the default axis instance.
        """
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        # allow callers to override the hold state by passing hold=True|False
        b = ax.ishold()
        h = popd(kwargs, 'hold', None)
        if h is not None:
            ax.hold(h)
        try:
            ret =  ax.quiver(x,y,u,v,*args,**kwargs)
            try:
                pylab.draw_if_interactive()
            except:
                pass
        except:
            ax.hold(b)
            raise
        ax.hold(b)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        return ret

    def drawlsmask(self,rgba_land,rgba_ocean,lsmask=None,
                   lsmask_lons=None,lsmask_lats=None,lakes=False,**kwargs):
        """
 draw land-sea mask image.

 land is colored with rgba integer tuple rgba_land.
 ocean is colored with rgba integer tuple rgba_ocean.

 For example, to color oceans blue and land green, use
 rgba_ocean = (0,0,255,255) and rgba_land  = (0,255,0,255).
 To make oceans transparent (useful if you just want to mask land
 regions over another image), use rgba_ocean = (0,0,255,0).

 If lakes=True, inland lakes are also colored with 
 rgba_ocean (default is lakes=False).

 Default land-sea mask is from
 http://www.ngdc.noaa.gov/seg/cdroms/graham/graham/graham.htm
 and has 5-minute resolution.

 To specify your own global land-sea mask, set the 
 lsmask keyword to a (nlats, nlons) array
 with 0's for ocean pixels, 1's for land pixels and
 optionally 2's for inland lake pixels.
 The lsmask_lons keyword should be a 1-d array
 with the longitudes of the mask, lsmask_lats should be
 a 1-d array with the latitudes.  Longitudes should be ordered
 from -180 W eastward, latitudes from -90 S northward.
 If any of the lsmask, lsmask_lons or lsmask_lats keywords are not
 set, the default land-sea mask is used.

 extra keyword 'ax' can be used to override the default axis instance.
        """
        # look for axes instance (as keyword, an instance variable
        # or from pylab.gca().
        if not kwargs.has_key('ax') and self.ax is None:
            try:
                ax = pylab.gca()
            except:
                import pylab
                ax = pylab.gca()
        elif not kwargs.has_key('ax') and self.ax is not None:
            ax = self.ax
        else:
            ax = popd(kwargs,'ax')
        # if lsmask,lsmask_lons,lsmask_lats keywords not given,
        # read default land-sea mask in from file.
        if lsmask is None or lsmask_lons is None or lsmask_lats is None:
            # if lsmask instance variable already set, data already
            # read in.
            if self.lsmask is None:
                # read in land/sea mask.
                lsmaskf = open(os.path.join(_datadir,'5minmask.bin'),'rb')
                nlons = 4320; nlats = nlons/2
                delta = 360./float(nlons)
                lsmask_lons = NX.arange(-180+0.5*delta,180.,delta)
                lsmask_lats = NX.arange(-90.+0.5*delta,90.,delta)
                lsmask = NX.reshape(NX.fromstring(lsmaskf.read(),NX.UInt8),(nlats,nlons))
                lsmaskf.close()
        # instance variable lsmask is set on first invocation,
        # it contains the land-sea mask interpolated to the native
        # projection grid.  Further calls to drawlsmask will not
        # redo the interpolation (unless a new land-sea mask is passed
        # in via the lsmask, lsmask_lons, lsmask_lats keywords).

        # transform mask to nx x ny regularly spaced native projection grid
        # nx and ny chosen to have roughly the same horizontal
        # resolution as mask.
        if self.lsmask is None:
            nlons = len(lsmask_lons)
            nlats = len(lsmask_lats)
            if self.projection == 'cyl':
                dx = lsmask_lons[1]-lsmask_lons[0]
            else:
                dx = 2.*math.pi*self.rmajor/float(nlons)
            nx = int((self.xmax-self.xmin)/dx)+1; ny = int((self.ymax-self.ymin)/dx)+1
        # interpolate rgba values from proj='cyl' (geographic coords)
        # to a rectangular map projection grid.
            mask = self.transform_scalar(lsmask,lsmask_lons,\
                                         lsmask_lats,nx,ny,order=0,masked=255)
            self.lsmask = mask
        # optionally, set lakes to ocean color.
        if lakes:
            mask = NX.where(self.lsmask==2,0,self.lsmask)
        else:
            mask = self.lsmask
        ny, nx = mask.shape
        rgba = NX.ones((ny,nx,4),NX.UInt8)
        rgba_land = NX.array(rgba_land,NX.UInt8)
        rgba_ocean = NX.array(rgba_ocean,NX.UInt8)
        for k in range(4):
            rgba[:,:,k] = NX.where(mask,rgba_land[k],rgba_ocean[k])
        # make points outside projection limb transparent.
        try:
            rgba[:,:,3] = NX.where(mask==255,0,rgba[:,:,3])
        except: # kluge for Numeric
            rgba[:,:,3] = NX.where(mask==255,0,rgba[:,:,3]).astype(NX.UInt8)
        # plot mask as rgba image.
        im = self.imshow(rgba,interpolation='nearest',ax=ax)

def _searchlist(a,x):
    """
 like bisect, but works for lists that are not sorted,
 and are not in increasing order.
 returns -1 if x does not fall between any two elements"""
    # make sure x is a float (and not an array scalar)
    x = float(x)
    itemprev = a[0]
    nslot = -1
    eps = 180.
    for n,item in enumerate(a[1:]):
        if item < itemprev:
            if itemprev-item>eps:
                if ((x>itemprev and x<=360.) or (x<item and x>=0.)):
                    nslot = n+1
                    break
            elif x <= itemprev and x > item and itemprev:
                nslot = n+1
                break
        else:
            if item-itemprev>eps:
                if ((x<itemprev and x>=0.) or (x>item and x<=360.)):
                    nslot = n+1
                    break
            elif x >= itemprev and x < item:
                nslot = n+1
                break
        itemprev = item
    return nslot

def interp_numeric(datain,xin,yin,xout,yout,checkbounds=False,masked=False,order=1):
    """
 dataout = interp(datain,xin,yin,xout,yout,order=1)

 This version uses 'take' instead of array indexing, and
 thus is compatible with Numeric.

 interpolate data (datain) on a rectilinear grid (with x=xin
 y=yin) to a grid with x=xout, y=yout.

 datain is a rank-2 array with 1st dimension corresponding to y,
 2nd dimension x.

 xin, yin are rank-1 arrays containing x and y of
 datain grid in increasing order.

 xout, yout are rank-2 arrays containing x and y of desired output grid.

 If checkbounds=True, values of xout and yout are checked to see that
 they lie within the range specified by xin and xin.  Default is False.
 If checkbounds=False, and xout,yout are outside xin,yin, interpolated
 values will be clipped to values on boundary of input grid (xin,yin)
 if masked=False. If masked=True, the return value will be a masked
 array with those points masked. If masked is set to a number, then
 points outside the range of xin and yin will be set to that number.

 The order keyword can be 0 for nearest-neighbor interpolation,
 or 1 for bilinear interpolation (default 1).

 If datain is a masked array and order=1 (bilinear interpolation) is 
 used, elements of dataout will be masked if any of the four surrounding
 points in datain are masked.  To avoid this, do the interpolation in two
 passes, first with order=1 (producing dataout1), then with order=0
 (producing dataout2).  Then replace all the masked values in dataout1
 with the corresponding elements in dataout2 (using numerix.where).
 This effectively uses nearest neighbor interpolation if any of the
 four surrounding points in datain are masked, and bilinear interpolation
 otherwise.
    """
    # xin and yin must be monotonically increasing.
    if xin[-1]-xin[0] < 0 or yin[-1]-yin[0] < 0:
        raise ValueError, 'xin and yin must be increasing!'
    if xout.shape != yout.shape:
        raise ValueError, 'xout and yout must have same shape!'
    # check that xout,yout are
    # within region defined by xin,yin.
    if checkbounds:
        if min(NX.ravel(xout)) < min(xin) or \
           max(NX.ravel(xout)) > max(xin) or \
           min(NX.ravel(yout)) < min(yin) or \
           max(NX.ravel(yout)) > max(yin):
            raise ValueError, 'yout or xout outside range of yin or xin'
    # compute grid coordinates of output grid.
    xoutflat = NX.ravel(xout)
    youtflat = NX.ravel(yout)
    datainflat = NX.ravel(datain)
    delx = xin[1:]-xin[0:-1]
    dely = yin[1:]-yin[0:-1]
    if max(delx)-min(delx) < 1.e-4 and max(dely)-min(dely) < 1.e-4:
        # regular input grid.
        xcoords = (len(xin)-1)*(xoutflat-xin[0])/(xin[-1]-xin[0])
        ycoords = (len(yin)-1)*(youtflat-yin[0])/(yin[-1]-yin[0])
    else:
        # irregular (but still rectilinear) input grid.
        ix = (NX.searchsorted(xin,xoutflat)-1).tolist()
        iy = (NX.searchsorted(yin,youtflat)-1).tolist()
        xoutflat = xoutflat.tolist(); xin = xin.tolist()
        youtflat = youtflat.tolist(); yin = yin.tolist()
        xcoords = []; ycoords = []
        for n,i in enumerate(ix):
            if i < 0:
                xcoords.append(-1) # outside of range on xin (lower end)
            elif i >= len(xin)-1:
                xcoords.append(len(xin)) # outside range on upper end.
            else:
                xcoords.append(float(i)+(xoutflat[n]-xin[i])/(xin[i+1]-xin[i]))
        xcoords = NX.array(xcoords,NX.Float32)
        for m,j in enumerate(iy):
            if j < 0:
                ycoords.append(-1) # outside of range of yin (on lower end)
            elif j >= len(yin)-1:
                ycoords.append(len(yin)) # outside range on upper end
            else:
                ycoords.append(float(j)+(youtflat[m]-yin[j])/(yin[j+1]-yin[j]))
        ycoords = NX.array(ycoords,NX.Float32)
    # data outside range xin,yin will be clipped to
    # values on boundary.
    if masked:
        xmask = NX.logical_or(NX.less(xcoords,0),NX.greater(xcoords,len(xin)-1))
        ymask = NX.logical_or(NX.less(ycoords,0),NX.greater(ycoords,len(yin)-1))
        xymask = NX.logical_or(xmask,ymask)   
        xymask = NX.reshape(xymask,xout.shape)
    xcoords = NX.clip(xcoords,0,len(xin)-1)
    ycoords = NX.clip(ycoords,0,len(yin)-1)
    # interpolate to output grid using bilinear interpolation.
    if order == 1:
        xi = xcoords.astype(NX.Int32)
        yi = ycoords.astype(NX.Int32)
        xip1 = xi+1
        yip1 = yi+1
        xip1 = NX.clip(xip1,0,len(xin)-1)
        yip1 = NX.clip(yip1,0,len(yin)-1)
        delx = xcoords-xi.astype(NX.Float32)
        dely = ycoords-yi.astype(NX.Float32)
        coords = yi*datain.shape[1]+xi
        data11 = NX.take(datainflat,coords)
        coords = yip1*datain.shape[1]+xip1
        data22 = NX.take(datainflat,coords)
        coords = yi*datain.shape[1]+xip1
        data12 = NX.take(datainflat,coords)
        coords = yip1*datain.shape[1]+xi
        data21 = NX.take(datainflat,coords)
        dataout = (1.-delx)*(1.-dely)*data11 + \
                  delx*dely*data22 + \
                  (1.-delx)*dely*data21 + \
                  delx*(1.-dely)*data12
        dataout = NX.reshape(dataout,xout.shape)
    elif order == 0:
        xcoordsi = NX.around(xcoords).astype(NX.Int32)
        ycoordsi = NX.around(ycoords).astype(NX.Int32)
        coords = ycoordsi*datain.shape[1]+xcoordsi
        dataout = NX.take(datainflat,coords)
        dataout = NX.reshape(dataout,xout.shape)
    else:
        raise ValueError,'order keyword must be 0 or 1'
    if masked and isinstance(masked,bool):
        dataout = ma.masked_array(dataout)
        newmask = ma.mask_or(ma.getmask(dataout), xymask)
        dataout = ma.masked_array(dataout,mask=newmask)
    elif masked and is_scalar(masked):
        dataout = NX.where(xymask,masked,dataout)
    return dataout

def interp(datain,xin,yin,xout,yout,checkbounds=False,masked=False,order=1):
    """
 dataout = interp(datain,xin,yin,xout,yout,order=1)

 This version uses array indexing and is not compatible with Numeric.

 interpolate data (datain) on a rectilinear grid (with x=xin
 y=yin) to a grid with x=xout, y=yout.

 datain is a rank-2 array with 1st dimension corresponding to y,
 2nd dimension x.

 xin, yin are rank-1 arrays containing x and y of
 datain grid in increasing order.

 xout, yout are rank-2 arrays containing x and y of desired output grid.

 If checkbounds=True, values of xout and yout are checked to see that
 they lie within the range specified by xin and xin.  Default is False.
 If checkbounds=False, and xout,yout are outside xin,yin, interpolated
 values will be clipped to values on boundary of input grid (xin,yin)
 if masked=False. If masked=True, the return value will be a masked
 array with those points masked. If masked is set to a number, then
 points outside the range of xin and yin will be set to that number.

 The order keyword can be 0 for nearest-neighbor interpolation,
 or 1 for bilinear interpolation (default 1).

 If datain is a masked array and order=1 (bilinear interpolation) is 
 used, elements of dataout will be masked if any of the four surrounding
 points in datain are masked.  To avoid this, do the interpolation in two
 passes, first with order=1 (producing dataout1), then with order=0
 (producing dataout2).  Then replace all the masked values in dataout1
 with the corresponding elements in dataout2 (using numerix.where).
 This effectively uses nearest neighbor interpolation if any of the
 four surrounding points in datain are masked, and bilinear interpolation
 otherwise.
    """
    # xin and yin must be monotonically increasing.
    if xin[-1]-xin[0] < 0 or yin[-1]-yin[0] < 0:
        raise ValueError, 'xin and yin must be increasing!'
    if xout.shape != yout.shape:
        raise ValueError, 'xout and yout must have same shape!'
    # check that xout,yout are
    # within region defined by xin,yin.
    if checkbounds:
        if min(NX.ravel(xout)) < min(xin) or \
           max(NX.ravel(xout)) > max(xin) or \
           min(NX.ravel(yout)) < min(yin) or \
           max(NX.ravel(yout)) > max(yin):
            raise ValueError, 'yout or xout outside range of yin or xin'
    # compute grid coordinates of output grid.
    delx = xin[1:]-xin[0:-1]
    dely = yin[1:]-yin[0:-1]
    if max(delx)-min(delx) < 1.e-4 and max(dely)-min(dely) < 1.e-4:
        # regular input grid.
        xcoords = (len(xin)-1)*(xout-xin[0])/(xin[-1]-xin[0])
        ycoords = (len(yin)-1)*(yout-yin[0])/(yin[-1]-yin[0])
    else:
        # irregular (but still rectilinear) input grid.
        xoutflat = NX.ravel(xout); youtflat = NX.ravel(yout)
        ix = (NX.searchsorted(xin,xoutflat)-1).tolist()
        iy = (NX.searchsorted(yin,youtflat)-1).tolist()
        xoutflat = xoutflat.tolist(); xin = xin.tolist()
        youtflat = youtflat.tolist(); yin = yin.tolist()
        xcoords = []; ycoords = []
        for n,i in enumerate(ix):
            if i < 0:
                xcoords.append(-1) # outside of range on xin (lower end)
            elif i >= len(xin)-1:
                xcoords.append(len(xin)) # outside range on upper end.
            else:
                xcoords.append(float(i)+(xoutflat[n]-xin[i])/(xin[i+1]-xin[i]))
        for m,j in enumerate(iy):
            if j < 0:
                ycoords.append(-1) # outside of range of yin (on lower end)
            elif j >= len(yin)-1:
                ycoords.append(len(yin)) # outside range on upper end
            else:
                ycoords.append(float(j)+(youtflat[m]-yin[j])/(yin[j+1]-yin[j]))
        xcoords = NX.reshape(xcoords,xout.shape)
        ycoords = NX.reshape(ycoords,yout.shape)
    # data outside range xin,yin will be clipped to
    # values on boundary.
    if masked:
        xmask = NX.logical_or(NX.less(xcoords,0),NX.greater(xcoords,len(xin)-1))
        ymask = NX.logical_or(NX.less(ycoords,0),NX.greater(ycoords,len(yin)-1))
        xymask = NX.logical_or(xmask,ymask)   
    xcoords = NX.clip(xcoords,0,len(xin)-1)
    ycoords = NX.clip(ycoords,0,len(yin)-1)
    # interpolate to output grid using bilinear interpolation.
    if order == 1:
        xi = xcoords.astype(NX.Int32)
        yi = ycoords.astype(NX.Int32)
        xip1 = xi+1
        yip1 = yi+1
        xip1 = NX.clip(xip1,0,len(xin)-1)
        yip1 = NX.clip(yip1,0,len(yin)-1)
        delx = xcoords-xi.astype(NX.Float32)
        dely = ycoords-yi.astype(NX.Float32)
        dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
                  delx*dely*datain[yip1,xip1] + \
                  (1.-delx)*dely*datain[yip1,xi] + \
                  delx*(1.-dely)*datain[yi,xip1]
    elif order == 0:
        xcoordsi = NX.around(xcoords).astype(NX.Int32)
        ycoordsi = NX.around(ycoords).astype(NX.Int32)
        dataout = datain[ycoordsi,xcoordsi]
    else:
        raise ValueError,'order keyword must be 0 or 1'
    if masked and isinstance(masked,bool):
        dataout = ma.masked_array(dataout)
        newmask = ma.mask_or(ma.getmask(dataout), xymask)
        dataout = ma.masked_array(dataout,mask=newmask)
    elif masked and is_scalar(masked):
        dataout = NX.where(xymask,masked,dataout)
    return dataout

def shiftgrid(lon0,datain,lonsin,start=True):
    """
 shift global lat/lon grid east or west.
 assumes wraparound (or cyclic point) is included.

 lon0:  starting longitude for shifted grid
        (ending longitude if start=False). lon0 must be on
        input grid (with the range of lonsin).
 datain:  original data.
 lonsin:  original longitudes.
 start[True]: if True, lon0 represents the starting longitude
 of the new grid. if False, lon0 is the ending longitude.

 returns dataout,lonsout (data and longitudes on shifted grid).
    """
    if NX.fabs(lonsin[-1]-lonsin[0]-360.) > 1.e-4:
        raise ValueError, 'cyclic point not included'
    if lon0 < lonsin[0] or lon0 > lonsin[-1]:
        raise ValueError, 'lon0 outside of range of lonsin'
    i0 = NX.argsort(NX.fabs(lonsin-lon0))[0]
    try:
        dataout = NX.zeros(datain.shape,datain.typecode())
        lonsout = NX.zeros(lonsin.shape,lonsin.typecode())
    except:
        dataout = NX.zeros(datain.shape,datain.dtype)
        lonsout = NX.zeros(lonsin.shape,lonsin.dtype)
    if start:
        lonsout[0:len(lonsin)-i0] = lonsin[i0:]
    else:
        lonsout[0:len(lonsin)-i0] = lonsin[i0:]-360.
    dataout[:,0:len(lonsin)-i0] = datain[:,i0:]
    if start:
        lonsout[len(lonsin)-i0:] = lonsin[1:i0+1]+360.
    else:
        lonsout[len(lonsin)-i0:] = lonsin[1:i0+1]
    dataout[:,len(lonsin)-i0:] = datain[:,1:i0+1]
    return dataout,lonsout

def addcyclic(arrin,lonsin):
    """
 arrout, lonsout = addcyclic(arrin, lonsin)

 Add cyclic (wraparound) point in longitude.
    """
    nlats = arrin.shape[0]
    nlons = arrin.shape[1]
    try:
        arrout  = NX.zeros((nlats,nlons+1),arrin.typecode())
    except:
        arrout  = NX.zeros((nlats,nlons+1),arrin.dtype)
    arrout[:,0:nlons] = arrin[:,:]
    arrout[:,nlons] = arrin[:,0]
    try:
        lonsout = NX.zeros(nlons+1,lonsin.typecode())
    except:
        lonsout = NX.zeros(nlons+1,lonsin.dtype)
    lonsout[0:nlons] = lonsin[:]
    lonsout[nlons]  = lonsin[-1] + lonsin[1]-lonsin[0]
    return arrout,lonsout

def _choosecorners(width,height,**kwargs):
    """
private function to determine lat/lon values of projection region corners, 
given width and height of projection region in meters."""
    p = pyproj.Proj(kwargs)
    urcrnrlon, urcrnrlat = p(0.5*width,0.5*height, inverse=True)
    llcrnrlon, llcrnrlat = p(-0.5*width,-0.5*height, inverse=True)
    corners = llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat
    # test for nans,infs in output
    if isnan(llcrnrlon) or isnan(llcrnrlon-llcrnrlon) or isnan(urcrnrlon) or isnan(urcrnrlon-urcrnrlon):
       raise ValueError, 'width and/or height too large for this projection, try smaller values'
    else:
       return corners