__init__.py

from __future__ import (absolute_import, division, print_function)

"""
Module for plotting data on maps with matplotlib.

Contains the :class:`Basemap` class (which does most of the
heavy lifting), and the following functions:

:func:`interp`:  bilinear interpolation between rectilinear grids.

:func:`maskoceans`:  mask 'wet' points of an input array.

:func:`shiftgrid`:  shifts global lat/lon grids east or west.

:func:`addcyclic`: Add cyclic (wraparound) point in longitude.
"""

import os
import sys
import math
import warnings
import functools
try:
    from urllib2 import urlopen
    from urllib import urlretrieve
except ImportError:
    from urllib.request import urlopen
    from urllib.request import urlretrieve

import numpy as np
import numpy.ma as ma

import matplotlib as mpl
from matplotlib.artist import Artist
from matplotlib.collections import LineCollection
from matplotlib.collections import PolyCollection
from matplotlib.image import imread
from matplotlib.lines import Line2D
from matplotlib.patches import Circle
from matplotlib.patches import Ellipse
from matplotlib.patches import FancyArrowPatch
from matplotlib.patches import Polygon
from matplotlib.transforms import Bbox
from mpl_toolkits.axes_grid1 import make_axes_locatable

import pyproj
import _geoslib
from . proj import Proj


__version__ = "2.1.0.dev0"

# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
# check to see if environment variable BASEMAPDATA set to a directory,
# and if so look for the data there.
if 'BASEMAPDATA' in os.environ:
    basemap_datadir = os.environ['BASEMAPDATA']
    if not os.path.isdir(basemap_datadir):
        raise RuntimeError('Path in environment BASEMAPDATA not a directory')
else:
    from mpl_toolkits import basemap_data
    basemap_datadir = os.path.abspath(list(basemap_data.__path__)[0])

# module variable that sets the default value for the 'latlon' kwarg.
# can be set to True by user so plotting functions can take lons,lats
# in degrees by default, instead of x,y (map projection coords in meters).
latlon_default = False

# supported map projections.
_projnames = {'cyl'      : 'Cylindrical Equidistant',
             'merc'     : 'Mercator',
             'tmerc'    : 'Transverse Mercator',
             'omerc'    : 'Oblique Mercator',
             'mill'     : 'Miller Cylindrical',
             'gall'     : 'Gall Stereographic Cylindrical',
             'cea'      : 'Cylindrical Equal Area',
             'lcc'      : 'Lambert Conformal',
             'laea'     : 'Lambert Azimuthal Equal Area',
             'nplaea'   : 'North-Polar Lambert Azimuthal',
             'splaea'   : 'South-Polar Lambert Azimuthal',
             'eqdc'     : 'Equidistant Conic',
             'aeqd'     : 'Azimuthal Equidistant',
             'npaeqd'   : 'North-Polar Azimuthal Equidistant',
             'spaeqd'   : 'South-Polar Azimuthal Equidistant',
             'aea'      : 'Albers Equal Area',
             'stere'    : 'Stereographic',
             'npstere'  : 'North-Polar Stereographic',
             'spstere'  : 'South-Polar Stereographic',
             'cass'     : 'Cassini-Soldner',
             'poly'     : 'Polyconic',
             'ortho'    : 'Orthographic',
             'geos'     : 'Geostationary',
             'nsper'    : 'Near-Sided Perspective',
             'sinu'     : 'Sinusoidal',
             'moll'     : 'Mollweide',
             'hammer'   : 'Hammer',
             'robin'    : 'Robinson',
             'kav7'     : 'Kavrayskiy VII',
             'eck4'     : 'Eckert IV',
             'vandg'    : 'van der Grinten',
             'mbtfpq'   : 'McBryde-Thomas Flat-Polar Quartic',
             'gnom'     : 'Gnomonic',
             'rotpole'  : 'Rotated Pole',
             }
supported_projections = []
for _items in _projnames.items():
    supported_projections.append(" %-17s%-40s\n" % (_items))
supported_projections = ''.join(supported_projections)

_cylproj = ['cyl','merc','mill','gall','cea']
_pseudocyl = ['moll','robin','eck4','kav7','sinu','mbtfpq','vandg','hammer']
_dg2rad = math.radians(1.)
_rad2dg = math.degrees(1.)

# projection specific parameters.
projection_params = {'cyl'      : 'corners only (no width/height)',
             'merc'     : 'corners plus lat_ts (no width/height)',
             'tmerc'    : 'lon_0,lat_0,k_0',
             'omerc'    : 'lon_0,lat_0,lat_1,lat_2,lon_1,lon_2,no_rot,k_0',
             'mill'     : 'corners only (no width/height)',
             'gall'     : 'corners only (no width/height)',
             'cea'      : 'corners only plus lat_ts (no width/height)',
             'lcc'      : 'lon_0,lat_0,lat_1,lat_2,k_0',
             'laea'     : 'lon_0,lat_0',
             'nplaea'   : 'bounding_lat,lon_0,lat_0,no corners or width/height',
             'splaea'   : 'bounding_lat,lon_0,lat_0,no corners or width/height',
             'eqdc'     : 'lon_0,lat_0,lat_1,lat_2',
             'aeqd'     : 'lon_0,lat_0',
             'npaeqd'   : 'bounding_lat,lon_0,lat_0,no corners or width/height',
             'spaeqd'   : 'bounding_lat,lon_0,lat_0,no corners or width/height',
             'aea'      : 'lon_0,lat_0,lat_1',
             'stere'    : 'lon_0,lat_0,lat_ts,k_0',
             'npstere'  : 'bounding_lat,lon_0,lat_0,no corners or width/height',
             'spstere'  : 'bounding_lat,lon_0,lat_0,no corners or width/height',
             'cass'     : 'lon_0,lat_0',
             'poly'     : 'lon_0,lat_0',
             'ortho'    : 'lon_0,lat_0,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
             'geos'     : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
             'nsper'    : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
             'sinu'     : 'lon_0,lat_0,no corners or width/height',
             'moll'     : 'lon_0,lat_0,no corners or width/height',
             'hammer'   : 'lon_0,lat_0,no corners or width/height',
             'robin'    : 'lon_0,lat_0,no corners or width/height',
             'eck4'    : 'lon_0,lat_0,no corners or width/height',
             'kav7'    : 'lon_0,lat_0,no corners or width/height',
             'vandg'    : 'lon_0,lat_0,no corners or width/height',
             'mbtfpq'   : 'lon_0,lat_0,no corners or width/height',
             'gnom'     : 'lon_0,lat_0',
             'rotpole'  : 'lon_0,o_lat_p,o_lon_p,corner lat/lon or corner x,y (no width/height)'
             }

# create dictionary that maps epsg codes to Basemap kwargs.
epsgf = open(os.path.join(basemap_datadir, 'epsg'))
epsg_dict={}
for line in epsgf:
    if line.startswith("#"):
        continue
    l = line.split()
    code = l[0].strip("<>")
    parms = ' '.join(l[1:-1])
    _kw_args={}
    for s in l[1:-1]:
        try:
            k,v = s.split('=')
        except:
            pass
        k = k.strip("+")
        if k=='proj':
            if v == 'longlat': v = 'cyl'
            if v not in _projnames:
                continue
            k='projection'
        if k=='k':
            k='k_0'
        if k in ['projection','lat_1','lat_2','lon_0','lat_0',\
                 'a','b','k_0','lat_ts','ellps','datum']:
            if k not in ['projection','ellps','datum']:
                v = float(v)
            _kw_args[k]=v
    if 'projection' in _kw_args:
        if 'a' in _kw_args:
            if 'b' in _kw_args:
                _kw_args['rsphere']=(_kw_args['a'],_kw_args['b'])
                del _kw_args['b']
            else:
                _kw_args['rsphere']=_kw_args['a']
            del _kw_args['a']
        if 'datum' in _kw_args:
            if _kw_args['datum'] == 'NAD83':
                _kw_args['ellps'] = 'GRS80'
            elif _kw_args['datum'] == 'NAD27':
                _kw_args['ellps'] = 'clrk66'
            elif _kw_args['datum'] == 'WGS84':
                _kw_args['ellps'] = 'WGS84'
            del _kw_args['datum']
        # supported epsg projections.
        # omerc not supported yet, since we can't handle
        # alpha,gamma and lonc keywords.
        if _kw_args['projection'] != 'omerc':
            epsg_dict[code]=_kw_args
epsgf.close()

# The __init__ docstring is pulled out here because it is so long;
# Having it in the usual place makes it hard to get from the
# __init__ argument list to the code that uses the arguments.
_Basemap_init_doc = """

 Sets up a basemap with specified map projection.
 and creates the coastline data structures in map projection
 coordinates.

 Calling a Basemap class instance with the arguments lon, lat will
 convert lon/lat (in degrees) to x/y map projection coordinates
 (in meters). The inverse transformation is done if the optional keyword
 ``inverse`` is set to True.

 The desired projection is set with the projection keyword. Default is ``cyl``.
 Supported values for the projection keyword are:

 ==============   ====================================================
 Value            Description
 ==============   ====================================================
%(supported_projections)s
 ==============   ====================================================

 For most map projections, the map projection region can either be
 specified by setting these keywords:

 .. tabularcolumns:: |l|L|

 ==============   ====================================================
 Keyword          Description
 ==============   ====================================================
 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

 .. tabularcolumns:: |l|L|

 ==============   ====================================================
 Keyword          Description
 ==============   ====================================================
 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 ``sinu``, ``moll``, ``hammer``, ``npstere``, ``spstere``, ``nplaea``, ``splaea``,
 ``npaeqd``, ``spaeqd``, ``robin``, ``eck4``, ``kav7``, or ``mbtfpq``, 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``, ``mill``, ``cea``  and ``gall``),
 the default is to use
 llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other
 projections except ``ortho``, ``geos`` and ``nsper``, either the lat/lon values of the
 corners or width and height must be specified by the user.

 For ``ortho``, ``geos`` and ``nsper``, the lat/lon values of the corners may be specified,
 or the x/y values of the corners (llcrnrx,llcrnry,urcrnrx,urcrnry) in the
 coordinate system of the global projection (with x=0,y=0 at the center
 of the global projection).  If the corners are not specified,
 the entire globe is plotted.

 For ``rotpole``, the lat/lon values of the corners on the unrotated sphere
 may be provided as llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat, or the lat/lon
 values of the corners on the rotated sphere can be given as
 llcrnrx,llcrnry,urcrnrx,urcrnry.

 Other keyword arguments:

 .. tabularcolumns:: |l|L|

 ==============   ====================================================
 Keyword          Description
 ==============   ====================================================
 resolution       resolution of boundary database to use. Can be ``c``
                  (crude), ``l`` (low), ``i`` (intermediate), ``h``
                  (high), ``f`` (full) or None.
                  If None, no boundary data will be read in (and
                  class methods such as drawcoastlines will raise an
                  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,1 for resolution
                  ``c``, ``l``, ``i``, ``h``, ``f``.
 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 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).
 ellps            string describing ellipsoid ('GRS80' or 'WGS84',
                  for example). If both rsphere and ellps are given,
                  rsphere is ignored. Default None. See pyproj.pj_ellps
                  for allowed values.
 suppress_ticks   suppress automatic drawing of axis ticks and labels
                  in map projection coordinates.  Default True,
                  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 if
                  suppress_ticks=False.  suppress_ticks=False
                  is useful if you want to use your own custom tick
                  formatter, or  if you want to let matplotlib label
                  the axes in meters using map projection
                  coordinates.
 fix_aspect       fix aspect ratio of plot to match aspect ratio
                  of map projection region (default True).
 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``, and ``W``.
 celestial        use astronomical conventions for longitude (i.e.
                  negative longitudes to the east of 0). Default False.
                  Implies resolution=None.
 ax               set default axes instance
                  (default None - matplotlib.pyplot.gca() may be used
                  to get the current axes instance).
                  If you do not want matplotlib.pyplot 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 keywords are map projection parameters which all default to
 None.  Not all parameters are used by all projections, some are ignored.
 The module variable ``projection_params`` is a dictionary which
 lists which parameters apply to which projections.

 .. tabularcolumns:: |l|L|

 ================ ====================================================
 Keyword          Description
 ================ ====================================================
 lat_ts           latitude of true scale. Optional for stereographic,
                  cylindrical equal area and mercator projections.
                  default is lat_0 for stereographic projection.
                  default is 0 for mercator and cylindrical equal area
                  projections.
 lat_1            first standard parallel for lambert conformal,
                  albers equal area and equidistant conic.
                  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 and equidistant conic.
                  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.
 k_0              Scale factor at natural origin (used
                  by 'tmerc', 'omerc', 'stere' and 'lcc').
 no_rot           only used by oblique mercator.
                  If set to True, the map projection coordinates will
                  not be rotated to true North.  Default is False
                  (projection coordinates are automatically rotated).
 lat_0            central latitude (y-axis origin) - used by all
                  projections.
 lon_0            central meridian (x-axis origin) - used by all
                  projections.
 o_lat_p          latitude of rotated pole (only used by 'rotpole')
 o_lon_p          longitude of rotated pole (only used by 'rotpole')
 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.
 round            cut off pole-centered projection at boundinglat
                  (so plot is a circle instead of a square). Only
                  relevant for npstere,spstere,nplaea,splaea,npaeqd
                  or spaeqd projections. Default False.
 satellite_height height of satellite (in m) above equator -
                  only relevant for geostationary
                  and near-sided perspective (``geos`` or ``nsper``)
                  projections. Default 35,786 km.
 ================ ====================================================

 Useful instance variables:

 .. tabularcolumns:: |l|L|

 ================ ====================================================
 Variable Name    Description
 ================ ====================================================
 projection       map projection. Print the module variable
                  ``supported_projections`` to see a list of allowed
                  values.
 epsg             EPSG code defining projection (see
                  http://spatialreference.org for a list of
                  EPSG codes and their definitions).
 aspect           map aspect ratio
                  (size of y dimension / size of x dimension).
 llcrnrlon        longitude of lower left hand corner of the
                  selected map domain.
 llcrnrlat        latitude of lower left hand corner of the
                  selected map domain.
 urcrnrlon        longitude of upper right hand corner of the
                  selected map domain.
 urcrnrlat        latitude of upper right hand corner of the
                  selected map domain.
 llcrnrx          x value of lower left hand corner of the
                  selected map domain in map projection coordinates.
 llcrnry          y value of lower left hand corner of the
                  selected map domain in map projection coordinates.
 urcrnrx          x value of upper right hand corner of the
                  selected map domain in map projection coordinates.
 urcrnry          y value of upper right hand corner of the
                  selected map domain in map projection coordinates.
 rmajor           equatorial radius of ellipsoid used (in meters).
 rminor           polar radius 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.
 proj4string      the string describing the map projection that is
                  used by PROJ.4.
 ================ ====================================================

 **Converting from Geographic (lon/lat) to Map Projection (x/y) Coordinates**

 Calling a Basemap class instance with the arguments lon, lat will
 convert lon/lat (in degrees) to x/y 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``,
 ``cea``, ``gall`` or ``merc``)
 the inverse transformation will return longitudes between
 self.llcrnrlon and self.llcrnrlat.

 Input arguments lon, lat can be either scalar floats, sequences
 or numpy arrays.

 **Example Usage:**

 >>> from mpl_toolkits.basemap import Basemap
 >>> import numpy as np
 >>> import matplotlib.pyplot as plt
 >>> # read in topo data (on a regular lat/lon grid)
 >>> etopo = np.loadtxt('etopo20data.gz')
 >>> lons  = np.loadtxt('etopo20lons.gz')
 >>> lats  = np.loadtxt('etopo20lats.gz')
 >>> # create Basemap instance for Robinson projection.
 >>> m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
 >>> # compute map projection coordinates for lat/lon grid.
 >>> x, y = m(*np.meshgrid(lons,lats))
 >>> # make filled contour plot.
 >>> cs = m.contourf(x,y,etopo,30,cmap=plt.cm.jet)
 >>> m.drawcoastlines() # draw coastlines
 >>> m.drawmapboundary() # draw a line around the map region
 >>> m.drawparallels(np.arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels
 >>> m.drawmeridians(np.arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
 >>> plt.title('Robinson Projection') # add a title
 >>> plt.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 matplotlib.pyplot) is called "simpletest_oo.py".]
""" % locals()

# unsupported projection error message.
_unsupported_projection = ["'%s' is an unsupported projection.\n"]
_unsupported_projection.append("The supported projections are:\n")
_unsupported_projection.append(supported_projections)
_unsupported_projection = ''.join(_unsupported_projection)

def _validated_ll(param, name, minval, maxval):
    item = np.squeeze(param)
    param = float(param if item.shape else item)
    if param > maxval or param < minval:
        raise ValueError('%s must be between %f and %f degrees' %
                                           (name, minval, maxval))
    return param


def _validated_or_none(param, name, minval, maxval):
    if param is None:
        return None
    return _validated_ll(param, name, minval, maxval)


def _insert_validated(d, param, name, minval, maxval):
    if param is not None:
        d[name] = _validated_ll(param, name, minval, maxval)

def _transform(plotfunc):
    # shift data and longitudes to map projection region, then compute
    # transformation to map projection coordinates.
    @functools.wraps(plotfunc)
    def with_transform(self,x,y,data,*args,**kwargs):
        # input coordinates are latitude/longitude, not map projection coords.
        if kwargs.pop('latlon', latlon_default):
            # shift data to map projection region for
            # cylindrical and pseudo-cylindrical projections.
            if self.projection in _cylproj or self.projection in _pseudocyl:
                x, data = self.shiftdata(x, data,
                                         fix_wrap_around=plotfunc.__name__ not in ["scatter"])
            # convert lat/lon coords to map projection coords.
            x, y = self(x,y)
        return plotfunc(self,x,y,data,*args,**kwargs)
    return with_transform

def _transform1d(plotfunc):
    # shift data and longitudes to map projection region, then compute
    # transformation to map projection coordinates.
    @functools.wraps(plotfunc)
    def with_transform(self,x,y,*args,**kwargs):
        x = np.asarray(x)
        # input coordinates are latitude/longitude, not map projection coords.
        if kwargs.pop('latlon', latlon_default):
            # shift data to map projection region for
            # cylindrical and pseudo-cylindrical projections.
            if self.projection in _cylproj or self.projection in _pseudocyl:
                if x.ndim == 1:
                    x = self.shiftdata(x, fix_wrap_around=plotfunc.__name__ not in ["scatter"])
                elif x.ndim == 0:
                    if x > 180:
                        x = x - 360.
            # convert lat/lon coords to map projection coords.
            x, y = self(x,y)
        return plotfunc(self,x,y,*args,**kwargs)
    return with_transform

def _transformuv(plotfunc):
    # shift data and longitudes to map projection region, then compute
    # transformation to map projection coordinates. Works when call
    # signature has two data arrays instead of one.
    @functools.wraps(plotfunc)
    def with_transform(self,x,y,u,v,*args,**kwargs):
        # input coordinates are latitude/longitude, not map projection coords.
        if kwargs.pop('latlon', latlon_default):
            # shift data to map projection region for
            # cylindrical and pseudo-cylindrical projections.
            if self.projection in _cylproj or self.projection in _pseudocyl:
                x1, u = self.shiftdata(x, u)
                x, v = self.shiftdata(x, v)
            # convert lat/lon coords to map projection coords.
            x, y = self(x,y)
        return plotfunc(self,x,y,u,v,*args,**kwargs)
    return with_transform

class Basemap(object):

    def __init__(self, llcrnrlon=None, llcrnrlat=None,
                       urcrnrlon=None, urcrnrlat=None,
                       llcrnrx=None, llcrnry=None,
                       urcrnrx=None, urcrnry=None,
                       width=None, height=None,
                       projection='cyl', resolution='c',
                       area_thresh=None, rsphere=6370997.0,
                       ellps=None, lat_ts=None,
                       lat_1=None, lat_2=None,
                       lat_0=None, lon_0=None,
                       lon_1=None, lon_2=None,
                       o_lon_p=None, o_lat_p=None,
                       k_0=None,
                       no_rot=False,
                       suppress_ticks=True,
                       satellite_height=35786000,
                       boundinglat=None,
                       fix_aspect=True,
                       anchor='C',
                       celestial=False,
                       round=False,
                       epsg=None,
                       ax=None):
        # docstring is added after __init__ method definition

        # set epsg code if given, set to 4326 for projection='cyl':
        if epsg is not None:
            self.epsg = epsg
        elif projection == 'cyl':
            self.epsg = 4326
        # replace kwarg values with those implied by epsg code,
        # if given.
        if hasattr(self,'epsg'):
            if str(self.epsg) not in epsg_dict:
                raise ValueError('%s is not a supported EPSG code' %
                        self.epsg)
            epsg_params = epsg_dict[str(self.epsg)]
            for k in epsg_params:
                if k == 'projection':
                    projection = epsg_params[k]
                elif k == 'rsphere':
                    rsphere = epsg_params[k]
                elif k == 'ellps':
                    ellps = epsg_params[k]
                elif k == 'lat_1':
                    lat_1 = epsg_params[k]
                elif k == 'lat_2':
                    lat_2 = epsg_params[k]
                elif k == 'lon_0':
                    lon_0 = epsg_params[k]
                elif k == 'lat_0':
                    lat_0 = epsg_params[k]
                elif k == 'lat_ts':
                    lat_ts = epsg_params[k]
                elif k == 'k_0':
                    k_0 = epsg_params[k]

        # fix aspect to ratio to match aspect ratio of map projection
        # region
        self.fix_aspect = fix_aspect
        # where to put plot in figure (default is 'C' or center)
        self.anchor = anchor
        # geographic or celestial coords?
        self.celestial = celestial
        # map projection.
        self.projection = projection
        # bounding lat (for pole-centered plots)
        self.boundinglat = boundinglat
        # is a round pole-centered plot desired?
        self.round = round
        # full disk projection?
        self._fulldisk = False # default value

        # set up projection parameter dict.
        projparams = {}
        projparams['proj'] = projection
        # if ellps keyword specified, it over-rides rsphere.
        if ellps is not None:
            try:
                elldict = pyproj.pj_ellps[ellps]
            except KeyError:
                raise ValueError(
                'illegal ellps definition, allowed values are %s' %
                pyproj.pj_ellps.keys())
            projparams['a'] = elldict['a']
            if 'b' in elldict:
                projparams['b'] = elldict['b']
            else:
                projparams['b'] = projparams['a']*(1.0-(1.0/elldict['rf']))
        else:
            try:
                if rsphere[0] > rsphere[1]:
                    projparams['a'] = rsphere[0]
                    projparams['b'] = rsphere[1]
                else:
                    projparams['a'] = rsphere[1]
                    projparams['b'] = rsphere[0]
            except:
                if projection == 'tmerc':
                # use bR_a instead of R because of obscure bug
                # in proj4 for tmerc projection.
                    projparams['bR_a'] = rsphere
                else:
                    projparams['R'] = rsphere
        # set units to meters.
        projparams['units']='m'
        # check for sane values of lon_0, lat_0, lat_ts, lat_1, lat_2
        lat_0 = _validated_or_none(lat_0, 'lat_0', -90, 90)
        lat_1 = _validated_or_none(lat_1, 'lat_1', -90, 90)
        lat_2 = _validated_or_none(lat_2, 'lat_2', -90, 90)
        lat_ts = _validated_or_none(lat_ts, 'lat_ts', -90, 90)
        lon_0 = _validated_or_none(lon_0, 'lon_0', -360, 720)
        lon_1 = _validated_or_none(lon_1, 'lon_1', -360, 720)
        lon_2 = _validated_or_none(lon_2, 'lon_2', -360, 720)
        llcrnrlon = _validated_or_none(llcrnrlon, 'llcrnrlon', -360, 720)
        urcrnrlon = _validated_or_none(urcrnrlon, 'urcrnrlon', -360, 720)
        llcrnrlat = _validated_or_none(llcrnrlat, 'llcrnrlat', -90, 90)
        urcrnrlat = _validated_or_none(urcrnrlat, 'urcrnrlat', -90, 90)

        _insert_validated(projparams, lat_0, 'lat_0', -90, 90)
        _insert_validated(projparams, lat_1, 'lat_1', -90, 90)
        _insert_validated(projparams, lat_2, 'lat_2', -90, 90)
        _insert_validated(projparams, lat_ts, 'lat_ts', -90, 90)
        _insert_validated(projparams, lon_0, 'lon_0', -360, 720)
        _insert_validated(projparams, lon_1, 'lon_1', -360, 720)
        _insert_validated(projparams, lon_2, 'lon_2', -360, 720)
        if projection in ['geos','nsper']:
            projparams['h'] = satellite_height
        # check for sane values of projection corners.
        using_corners = (None not in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat])
        if using_corners:
            self.llcrnrlon = _validated_ll(llcrnrlon, 'llcrnrlon', -360, 720)
            self.urcrnrlon = _validated_ll(urcrnrlon, 'urcrnrlon', -360, 720)
            self.llcrnrlat = _validated_ll(llcrnrlat, 'llcrnrlat', -90, 90)
            self.urcrnrlat = _validated_ll(urcrnrlat, 'urcrnrlat', -90, 90)

        # for each of the supported projections,
        # compute lat/lon of domain corners
        # and set values in projparams dict as needed.

        if projection in ['lcc', 'eqdc', 'aea']:
            if projection == 'lcc' and k_0 is not None:
                projparams['k_0']=k_0
            # 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 %s basemap (lat_2 is optional)' % _projnames[projection])
            if lat_2 is None:
                projparams['lat_2'] = lat_1
            if not using_corners:
                using_cornersxy = (None not in [llcrnrx,llcrnry,urcrnrx,urcrnry])
                if using_cornersxy:
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecornersllur(llcrnrx,llcrnry,urcrnrx,urcrnry,**projparams)
                    self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                    self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
                else:
                    if width is None or height is None:
                        raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
                    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 k_0 is not None:
                projparams['k_0']=k_0
            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 not using_corners:
                if width is None or height is None:
                    raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
                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 in ['spstere', 'npstere',
                            'splaea', 'nplaea',
                            'spaeqd', 'npaeqd']:
            if (projection == 'splaea' and boundinglat >= 0) or\
               (projection == 'nplaea' and boundinglat <= 0):
                raise ValueError('boundinglat cannot extend into opposite hemisphere')
            if boundinglat is None or lon_0 is None:
                raise ValueError('must specify boundinglat and lon_0 for %s basemap' % _projnames[projection])
            if projection[0] == 's':
                sgn = -1
            else:
                sgn = 1
            rootproj = projection[2:]
            projparams['proj'] = rootproj
            if rootproj == 'stere':
                projparams['lat_ts'] = sgn * 90.
            projparams['lat_0'] = sgn * 90.
            self.llcrnrlon = lon_0 - sgn*45.
            self.urcrnrlon = lon_0 + sgn*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:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[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 not using_corners:
                if width is None or height is None:
                    raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
                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 in ['tmerc','gnom','cass','poly'] :
            if projection == 'tmerc' and k_0 is not None:
                projparams['k_0']=k_0
            if projection == 'gnom' and 'R' not in projparams:
                raise ValueError('gnomonic projection only works for perfect spheres - not ellipsoids')
            if lat_0 is None or lon_0 is None:
                raise ValueError('must specify lat_0 and lon_0 for Transverse Mercator, Gnomonic, Cassini-Soldnerr and Polyconic basemap')
            if not using_corners:
                if width is None or height is None:
                    raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
                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 'R' not in projparams:
                raise ValueError('orthographic projection only works for perfect spheres - not ellipsoids')
            if lat_0 is None or lon_0 is None:
                raise ValueError('must specify lat_0 and lon_0 for Orthographic basemap')
            if (lat_0 == 90 or lat_0 == -90) and\
               None in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
                # for ortho plot centered on pole, set boundinglat to equator.
                # (so meridian labels can be drawn in this special case).
                self.boundinglat = 0
                self.round = True
            if width is not None or height is not None:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
            if not using_corners:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self._fulldisk = True
            else:
                self._fulldisk = False
            self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
            self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
            # FIXME: won't work for points exactly on equator??
            if np.abs(lat_0) < 1.e-2: lat_0 = 1.e-2
            projparams['lat_0'] = lat_0
        elif projection == 'geos':
            if lat_0 is not None and lat_0 != 0:
                raise ValueError('lat_0 must be zero for Geostationary basemap')
            if lon_0 is None:
                raise ValueError('must specify lon_0 for Geostationary basemap')
            if width is not None or height is not None:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
            if not using_corners:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self._fulldisk = True
            else:
                self._fulldisk = False
            self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
            self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'nsper':
            if 'R' not in projparams:
                raise ValueError('near-sided perspective projection only works for perfect spheres - not ellipsoids')
            if lat_0 is None or lon_0 is None:
                raise ValueError('must specify lon_0 and lat_0 for near-sided perspective Basemap')
            if width is not None or height is not None:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
            if not using_corners:
                llcrnrlon = -180.
                llcrnrlat = -90.
                urcrnrlon = 180
                urcrnrlat = 90.
                self._fulldisk = True
            else:
                self._fulldisk = False
            self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
            self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection in _pseudocyl:
            if lon_0 is None:
                raise ValueError('must specify lon_0 for %s projection' % _projnames[self.projection])
            if width is not None or height is not None:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
            llcrnrlon = lon_0-180.
            llcrnrlat = -90.
            urcrnrlon = lon_0+180
            urcrnrlat = 90.
            self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
            self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection == 'omerc':
            if k_0 is not None:
                projparams['k_0']=k_0
            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
            projparams['lat_0'] = lat_0
            if no_rot:
                projparams['no_rot']=''
            #if not using_corners:
            #    raise ValueError, 'cannot specify map region with width and height keywords for this projection, please specify lat/lon values of corners'
            if not using_corners:
                if width is None or height is None:
                    raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
                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 not using_corners:
                if width is None or height is None:
                    self._fulldisk = True
                    llcrnrlon = -180.
                    llcrnrlat = -90.
                    urcrnrlon = 180
                    urcrnrlat = 90.
                else:
                    self._fulldisk = False
                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')
                if not self._fulldisk:
                    llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        elif projection in _cylproj:
            if projection == 'merc' or projection == 'cea':
                if lat_ts is None:
                    lat_ts = 0.
                    projparams['lat_ts'] = lat_ts
            if not using_corners:
                llcrnrlat = -90.
                urcrnrlat = 90.
                if lon_0 is not None:
                    llcrnrlon = lon_0-180.
                    urcrnrlon = lon_0+180.
                else:
                    llcrnrlon = -180.
                    urcrnrlon = 180
                if projection == 'merc':
                    # clip plot region to be within -89.99S to 89.99N
                    # (mercator is singular at poles)
                    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:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
            if lon_0 is not None:
                projparams['lon_0'] = lon_0
            else:
                projparams['lon_0']=0.5*(llcrnrlon+urcrnrlon)
        elif projection == 'rotpole':
            if lon_0 is None or o_lon_p is None or o_lat_p is None:
                raise ValueError('must specify lon_0,o_lat_p,o_lon_p for rotated pole Basemap')
            if width is not None or height is not None:
                sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
            projparams['lon_0']=lon_0
            projparams['o_lon_p']=o_lon_p
            projparams['o_lat_p']=o_lat_p
            projparams['o_proj']='longlat'
            projparams['proj']='ob_tran'
            if not using_corners and None in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
                raise ValueError('must specify lat/lon values of corners in degrees')
            if None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
                p = pyproj.Proj(projparams)
                llcrnrx = _dg2rad*llcrnrx; llcrnry = _dg2rad*llcrnry
                urcrnrx = _dg2rad*urcrnrx; urcrnry = _dg2rad*urcrnry
                llcrnrlon, llcrnrlat = p(llcrnrx,llcrnry,inverse=True)
                urcrnrlon, urcrnrlat = p(urcrnrx,urcrnry,inverse=True)
                self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
                self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
        else:
            raise ValueError(_unsupported_projection % projection)

        # initialize proj4
        proj = Proj(projparams,self.llcrnrlon,self.llcrnrlat,self.urcrnrlon,self.urcrnrlat)

        # make sure axis ticks are suppressed.
        self.noticks = suppress_ticks
        # map boundary not yet drawn.
        self._mapboundarydrawn = False

        # 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]
        # these only exist for geostationary projection.
        if hasattr(proj,'_width'):
            self.__dict__['_width'] = proj.__dict__['_width']
        if hasattr(proj,'_height'):
            self.__dict__['_height'] = proj.__dict__['_height']
        # spatial reference string (useful for georeferencing output
        # images with gdal_translate).
        if hasattr(self,'_proj4'):
            #self.srs = proj._proj4.srs
            self.srs = proj._proj4.pjinitstring
        else:
            pjargs = []
            for key,value in self.projparams.items():
                # 'cyl' projection translates to 'eqc' in PROJ.4
                if projection == 'cyl' and key == 'proj':
                    value = 'eqc'
                # ignore x_0 and y_0 settings for 'cyl' projection
                # (they are not consistent with what PROJ.4 uses)
                elif projection == 'cyl' and key in ['x_0','y_0']:
                    continue
                pjargs.append('+'+key+"="+str(value)+' ')
            self.srs = ''.join(pjargs)
        self.proj4string = self.srs
        # 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)
        if projection in ['geos','ortho','nsper'] and \
           None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
            self.llcrnrx = llcrnrx+0.5*proj.xmax
            self.llcrnry = llcrnry+0.5*proj.ymax
            self.urcrnrx = urcrnrx+0.5*proj.xmax
            self.urcrnry = urcrnry+0.5*proj.ymax
            self._fulldisk = False
        else:
            self.llcrnrx = proj.llcrnrx
            self.llcrnry = proj.llcrnry
            self.urcrnrx = proj.urcrnrx
            self.urcrnry = proj.urcrnry

        if self.projection == 'rotpole':
            lon0,lat0 = self(0.5*(self.llcrnrx + self.urcrnrx),\
                             0.5*(self.llcrnry + self.urcrnry),\
                             inverse=True)
            self.projparams['lat_0']=lat0

        # if ax == None, pyplot.gca may be used.
        self.ax = ax
        self.lsmask = None
        # This will record hashs of Axes instances.
        self._initialized_axes = set()

        # set defaults for area_thresh.
        self.resolution = resolution
        # celestial=True implies resolution=None (no coastlines).
        if self.celestial:
            self.resolution=None
        if area_thresh is None and self.resolution is not 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.
            elif resolution == 'f':
                area_thresh = 1.
            else:
                raise ValueError("boundary resolution must be one of 'c','l','i','h' or 'f'")
        self.area_thresh = area_thresh
        # define map boundary polygon (in lat/lon coordinates)
        blons, blats, self._boundarypolyll, self._boundarypolyxy = self._getmapboundary()
        self.boundarylats = blats
        self.boundarylons = blons
        # set min/max lats for projection domain.
        if self.projection in _cylproj:
            self.latmin = self.llcrnrlat
            self.latmax = self.urcrnrlat
            self.lonmin = self.llcrnrlon
            self.lonmax = self.urcrnrlon
        elif self.projection in ['ortho','geos','nsper'] + _pseudocyl:
            self.latmin = -90.
            self.latmax = 90.
            self.lonmin = self.llcrnrlon
            self.lonmax = self.urcrnrlon
        else:
            lons, lats = self.makegrid(1001,1001)
            lats = ma.masked_where(lats > 1.e20,lats)
            lons = ma.masked_where(lons > 1.e20,lons)
            self.latmin = lats.min()
            self.latmax = lats.max()
            self.lonmin = lons.min()
            self.lonmax = lons.max()
            NPole = _geoslib.Point(self(0.,90.))
            SPole = _geoslib.Point(self(0.,-90.))
            if lat_0 is None:
                lon_0, lat_0 =\
                self(0.5*(self.xmin+self.xmax),
                     0.5*(self.ymin+self.ymax),inverse=True)
            Dateline = _geoslib.Point(self(180.,lat_0))
            Greenwich = _geoslib.Point(self(0.,lat_0))
            hasNP = NPole.within(self._boundarypolyxy)
            hasSP = SPole.within(self._boundarypolyxy)
            hasPole = hasNP or hasSP
            hasDateline = Dateline.within(self._boundarypolyxy)
            hasGreenwich = Greenwich.within(self._boundarypolyxy)
            # projection crosses dateline (and not Greenwich or pole).
            if not hasPole and hasDateline and not hasGreenwich:
                if self.lonmin < 0 and self.lonmax > 0.:
                    lons = np.where(lons < 0, lons+360, lons)
                    self.lonmin = lons.min()
                    self.lonmax = lons.max()
        # read in coastline polygons, only keeping those that
        # intersect map boundary polygon.
        if self.resolution is not None:
            self.coastsegs, self.coastpolygontypes =\
            self._readboundarydata('gshhs',as_polygons=True)
            # reformat for use in matplotlib.patches.Polygon.
            self.coastpolygons = []
            for seg in self.coastsegs:
                x, y = list(zip(*seg))
                self.coastpolygons.append((x,y))
            # replace coastsegs with line segments (instead of polygons)
            self.coastsegs, types =\
            self._readboundarydata('gshhs',as_polygons=False)
            self.coastsegs = [sg for sg in self.coastsegs if len(sg) > 0]
        # create geos Polygon structures for land areas.
        # currently only used in is_land method.
        self.landpolygons=[]
        self.lakepolygons=[]
        if self.resolution is not None and len(self.coastpolygons) > 0:
            #self.islandinlakepolygons=[]
            #self.lakeinislandinlakepolygons=[]
            x, y = list(zip(*self.coastpolygons))
            for x,y,typ in zip(x,y,self.coastpolygontypes):
                b = np.asarray([x,y]).T
                if typ == 1: self.landpolygons.append(_geoslib.Polygon(b))
                if typ == 2: self.lakepolygons.append(_geoslib.Polygon(b))
                #if typ == 3: self.islandinlakepolygons.append(_geoslib.Polygon(b))
                #if typ == 4: self.lakeinislandinlakepolygons.append(_geoslib.Polygon(b))

    # set __init__'s docstring
    __init__.__doc__ = _Basemap_init_doc

    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 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``,
        ``cea``, ``mill``, ``gall`` or ``merc``)
        the inverse transformation will return longitudes between
        self.llcrnrlon and self.llcrnrlat.

        Input arguments lon, lat can be either scalar floats,
        sequences, or numpy arrays.
        """
        if self.celestial:
            # don't assume center of map is at greenwich
            # (only relevant for cyl or pseudo-cyl projections)
            if self.projection in _pseudocyl or self.projection in _cylproj:
                lon_0=self.projparams['lon_0']
            else:
                lon_0 = 0.
        if self.celestial and not inverse:
            try:
                x = 2.*lon_0-x
            except TypeError:
                x = [2*lon_0-xx for xx in x]
        if self.projection == 'rotpole' and inverse:
            try:
                x = _dg2rad*x
            except TypeError:
                x = [_dg2rad*xx for xx in x]
            try:
                y = _dg2rad*y
            except TypeError:
                y = [_dg2rad*yy for yy in y]
        xout,yout = self.projtran(x,y,inverse=inverse)
        if self.celestial and inverse:
            try:
                xout = -2.*lon_0-xout
            except:
                xout = [-2.*lon_0-xx for xx in xout]
        if self.projection == 'rotpole' and not inverse:
            try:
                xout = _rad2dg*xout
                xout = np.where(xout < 0., xout+360, xout)
            except TypeError:
                xout = [_rad2dg*xx for xx in xout]
                xout = [xx+360. if xx < 0 else xx for xx in xout]
            try:
                yout = _rad2dg*yout
            except TypeError:
                yout = [_rad2dg*yy for yy in yout]
        return xout,yout

    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 _readboundarydata(self,name,as_polygons=False):
        """
        read boundary data, clip to map projection region.
        """

        # only gshhs coastlines can be polygons.
        if name != 'gshhs': as_polygons=False
        try:
            bdatfile = open(os.path.join(basemap_datadir,name+'_'+self.resolution+'.dat'),'rb')
            bdatmetafile = open(os.path.join(basemap_datadir,name+'meta_'+self.resolution+'.dat'),'r')
        except:
            raise IOError(
                "Unable to open boundary dataset file. Only the 'crude', 'low' "
                "and 'intermediate' resolution datasets are installed by default. "
                "If you are requesting a 'high' or 'full' resolution dataset, "
                "you need to install the `basemap-data-hires` package")
        polygons = []
        polygon_types = []
        # coastlines are polygons, other boundaries are line segments.
        if name == 'gshhs':
            Shape = _geoslib.Polygon
        else:
            Shape = _geoslib.LineString
        # see if map projection region polygon contains a pole.
        NPole = _geoslib.Point(self(0.,90.))
        SPole = _geoslib.Point(self(0.,-90.))
        boundarypolyxy = self._boundarypolyxy
        boundarypolyll = self._boundarypolyll
        hasNP = NPole.within(boundarypolyxy)
        hasSP = SPole.within(boundarypolyxy)
        containsPole = hasNP or hasSP
        # these projections cannot cross pole.
        if containsPole and\
            self.projection in _cylproj + _pseudocyl + ['geos']:
            raise ValueError('%s projection cannot cross pole'%(self.projection))
        # make sure some projections have has containsPole=True
        # we will compute the intersections in stereographic
        # coordinates, then transform back. This is
        # because these projections are only defined on a hemisphere, and
        # some boundary features (like Eurasia) would be undefined otherwise.
        tostere =\
        ['omerc','ortho','gnom','nsper','nplaea','npaeqd','splaea','spaeqd']
        if self.projection in tostere and name == 'gshhs':
            containsPole = True
            lon_0=self.projparams['lon_0']
            lat_0=self.projparams['lat_0']
            re = self.projparams['R']
            # center of stereographic projection restricted to be
            # nearest one of 6 points on the sphere (every 90 deg lat/lon).
            lon0 = 90.*(np.around(lon_0/90.))
            lat0 = 90.*(np.around(lat_0/90.))
            if np.abs(int(lat0)) == 90: lon0=0.
            maptran = pyproj.Proj(proj='stere',lon_0=lon0,lat_0=lat0,R=re)
            # boundary polygon for ortho/gnom/nsper projection
            # in stereographic coordinates.
            b = self._boundarypolyll.boundary
            blons = b[:,0]; blats = b[:,1]
            b[:,0], b[:,1] = maptran(blons, blats)
            boundarypolyxy = _geoslib.Polygon(b)
        for line in bdatmetafile:
            linesplit = line.split()
            area = float(linesplit[1])
            south = float(linesplit[3])
            north = float(linesplit[4])
            crossdatelineE=False; crossdatelineW=False
            if name == 'gshhs':
                id = linesplit[7]
                if id.endswith('E'):
                    crossdatelineE = True
                elif id.endswith('W'):
                    crossdatelineW = True
            # make sure south/north limits of dateline crossing polygons
            # (Eurasia) are the same, since they will be merged into one.
            # (this avoids having one filtered out and not the other).
            if crossdatelineE:
                south_save=south
                north_save=north
            if crossdatelineW:
                south=south_save
                north=north_save
            if area < 0.: area = 1.e30
            useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh
            if useit:
                typ = int(linesplit[0])
                npts = int(linesplit[2])
                offsetbytes = int(linesplit[5])
                bytecount = int(linesplit[6])
                bdatfile.seek(offsetbytes,0)
                # read in binary string convert into an npts by 2
                # numpy array (first column is lons, second is lats).
                polystring = bdatfile.read(bytecount)
                # binary data is little endian.
                b = np.array(np.frombuffer(polystring,dtype='<f4'),'f8').reshape(npts, 2)
                b2 = b.copy()
                # merge polygons that cross dateline.
                poly = Shape(b)
                # hack to try to avoid having Antartica filled polygon
                # covering entire map (if skipAnart = False, this happens
                # for ortho lon_0=-120, lat_0=60, for example).
                skipAntart = self.projection in tostere and south < -89 and \
                 not hasSP
                if crossdatelineE and not skipAntart:
                    if not poly.is_valid(): poly=poly.fix()
                    polyE = poly
                    continue
                elif crossdatelineW and not skipAntart:
                    if not poly.is_valid(): poly=poly.fix()
                    b = poly.boundary
                    b[:,0] = b[:,0]+360.
                    poly = Shape(b)
                    poly = poly.union(polyE)
                    if not poly.is_valid(): poly=poly.fix()
                    b = poly.boundary
                    b2 = b.copy()
                    # fix Antartica.
                    if name == 'gshhs' and south < -89:
                        b = b[3:,:]
                        b2 = b.copy()
                        poly = Shape(b)
                # if map boundary polygon is a valid one in lat/lon
                # coordinates (i.e. it does not contain either pole),
                # the intersections of the boundary geometries
                # and the map projection region can be computed before
                # transforming the boundary geometry to map projection
                # coordinates (this saves time, especially for small map
                # regions and high-resolution boundary geometries).
                if not containsPole:
                    # close Antarctica.
                    if name == 'gshhs' and south < -89:
                        lons2 = b[:,0]
                        lats = b[:,1]
                        lons1 = lons2 - 360.
                        lons3 = lons2 + 360.
                        lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
                        lats = lats.tolist()+lats.tolist()+lats.tolist()
                        lonstart,latstart = lons[0], lats[0]
                        lonend,latend = lons[-1], lats[-1]
                        lons.insert(0,lonstart)
                        lats.insert(0,-90.)
                        lons.append(lonend)
                        lats.append(-90.)
                        b = np.empty((len(lons),2),np.float64)
                        b[:,0] = lons; b[:,1] = lats
                        poly = Shape(b)
                        if not poly.is_valid(): poly=poly.fix()
                        # if polygon instersects map projection
                        # region, process it.
                        if poly.intersects(boundarypolyll):
                            if name != 'gshhs' or as_polygons:
                                geoms = poly.intersection(boundarypolyll)
                            else:
                                # convert polygons to line segments
                                poly = _geoslib.LineString(poly.boundary)
                                geoms = poly.intersection(boundarypolyll)
                            # iterate over geometries in intersection.
                            for psub in geoms:
                                b = psub.boundary
                                blons = b[:,0]; blats = b[:,1]
                                bx, by = self(blons, blats)
                                polygons.append(list(zip(bx,by)))
                                polygon_types.append(typ)
                    else:
                        # create duplicate polygons shifted by -360 and +360
                        # (so as to properly treat polygons that cross
                        # Greenwich meridian).
                        b2[:,0] = b[:,0]-360
                        poly1 = Shape(b2)
                        b2[:,0] = b[:,0]+360
                        poly2 = Shape(b2)
                        polys = [poly1,poly,poly2]
                        for poly in polys:
                            # try to fix "non-noded intersection" errors.
                            if not poly.is_valid(): poly=poly.fix()
                            # if polygon instersects map projection
                            # region, process it.
                            if poly.intersects(boundarypolyll):
                                if name != 'gshhs' or as_polygons:
                                    geoms = poly.intersection(boundarypolyll)
                                else:
                                    # convert polygons to line segments
                                    # note: use fix method here or Eurasia
                                    # line segments sometimes disappear.
                                    poly = _geoslib.LineString(poly.fix().boundary)
                                    geoms = poly.intersection(boundarypolyll)
                                # iterate over geometries in intersection.
                                for psub in geoms:
                                    b = psub.boundary
                                    blons = b[:,0]; blats = b[:,1]
                                    # transformation from lat/lon to
                                    # map projection coordinates.
                                    bx, by = self(blons, blats)
                                    if not as_polygons or len(bx) > 4:
                                        polygons.append(list(zip(bx,by)))
                                        polygon_types.append(typ)
                # if map boundary polygon is not valid in lat/lon
                # coordinates, compute intersection between map
                # projection region and boundary geometries in map
                # projection coordinates.
                else:
                    # transform coordinates from lat/lon
                    # to map projection coordinates.
                    # special case for ortho/gnom/nsper, compute coastline polygon
                    # vertices in stereographic coords.
                    if name == 'gshhs' and as_polygons and self.projection in tostere:
                        b[:,0], b[:,1] = maptran(b[:,0], b[:,1])
                    else:
                        b[:,0], b[:,1] = self(b[:,0], b[:,1])
                    goodmask = np.logical_and(b[:,0]<1.e20,b[:,1]<1.e20)
                    # if less than two points are valid in
                    # map proj coords, skip this geometry.
                    if np.sum(goodmask) <= 1: continue
                    if name != 'gshhs' or (name == 'gshhs' and not as_polygons):
                        # if not a polygon,
                        # just remove parts of geometry that are undefined
                        # in this map projection.
                        bx = np.compress(goodmask, b[:,0])
                        by = np.compress(goodmask, b[:,1])
                        # split coastline segments that jump across entire plot.
                        xd = (bx[1:]-bx[0:-1])**2
                        yd = (by[1:]-by[0:-1])**2
                        dist = np.sqrt(xd+yd)
                        split = dist > 0.1*(self.xmax-self.xmin)
                        if np.sum(split) and self.projection not in _cylproj:
                            ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist()
                            iprev = 0
                            ind.append(len(xd))
                            for i in ind:
                                # don't add empty lists.
                                if len(list(range(iprev,i))):
                                    polygons.append(list(zip(bx[iprev:i],by[iprev:i])))
                                iprev = i
                        else:
                            polygons.append(list(zip(bx,by)))
                        polygon_types.append(typ)
                        continue
                    # create a GEOS geometry object.
                    if name == 'gshhs' and not as_polygons:
                        # convert polygons to line segments
                        poly = _geoslib.LineString(poly.boundary)
                    else:
                        # this is a workaround to avoid
                        # GEOS_ERROR: CGAlgorithmsDD::orientationIndex encountered NaN/Inf numbers
                        b[np.isposinf(b)] = 1e20
                        b[np.isneginf(b)] = -1e20
                        poly = Shape(b)
                    # this is a workaround to avoid
                    # "GEOS_ERROR: TopologyException:
                    # found non-noded intersection between ..."
                    if not poly.is_valid(): poly=poly.fix()
                    # if geometry instersects map projection
                    # region, and doesn't have any invalid points, process it.
                    if goodmask.all() and poly.intersects(boundarypolyxy):
                        # if geometry intersection calculation fails,
                        # just move on.
                        try:
                            geoms = poly.intersection(boundarypolyxy)
                        except:
                            continue
                        # iterate over geometries in intersection.
                        for psub in geoms:
                            b = psub.boundary
                            # if projection in ['ortho','gnom','nsper'],
                            # transform polygon from stereographic
                            # to ortho/gnom/nsper coordinates.
                            if self.projection in tostere:
                                # if coastline polygon covers more than 99%
                                # of map region for fulldisk projection,
                                # it's probably bogus, so skip it.
                                #areafrac = psub.area()/boundarypolyxy.area()
                                #if self.projection == ['ortho','nsper']:
                                #    if name == 'gshhs' and\
                                #       self._fulldisk and\
                                #       areafrac > 0.99: continue
                                # inverse transform from stereographic
                                # to lat/lon.
                                b[:,0], b[:,1] = maptran(b[:,0], b[:,1], inverse=True)
                                # orthographic/gnomonic/nsper.
                                b[:,0], b[:,1]= self(b[:,0], b[:,1])
                            if not as_polygons or len(b) > 4:
                                polygons.append(list(zip(b[:,0],b[:,1])))
                                polygon_types.append(typ)
        bdatfile.close()
        bdatmetafile.close()
        return polygons, polygon_types

    def _getmapboundary(self):
        """
        create map boundary polygon (in lat/lon and x/y coordinates)
        """
        nx = 100; ny = 100
        maptran = self
        if self.projection in ['ortho','geos','nsper']:
            # circular region.
            thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
            rminor = self._height
            rmajor = self._width
            x = rmajor*np.cos(thetas) + rmajor
            y = rminor*np.sin(thetas) + rminor
            b = np.empty((len(x),2),np.float64)
            b[:,0]=x; b[:,1]=y
            boundaryxy = _geoslib.Polygon(b)
            # compute proj instance for full disk, if necessary.
            if not self._fulldisk:
                projparms = self.projparams.copy()
                del projparms['x_0']
                del projparms['y_0']
                if self.projection == 'ortho':
                    llcrnrx = -self.rmajor
                    llcrnry = -self.rmajor
                    urcrnrx = -llcrnrx
                    urcrnry = -llcrnry
                else:
                    llcrnrx = -self._width
                    llcrnry = -self._height
                    urcrnrx = -llcrnrx
                    urcrnry = -llcrnry
                projparms['x_0']=-llcrnrx
                projparms['y_0']=-llcrnry
                maptran = pyproj.Proj(projparms)
        elif self.projection == 'aeqd' and self._fulldisk:
            # circular region.
            thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
            rminor = self._height
            rmajor = self._width
            x = rmajor*np.cos(thetas) + rmajor
            y = rminor*np.sin(thetas) + rminor
            b = np.empty((len(x),2),np.float64)
            b[:,0]=x; b[:,1]=y
            boundaryxy = _geoslib.Polygon(b)
        elif self.projection in _pseudocyl:
            nx = 10*nx; ny = 10*ny
            # quasi-elliptical region.
            lon_0 = self.projparams['lon_0']
            # left side
            lats1 = np.linspace(-89.9999,89.9999,ny).tolist()
            lons1 = len(lats1)*[lon_0-179.9]
            # top.
            lons2 = np.linspace(lon_0-179.9,lon_0+179.9,nx).tolist()
            lats2 = len(lons2)*[89.9999]
            # right side
            lats3 = np.linspace(89.9999,-89.9999,ny).tolist()
            lons3 = len(lats3)*[lon_0+179.9]
            # bottom.
            lons4 = np.linspace(lon_0+179.9,lon_0-179.9,nx).tolist()
            lats4 = len(lons4)*[-89.9999]
            lons = np.array(lons1+lons2+lons3+lons4,np.float64)
            lats = np.array(lats1+lats2+lats3+lats4,np.float64)
            x, y = maptran(lons,lats)
            b = np.empty((len(x),2),np.float64)
            b[:,0]=x; b[:,1]=y
            boundaryxy = _geoslib.Polygon(b)
        else: # all other projections are rectangular.
            nx = 100*nx; ny = 100*ny
            # left side (x = xmin, ymin <= y <= ymax)
            yy = np.linspace(self.ymin, self.ymax, ny)[:-1]
            x = len(yy)*[self.xmin]; y = yy.tolist()
            # top (y = ymax, xmin <= x <= xmax)
            xx = np.linspace(self.xmin, self.xmax, nx)[:-1]
            x = x + xx.tolist()
            y = y + len(xx)*[self.ymax]
            # right side (x = xmax, ymin <= y <= ymax)
            yy = np.linspace(self.ymax, self.ymin, ny)[:-1]
            x = x + len(yy)*[self.xmax]; y = y + yy.tolist()
            # bottom (y = ymin, xmin <= x <= xmax)
            xx = np.linspace(self.xmax, self.xmin, nx)[:-1]
            x = x + xx.tolist()
            y = y + len(xx)*[self.ymin]
            x = np.array(x,np.float64)
            y = np.array(y,np.float64)
            b = np.empty((4,2),np.float64)
            b[:,0]=[self.xmin,self.xmin,self.xmax,self.xmax]
            b[:,1]=[self.ymin,self.ymax,self.ymax,self.ymin]
            boundaryxy = _geoslib.Polygon(b)
        if self.projection in _cylproj:
            # make sure map boundary doesn't quite include pole.
            if self.urcrnrlat > 89.9999:
                urcrnrlat = 89.9999
            else:
                urcrnrlat = self.urcrnrlat
            if self.llcrnrlat < -89.9999:
                llcrnrlat = -89.9999
            else:
                llcrnrlat = self.llcrnrlat
            lons = [self.llcrnrlon, self.llcrnrlon, self.urcrnrlon, self.urcrnrlon]
            lats = [llcrnrlat, urcrnrlat, urcrnrlat, llcrnrlat]
            self.boundarylonmin = min(lons)
            self.boundarylonmax = max(lons)
            x, y = self(lons, lats)
            b = np.empty((len(x),2),np.float64)
            b[:,0]=x; b[:,1]=y
            boundaryxy = _geoslib.Polygon(b)
        else:
            if self.projection not in _pseudocyl:
                lons, lats = maptran(x,y,inverse=True)
                # fix lons so there are no jumps.
                n = 1
                lonprev = lons[0]
                for lon,lat in zip(lons[1:],lats[1:]):
                    if np.abs(lon-lonprev) > 90.:
                        if lonprev < 0:
                            lon = lon - 360.
                        else:
                            lon = lon + 360
                        lons[n] = lon
                    lonprev = lon
                    n = n + 1
                self.boundarylonmin = lons.min()
                self.boundarylonmax = lons.max()
                # for circular full disk projections where boundary is
                # a latitude circle, set boundarylonmax and boundarylonmin
                # to cover entire world (so parallels will be drawn).
                if self._fulldisk and \
                   np.abs(self.boundarylonmax-self.boundarylonmin) < 1.:
                   self.boundarylonmin = -180.
                   self.boundarylonmax = 180.
        b = np.empty((len(lons),2),np.float64)
        b[:,0] = lons; b[:,1] = lats
        boundaryll = _geoslib.Polygon(b)
        return lons, lats, boundaryll, boundaryxy

    def drawmapboundary(self,color='k',linewidth=1.0,fill_color=None,\
                        zorder=None,ax=None):
        """
        draw boundary around map projection region, optionally
        filling interior of region.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        linewidth        line width for boundary (default 1.)
        color            color of boundary line (default black)
        fill_color       fill the map region background with this
                         color (default is to fill with axis
                         background color). If set to the string
                         'none', no filling is done.
        zorder           sets the zorder for filling map background
                         (default 0).
        ax               axes instance to use
                         (default None, use default axes instance).
        ==============   ====================================================

        returns matplotlib.collections.PatchCollection representing map boundary.
        """
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # if no fill_color given, use axes background color.
        # if fill_color is string 'none', really don't fill.
        if fill_color is None:
            mpl_version = tuple(map(int, mpl.__version__.split(".")[:2]))
            if mpl_version >= (2, 0):
                fill_color = ax.get_facecolor()
            else:
                fill_color = ax.get_axis_bgcolor()
        elif fill_color == 'none' or fill_color == 'None':
            fill_color = None
        limb = None
        if self.projection in ['ortho','geos','nsper'] or (self.projection=='aeqd' and\
           self._fulldisk):
            limb = Ellipse((self._width,self._height),2.*self._width,2.*self._height)
        if self.projection in ['ortho','geos','nsper','aeqd'] and self._fulldisk:
            # elliptical region.
            ax.set_frame_on(False)
        elif self.projection in _pseudocyl:  # elliptical region.
            ax.set_frame_on(False)
            nx = 100; ny = 100
            if self.projection == 'vandg':
                nx = 10*nx; ny = 10*ny
            # quasi-elliptical region.
            lon_0 = self.projparams['lon_0']
            # left side
            lats1 = np.linspace(-89.9999,89.99999,ny).tolist()
            lons1 = len(lats1)*[lon_0-179.9]
            # top.
            lons2 = np.linspace(lon_0-179.9999,lon_0+179.9999,nx).tolist()
            lats2 = len(lons2)*[89.9999]
            # right side
            lats3 = np.linspace(89.9999,-89.9999,ny).tolist()
            lons3 = len(lats3)*[lon_0+179.9999]
            # bottom.
            lons4 = np.linspace(lon_0+179.9999,lon_0-179.9999,nx).tolist()
            lats4 = len(lons4)*[-89.9999]
            lons = np.array(lons1+lons2+lons3+lons4,np.float64)
            lats = np.array(lats1+lats2+lats3+lats4,np.float64)
            x, y = self(lons,lats)
            xy = list(zip(x,y))
            limb = Polygon(xy)
        elif self.round:
            ax.set_frame_on(False)
            limb = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
                    radius=0.5*(self.xmax-self.xmin),fc='none')
        else: # all other projections are rectangular.
            ax.set_frame_on(True)
            for spine in ax.spines.values():
                spine.set_linewidth(linewidth)
                spine.set_edgecolor(color)
                if zorder is not None:
                    spine.set_zorder(zorder)
            if self.projection not in ['geos','ortho','nsper']:
                limb = ax.patch

        if limb is not None:
            if limb is not ax.patch:
                ax.add_patch(limb)
            self._mapboundarydrawn = limb
            if fill_color is None:
                limb.set_fill(False)
            else:
                limb.set_facecolor(fill_color)
                limb.set_zorder(0)
            limb.set_edgecolor(color)
            limb.set_linewidth(linewidth)
            if zorder is not None:
                limb.set_zorder(zorder)
            limb.set_clip_on(True)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        return limb

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

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        color            color to fill continents (default gray).
        lake_color       color to fill inland lakes (default axes background).
        ax               axes instance (overrides default axes instance).
        zorder           sets the zorder for the continent polygons (if not
                         specified, uses default zorder for a Polygon patch).
                         Set to zero if you want to paint over the filled
                         continents).
        alpha            sets alpha transparency for continent polygons
        ==============   ====================================================

        After filling continents, lakes are re-filled with
        axis background color.

        returns a list of matplotlib.patches.Polygon objects.
        """
        if self.resolution is None:
            raise AttributeError('there are no boundary datasets associated with this Basemap instance')
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # get axis background color.
        mpl_version = tuple(map(int, mpl.__version__.split(".")[:2]))
        if mpl_version >= (2, 0):
            axisbgc = ax.get_facecolor()
        else:
            axisbgc = ax.get_axis_bgcolor()
        npoly = 0
        polys = []
        for x,y in self.coastpolygons:
            xa = np.array(x,np.float32)
            ya = np.array(y,np.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).
        # ** turn this off for now since it prevents continents that
        # fill the whole map from being filled **
            #delx = 10; dely = 10
            #if self.projection in ['cyl']:
            #    delx = 0.1
            #    dely = 0.1
            #test1 = np.fabs(xa-self.urcrnrx) < delx
            #test2 = np.fabs(xa-self.llcrnrx) < delx
            #test3 = np.fabs(ya-self.urcrnry) < dely
            #test4 = np.fabs(ya-self.llcrnry) < dely
            #hasp1 = np.sum(test1*test3)
            #hasp2 = np.sum(test2*test3)
            #hasp4 = np.sum(test2*test4)
            #hasp3 = np.sum(test1*test4)
            #if not hasp1 or not hasp2 or not hasp3 or not hasp4:
            if 1:
                xy = list(zip(xa.tolist(),ya.tolist()))
                if self.coastpolygontypes[npoly] not in [2,4]:
                    poly = Polygon(xy,facecolor=color,edgecolor=color,linewidth=0)
                else: # lakes filled with background color by default
                    if lake_color is None:
                        poly = Polygon(xy,facecolor=axisbgc,edgecolor=axisbgc,linewidth=0)
                    else:
                        poly = Polygon(xy,facecolor=lake_color,edgecolor=lake_color,linewidth=0)
                if zorder is not None:
                    poly.set_zorder(zorder)
                if alpha is not None:
                    poly.set_alpha(alpha)
                ax.add_patch(poly)
                polys.append(poly)
            npoly = npoly + 1
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip continent polygons to map limbs
        polys,c = self._cliplimb(ax,polys)
        return polys

    def _cliplimb(self,ax,coll):
        if not self._mapboundarydrawn:
            return coll, None
        c = self._mapboundarydrawn
        if c not in ax.patches:
            p = ax.add_patch(c)
            #p.set_clip_on(False)
        try:
            coll.set_clip_path(c)
        except:
            for item in coll:
                item.set_clip_path(c)
        return coll,c

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

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        linewidth        coastline width (default 1.)
        linestyle        coastline linestyle (default solid)
        color            coastline color (default black)
        antialiased      antialiasing switch for coastlines (default True).
        ax               axes instance (overrides default axes instance)
        zorder           sets the zorder for the coastlines (if not specified,
                         uses default zorder for
                         matplotlib.patches.LineCollections).
        ==============   ====================================================

        returns a matplotlib.patches.LineCollection object.
        """
        if self.resolution is None:
            raise AttributeError('there are no boundary datasets associated with this Basemap instance')
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        coastlines = LineCollection(self.coastsegs,antialiaseds=(antialiased,))
        coastlines.set_color(color)
        coastlines.set_linestyle(linestyle)
        coastlines.set_linewidth(linewidth)
        coastlines.set_label('_nolabel_')
        if zorder is not None:
            coastlines.set_zorder(zorder)
        ax.add_collection(coastlines)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        coastlines,c = self._cliplimb(ax,coastlines)
        return coastlines

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

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        linewidth        country boundary line width (default 0.5)
        linestyle        coastline linestyle (default solid)
        color            country boundary line color (default black)
        antialiased      antialiasing switch for country boundaries (default
                         True).
        ax               axes instance (overrides default axes instance)
        zorder           sets the zorder for the country boundaries (if not
                         specified uses default zorder for
                         matplotlib.patches.LineCollections).
        ==============   ====================================================

        returns a matplotlib.patches.LineCollection object.
        """
        if self.resolution is None:
            raise AttributeError('there are no boundary datasets associated with this Basemap instance')
        # read in country line segments, only keeping those that
        # intersect map boundary polygon.
        if not hasattr(self,'cntrysegs'):
            self.cntrysegs, types = self._readboundarydata('countries')
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        countries = LineCollection(self.cntrysegs,antialiaseds=(antialiased,))
        countries.set_color(color)
        countries.set_linestyle(linestyle)
        countries.set_linewidth(linewidth)
        countries.set_label('_nolabel_')
        if zorder is not None:
            countries.set_zorder(zorder)
        ax.add_collection(countries)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip countries to map limbs
        countries,c = self._cliplimb(ax,countries)
        return countries

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

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        linewidth        state boundary line width (default 0.5)
        linestyle        coastline linestyle (default solid)
        color            state boundary line color (default black)
        antialiased      antialiasing switch for state boundaries
                         (default True).
        ax               axes instance (overrides default axes instance)
        zorder           sets the zorder for the state boundaries (if not
                         specified, uses default zorder for
                         matplotlib.patches.LineCollections).
        ==============   ====================================================

        returns a matplotlib.patches.LineCollection object.
        """
        if self.resolution is None:
            raise AttributeError('there are no boundary datasets associated with this Basemap instance')
        # read in state line segments, only keeping those that
        # intersect map boundary polygon.
        if not hasattr(self,'statesegs'):
            self.statesegs, types = self._readboundarydata('states')
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        states = LineCollection(self.statesegs,antialiaseds=(antialiased,))
        states.set_color(color)
        states.set_linestyle(linestyle)
        states.set_linewidth(linewidth)
        states.set_label('_nolabel_')
        if zorder is not None:
            states.set_zorder(zorder)
        ax.add_collection(states)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip states to map limbs
        states,c = self._cliplimb(ax,states)
        return states

    def drawcounties(self,linewidth=0.1,linestyle='solid',color='k',antialiased=1,
                     facecolor='none',ax=None,zorder=None,drawbounds=False):
        """
        Draw county boundaries in US. The county boundary shapefile
        originates with the NOAA Coastal Geospatial Data Project
        (http://coastalgeospatial.noaa.gov/data_gis.html).

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        linewidth        county boundary line width (default 0.1)
        linestyle        coastline linestyle (default solid)
        color            county boundary line color (default black)
        facecolor        fill color of county (default is no fill)
        antialiased      antialiasing switch for county boundaries
                         (default True).
        ax               axes instance (overrides default axes instance)
        zorder           sets the zorder for the county boundaries (if not
                         specified, uses default zorder for
                         matplotlib.patches.LineCollections).
        ==============   ====================================================

        returns a matplotlib.patches.LineCollection object.
        """
        ax = ax or self._check_ax()
        gis_file = os.path.join(basemap_datadir,'UScounties')
        county_info = self.readshapefile(gis_file,'counties',\
                      default_encoding='latin-1',drawbounds=drawbounds)
        counties = [coords for coords in self.counties]
        counties = PolyCollection(counties, antialiaseds=(antialiased,))
        counties.set_linestyle(linestyle)
        counties.set_linewidth(linewidth)
        counties.set_edgecolor(color)
        counties.set_facecolor(facecolor)
        counties.set_label('counties')
        if zorder:
            counties.set_zorder(zorder)
        ax.add_collection(counties)
        return counties

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

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        linewidth        river boundary line width (default 0.5)
        linestyle        coastline linestyle (default solid)
        color            river boundary line color (default black)
        antialiased      antialiasing switch for river boundaries (default
                         True).
        ax               axes instance (overrides default axes instance)
        zorder           sets the zorder for the rivers (if not
                         specified uses default zorder for
                         matplotlib.patches.LineCollections).
        ==============   ====================================================

        returns a matplotlib.patches.LineCollection object.
        """
        if self.resolution is None:
            raise AttributeError('there are no boundary datasets associated with this Basemap instance')
        # read in river line segments, only keeping those that
        # intersect map boundary polygon.
        if not hasattr(self,'riversegs'):
            self.riversegs, types = self._readboundarydata('rivers')
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        rivers = LineCollection(self.riversegs,antialiaseds=(antialiased,))
        rivers.set_color(color)
        rivers.set_linestyle(linestyle)
        rivers.set_linewidth(linewidth)
        rivers.set_label('_nolabel_')
        if zorder is not None:
            rivers.set_zorder(zorder)
        ax.add_collection(rivers)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip rivers to map limbs
        rivers,c = self._cliplimb(ax,rivers)
        return rivers

    def is_land(self,xpt,ypt):
        """
        Returns True if the given x,y point (in projection coordinates) is
        over land, False otherwise.  The definition of land is based upon
        the GSHHS coastline polygons associated with the class instance.
        Points over lakes inside land regions are not counted as land points.
        """
        if self.resolution is None: return None
        landpt = False
        for poly in self.landpolygons:
            landpt = _geoslib.Point((xpt,ypt)).within(poly)
            if landpt: break
        lakept = False
        for poly in self.lakepolygons:
            lakept = _geoslib.Point((xpt,ypt)).within(poly)
            if lakept: break
        return landpt and not lakept

    def readshapefile(self,shapefile,name,drawbounds=True,zorder=None,
                      linewidth=0.5,color='k',antialiased=1,ax=None,
                      default_encoding='utf-8',encoding_errors='strict'):
        """
        Read in shape file, optionally draw boundaries on map.

        .. note::
          - Assumes shapes are 2D
          - only works for Point, MultiPoint, Polyline and Polygon shapes.
          - vertices/points must be in geographic (lat/lon) coordinates.

        Mandatory Arguments:

        .. tabularcolumns:: |l|L|

        ================ ====================================================
        Argument         Description
        ================ ====================================================
        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 or points in 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 x,y vertices for each shape in
                         map projection  coordinates and self.counties_info
                         will be a list of dictionaries with shape
                         attributes.  Rings in individual Polygon
                         shapes are split out into separate polygons, and
                         additional keys 'RINGNUM' and 'SHAPENUM' are added
                         to the shape attribute dictionary.
        ================ ====================================================

        The following optional keyword arguments are only relevant for Polyline
        and Polygon shape types, for Point and MultiPoint shapes they are
        ignored.

        .. tabularcolumns:: |l|L|

        ================ =====================================================
        Keyword          Description
        ================ =====================================================
        drawbounds       draw boundaries of shapes (default True).
        zorder           shape boundary zorder (if not specified,
                         default for matplotlib.lines.LineCollection
                         is used).
        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)
        default_encoding encoding used to parse properties from .dbf files
                         (default utf-8)
        encoding_errors  encoding error handling (default strict), other
                         possible values: ignore, replace and backslashreplace
        ================ =====================================================

        A tuple (num_shapes, type, min, max) containing shape file info
        is returned.
        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. If ``drawbounds=True`` a
        matplotlib.patches.LineCollection object is appended to the tuple.
        """
        import shapefile as shp
        from shapefile import Reader
        shp.default_encoding = default_encoding
        if not os.path.exists('%s.shp'%shapefile):
            raise IOError('cannot locate %s.shp'%shapefile)
        if not os.path.exists('%s.shx'%shapefile):
            raise IOError('cannot locate %s.shx'%shapefile)
        if not os.path.exists('%s.dbf'%shapefile):
            raise IOError('cannot locate %s.dbf'%shapefile)
        # open shapefile, read vertices for each object, convert
        # to map projection coordinates (only works for 2D shape types).
        try:
            shf = Reader(shapefile, encoding=default_encoding,
                         encodingErrors=encoding_errors)
        except:
            raise IOError('error reading shapefile %s.shp' % shapefile)
        fields = shf.fields
        coords = []; attributes = []
        msg = " ".join([
            "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)"])
        shptype = shf.shapes()[0].shapeType
        bbox = shf.bbox.tolist()
        info = (shf.numRecords,shptype,bbox[0:2]+[0.,0.],bbox[2:]+[0.,0.])
        npoly = 0
        for shprec in shf.shapeRecords():
            shp = shprec.shape; rec = shprec.record
            npoly = npoly + 1
            if shptype != shp.shapeType:
                raise ValueError('readshapefile can only handle a single shape type per file')
            if shptype not in [1,3,5,8]:
                raise ValueError('readshapefile can only handle 2D shape types')
            verts = shp.points
            if shptype in [1,8]: # a Point or MultiPoint shape.
                lons, lats = list(zip(*verts))
                if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
                    raise ValueError(msg)
                # if latitude is slightly greater than 90, truncate to 90
                lats = [max(min(lat, 90.0), -90.0) for lat in lats]
                if len(verts) > 1: # MultiPoint
                    x,y = self(lons, lats)
                    coords.append(list(zip(x,y)))
                else: # single Point
                    x,y = self(lons[0], lats[0])
                    coords.append((x,y))
                attdict={}
                for r,key in zip(rec,fields[1:]):
                    attdict[key[0]]=r
                attributes.append(attdict)
            else: # a Polyline or Polygon shape.
                parts = shp.parts.tolist()
                ringnum = 0
                for indx1,indx2 in zip(parts,parts[1:]+[len(verts)]):
                    ringnum = ringnum + 1
                    lons, lats = list(zip(*verts[indx1:indx2]))
                    if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
                        raise ValueError(msg)
                    # if latitude is slightly greater than 90, truncate to 90
                    lats = [max(min(lat, 90.0), -90.0) for lat in lats]
                    x, y = self(lons, lats)
                    coords.append(list(zip(x,y)))
                    attdict={}
                    for r,key in zip(rec,fields[1:]):
                        attdict[key[0]]=r
                    # add information about ring number to dictionary.
                    attdict['RINGNUM'] = ringnum
                    attdict['SHAPENUM'] = npoly
                    attributes.append(attdict)
        # draw shape boundaries for polylines, polygons  using LineCollection.
        if shptype not in [1,8] and drawbounds:
            # get current axes instance (if none specified).
            ax = ax or self._check_ax()
            # make LineCollections for each polygon.
            lines = LineCollection(coords,antialiaseds=(antialiased,))
            lines.set_color(color)
            lines.set_linewidth(linewidth)
            lines.set_label('_nolabel_')
            if zorder is not None:
               lines.set_zorder(zorder)
            ax.add_collection(lines)
            # set axes limits to fit map region.
            self.set_axes_limits(ax=ax)
            # clip boundaries to map limbs
            lines,c = self._cliplimb(ax,lines)
            info = info + (lines,)
        self.__dict__[name]=coords
        self.__dict__[name+'_info']=attributes
        return info

    def drawparallels(self,circles,color='k',textcolor='k',linewidth=1.,zorder=None, \
                      dashes=[1,1],labels=[0,0,0,0],labelstyle=None, \
                      fmt='%g',xoffset=None,yoffset=None,ax=None,latmax=None,
                      **text_kwargs):
        r"""
        Draw and label parallels (latitude lines) for values (in degrees)
        given in the sequence ``circles``.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        color            color to draw parallels (default black).
        textcolor        color to draw labels (default black).
        linewidth        line width for parallels (default 1.)
        zorder           sets the zorder for parallels (if not specified,
                         uses default zorder for matplotlib.lines.Line2D
                         objects).
        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
                         and bottom of the plot, but not the right and top.
        labelstyle       if set to "+/-", north and south latitudes are
                         labelled with "+" and "-", otherwise they are
                         labelled with "N" and "S".
        fmt              a format string to format the parallel labels
                         (default '%g') **or** a function that takes a
                         latitude value in degrees as it's only argument
                         and returns a formatted string.
        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)
        latmax           absolute value of latitude to which meridians are drawn
                         (default is 80).
        \**text_kwargs   additional keyword arguments controlling text
                         for labels that are passed on to
                         the text method of the axes instance (see
                         matplotlib.pyplot.text documentation).
        ==============   ====================================================

        returns a dictionary whose keys are the parallel values, and
        whose values are tuples containing lists of the
        matplotlib.lines.Line2D and matplotlib.text.Text instances
        associated with each parallel. Deleting an item from the
        dictionary removes the corresponding parallel from the plot.
        """
        text_kwargs['color']=textcolor # pass textcolor kwarg on to ax.text
        # if celestial=True, don't use "N" and "S" labels.
        if labelstyle is None and self.celestial:
            labelstyle="+/-"
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # don't draw meridians past latmax, always draw parallel at latmax.
        if latmax is None: 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 _cylproj + _pseudocyl:
            lons = np.linspace(self.llcrnrlon, self.urcrnrlon, 10001)
        elif self.projection in ['tmerc']:
            lon_0 = self.projparams['lon_0']
            # tmerc only defined within +/- 90 degrees of lon_0
            lons = np.linspace(lon_0-90,lon_0+90,100001)
        else:
            lonmin = self.boundarylonmin; lonmax = self.boundarylonmax
            lons = np.linspace(lonmin, lonmax, 10001)
        # make sure latmax degree parallel is drawn if projection not merc or cyl or miller
        try:
            circlesl = list(circles)
        except:
            circlesl = circles
        if self.projection not in _cylproj + _pseudocyl:
            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)
        linecolls = {}
        for circ in circlesl:
            lats = circ*np.ones(len(lons),np.float32)
            x,y = self(lons,lats)
            # remove points outside domain.
            # leave a little slop around edges (3*xdelta)
            # don't really know why, but this appears to be needed to
            # or lines sometimes don't reach edge of plot.
            testx = np.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta)
            x = np.compress(testx, x)
            y = np.compress(testx, y)
            testy = np.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta)
            x = np.compress(testy, x)
            y = np.compress(testy, y)
            lines = []
            if len(x) > 1 and len(y) > 1:
                # split into separate line segments if necessary.
                # (not necessary for cylindrical or pseudocylindricl projections)
                xd = (x[1:]-x[0:-1])**2
                yd = (y[1:]-y[0:-1])**2
                dist = np.sqrt(xd+yd)
                if self.projection not in ['cyl','rotpole']:
                    split = dist > self.rmajor/10.
                else:
                    split = dist > 1.
                if np.sum(split) and self.projection not in _cylproj:
                    ind = (np.compress(split,np.squeeze(split*np.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)
                        l.set_color(color)
                        l.set_dashes(dashes)
                        l.set_label('_nolabel_')
                        if zorder is not None:
                            l.set_zorder(zorder)
                        ax.add_line(l)
                        lines.append(l)
            linecolls[circ] = (lines,[])
        # draw labels for parallels
        # parallels not labelled for fulldisk orthographic or geostationary
        if self.projection in ['ortho','geos','nsper','vandg','aeqd'] and max(labels):
            if self.projection == 'vandg' or self._fulldisk:
                sys.stdout.write('Warning: Cannot label parallels on %s basemap' % _projnames[self.projection])
                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.
        dx = (self.xmax-self.xmin)/1000.
        dy = (self.ymax-self.ymin)/1000.
        if self.projection in _pseudocyl:
            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 _cylproj + _pseudocyl and side in ['t','b']: continue
            if side in ['l','r']:
                nmax = int((self.ymax-self.ymin)/dy+1)
                yy = np.linspace(self.llcrnry,self.urcrnry,nmax)
                if side == 'l':
                    if self.projection in _pseudocyl:
                        lats = np.linspace(-89.99,89,99,nmax)
                        if self.celestial:
                            lons = (self.projparams['lon_0']+180.)*np.ones(len(lats),lats.dtype)
                        else:
                            lons = (self.projparams['lon_0']-180.)*np.ones(len(lats),lats.dtype)
                        xx, yy = self(lons, lats)
                    else:
                        xx = self.llcrnrx*np.ones(yy.shape,yy.dtype)
                        lons,lats = self(xx,yy,inverse=True)
                        lons = lons.tolist(); lats = lats.tolist()
                else:
                    if self.projection in _pseudocyl:
                        lats = np.linspace(-89.99,89,99,nmax)
                        if self.celestial:
                           lons = (self.projparams['lon_0']-180.)*np.ones(len(lats),lats.dtype)
                        else:
                           lons = (self.projparams['lon_0']+180.)*np.ones(len(lats),lats.dtype)
                        xx, yy = self(lons, lats)
                    else:
                        xx = self.urcrnrx*np.ones(yy.shape,yy.dtype)
                        lons,lats = self(xx,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 = np.linspace(self.llcrnrx,self.urcrnrx,nmax)
                if side == 'b':
                    lons,lats = self(xx,self.llcrnry*np.ones(xx.shape,np.float32),inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(xx,self.urcrnry*np.ones(xx.shape,np.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:
                # don't label parallels for round polar plots
                if self.round: continue
                # 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
                latlab = _setlatlab(fmt,lat,labelstyle)
                # 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:
                        t = None
                        if side == 'l':
                            if self.projection in _pseudocyl:
                                if self.celestial:
                                    xlab,ylab = self(lon_0+179.9,lat)
                                else:
                                    xlab,ylab = self(lon_0-179.9,lat)
                            else:
                                xlab = self.llcrnrx
                            xlab = xlab-xoffset
                            if self.projection in _pseudocyl:
                                if lat>0:
                                   t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='bottom',**text_kwargs)
                                elif lat<0:
                                   t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='top',**text_kwargs)
                                else:
                                   t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='center',**text_kwargs)
                            else:
                               t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='center',**text_kwargs)
                        elif side == 'r':
                            if self.projection in _pseudocyl:
                                if self.celestial:
                                   xlab,ylab = self(lon_0-179.9,lat)
                                else:
                                   xlab,ylab = self(lon_0+179.9,lat)
                            else:
                                xlab = self.urcrnrx
                            xlab = xlab+xoffset
                            if self.projection in _pseudocyl:
                                if lat>0:
                                   t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='bottom',**text_kwargs)
                                elif lat<0:
                                   t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='top',**text_kwargs)
                                else:
                                   t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='center',**text_kwargs)
                            else:
                               t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='center',**text_kwargs)
                        elif side == 'b':
                            t = ax.text(xx[n],self.llcrnry-yoffset,latlab,horizontalalignment='center',verticalalignment='top',**text_kwargs)
                        else:
                            t = ax.text(xx[n],self.urcrnry+yoffset,latlab,horizontalalignment='center',verticalalignment='bottom',**text_kwargs)
                        if t is not None: linecolls[lat][1].append(t)

        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        keys = list(linecolls.keys()); vals = list(linecolls.values())
        for k,v in zip(keys,vals):
            if v == ([], []):
                del linecolls[k]
            # add a remove method to each tuple.
            else:
                linecolls[k] = _tup(linecolls[k])
        # override __delitem__ in dict to call remove() on values.
        pardict = _dict(linecolls)
        # clip parallels for round polar plots (and delete labels).
        for lines, _ in pardict.values():
            self._cliplimb(ax, lines)
        return pardict

    def drawmeridians(self,meridians,color='k',textcolor='k',linewidth=1., zorder=None,\
                      dashes=[1,1],labels=[0,0,0,0],labelstyle=None,\
                      fmt='%g',xoffset=None,yoffset=None,ax=None,latmax=None,
                      **text_kwargs):
        r"""
        Draw and label meridians (longitude lines) for values (in degrees)
        given in the sequence ``meridians``.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        color            color to draw meridians (default black).
        textcolor        color to draw labels (default black).
        linewidth        line width for meridians (default 1.)
        zorder           sets the zorder for meridians (if not specified,
                         uses default zorder for matplotlib.lines.Line2D
                         objects).
        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
                         and bottom of the plot, but not the right and top.
        labelstyle       if set to "+/-", east and west longitudes are
                         labelled with "+" and "-", otherwise they are
                         labelled with "E" and "W".
        fmt              a format string to format the meridian labels
                         (default '%g') **or** a function that takes a
                         longitude value in degrees as it's only argument
                         and returns a formatted string.
        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)
        latmax           absolute value of latitude to which meridians are drawn
                         (default is 80).
        \**text_kwargs   additional keyword arguments controlling text
                         for labels that are passed on to
                         the text method of the axes instance (see
                         matplotlib.pyplot.text documentation).
        ==============   ====================================================

        returns a dictionary whose keys are the meridian values, and
        whose values are tuples containing lists of the
        matplotlib.lines.Line2D and matplotlib.text.Text instances
        associated with each meridian. Deleting an item from the
        dictionary removes the correpsonding meridian from the plot.
        """
        text_kwargs['color']=textcolor # pass textcolor kwarg on to ax.text
        # for cylindrical projections, try to handle wraparound (i.e. if
        # projection is defined in -180 to 0 and user asks for meridians from
        # 180 to 360 to be drawn, it should work)
        if self.projection in _cylproj or self.projection in _pseudocyl:
            def addlon(meridians,madd):
                minside = (madd >= self.llcrnrlon and madd <= self.urcrnrlon)
                if minside and madd not in meridians: meridians.append(madd)
                return meridians
            merids = list(meridians)
            meridians = []
            for m in merids:
                meridians = addlon(meridians,m)
                meridians = addlon(meridians,m+360)
                meridians = addlon(meridians,m-360)
            meridians.sort()
        # if celestial=True, don't use "E" and "W" labels.
        if labelstyle is None and self.celestial:
            labelstyle="+/-"
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # don't draw meridians past latmax, always draw parallel at latmax.
        if latmax is None: latmax = 80. # unused w/ cyl, merc or miller proj.
        # 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.

        lats = np.linspace(self.latmin,self.latmax,10001)
        if self.projection not in _cylproj + _pseudocyl:
            testlat = np.logical_and(lats>-latmax,lats<latmax)
            lats = np.compress(testlat,lats)

        xdelta = 0.01*(self.xmax-self.xmin)
        ydelta = 0.01*(self.ymax-self.ymin)
        linecolls = {}
        for merid in meridians:
            lons = merid*np.ones(len(lats),np.float32)
            x,y = self(lons,lats)
            # remove points outside domain.
            # leave a little slop around edges (3*xdelta)
            # don't really know why, but this appears to be needed to
            # or lines sometimes don't reach edge of plot.
            testx = np.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta)
            x = np.compress(testx, x)
            y = np.compress(testx, y)
            testy = np.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta)
            x = np.compress(testy, x)
            y = np.compress(testy, y)
            lines = []
            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 = np.sqrt(xd+yd)
                if self.projection not in ['cyl','rotpole']:
                    split = dist > self.rmajor/10.
                else:
                    split = dist > 1.
                if np.sum(split) and self.projection not in _cylproj:
                    ind = (np.compress(split,np.squeeze(split*np.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)
                        l.set_color(color)
                        l.set_dashes(dashes)
                        l.set_label('_nolabel_')
                        if zorder is not None:
                            l.set_zorder(zorder)
                        ax.add_line(l)
                        lines.append(l)
            linecolls[merid] = (lines,[])
        # draw labels for meridians.
        # meridians not labelled for sinusoidal, hammer, mollweide,
        # VanDerGrinten or full-disk orthographic/geostationary.
        if self.projection in ['sinu','moll','hammer','vandg'] and max(labels):
            sys.stdout.write('Warning: Cannot label meridians on %s basemap' % _projnames[self.projection])
            labels = [0,0,0,0]
        if self.projection in ['ortho','geos','nsper','aeqd'] and max(labels):
            if self._fulldisk and self.boundinglat is None:
                sys.stdout.write(" ".join([
                    "Warning: Cannot label meridians on full-disk Geostationary,"
                    "Orthographic or Azimuthal equidistant 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.
        dx = (self.xmax-self.xmin)/1000.
        dy = (self.ymax-self.ymin)/1000.
        if self.projection in _pseudocyl:
            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 or self.round: continue
            # for cylindrical projections, don't draw meridians on left or right.
            if self.projection in _cylproj + _pseudocyl and side in ['l','r']: continue
            if side in ['l','r']:
                nmax = int((self.ymax-self.ymin)/dy+1)
                yy = np.linspace(self.llcrnry,self.urcrnry,nmax)
                if side == 'l':
                    lons,lats = self(self.llcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(self.urcrnrx*np.ones(yy.shape,np.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)
                if self.projection in _pseudocyl:
                    xx = np.linspace(xmin,xmax,nmax)
                else:
                    xx = np.linspace(self.llcrnrx,self.urcrnrx,nmax)
                if side == 'b':
                    lons,lats = self(xx,self.llcrnry*np.ones(xx.shape,np.float32),inverse=True)
                    lons = lons.tolist(); lats = lats.tolist()
                else:
                    lons,lats = self(xx,self.urcrnry*np.ones(xx.shape,np.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
                lon2 = (lon+360) % 360
                # find index of meridian (there may be two, so
                # search from left and right).
                nl = _searchlist(lons,lon2)
                nr = _searchlist(lons[::-1],lon2)
                if nr != -1: nr = len(lons)-nr-1
                lonlab = _setlonlab(fmt,lon2,labelstyle)
                # 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 _cylproj and lat > latmax: continue
                    # don't bother if close to the first label.
                    if i and abs(nr-nl) < 100: continue
                    if n >= 0:
                        t = None
                        if side == 'l':
                            t = ax.text(self.llcrnrx-xoffset,yy[n],lonlab,horizontalalignment='right',verticalalignment='center',**text_kwargs)
                        elif side == 'r':
                            t = ax.text(self.urcrnrx+xoffset,yy[n],lonlab,horizontalalignment='left',verticalalignment='center',**text_kwargs)
                        elif side == 'b':
                            t = ax.text(xx[n],self.llcrnry-yoffset,lonlab,horizontalalignment='center',verticalalignment='top',**text_kwargs)
                        else:
                            t = ax.text(xx[n],self.urcrnry+yoffset,lonlab,horizontalalignment='center',verticalalignment='bottom',**text_kwargs)

                        if t is not None: linecolls[lon][1].append(t)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # remove empty values from linecolls dictionary
        keys = list(linecolls.keys()); vals = list(linecolls.values())
        for k,v in zip(keys,vals):
            if v == ([], []):
                del linecolls[k]
            else:
            # add a remove method to each tuple.
                linecolls[k] = _tup(linecolls[k])
        # override __delitem__ in dict to call remove() on values.
        meridict = _dict(linecolls)
        # for round polar plots, clip meridian lines and label them.
        if self.round:
            # label desired?
            label = False
            for lab in labels:
                if lab: label = True
            for merid in meridict:
                if not label: continue
                # label
                lonlab = _setlonlab(fmt,merid,labelstyle)
                x,y = self(merid,self.boundinglat)
                r = np.sqrt((x-0.5*(self.xmin+self.xmax))**2+
                            (y-0.5*(self.ymin+self.ymax))**2)
                r = r + np.sqrt(xoffset**2+yoffset**2)
                if self.projection.startswith('np'):
                    pole = 1
                elif self.projection.startswith('sp'):
                    pole = -1
                elif self.projection == 'ortho' and self.round:
                    pole = 1
                if pole == 1:
                    theta = (np.pi/180.)*(merid-self.projparams['lon_0']-90)
                    if self.projection == 'ortho' and\
                       self.projparams['lat_0'] == -90:
                        theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90)
                    x = r*np.cos(theta)+0.5*(self.xmin+self.xmax)
                    y = r*np.sin(theta)+0.5*(self.ymin+self.ymax)
                    if x > 0.5*(self.xmin+self.xmax)+xoffset:
                        horizalign = 'left'
                    elif x < 0.5*(self.xmin+self.xmax)-xoffset:
                        horizalign = 'right'
                    else:
                        horizalign = 'center'
                    if y > 0.5*(self.ymin+self.ymax)+yoffset:
                        vertalign = 'bottom'
                    elif y < 0.5*(self.ymin+self.ymax)-yoffset:
                        vertalign = 'top'
                    else:
                        vertalign = 'center'
                    # labels [l,r,t,b]
                    if labels[0] and not labels[1] and x >= 0.5*(self.xmin+self.xmax)+xoffset: continue
                    if labels[1] and not labels[0] and x <= 0.5*(self.xmin+self.xmax)-xoffset: continue
                    if labels[2] and not labels[3] and y <= 0.5*(self.ymin+self.ymax)-yoffset: continue
                    if labels[3] and not labels[2]and y >= 0.5*(self.ymin+self.ymax)+yoffset: continue
                elif pole == -1:
                    theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90)
                    x = r*np.cos(theta)+0.5*(self.xmin+self.xmax)
                    y = r*np.sin(theta)+0.5*(self.ymin+self.ymax)
                    if x > 0.5*(self.xmin+self.xmax)-xoffset:
                        horizalign = 'right'
                    elif x < 0.5*(self.xmin+self.xmax)+xoffset:
                        horizalign = 'left'
                    else:
                        horizalign = 'center'
                    if y > 0.5*(self.ymin+self.ymax)-yoffset:
                        vertalign = 'top'
                    elif y < 0.5*(self.ymin+self.ymax)+yoffset:
                        vertalign = 'bottom'
                    else:
                        vertalign = 'center'
                    # labels [l,r,t,b]
                    if labels[0] and not labels[1] and x <=  0.5*(self.xmin+self.xmax)+xoffset: continue
                    if labels[1] and not labels[0] and x >=  0.5*(self.xmin+self.xmax)-xoffset: continue
                    if labels[2] and not labels[3] and y >=  0.5*(self.ymin+self.ymax)-yoffset: continue
                    if labels[3] and not labels[2] and y <=  0.5*(self.ymin+self.ymax)+yoffset: continue
                t=ax.text(x,y,lonlab,horizontalalignment=horizalign,verticalalignment=vertalign,**text_kwargs)
                meridict[merid][1].append(t)
        for lines, _ in meridict.values():
            self._cliplimb(ax, lines)
        return meridict

    def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs):
        r"""
        Draw a polygon centered at ``lon_0,lat_0``.  The polygon
        approximates a circle on the surface of the earth with radius
        ``radius_deg`` degrees latitude along longitude ``lon_0``,
        made up of ``npts`` vertices.
        The polygon represents a Tissot's indicatrix
        (http://en.wikipedia.org/wiki/Tissot's_Indicatrix),
        which when drawn on a map shows the distortion
        inherent in the map projection.

        .. note::
         Cannot handle situations in which the polygon intersects
         the edge of the map projection domain, and then re-enters the domain.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.patches.Polygon.

        returns a matplotlib.patches.Polygon object."""
        ax = kwargs.pop('ax', None) or self._check_ax()
        g = pyproj.Geod(a=self.rmajor,b=self.rminor)
        az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
        seg = [self(lon_0,lat_0+radius_deg)]
        delaz = 360./npts
        az = az12
        for n in range(npts):
            az = az+delaz
            lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
            x,y = self(lon,lat)
            # add segment if it is in the map projection region.
            if x < 1.e20 and y < 1.e20:
                seg.append((x,y))
        poly = Polygon(seg,**kwargs)
        ax.add_patch(poly)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip polygons to map limbs
        poly,c = self._cliplimb(ax,poly)
        return poly

    def gcpoints(self,lon1,lat1,lon2,lat2,npoints):
        """
        compute ``points`` points along a great circle with endpoints
        ``(lon1,lat1)`` and ``(lon2,lat2)``.

        Returns arrays x,y with map projection coordinates.
        """
        gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
        lonlats = gc.npts(lon1,lat1,lon2,lat2,npoints-2)
        lons=[lon1];lats=[lat1]
        for lon,lat in lonlats:
            lons.append(lon); lats.append(lat)
        lons.append(lon2); lats.append(lat2)
        x, y = self(lons, lats)
        return x,y

    def drawgreatcircle(self,lon1,lat1,lon2,lat2,del_s=100.,**kwargs):
        r"""
        Draw a great circle on the map from the longitude-latitude
        pair ``lon1,lat1`` to ``lon2,lat2``

        .. tabularcolumns:: |l|L|

        ==============   =======================================================
        Keyword          Description
        ==============   =======================================================
        del_s            points on great circle computed every del_s kilometers
                         (default 100).
        \**kwargs        other keyword arguments are passed on to :meth:`plot`
                         method of Basemap instance.
        ==============   =======================================================

        Returns a list with a single ``matplotlib.lines.Line2D`` object like a
        call to ``pyplot.plot()``.
        """
        # use great circle formula for a perfect sphere.
        gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
        az12,az21,dist = gc.inv(lon1,lat1,lon2,lat2)
        npoints = int((dist+0.5*1000.*del_s)/(1000.*del_s))
        lonlats = gc.npts(lon1,lat1,lon2,lat2,npoints)
        lons = [lon1]; lats = [lat1]
        for lon, lat in lonlats:
            lons.append(lon)
            lats.append(lat)
        lons.append(lon2); lats.append(lat2)
        x, y = self(lons, lats)

        # Correct wrap around effect of great circles

        # get points
        _p = self.plot(x,y,**kwargs)
        p = _p[0].get_path()

        # since we know the difference between any two points, we can use this to find wrap arounds on the plot
        max_dist = 1000*del_s*2

        # calculate distances and compare with max allowable distance
        dists = np.abs(np.diff(p.vertices[:,0]))
        cuts = np.where( dists > max_dist )[0]

        # if there are any cut points, cut them and begin again at the next point
        for i,k in enumerate(cuts):
            # vertex to cut at
            cut_point = cuts[i]

            # create new vertices with a nan inbetween and set those as the path's vertices
            verts = np.concatenate(
                                       [p.vertices[:cut_point, :],
                                        [[np.nan, np.nan]],
                                        p.vertices[cut_point+1:, :]]
                                       )
            p.codes = None
            p.vertices = verts

        return _p

    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`` by ``nx``
        map projection grid.  Typically used to transform data to
        map projection coordinates for plotting on a map with
        the :meth:`imshow`.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Argument         Description
        ==============   ====================================================
        datin            input data on a lat/lon grid.
        lons, lats       rank-1 arrays containing longitudes and latitudes
                         (in degrees) of input data in increasing order.
                         For non-cylindrical projections (those other than
                         ``cyl``, ``merc``, ``cea``, ``gall`` and ``mill``) lons
                         must fit within range -180 to 180.
        nx, ny           The size of the output regular grid in map
                         projection coordinates
        ==============   ====================================================

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        returnxy         If True, the x and y values of the map
                         projection grid are also returned (Default False).
        checkbounds      If True, values of lons and lats are checked to see
                         that they lie within the map projection region.
                         Default is False, and data outside map projection
                         region is clipped to values on boundary.
        masked           If True, interpolated data is returned as a masked
                         array with values outside map projection region
                         masked (Default False).
        order            0 for nearest-neighbor interpolation, 1 for
                         bilinear, 3 for cubic spline (Default 1).
                         Cubic spline interpolation requires scipy.ndimage.
        ==============   ====================================================

        Returns ``datout`` (data on map projection grid).
        If returnxy=True, returns ``data,x,y``.
        """
        # 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 _cylproj:
            lonsa = np.array(lons)
            count = np.sum(lonsa < -180.00001) + np.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)
        datout = interp(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`` by ``nx`` map projection grid.

        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.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Arguments        Description
        ==============   ====================================================
        uin, vin         input vector field on a lat/lon grid.
        lons, lats       rank-1 arrays containing longitudes and latitudes
                         (in degrees) of input data in increasing order.
                         For non-cylindrical projections (those other than
                         ``cyl``, ``merc``, ``cea``, ``gall`` and ``mill``) lons
                         must fit within range -180 to 180.
        nx, ny           The size of the output regular grid in map
                         projection coordinates
        ==============   ====================================================

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keyword          Description
        ==============   ====================================================
        returnxy         If True, the x and y values of the map
                         projection grid are also returned (Default False).
        checkbounds      If True, values of lons and lats are checked to see
                         that they lie within the map projection region.
                         Default is False, and data outside map projection
                         region is clipped to values on boundary.
        masked           If True, interpolated data is returned as a masked
                         array with values outside map projection region
                         masked (Default False).
        order            0 for nearest-neighbor interpolation, 1 for
                         bilinear, 3 for cubic spline (Default 1).
                         Cubic spline interpolation requires scipy.ndimage.
        ==============   ====================================================

        Returns ``uout, vout`` (vector field on map projection grid).
        If returnxy=True, returns ``uout,vout,x,y``.
        """
        # 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 _cylproj:
            lonsa = np.array(lons)
            count = np.sum(lonsa < -180.00001) + np.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.
        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)
        # rotate from geographic to map coordinates.
        return self.rotate_vector(uin,vin,lonsout,latsout,returnxy=returnxy)

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

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

        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.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Arguments        Description
        ==============   ====================================================
        uin, vin         input vector field on a lat/lon grid.
        lons, lats       Arrays containing longitudes and latitudes
                         (in degrees) of input data in increasing order.
                         For non-cylindrical projections (those other than
                         ``cyl``, ``merc``, ``cyl``, ``gall`` and ``mill``) lons
                         must fit within range -180 to 180.
        ==============   ====================================================

        Returns ``uout, vout`` (rotated vector field).
        If the optional keyword argument
        ``returnxy`` is True (default is False),
        returns ``uout,vout,x,y`` (where ``x,y`` are the map projection
        coordinates of the grid defined by ``lons,lats``).
        """
        # if lons,lats are 1d and uin,vin are 2d, and
        # lats describes 1st dim of uin,vin, and
        # lons describes 2nd dim of uin,vin, make lons,lats 2d
        # with meshgrid.
        if lons.ndim == lats.ndim == 1 and uin.ndim == vin.ndim == 2 and\
           uin.shape[1] == vin.shape[1] == lons.shape[0] and\
           uin.shape[0] == vin.shape[0] == lats.shape[0]:
            lons, lats = np.meshgrid(lons, lats)
        else:
            if not lons.shape == lats.shape == uin.shape == vin.shape:
                raise TypeError("shapes of lons,lats and uin,vin don't match")
        x, y = self(lons, lats)
        # rotate from geographic to map coordinates.
        if ma.isMaskedArray(uin):
            mask = ma.getmaskarray(uin)
            masked = True
            uin = uin.filled(1)
            vin = vin.filled(1)
        else:
            masked = False

        # Map the (lon, lat) vector in the complex plane.
        uvc = uin + 1j*vin
        uvmag = np.abs(uvc)
        theta = np.angle(uvc)

        # Define a displacement (dlon, dlat) that moves all
        # positions (lons, lats) a small distance in the
        # direction of the original vector.
        dc = 1E-5 * np.exp(theta*1j)
        dlat = dc.imag * np.cos(np.radians(lats))
        dlon = dc.real

        # Deal with displacements that overshoot the North or South Pole.
        farnorth = np.abs(lats+dlat) >= 90.0
        somenorth = farnorth.any()
        if somenorth:
            dlon[farnorth] *= -1.0
            dlat[farnorth] *= -1.0

        # Add displacement to original location and find the native coordinates.
        lon1 = lons + dlon
        lat1 = lats + dlat
        xn, yn = self(lon1, lat1)

        # Determine the angle of the displacement in the native coordinates.
        vecangle = np.arctan2(yn-y, xn-x)
        if somenorth:
            vecangle[farnorth] += np.pi

        # Compute the x-y components of the original vector.
        uvcout = uvmag * np.exp(1j*vecangle)
        uout = uvcout.real
        vout = uvcout.imag

        if masked:
            uout = ma.array(uout, mask=mask)
            vout = ma.array(vout, mask=mask)
        if returnxy:
            return uout,vout,x,y
        else:
            return uout,vout

    def set_axes_limits(self,ax=None):
        """
        Final step in Basemap method wrappers of Axes plotting methods:

        Set axis limits, fix aspect ratio for map domain using current
        or specified axes instance.  This is done only once per axes
        instance.

        In interactive mode, this method always calls draw_if_interactive
        before returning.

        """
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()

        # If we have already set the axes limits, and if the user
        # has not defeated this by turning autoscaling back on,
        # then all we need to do is plot if interactive.
        if (hash(ax) in self._initialized_axes
                                 and not ax.get_autoscalex_on()
                                 and not ax.get_autoscaley_on()):
            if mpl.is_interactive():
                import matplotlib.pyplot as plt
                plt.draw_if_interactive()
            return

        self._initialized_axes.add(hash(ax))
        # Take control of axis scaling:
        ax.set_autoscale_on(False)
        # 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))
        # if map boundary not yet drawn for elliptical maps, draw it with default values.
        if not self._mapboundarydrawn or self._mapboundarydrawn not in ax.patches:
            # elliptical map, draw boundary manually.
            if ((self.projection in ['ortho', 'geos', 'nsper', 'aeqd'] and
                    self._fulldisk) or self.round or
                    self.projection in _pseudocyl):
                # first draw boundary, no fill
                limb1 = self.drawmapboundary(fill_color='none', ax=ax)
                # draw another filled patch, with no boundary.
                limb2 = self.drawmapboundary(fill_color='none', linewidth=0, ax=ax)
                self._mapboundarydrawn = limb2
        # for elliptical map, always turn off axis_frame.
        if ((self.projection in ['ortho', 'geos', 'nsper', 'aeqd'] and
                self._fulldisk) or self.round or
                self.projection in _pseudocyl):
            # turn off axes frame.
            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)
        if self.fix_aspect:
            ax.set_aspect('equal',anchor=self.anchor)
        else:
            ax.set_aspect('auto',anchor=self.anchor)
        # make sure axis ticks are turned off.
        if self.noticks:
            ax.set_xticks([])
            ax.set_yticks([])
        # force draw if in interactive mode.
        if mpl.is_interactive():
            import matplotlib.pyplot as plt
            plt.draw_if_interactive()

    def _save_use_hold(self, ax, kwargs):
        h = kwargs.pop('hold', None)
        if hasattr(ax, '_hold'):
            self._tmp_hold = ax._hold
            if h is not None:
                ax._hold = h

    def _restore_hold(self, ax):
        if hasattr(ax, '_hold'):
            ax._hold = self._tmp_hold

    @_transform1d
    def scatter(self, *args, **kwargs):
        r"""
        Plot points with markers on the map
        (see matplotlib.pyplot.scatter documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axes instance.

        Other \**kwargs passed on to matplotlib.pyplot.scatter.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.scatter(*args, **kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        return ret

    @_transform1d
    def plot(self, *args, **kwargs):
        r"""
        Draw lines and/or markers on the map
        (see matplotlib.pyplot.plot documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.pyplot.plot.
        """
        ax = kwargs.pop('ax', None) or self._check_ax()
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.plot(*args, **kwargs)
        finally:
            self._restore_hold(ax)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        return ret

    def imshow(self, *args, **kwargs):
        r"""
        Display an image over the map
        (see matplotlib.pyplot.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.

        Other \**kwargs passed on to matplotlib.pyplot.plot.

        returns an matplotlib.image.AxesImage instance.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        kwargs['extent']=(self.llcrnrx,self.urcrnrx,self.llcrnry,self.urcrnry)
        # use origin='lower', unless overridden.
        if 'origin' not in kwargs:
            kwargs['origin']='lower'
        self._save_use_hold(ax, kwargs)
        try:
            with warnings.catch_warnings():
                warnings.simplefilter("ignore", category=DeprecationWarning)
                ret = ax.imshow(*args, **kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip image to map limbs
        ret,c = self._cliplimb(ax,ret)
        return ret

    @_transform
    def pcolor(self,x,y,data,**kwargs):
        r"""
        Make a pseudo-color plot over the map
        (see matplotlib.pyplot.pcolor documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        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.

        If ``tri`` is set to ``True``, an unstructured grid is assumed
        (x,y,data must be 1-d) and matplotlib.pyplot.tripcolor is used.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.pyplot.pcolor (or tripcolor if
        ``tri=True``).

        Note: (taken from matplotlib.pyplot.pcolor documentation)
        Ideally the dimensions of x and y should be one greater than those of data;
        if the dimensions are the same, then the last row and column of data will be ignored.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            if kwargs.pop('tri', False):
                import matplotlib.tri as tri
                # for unstructured grids, toss out points outside
                # projection limb (don't use those points in triangulation).
                if ma.isMA(data):
                    data = data.filled(fill_value=1.e30)
                    masked=True
                else:
                    masked=False
                mask = np.logical_or(x<1.e20,y<1.e20)
                x = np.compress(mask,x)
                y = np.compress(mask,y)
                data = np.compress(mask,data)
                if masked:
                    triang = tri.Triangulation(x, y)
                    z = data[triang.triangles]
                    mask = (z > 1.e20).sum(axis=-1)
                    triang.set_mask(mask)
                    ret = ax.tripcolor(triang,data,**kwargs)
                else:
                    ret = ax.tripcolor(x,y,data,**kwargs)
            else:
                # make x,y masked arrays
                # (masked where data is outside of projection limb)
                x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
                y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
                ret = ax.pcolor(x,y,data,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        if self.round:
            # for some reason, frame gets turned on.
            ax.set_frame_on(False)
        return ret

    @_transform
    def pcolormesh(self,x,y,data,**kwargs):
        r"""
        Make a pseudo-color plot over the map
        (see matplotlib.pyplot.pcolormesh documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.pyplot.pcolormesh.

        Note: (taken from matplotlib.pyplot.pcolor documentation)
        Ideally the dimensions of x and y should be one greater than those of data;
        if the dimensions are the same, then the last row and column of data will be ignored.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        # fix for invalid grid points
        if ((np.any(x > 1e20) or np.any(y > 1e20)) and
            x.ndim == 2 and y.ndim == 2):
            if x.shape != y.shape:
                raise ValueError('pcolormesh: x and y need same dimension')
            nx,ny = x.shape
            if nx < data.shape[0] or ny < data.shape[1]:
                raise ValueError('pcolormesh: data dimension needs to be at least that of x and y.')
            mask = (
                (x[:-1,:-1] > 1e20) |
                (x[1:,:-1] > 1e20) |
                (x[:-1,1:] > 1e20) |
                (x[1:,1:] > 1e20) |
                (y[:-1,:-1] > 1e20) |
                (y[1:,:-1] > 1e20) |
                (y[:-1,1:] > 1e20) |
                (y[1:,1:] > 1e20)
                )
            # we do not want to overwrite original array
            data = data[:nx-1,:ny-1].copy()
            data[mask] = np.nan
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.pcolormesh(x,y,data,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        if self.round:
            # for some reason, frame gets turned on.
            ax.set_frame_on(False)
        return ret

    def hexbin(self,x,y,**kwargs):
        r"""
        Make a hexagonal binning plot of x versus y, where x, y are 1-D
        sequences of the same length, N. If C is None (the default), this is a
        histogram of the number of occurences of the observations at
        (x[i],y[i]).

        If C is specified, it specifies values at the coordinate (x[i],y[i]).
        These values are accumulated for each hexagonal bin and then reduced
        according to reduce_C_function, which defaults to the numpy mean function
        (np.mean). (If C is specified, it must also be a 1-D sequence of the
        same length as x and y.)

        x, y and/or C may be masked arrays, in which case only unmasked points
        will be plotted.

        (see matplotlib.pyplot.hexbin documentation).

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.pyplot.hexbin
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            # make x,y masked arrays
            # (masked where data is outside of projection limb)
            x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
            y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
            ret = ax.hexbin(x,y,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        return ret

    @_transform
    def contour(self,x,y,data,*args,**kwargs):
        r"""
        Make a contour plot over the map
        (see matplotlib.pyplot.contour documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        If ``tri`` is set to ``True``, an unstructured grid is assumed
        (x,y,data must be 1-d) and matplotlib.pyplot.tricontour is used.

        Other \*args and \**kwargs passed on to matplotlib.pyplot.contour
        (or tricontour if ``tri=True``).
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            if kwargs.pop('tri', False):
                import matplotlib.tri as tri
                # for unstructured grids, toss out points outside
                # projection limb (don't use those points in triangulation).
                if ma.isMA(data):
                    data = data.filled(fill_value=1.e30)
                    masked=True
                else:
                    masked=False
                mask = np.logical_or(x<self.xmin,y<self.xmin) +\
                       np.logical_or(x>self.xmax,y>self.xmax)
                x = np.compress(mask,x)
                y = np.compress(mask,y)
                data = np.compress(mask,data)
                if masked:
                    triang = tri.Triangulation(x, y)
                    z = data[triang.triangles]
                    mask = (z > 1.e20).sum(axis=-1)
                    triang.set_mask(mask)
                    CS = ax.tricontour(triang,data,*args,**kwargs)
                else:
                    CS = ax.tricontour(x,y,data,*args,**kwargs)
            else:
                # make sure x is monotonically increasing - if not,
                # print warning suggesting that the data be shifted in longitude
                # with the shiftgrid function.
                # only do this check for global projections.
                if self.projection in _cylproj + _pseudocyl:
                    xx = x[x.shape[0]//2,:]
                    condition = (xx >= self.xmin) & (xx <= self.xmax)
                    xl = xx.compress(condition).tolist()
                    xs = xl[:]
                    xs.sort()
                    if xl != xs:
                        sys.stdout.write(" ".join([
                            "WARNING: x coordinate not montonically increasing",
                            "- contour plot may not be what you expect. If it",
                            "looks odd, you can either adjust the map projection",
                            "region to be consistent with your data, or (if your",
                            "data is on a global lat/lon grid) use the shiftdata",
                            "method to adjust the data to be consistent with the",
                            "map projection region (see doc/examples/shiftdata.py)"]))
                # mask for points more than one grid length outside projection limb.
                xx = ma.masked_where(x > 1.e20, x)
                yy = ma.masked_where(y > 1.e20, y)
                epsx = np.abs(xx[:,1:]-xx[:,0:-1]).max()
                epsy = np.abs(yy[1:,:]-yy[0:-1,:]).max()
                xymask = \
                np.logical_or(np.greater(x,self.xmax+epsx),np.greater(y,self.ymax+epsy))
                xymask = xymask + \
                np.logical_or(np.less(x,self.xmin-epsx),np.less(y,self.ymin-epsy))
                data = ma.asarray(data)
                # combine with data mask.
                mask = np.logical_or(ma.getmaskarray(data),xymask)
                data = ma.masked_array(data,mask=mask)
                CS = ax.contour(x,y,data,*args,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt and CS.get_array() is not None:
            plt.sci(CS)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        if isinstance(CS, Artist):
            # Since MPL 3.8, `QuadContourSet` objects are `Artist` objects too.
            CS, c = self._cliplimb(ax, CS)
        else:
            CS.collections, c = self._cliplimb(ax, CS.collections)
        return CS

    @_transform
    def contourf(self,x,y,data,*args,**kwargs):
        r"""
        Make a filled contour plot over the map
        (see matplotlib.pyplot.contourf documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        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 ``tri`` is set to ``True``, an unstructured grid is assumed
        (x,y,data must be 1-d) and matplotlib.pyplot.tricontourf is used.

        Other \*args and \**kwargs passed on to matplotlib.pyplot.contourf
        (or tricontourf if ``tri=True``).
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            if kwargs.get('tri', False):
                import matplotlib.tri as tri
                # for unstructured grids, toss out points outside
                # projection limb (don't use those points in triangulation).
                if ma.isMA(data):
                    data = data.filled(fill_value=1.e30)
                    masked=True
                else:
                    masked=False
                mask = np.logical_or(x<1.e20,y<1.e20)
                x = np.compress(mask,x)
                y = np.compress(mask,y)
                data = np.compress(mask,data)
                if masked:
                    triang = tri.Triangulation(x, y)
                    z = data[triang.triangles]
                    mask = (z > 1.e20).sum(axis=-1)
                    triang.set_mask(mask)
                    CS = ax.tricontourf(triang,data,*args,**kwargs)
                else:
                    CS = ax.tricontourf(x,y,data,*args,**kwargs)
            else:
                # make sure x is monotonically increasing - if not,
                # print warning suggesting that the data be shifted in longitude
                # with the shiftgrid function.
                # only do this check for global projections.
                if self.projection in _cylproj + _pseudocyl:
                    xx = x[x.shape[0]//2,:]
                    condition = (xx >= self.xmin) & (xx <= self.xmax)
                    xl = xx.compress(condition).tolist()
                    xs = xl[:]
                    xs.sort()
                    if xl != xs:
                        sys.stdout.write(" ".join([
                            "WARNING: x coordinate not montonically increasing",
                            "- contour plot may not be what you expect. If it",
                            "looks odd, you can either adjust the map projection",
                            "region to be consistent with your data, or (if your",
                            "data is on a global lat/lon grid) use the shiftgrid",
                            "function to adjust the data to be consistent with the",
                            "map projection region (see doc/examples/contour_demo.py)"]))
                # mask for points more than one grid length outside projection limb.
                xx = ma.masked_where(x > 1.e20, x)
                yy = ma.masked_where(y > 1.e20, y)
                if self.projection != 'omerc':
                    epsx = np.abs(xx[:,1:]-xx[:,0:-1]).max()
                    epsy = np.abs(yy[1:,:]-yy[0:-1,:]).max()
                else: # doesn't work for omerc (FIXME)
                    epsx = 0.; epsy = 0
                xymask = \
                np.logical_or(np.greater(x,self.xmax+epsx),np.greater(y,self.ymax+epsy))
                xymask = xymask + \
                np.logical_or(np.less(x,self.xmin-epsx),np.less(y,self.ymin-epsy))
                data = ma.asarray(data)
                # combine with data mask.
                mask = np.logical_or(ma.getmaskarray(data),xymask)
                data = ma.masked_array(data,mask=mask)
                CS = ax.contourf(x,y,data,*args,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt and CS.get_array() is not None:
            plt.sci(CS)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        if isinstance(CS, Artist):
            # Since MPL 3.8, `QuadContourSet` objects are `Artist` objects too.
            CS, c = self._cliplimb(ax, CS)
        else:
            CS.collections, c = self._cliplimb(ax, CS.collections)
        return CS

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

        Arguments may be 1-D or 2-D arrays or sequences
        (see matplotlib.pyplot.quiver documentation for details).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \*args and \**kwargs passed on to matplotlib.pyplot.quiver.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.quiver(x,y,u,v,*args,**kwargs)
        finally:
            self._restore_hold(ax)
        if plt is not None and ret.get_array() is not None:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        return ret

    @_transformuv
    def streamplot(self, x, y, u, v, *args, **kwargs):
        r"""
        Draws streamlines of a vector flow.
        (see matplotlib.pyplot.streamplot documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \*args and \**kwargs passed on to matplotlib.pyplot.streamplot.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.streamplot(x,y,u,v,*args,**kwargs)
        finally:
            self._restore_hold(ax)
        if plt is not None and ret.lines.get_array() is not None:
            plt.sci(ret.lines)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret.lines,c = self._cliplimb(ax,ret.lines)
        ret.arrows,c = self._cliplimb(ax,ret.arrows)
        # streamplot arrows not returned in matplotlib 1.1.1, so clip all
        # FancyArrow patches attached to axes instance.
        if c is not None:
            for p in ax.patches:
                if isinstance(p,FancyArrowPatch): p.set_clip_path(c)
        return ret

    @_transformuv
    def barbs(self, x, y, u, v, *args, **kwargs):
        r"""
        Make a wind barb plot (u, v) with on the map.
        (see matplotlib.pyplot.barbs documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \*args and \**kwargs passed on to matplotlib.pyplot.barbs

        Returns two matplotlib.axes.Barbs instances, one for the Northern
        Hemisphere and one for the Southern Hemisphere.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        lons, lats = self(x, y, inverse=True)
        unh = ma.masked_where(lats <= 0, u)
        vnh = ma.masked_where(lats <= 0, v)
        ush = ma.masked_where(lats > 0, u)
        vsh = ma.masked_where(lats > 0, v)
        self._save_use_hold(ax, kwargs)
        try:
            retnh =  ax.barbs(x,y,unh,vnh,*args,**kwargs)
            kwargs['flip_barb']=True
            retsh =  ax.barbs(x,y,ush,vsh,*args,**kwargs)
        finally:
            self._restore_hold(ax)
        # Because there are two collections returned in general,
        # we can't set the current image...
        #if plt is not None and ret.get_array() is not None:
        #    plt.sci(retnh)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        retnh,c = self._cliplimb(ax,retnh)
        retsh,c = self._cliplimb(ax,retsh)

        return retnh,retsh

    def drawlsmask(self,land_color="0.8",ocean_color="w",lsmask=None,
                   lsmask_lons=None,lsmask_lats=None,lakes=True,resolution='l',grid=5,**kwargs):
        r"""
        Draw land-sea mask image.

        .. note::
         The land-sea mask image cannot be overlaid on top
         of other images, due to limitations in matplotlib image handling
         (you can't specify the zorder of an image).

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keywords         Description
        ==============   ====================================================
        land_color       desired land color (color name or rgba tuple).
                         Default gray ("0.8").
        ocean_color      desired water color (color name or rgba tuple).
                         Default white.
        lsmask           An array of 0's for ocean pixels, 1's for
                         land pixels and 2's for lake/pond pixels.
                         Default is None
                         (default 5-minute resolution land-sea mask is used).
        lakes            Plot lakes and ponds (Default True)
        lsmask_lons      1d array of longitudes for lsmask (ignored
                         if lsmask is None). Longitudes must be ordered
                         from -180 W eastward.
        lsmask_lats      1d array of latitudes for lsmask (ignored
                         if lsmask is None). Latitudes must be ordered
                         from -90 S northward.
        resolution       gshhs coastline resolution used to define land/sea
                         mask (default 'l', available 'c','l','i','h' or 'f')
        grid             land/sea mask grid spacing in minutes (Default 5;
                         10, 2.5 and 1.25 are also available).
        \**kwargs        extra keyword arguments passed on to
                         :meth:`imshow`
        ==============   ====================================================

        If any of the lsmask, lsmask_lons or lsmask_lats keywords are not
        set, the built in GSHHS land-sea mask datasets are used.

        Extra keyword ``ax`` can be used to override the default axis instance.

        returns a matplotlib.image.AxesImage instance.
        """
        # convert land and water colors to integer rgba tuples with
        # values between 0 and 255.
        from matplotlib.colors import ColorConverter
        c = ColorConverter()
        # if conversion fails, assume it's because the color
        # given is already an rgba tuple with values between 0 and 255.
        try:
            cl = c.to_rgba(land_color)
            rgba_land = tuple([int(255*x) for x in cl])
        except:
            rgba_land = land_color
        try:
            co = c.to_rgba(ocean_color)
            rgba_ocean = tuple([int(255*x) for x in co])
        except:
            rgba_ocean = ocean_color
        # look for axes instance (as keyword, an instance variable
        # or from plt.gca().
        ax = kwargs.pop('ax', None) or self._check_ax()
        # Clear saved lsmask if new lsmask is passed
        if lsmask is not None or lsmask_lons is not None \
                or lsmask_lats is not None:
            # Make sure passed lsmask is not the same as cached mask
            if lsmask is not self.lsmask:
                self.lsmask = None
        # 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.
                lsmask_lons, lsmask_lats, lsmask =\
                _readlsmask(lakes=lakes,resolution=resolution,grid=grid)
                # 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).

                # is it a cylindrical projection whose limits lie
                # outside the limits of the image?
                cylproj =  self.projection in _cylproj and \
                          (self.urcrnrlon > lsmask_lons[-1] or \
                           self.llcrnrlon < lsmask_lons[0])
                if cylproj:
                    # stack grids side-by-side (in longitiudinal direction), so
                    # any range of longitudes may be plotted on a world map.
                    # in versions of NumPy later than 1.10.0, concatenate will
                    # not stack these arrays as expected. If axis 1 is outside
                    # the dimensions of the array, concatenate will now raise
                    # an IndexError. Using hstack instead.
                    lsmask_lons = \
                            np.hstack((lsmask_lons,lsmask_lons[1:] + 360))
                    lsmask = \
                            np.hstack((lsmask,lsmask[:,1:]))
        else:
            if lakes: lsmask = np.where(lsmask==2,np.array(0,np.uint8),lsmask)

        # 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 = (np.pi/180.)*(lsmask_lons[1]-lsmask_lons[0])*self.rmajor
            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,x,y = self.transform_scalar(lsmask,lsmask_lons,\
                       lsmask_lats,nx,ny,returnxy=True,order=0,masked=255)
            lsmask_lats.dtype
            # for these projections, points outside the projection
            # limb have to be set to transparent manually.
            if self.projection in _pseudocyl:
                lons, lats = self(x, y, inverse=True)
                lon_0 = self.projparams['lon_0']
                lats = lats[:,nx//2]
                lons1 = (lon_0+180.)*np.ones(lons.shape[0],np.float64)
                lons2 = (lon_0-180.)*np.ones(lons.shape[0],np.float64)
                xmax,ytmp = self(lons1,lats)
                xmin,ytmp = self(lons2,lats)
                for j in range(lats.shape[0]):
                    xx = x[j,:]
                    mask[j,:]=np.where(np.logical_or(xx<xmin[j],xx>xmax[j]),\
                                        255,mask[j,:])
            self.lsmask = mask
        ny, nx = self.lsmask.shape
        rgba = np.ones((ny,nx,4),np.uint8)
        rgba_land = np.array(rgba_land,np.uint8)
        rgba_ocean = np.array(rgba_ocean,np.uint8)
        for k in range(4):
            rgba[:,:,k] = np.where(self.lsmask,rgba_land[k],rgba_ocean[k])
        # make points outside projection limb transparent.
        rgba[:,:,3] = np.where(self.lsmask==255,0,rgba[:,:,3])
        # plot mask as rgba image.
        im = self.imshow(rgba,interpolation='nearest',ax=ax,**kwargs)
        # clip to map limbs.
        im,c = self._cliplimb(ax,im)
        return im

    def bluemarble(self, ax=None, scale=None, **kwargs):
        r"""Display blue marble image as map background.

        The original image is available at::

            https://visibleearth.nasa.gov/

        The default image size is 5400x2700, which can be a bit slow
        memory-consuming. The ``scale`` keyword allows to downsample
        the image (e.g. ``scale=0.5`` downsamples to 2700x1350).

        Parameters
        ----------

        ax : matplotlib.image.AxesImage, optional
            if given, alternative axes instance where the image is drawn

        scale : float, optional
            if given, rescale the image by the given factor

        \**kwargs : dict, optional
            keyword arguments passed on to :meth:`Basemap.imshow`.

        Returns
        -------

        ax : matplotlib.image.AxesImage
            axes instance
        """

        image = "bluemarble"
        if ax is not None:
            return self.warpimage(image=image, ax=ax, scale=scale, **kwargs)
        return self.warpimage(image=image, scale=scale, **kwargs)

    def shadedrelief(self, ax=None, scale=None, **kwargs):
        r"""Display shaded relief image as map background.

        The original image is available at::

            https://www.shadedrelief.com/

        The default image size is 5400x2700, which can be a bit slow
        memory-consuming. The ``scale`` keyword allows to downsample
        the image (e.g. ``scale=0.5`` downsamples to 2700x1350).

        Parameters
        ----------

        ax : matplotlib.image.AxesImage, optional
            if given, alternative axes instance where the image is drawn

        scale : float, optional
            if given, rescale the image by the given factor

        \**kwargs : dict, optional
            keyword arguments passed on to :meth:`Basemap.imshow`.

        Returns
        -------

        ax : matplotlib.image.AxesImage
            axes instance
        """

        image = "shadedrelief"
        if ax is not None:
            return self.warpimage(image=image, ax=ax, scale=scale, **kwargs)
        return self.warpimage(image=image, scale=scale, **kwargs)

    def etopo(self, ax=None, scale=None, **kwargs):
        r"""Display ETOPO relief image as map background.

        The original image is available at::

            http://www.ngdc.noaa.gov/mgg/global/global.html

        The default image size is 5400x2700, which can be a bit slow
        memory-consuming. The ``scale`` keyword allows to downsample
        the image (e.g. ``scale=0.5`` downsamples to 2700x1350).

        Parameters
        ----------

        ax : matplotlib.image.AxesImage, optional
            if given, alternative axes instance where the image is drawn

        scale : float, optional
            if given, rescale the image by the given factor

        \**kwargs : dict, optional
            keyword arguments passed on to :meth:`Basemap.imshow`.

        Returns
        -------

        ax : matplotlib.image.AxesImage
            axes instance
        """

        image = "etopo"
        if ax is not None:
            return self.warpimage(image=image, ax=ax, scale=scale, **kwargs)
        return self.warpimage(image=image, scale=scale, **kwargs)

    def warpimage(self,image="bluemarble",scale=None,**kwargs):
        r"""
        Display an image (filename given by ``image`` keyword) as a map background.
        If image is a URL (starts with 'http'), it is downloaded to a temp
        file using urllib.urlretrieve.

        Default (if ``image`` not specified) is to display
        'blue marble next generation' image from http://visibleearth.nasa.gov/.

        Specified image must have pixels covering the whole globe in a regular
        lat/lon grid, starting and -180W and the South Pole.
        Works with the global images from
        http://earthobservatory.nasa.gov/Features/BlueMarble/BlueMarble_monthlies.php.

        The ``scale`` keyword can be used to downsample (rescale) the image.
        Values less than 1.0 will speed things up at the expense of image
        resolution.

        Extra keyword ``ax`` can be used to override the default axis instance.

        \**kwargs passed on to :meth:`imshow`.

        returns a matplotlib.image.AxesImage instance.
        """

        # fix PIL import on some versions of OSX and scipy
        try:
            from PIL import Image
        except ImportError:
            try:
                import Image
            except ImportError:
                raise ImportError("warpimage method requires PIL "
                                  "(http://pillow.readthedocs.io)")

        import matplotlib.image as mpimg

        def pil_to_array(*args, **kwargs):
            """Call :func:`~mpimg.pil_to_array` ignoring deprecation warnings."""

            with warnings.catch_warnings():
                warnings.simplefilter("ignore", category=DeprecationWarning)
                return mpimg.pil_to_array(*args, **kwargs)

        if self.celestial:
            raise ValueError("warpimage does not work in celestial coordinates")
        ax = kwargs.pop('ax', None) or self._check_ax()
        # default image file is blue marble next generation
        # from NASA (http://visibleearth.nasa.gov).
        if image == "bluemarble":
            file = os.path.join(basemap_datadir,'bmng.jpg')
        # display shaded relief image (from
        # http://www.shadedreliefdata.com)
        elif image == "shadedrelief":
            file = os.path.join(basemap_datadir,'shadedrelief.jpg')
        # display etopo image (from
        # http://www.ngdc.noaa.gov/mgg/image/globalimages.html)
        elif image == "etopo":
            file = os.path.join(basemap_datadir,'etopo1.jpg')
        else:
            file = image
        # if image is same as previous invocation, used cached data.
        # if not, regenerate rgba data.
        if not hasattr(self,'_bm_file') or self._bm_file != file:
            newfile = True
        else:
            newfile = False
        if file.startswith('http'):
            self._bm_file, headers = urlretrieve(file)
        else:
            self._bm_file = file
        # bmproj is True if map projection region is same as
        # image region.
        bmproj = self.projection == 'cyl' and \
                 self.llcrnrlon == -180 and self.urcrnrlon == 180 and \
                 self.llcrnrlat == -90 and self.urcrnrlat == 90
        # read in jpeg image to rgba array of normalized floats.
        if not hasattr(self,'_bm_rgba') or newfile:
            pilImage = Image.open(self._bm_file)
            if scale is not None:
                w, h = pilImage.size
                width = int(np.round(w*scale))
                height = int(np.round(h*scale))
                pilImage = pilImage.resize((width,height),Image.LANCZOS)
            # orientation of arrays returned by pil_to_array changed in
            # matplotlib 1.2 (https://github.com/matplotlib/matplotlib/pull/616)
            self._bm_rgba = pil_to_array(pilImage)[::-1,:]
            # define lat/lon grid that image spans.
            nlons = self._bm_rgba.shape[1]; nlats = self._bm_rgba.shape[0]
            delta = 360./float(nlons)
            self._bm_lons = np.arange(-180.+0.5*delta,180.,delta)
            self._bm_lats = np.arange(-90.+0.5*delta,90.,delta)
            # is it a cylindrical projection whose limits lie
            # outside the limits of the image?
            cylproj =  self.projection in _cylproj and \
                      (self.urcrnrlon > self._bm_lons[-1] or \
                       self.llcrnrlon < self._bm_lons[0])
            # if pil_to_array returns a 2D array, it's a grayscale image.
            # create an RGB image, with R==G==B.
            if self._bm_rgba.ndim == 2:
                tmp = np.empty(self._bm_rgba.shape+(3,),np.uint8)
                for k in range(3):
                    tmp[:,:,k] = self._bm_rgba
                self._bm_rgba = tmp
            if cylproj and not bmproj:
                # stack grids side-by-side (in longitiudinal direction), so
                # any range of longitudes may be plotted on a world map.
                self._bm_lons = \
                np.concatenate((self._bm_lons,self._bm_lons+360),0)
                self._bm_rgba = \
                np.concatenate((self._bm_rgba,self._bm_rgba),1)
            # convert to normalized floats.
            self._bm_rgba = self._bm_rgba.astype(np.float32)/255.
        if not bmproj: # interpolation necessary.
            if newfile or not hasattr(self,'_bm_rgba_warped'):
                # transform to nx x ny regularly spaced native
                # projection grid.
                # nx and ny chosen to have roughly the
                # same horizontal res as original image.
                if self.projection != 'cyl':
                    dx = 2.*np.pi*self.rmajor/float(nlons)
                    nx = int((self.xmax-self.xmin)/dx)+1
                    ny = int((self.ymax-self.ymin)/dx)+1
                else:
                    dx = 360./float(nlons)
                    nx = int((self.urcrnrlon-self.llcrnrlon)/dx)+1
                    ny = int((self.urcrnrlat-self.llcrnrlat)/dx)+1
                self._bm_rgba_warped = np.ones((ny,nx,4),np.float64)
                # interpolate rgba values from geographic coords (proj='cyl')
                # to map projection coords.
                # if masked=True, values outside of
                # projection limb will be masked.
                for k in range(self._bm_rgba.shape[2]):
                    self._bm_rgba_warped[:,:,k],x,y = \
                    self.transform_scalar(self._bm_rgba[:,:,k],\
                    self._bm_lons,self._bm_lats,nx,ny,returnxy=True)
                # for ortho,geos mask pixels outside projection limb.
                if self.projection in ['geos','ortho','nsper'] or \
                   (self.projection == 'aeqd' and self._fulldisk):
                    lonsr,latsr = self(x,y,inverse=True)
                    mask = ma.zeros((ny,nx,4),np.int8)
                    mask[:,:,0] = np.logical_or(lonsr>1.e20,latsr>1.e30)
                    for k in range(1,4):
                        mask[:,:,k] = mask[:,:,0]
                    self._bm_rgba_warped = \
                    ma.masked_array(self._bm_rgba_warped,mask=mask)
                    # make points outside projection limb transparent.
                    self._bm_rgba_warped = self._bm_rgba_warped.filled(0.)
                # treat pseudo-cyl projections such as mollweide, robinson and sinusoidal.
                elif self.projection in _pseudocyl and \
                     self.projection != 'hammer':
                    lonsr,latsr = self(x,y,inverse=True)
                    mask = ma.zeros((ny,nx,4),np.int8)
                    lon_0 = self.projparams['lon_0']
                    lonright = lon_0 + 180. - 1E-10
                    lonleft = lon_0 - 180. + 1E-10
                    x1 = np.array(ny*[0.5*(self.xmax + self.xmin)],np.float64)
                    y1 = np.linspace(self.ymin, self.ymax, ny)
                    lons1, lats1 = self(x1,y1,inverse=True)
                    lats1 = np.where(lats1 < -89.999999, -89.999999, lats1)
                    lats1 = np.where(lats1 > 89.999999, 89.999999, lats1)
                    for j,lat in enumerate(lats1):
                        xmax,ymax = self(lonright,lat)
                        xmin,ymin = self(lonleft,lat)
                        mask[j,:,0] = np.logical_or(x[j,:]>xmax,x[j,:]<xmin)
                    for k in range(1,4):
                        mask[:,:,k] = mask[:,:,0]
                    self._bm_rgba_warped = \
                    ma.masked_array(self._bm_rgba_warped,mask=mask)
                    # make points outside projection limb transparent.
                    # FIXME: Probably not needed anymore
                    self._bm_rgba_warped = self._bm_rgba_warped.filled(0.)
            # plot warped rgba image.
            im = self.imshow(self._bm_rgba_warped,ax=ax,**kwargs)
            # for hammer projection, use clip path defined by
            # projection limb (patch created in drawmapboundary).
            # FIXME: Is this now redundant?
            if self.projection == 'hammer':
                if not self._mapboundarydrawn:
                    self.drawmapboundary(color='none',linewidth=None)
                im.set_clip_path(self._mapboundarydrawn)
        else:
            # bmproj True, no interpolation necessary.
            im = self.imshow(self._bm_rgba,ax=ax,**kwargs)
        # clip to map limbs
        im,c = self._cliplimb(ax,im)
        return im

    def arcgisimage(self, server="http://server.arcgisonline.com/ArcGIS",
                    service="World_Imagery", xpixels=400, ypixels=None,
                    dpi=96, cachedir=None, verbose=False, **kwargs):
        r"""Display background image using ArcGIS Server REST API.

        In order to use this method, the :class:`Basemap` instance
        must be created using the ``epsg`` keyword to define the
        map projection, unless the "cyl" projection is used (in
        which case the EPSG code 4326 is assumed).

        Parameters
        ----------

        server : str, optional
            base URL of the web map server

        service : str, optional
            service (image type) hosted by the server

        xpixels : int, optional
            requested number of image pixels in the `x`-direction

        ypixels : int, optional
            requested number of image pixels in the `y`-direction;
            if not given, it is inferred from ``xpixels`` and the
            aspect ratio of the map projection region

        dpi : int, optional
            device resolution of the exported image

        cachedir : str, optional
            if given, directory to use as cache folder for the images
            retrieved from the server

        verbose : bool, optional
            if True, print debugging information

        \**kwargs : dict, optional
            keyword-only arguments; currently, only ``ax`` is supported
            to override the default :class:`matplotlib.axes.Axes`
            instance

        Returns
        -------

        aximg : matplotlib.image.AxesImage
            image axes instance
        """

        # Fix PIL import on some versions of OSX and scipy.
        try:
            from PIL import Image
        except ImportError:
            try:
                import Image
            except ImportError:
                raise ImportError("arcgisimage method requires PIL "
                                  "(http://pillow.readthedocs.io)")

        if not hasattr(self, "epsg"):
            raise ValueError("the Basemap instance must be created using "
                             "an EPSG code (http://spatialreference.org) "
                             "in order to use the wmsmap method")

        ax = kwargs.pop("ax", None) or self._check_ax()

        # Find the `(x, y)` values at the corner points.
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=FutureWarning)
            p = pyproj.Proj(init="epsg:%s" % self.epsg, preserve_units=True)
        xmin, ymin = p(self.llcrnrlon, self.llcrnrlat)
        xmax, ymax = p(self.urcrnrlon, self.urcrnrlat)
        if self.projection in _cylproj:
            dateline = _geoslib.Point(self(180., 0.5 * (self.llcrnrlat + self.urcrnrlat)))
            if dateline.within(self._boundarypolyxy):
                raise ValueError("arcgisimage cannot handle images that cross "
                                 "the dateline for cylindrical projections")

        # If ypixels is not given, compute it with xpixels and aspect ratio.
        if ypixels is None:
            ypixels = int(self.aspect * xpixels)

        # Construct a URL using the ArcGIS Server REST API.
        basemap_url = "".join([
            "%s/rest/services/%s/MapServer/export?",
            "bbox=%s,%s,%s,%s&",
            "bboxSR=%s&",
            "imageSR=%s&",
            "size=%s,%s&",
            "dpi=%s&",
            "format=png32&",
            "transparent=true&",
            "f=image",
        ]) % (server, service, xmin, ymin, xmax, ymax, self.epsg, self.epsg, xpixels, ypixels, dpi)

        # Print URL in verbose mode.
        if verbose:  # pragma: no cover
            print(basemap_url)

        # Try to return fast if cache is enabled.
        if cachedir is not None:
            # Generate a filename for the cached file.
            filename = "%s-bbox-%s-%s-%s-%s-bboxsr%s-imagesr%s-size-%s-%s-dpi%s.png" % \
                (service, xmin, ymin, xmax, ymax, self.epsg, self.epsg, xpixels, ypixels, dpi)
            # Return fast if the image is already in the cache.
            cache_path = os.path.join(cachedir, filename)
            if os.path.isfile(cache_path):
                if verbose:  # pragma: no cover
                    print("Image already in cache")
                img = Image.open(cache_path)
                return self.imshow(img, ax=ax, origin="upper")

        # Retrieve image from the remote server.
        import contextlib
        conn = urlopen(basemap_url)
        with contextlib.closing(conn):
            img = Image.open(conn)
            # Save to cache if requested.
            if cachedir is not None:
                # Check if the cache directory exists, if not create it.
                if not os.path.exists(cachedir):
                    os.makedirs(cachedir)
                img.save(cache_path)
        return self.imshow(img, ax=ax, origin="upper")

    def wmsimage(self,server,\
                 xpixels=400,ypixels=None,\
                 format='png',alpha=None,verbose=False,**kwargs):
        r"""
        Retrieve an image using from a WMS server using the
        Open Geospatial Consortium (OGC) standard interface
        and display on the map. Requires OWSLib
        (http://pypi.python.org/pypi/OWSLib).
        In order to use this method, the Basemap instance must be
        created using the ``epsg`` keyword to define the map projection, unless
        the ``cyl`` projection is used (in which case the epsg code 4326 is
        assumed).

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keywords         Description
        ==============   ====================================================
        server           WMS server URL.
        xpixels          requested number of image pixels in x-direction
                         (default 400).
        ypixels          requested number of image pixels in y-direction.
                         Default (None) is to infer the number from
                         from xpixels and the aspect ratio of the
                         map projection region.
        format           desired image format (default 'png')
        alpha            The alpha blending value,
                         between 0 (transparent) and 1 (opaque) (default None)
        verbose          if True, print WMS server info (default
                         False).
        \**kwargs        extra keyword arguments passed on to
                         OWSLib.wms.WebMapService and
                         OWSLib.wms.WebMapService.getmap.
        ==============   ====================================================

        Extra keyword ``ax`` can be used to override the default axis instance.

        returns a matplotlib.image.AxesImage instance.
        """
        try:
            from owslib.wms import WebMapService
        except ImportError:
            raise ImportError('OWSLib required to use wmsimage method')
        import io
        ax = kwargs.pop('ax', None) or self._check_ax()
        if not hasattr(self,'epsg'):
            raise ValueError("the Basemap instance must be created using "
                             "an EPSG code (http://spatialreference.org) "
                             "in order to use the wmsmap method")
        if 'layers' not in kwargs:
            raise ValueError('no layers specified')
        # find the x,y values at the corner points.
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=FutureWarning)
            p = pyproj.Proj(init="epsg:%s" % self.epsg, preserve_units=True)
        xmin,ymin = p(self.llcrnrlon,self.llcrnrlat)
        xmax,ymax = p(self.urcrnrlon,self.urcrnrlat)
        if self.projection in _cylproj:
            Dateline =\
            _geoslib.Point(self(180.,0.5*(self.llcrnrlat+self.urcrnrlat)))
            hasDateline = Dateline.within(self._boundarypolyxy)
            if hasDateline:
                raise ValueError("wmsimage cannot handle images that cross "
                                 "the dateline for cylindrical projections")
        if self.projection == 'cyl':
            xmin = (180./np.pi)*xmin; xmax = (180./np.pi)*xmax
            ymin = (180./np.pi)*ymin; ymax = (180./np.pi)*ymax
        # ypixels not given, find by scaling xpixels by the map aspect ratio.
        if ypixels is None:
            ypixels = int(self.aspect*xpixels)
        if verbose: print(server)
        wms_keys = ["version", "xml", "username", "password",
                    "parse_remote_metadata", "timeout", "headers", "auth"]
        wms_options = {k: kwargs[k] for k in wms_keys if k in kwargs}
        kwargs = {k: kwargs[k] for k in kwargs if k not in wms_keys}
        wms = WebMapService(server, **wms_options)
        if verbose:
            print('id: %s, version: %s' %
            (wms.identification.type,wms.identification.version))
            print('title: %s, abstract: %s' %
            (wms.identification.title,wms.identification.abstract))
            print('available layers:')
            print(list(wms.contents))
            print('projection options:')
            print(wms[kwargs['layers'][0]].crsOptions)
        # remove keys from kwargs that are over-ridden
        for k in ['format','bbox','service','size','srs']:
            if 'format' in kwargs: del kwargs['format']
        img = wms.getmap(service='wms',bbox=(xmin,ymin,xmax,ymax),
                         size=(xpixels,ypixels),format='image/%s'%format,
                         srs='EPSG:%s' % self.epsg, **kwargs)
        # return AxesImage instance.
        # this works for png and jpeg.
        return self.imshow(imread(io.BytesIO(img.read()),
                           format=format),origin='upper',alpha=alpha,ax=ax)

    def drawmapscale(self,lon,lat,lon0,lat0,length,barstyle='simple',\
                     units='km',fontsize=9,yoffset=None,labelstyle='simple',\
                     fontcolor='k',fillcolor1='w',fillcolor2='k',ax=None,\
                     format='%d',zorder=None,linecolor=None,linewidth=None):
        """
        Draw a map scale at ``lon,lat`` of length ``length``
        representing distance in the map
        projection coordinates at ``lon0,lat0``.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keywords         Description
        ==============   ====================================================
        units            the units of the length argument (Default km).
        barstyle         ``simple`` or ``fancy`` (roughly corresponding
                         to the styles provided by Generic Mapping Tools).
                         Default ``simple``.
        fontsize         for map scale annotations, default 9.
        fontcolor            for map scale annotations, default black.
        labelstyle       ``simple`` (default) or ``fancy``.  For
                         ``fancy`` the map scale factor (ratio betwee
                         the actual distance and map projection distance
                         at lon0,lat0) and the value of lon0,lat0 are also
                         displayed on the top of the scale bar. For
                         ``simple``, just the units are display on top
                         and the distance below the scale bar.
                         If equal to False, plot an empty label.
        format           a string formatter to format numeric values
        yoffset          yoffset controls how tall the scale bar is,
                         and how far the annotations are offset from the
                         scale bar.  Default is 0.02 times the height of
                         the map (0.02*(self.ymax-self.ymin)).
        fillcolor1(2)    colors of the alternating filled regions
                         (default white and black).  Only relevant for
                         'fancy' barstyle.
        zorder           sets the zorder for the map scale.
        linecolor        sets the color of the scale, by default, fontcolor
                         is used
        linewidth        linewidth for scale and ticks
        ==============   ====================================================

        Extra keyword ``ax`` can be used to override the default axis instance.
        """
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # not valid for cylindrical projection
        if self.projection == 'cyl':
            raise ValueError("cannot draw map scale for projection='cyl'")
        # convert length to meters
        lenlab = length
        if units == 'km':
            length = length*1000
        elif units == 'mi':
            length = length*1609.344
        elif units == 'nmi':
            length = length*1852
        elif units == 'ft':
            length = length*0.3048
        elif units != 'm':
            raise KeyError("units must be 'm' (meters), 'km' (kilometers), "
                           "'mi' (miles), 'nmi' (nautical miles), or 'ft' (feet)")
        # reference point and center of scale.
        x0,y0 = self(lon0,lat0)
        xc,yc = self(lon,lat)
        # make sure lon_0 between -180 and 180
        lon_0 = ((lon0+360) % 360) - 360
        degchar = b"\xc2\xb0".decode("utf-8")
        if lat0>0:
            if lon>0:
                lonlatstr = '%g%sN, %g%sE' % (lat0, degchar, lon_0, degchar)
            elif lon<0:
                lonlatstr = '%g%sN, %g%sW' % (lat0, degchar, lon_0, degchar)
            else:
                lonlatstr = '%g%s, %g%sW' % (lat0, degchar, lon_0, degchar)
        else:
            if lon>0:
                lonlatstr = '%g%sS, %g%sE' % (lat0, degchar, lon_0, degchar)
            elif lon<0:
                lonlatstr = '%g%sS, %g%sW' % (lat0, degchar, lon_0, degchar)
            else:
                lonlatstr = '%g%sS, %g%s' % (lat0, degchar, lon_0, degchar)
        # left edge of scale
        lon1,lat1 = self(x0-length/2,y0,inverse=True)
        x1,y1 = self(lon1,lat1)
        # right edge of scale
        lon4,lat4 = self(x0+length/2,y0,inverse=True)
        x4,y4 = self(lon4,lat4)
        x1 = x1-x0+xc; y1 = y1-y0+yc
        x4 = x4-x0+xc; y4 = y4-y0+yc
        if x1 > 1.e20 or x4 > 1.e20 or y1 > 1.e20 or y4 > 1.e20:
            raise ValueError("scale bar positioned outside projection limb")
        # scale factor for true distance
        gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
        az12,az21,dist = gc.inv(lon1,lat1,lon4,lat4)
        scalefact = dist/length
        # label to put on top of scale bar.
        if labelstyle=='simple':
            labelstr = units
        elif labelstyle == 'fancy':
            labelstr = units+" (scale factor %4.2f at %s)"%(scalefact,lonlatstr)
        elif labelstyle == False:
            labelstr = ''
        else:
            raise KeyError("labelstyle must be 'simple' or 'fancy'")
        # default y offset is 2 percent of map height.
        if yoffset is None: yoffset = 0.02*(self.ymax-self.ymin)
        rets = [] # will hold all plot objects generated.
        # set linecolor
        if linecolor is None:
            linecolor = fontcolor
        # 'fancy' style
        if barstyle == 'fancy':
            #we need 5 sets of x coordinates (in map units)
            #quarter scale
            lon2,lat2 = self(x0-length/4,y0,inverse=True)
            x2,y2 = self(lon2,lat2)
            x2 = x2-x0+xc; y2 = y2-y0+yc
            #three quarter scale
            lon3,lat3 = self(x0+length/4,y0,inverse=True)
            x3,y3 = self(lon3,lat3)
            x3 = x3-x0+xc; y3 = y3-y0+yc
            #plot top line
            ytop = yc+yoffset/2
            ybottom = yc-yoffset/2
            ytick = ybottom - yoffset/2
            ytext = ytick - yoffset/2
            rets.append(self.plot([x1,x4],[ytop,ytop],color=linecolor, linewidth=linewidth)[0])
            #plot bottom line
            rets.append(self.plot([x1,x4],[ybottom,ybottom],color=linecolor, linewidth=linewidth)[0])
            #plot left edge
            rets.append(self.plot([x1,x1],[ybottom,ytop],color=linecolor, linewidth=linewidth)[0])
            #plot right edge
            rets.append(self.plot([x4,x4],[ybottom,ytop],color=linecolor, linewidth=linewidth)[0])
            #make a filled black box from left edge to 1/4 way across
            rets.append(ax.fill([x1,x2,x2,x1,x1],[ytop,ytop,ybottom,ybottom,ytop],\
                        ec=fontcolor,fc=fillcolor1)[0])
            #make a filled white box from 1/4 way across to 1/2 way across
            rets.append(ax.fill([x2,xc,xc,x2,x2],[ytop,ytop,ybottom,ybottom,ytop],\
                        ec=fontcolor,fc=fillcolor2)[0])
            #make a filled white box from 1/2 way across to 3/4 way across
            rets.append(ax.fill([xc,x3,x3,xc,xc],[ytop,ytop,ybottom,ybottom,ytop],\
                        ec=fontcolor,fc=fillcolor1)[0])
            #make a filled white box from 3/4 way across to end
            rets.append(ax.fill([x3,x4,x4,x3,x3],[ytop,ytop,ybottom,ybottom,ytop],\
                        ec=fontcolor,fc=fillcolor2)[0])
            #plot 3 tick marks at left edge, center, and right edge
            rets.append(self.plot([x1,x1],[ytick,ybottom],color=linecolor, linewidth=linewidth)[0])
            rets.append(self.plot([xc,xc],[ytick,ybottom],color=linecolor, linewidth=linewidth)[0])
            rets.append(self.plot([x4,x4],[ytick,ybottom],color=linecolor, linewidth=linewidth)[0])
            #label 3 tick marks
            rets.append(ax.text(x1,ytext,format % (0),\
            horizontalalignment='center',\
            verticalalignment='top',\
            fontsize=fontsize,color=fontcolor))
            rets.append(ax.text(xc,ytext,format % (0.5*lenlab),\
            horizontalalignment='center',\
            verticalalignment='top',\
            fontsize=fontsize,color=fontcolor))
            rets.append(ax.text(x4,ytext,format % (lenlab),\
            horizontalalignment='center',\
            verticalalignment='top',\
            fontsize=fontsize,color=fontcolor))
            #put units, scale factor on top
            rets.append(ax.text(xc,ytop+yoffset/2,labelstr,\
            horizontalalignment='center',\
            verticalalignment='bottom',\
            fontsize=fontsize,color=fontcolor))
        # 'simple' style
        elif barstyle == 'simple':
            rets.append(self.plot([x1,x4],[yc,yc],color=linecolor, linewidth=linewidth)[0])
            rets.append(self.plot([x1,x1],[yc-yoffset,yc+yoffset],color=linecolor, linewidth=linewidth)[0])
            rets.append(self.plot([x4,x4],[yc-yoffset,yc+yoffset],color=linecolor, linewidth=linewidth)[0])
            rets.append(ax.text(xc,yc-yoffset,format % lenlab,\
            verticalalignment='top',horizontalalignment='center',\
            fontsize=fontsize,color=fontcolor))
            #put units, scale factor on top
            rets.append(ax.text(xc,yc+yoffset,labelstr,\
            horizontalalignment='center',\
            verticalalignment='bottom',\
            fontsize=fontsize,color=fontcolor))
        else:
            raise KeyError("barstyle must be 'simple' or 'fancy'")
        if zorder is not None:
            for ret in rets:
                try:
                    ret.set_zorder(zorder)
                except:
                    pass
        return rets

    def colorbar(self,mappable=None,location='right',size="5%",pad='2%',fig=None,ax=None,**kwargs):
        r"""
        Add colorbar to axes associated with a map.
        The colorbar axes instance is created using the axes_grid toolkit.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keywords         Description
        ==============   ====================================================
        mappable         the Image, ContourSet, etc. to which the colorbar
                         applies.  Default None, matplotlib.pyplot.gci() is
                         used to retrieve the current image mappable.
        location         where to put colorbar ('top','bottom','left','right')
                         Default 'right'.
        size             width of colorbar axes (string 'N%', where N is
                         an integer describing the fractional width of the parent
                         axes). Default '5%'.
        pad              Padding between parent axes and colorbar axes in
                         same units as size. Default '2%'.
        fig              Figure instance the map axes instance is associated
                         with. Default None, and matplotlib.pyplot.gcf() is used
                         to retrieve the current active figure instance.
        ax               The axes instance which the colorbar will be
                         associated with.  Default None, searches for self.ax,
                         and if None uses matplotlib.pyplot.gca().
        \**kwargs        extra keyword arguments passed on to
                         colorbar method of the figure instance.
        ==============   ====================================================

        Returns a matplotlib colorbar instance.
        """
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # get current figure instance (if none specified).
        if fig is None or mappable is None:
            import matplotlib.pyplot as plt
        if fig is None:
            fig = plt.gcf()
        # get current mappable if none specified.
        if mappable is None:
            mappable = plt.gci()
        # create colorbar axes uses axes_grid toolkit.
        divider = make_axes_locatable(ax)
        if location in ['left','right']:
            orientation = 'vertical'
        elif location in ['top','bottom']:
            orientation = 'horizontal'
        else:
            raise ValueError('location must be top,bottom,left or right')
        cax = divider.append_axes(location, size=size, pad=pad)
        # create colorbar.
        cb = fig.colorbar(mappable,orientation=orientation,cax=cax,**kwargs)
        fig.sca(ax) # reset parent axes as current axes.
        return cb

    def nightshade(self,date,color="k",delta=0.25,alpha=0.5,ax=None,zorder=2):
        """
        Shade the regions of the map that are in darkness at the time
        specifed by ``date``.  ``date`` is a datetime instance,
        assumed to be UTC.

        .. tabularcolumns:: |l|L|

        ==============   ====================================================
        Keywords         Description
        ==============   ====================================================
        color            color to shade night regions (default black).
        delta            day/night terminator is computed with a
                         a resolution of ``delta`` degrees (default 0.25).
        alpha            alpha transparency for shading (default 0.5, so
                         map background shows through).
        zorder           zorder for shading (default 2).
        ==============   ====================================================

        Extra keyword ``ax`` can be used to override the default axis instance.

        returns a matplotlib.contour.ContourSet instance.
        """
        from .solar import daynight_grid
        # make sure date is UTC, or naive with repect to time zones
        if date.utcoffset():
            raise ValueError('datetime instance must be UTC, not {0}'.format(date.tzname()))
        # get current axes instance (if none specified).
        ax = ax or self._check_ax()
        # create grid of day=0, night=1
        lons,lats,daynight = daynight_grid(date,delta,self.lonmin,self.lonmax)
        x,y = self(lons,lats)
        # contour the day-night grid, coloring the night area
        # with the specified color and transparency.
        CS = self.contourf(x,y,daynight,1,colors=[color],alpha=alpha,ax=ax)
        if isinstance(CS, Artist):
            # Since MPL 3.8, `QuadContourSet` objects are `Artist` objects too.
            CS.set_zorder(zorder)
            # clip to map limbs
            CS, c = self._cliplimb(ax, CS)
        else:
            # set zorder on ContourSet collections show night shading
            # is on top.
            for c in CS.collections:
                c.set_zorder(zorder)
            # clip to map limbs
            CS.collections, c = self._cliplimb(ax, CS.collections)
        return CS

    def _check_ax(self):
        """
        Returns the axis on which to draw.
        Returns self.ax, or if self.ax=None returns plt.gca().
        """
        if self.ax is None:
            try:
                ax = plt.gca()
            except:
                import matplotlib.pyplot as plt
                ax = plt.gca()
            # associate an axes instance with this Basemap instance
            # the first time this method is called.
            #self.ax = ax
        else:
            ax = self.ax
        return ax

    def _ax_plt_from_kw(self, kw):
        """
        Return (ax, plt), where ax is the current axes, and plt is
        None or a reference to the pyplot module.

        plt will be None if ax was popped from kw or taken from self.ax;
        otherwise, pyplot was used and is returned.
        """
        plt = None
        _ax = kw.pop('ax', None)
        if _ax is None:
            _ax = self.ax
            if _ax is None:
                import matplotlib.pyplot as plt
                _ax = plt.gca()
        return _ax, plt

    def shiftdata(self,lonsin,datain=None,lon_0=None,fix_wrap_around=True):
        """
        Shift longitudes (and optionally data) so that they match map projection region.
        Only valid for cylindrical/pseudo-cylindrical global projections and data
        on regular lat/lon grids. longitudes and data can be 1-d or 2-d, if 2-d
        it is assumed longitudes are 2nd (rightmost) dimension.

        .. tabularcolumns:: |l|L|

        ================ ======================================================
        Arguments        Description
        ================ ======================================================
        lonsin           original 1-d or 2-d longitudes.
        ================ ======================================================

        .. tabularcolumns:: |l|L|

        ================ ======================================================
        Keywords         Description
        ================ ======================================================
        datain           original 1-d or 2-d data. Default None.
        lon_0            center of map projection region. Defaut None,
                         given by current map projection.
        fix_wrap_around  if True reindex (if required) longitudes (and data) to
                         avoid jumps caused by remapping of longitudes of
                         points from outside of the [lon_0-180, lon_0+180]
                         interval back into the interval.
                         If False do not reindex longitudes and data, but do
                         make sure that longitudes are in the
                         [lon_0-180, lon_0+180] range.
        ================ ======================================================

        If datain is given, returns ``lonsout, dataout`` (longitudes and data
        shifted to fit in the interval [lon_0-180, lon_0+180]); otherwise,
        returns just the shifted longitudes. If transformed longitudes lie
        outside the map projection region, data is masked and longitudes are
        set to 1.e30.
        """
        if lon_0 is None and 'lon_0' not in self.projparams:
            raise ValueError('lon_0 keyword must be provided')
        lonsin = np.asarray(lonsin)
        if lonsin.ndim not in [1,2]:
            raise ValueError('1-d or 2-d longitudes required')
        if datain is not None:
            # if it's a masked array, leave it alone.
            if not ma.isMA(datain): datain = np.asarray(datain)
            if datain.ndim not in [1,2]:
                raise ValueError('1-d or 2-d data required')
        if lon_0 is None:
            lon_0 = self.projparams['lon_0']
        # 2-d data.
        if lonsin.ndim == 2:
            nlats = lonsin.shape[0]
            nlons = lonsin.shape[1]
            lonsin1 = lonsin[0,:]
            lonsin1 = np.where(lonsin1 > lon_0+180, lonsin1-360 ,lonsin1)
            lonsin1 = np.where(lonsin1 < lon_0-180, lonsin1+360 ,lonsin1)
            if nlons > 1:
                londiff = np.abs(lonsin1[0:-1]-lonsin1[1:])
                londiff_sort = np.sort(londiff)
                thresh = 360.-londiff_sort[-2] if nlons > 2 else 360.-londiff_sort[-1]
                itemindex = nlons-np.where(londiff>=thresh)[0]
            else:
                lonsin[0, :] = lonsin1
                itemindex = 0

            # if no shift necessary, itemindex will be
            # empty, so don't do anything
            if fix_wrap_around and np.all(itemindex != 0) and itemindex.size:
                # check to see if cyclic (wraparound) point included
                # if so, remove it.
                if np.abs(lonsin1[0]-lonsin1[-1]) < 1.e-4:
                    hascyclic = True
                    lonsin_save = lonsin.copy()
                    lonsin = lonsin[:,1:]
                    if datain is not None:
                       datain_save = datain.copy()
                       datain = datain[:,1:]
                else:
                    hascyclic = False
                lonsin = np.where(lonsin > lon_0+180, lonsin-360 ,lonsin)
                lonsin = np.where(lonsin < lon_0-180, lonsin+360 ,lonsin)
                lonsin = np.roll(lonsin,itemindex-1,axis=1)
                if datain is not None:
                    # np.roll works on ndarrays and on masked arrays
                    datain = np.roll(datain,itemindex-1,axis=1)
                # add cyclic point back at beginning.
                if hascyclic:
                    lonsin_save[:,1:] = lonsin
                    lonsin_save[:,0] = lonsin[:,-1]-360.
                    lonsin = lonsin_save
                    if datain is not None:
                        datain_save[:,1:] = datain
                        datain_save[:,0] = datain[:,-1]
                        datain = datain_save

        # 1-d data.
        elif lonsin.ndim == 1:
            nlons = len(lonsin)
            lonsin = np.where(lonsin > lon_0+180, lonsin-360 ,lonsin)
            lonsin = np.where(lonsin < lon_0-180, lonsin+360 ,lonsin)

            if nlons > 1:
                londiff = np.abs(lonsin[0:-1]-lonsin[1:])
                londiff_sort = np.sort(londiff)
                thresh = 360.-londiff_sort[-2] if nlons > 2 else 360.0 - londiff_sort[-1]
                itemindex = len(lonsin)-np.where(londiff>=thresh)[0]
            else:
                itemindex = 0

            if fix_wrap_around and np.all(itemindex != 0) and itemindex.size:
                # check to see if cyclic (wraparound) point included
                # if so, remove it.
                if np.abs(lonsin[0]-lonsin[-1]) < 1.e-4:
                    hascyclic = True
                    lonsin_save = lonsin.copy()
                    lonsin = lonsin[1:]
                    if datain is not None:
                        datain_save = datain.copy()
                        datain = datain[1:]
                else:
                    hascyclic = False
                lonsin = np.roll(lonsin,itemindex-1)
                if datain is not None:
                    datain = np.roll(datain,itemindex-1)
                # add cyclic point back at beginning.
                if hascyclic:
                    lonsin_save[1:] = lonsin
                    lonsin_save[0] = lonsin[-1]-360.
                    lonsin = lonsin_save
                    if datain is not None:
                        datain_save[1:] = datain
                        datain_save[0] = datain[-1]
                        datain = datain_save

        # mask points outside
        # map region so they don't wrap back in the domain.
        mask = np.logical_or(lonsin<lon_0-180,lonsin>lon_0+180)
        lonsin = np.where(mask,1.e30,lonsin)
        if datain is not None and mask.any():
            datain = ma.masked_where(mask, datain)

        if datain is not None:
            return lonsin, datain
        else:
            return lonsin

### End of Basemap class

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(datain,xin,yin,xout,yout,checkbounds=False,masked=False,order=1):
    """
    Interpolate data (``datain``) on a rectilinear grid (with x = ``xin``
    y = ``yin``) to a grid with x = ``xout``, y= ``yout``.

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Arguments        Description
    ==============   ====================================================
    datain           a rank-2 array with 1st dimension corresponding to
                     y, 2nd dimension x.
    xin, yin         rank-1 arrays containing x and y of
                     datain grid in increasing order.
    xout, yout       rank-2 arrays containing x and y of desired output grid.
    ==============   ====================================================

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Keywords         Description
    ==============   ====================================================
    checkbounds      If True, values of xout and yout are checked to see
                     that they lie within the range specified by xin
                     and xin.
                     If False, and xout,yout are outside xin,yin,
                     interpolated values will be clipped to values on
                     boundary of input grid (xin,yin)
                     Default is False.
    masked           If True, points outside the range of xin and yin
                     are masked (in a masked array).
                     If masked is set to a number, then
                     points outside the range of xin and yin will be
                     set to that number. Default False.
    order            0 for nearest-neighbor interpolation, 1 for
                     bilinear interpolation, 3 for cublic spline
                     (default 1). order=3 requires scipy.ndimage.
    ==============   ====================================================

    .. note::
     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 numpy.where).
     This effectively uses nearest neighbor interpolation if any of the
     four surrounding points in datain are masked, and bilinear interpolation
     otherwise.

    Returns ``dataout``, the interpolated data on the grid ``xout, yout``.
    """
    # 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 xout.min() < xin.min() or \
           xout.max() > xin.max() or \
           yout.min() < yin.min() or \
           yout.max() > yin.max():
            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 = xout.flatten(); youtflat = yout.flatten()
        ix = (np.searchsorted(xin,xoutflat)-1).tolist()
        iy = (np.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 = np.reshape(xcoords,xout.shape)
        ycoords = np.reshape(ycoords,yout.shape)
    # data outside range xin,yin will be clipped to
    # values on boundary.
    if masked:
        xmask = np.logical_or(np.less(xcoords,0),np.greater(xcoords,len(xin)-1))
        ymask = np.logical_or(np.less(ycoords,0),np.greater(ycoords,len(yin)-1))
        xymask = np.logical_or(xmask,ymask)
    xcoords = np.clip(xcoords,0,len(xin)-1)
    ycoords = np.clip(ycoords,0,len(yin)-1)
    # interpolate to output grid using bilinear interpolation.
    if order == 1:
        xi = xcoords.astype(np.int32)
        yi = ycoords.astype(np.int32)
        xip1 = xi+1
        yip1 = yi+1
        xip1 = np.clip(xip1,0,len(xin)-1)
        yip1 = np.clip(yip1,0,len(yin)-1)
        delx = xcoords-xi.astype(np.float32)
        dely = ycoords-yi.astype(np.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 = np.around(xcoords).astype(np.int32)
        ycoordsi = np.around(ycoords).astype(np.int32)
        dataout = datain[ycoordsi,xcoordsi]
    elif order == 3:
        try:
            from scipy.ndimage import map_coordinates
        except ImportError:
            raise ValueError('scipy.ndimage must be installed if order=3')
        coords = [ycoords,xcoords]
        dataout = map_coordinates(datain,coords,order=3,mode='nearest')
    else:
        raise ValueError('order keyword must be 0, 1 or 3')
    if masked:
        newmask = ma.mask_or(ma.getmask(dataout), xymask)
        dataout = ma.masked_array(dataout, mask=newmask)
        if not isinstance(masked, bool):
            dataout = dataout.filled(masked)
    return dataout

def shiftgrid(lon0,datain,lonsin,start=True,cyclic=360.0):
    """
    Shift global lat/lon grid east or west.

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Arguments        Description
    ==============   ====================================================
    lon0             starting longitude for shifted grid
                     (ending longitude if start=False). lon0 must be on
                     input grid (within the range of lonsin).
    datain           original data with longitude the right-most
                     dimension.
    lonsin           original longitudes.
    ==============   ====================================================

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Keywords         Description
    ==============   ====================================================
    start            if True, lon0 represents the starting longitude
                     of the new grid. if False, lon0 is the ending
                     longitude. Default True.
    cyclic           width of periodic domain (default 360)
    ==============   ====================================================

    returns ``dataout,lonsout`` (data and longitudes on shifted grid).
    """
    if np.fabs(lonsin[-1]-lonsin[0]-cyclic) > 1.e-4:
        # Use all data instead of raise ValueError, 'cyclic point not included'
        start_idx = 0
    else:
        # If cyclic, remove the duplicate point
        start_idx = 1
    if lon0 < lonsin[0] or lon0 > lonsin[-1]:
        raise ValueError('lon0 outside of range of lonsin')
    i0 = np.argmin(np.fabs(lonsin-lon0))
    i0_shift = len(lonsin)-i0
    if ma.isMA(datain):
        dataout  = ma.zeros(datain.shape,datain.dtype)
    else:
        dataout  = np.zeros(datain.shape,datain.dtype)
    if ma.isMA(lonsin):
        lonsout = ma.zeros(lonsin.shape,lonsin.dtype)
    else:
        lonsout = np.zeros(lonsin.shape,lonsin.dtype)
    if start:
        lonsout[0:i0_shift] = lonsin[i0:]
    else:
        lonsout[0:i0_shift] = lonsin[i0:]-cyclic
    dataout[...,0:i0_shift] = datain[...,i0:]
    if start:
        lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]+cyclic
    else:
        lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]
    dataout[...,i0_shift:] = datain[...,start_idx:i0+start_idx]
    return dataout,lonsout

def addcyclic(*arr,**kwargs):
    """
    Adds cyclic (wraparound) points in longitude to one or several arrays,
    the last array being longitudes in degrees. e.g.

   ``data1out, data2out, lonsout = addcyclic(data1,data2,lons)``

    ==============   ====================================================
    Keywords         Description
    ==============   ====================================================
    axis             the dimension representing longitude (default -1,
                     or right-most)
    cyclic           width of periodic domain (default 360)
    ==============   ====================================================
    """
    # get (default) keyword arguments
    axis = kwargs.get('axis',-1)
    cyclic = kwargs.get('cyclic',360)
    # define functions
    def _addcyclic(a):
        """addcyclic function for a single data array"""
        npsel = np.ma if np.ma.is_masked(a) else np
        slicer = [slice(None)] * np.ndim(a)
        try:
            slicer[axis] = slice(0, 1)
        except IndexError:
            raise ValueError('The specified axis does not correspond to an '
                    'array dimension.')
        return npsel.concatenate((a,a[tuple(slicer)]),axis=axis)
    def _addcyclic_lon(a):
        """addcyclic function for a single longitude array"""
        # select the right numpy functions
        npsel = np.ma if np.ma.is_masked(a) else np
        # get cyclic longitudes
        clon = (np.take(a,[0],axis=axis)
                + cyclic * np.sign(np.diff(np.take(a,[0,-1],axis=axis),axis=axis)))
        # ensure the values do not exceed cyclic
        clonmod = npsel.where(clon<=cyclic,clon,np.mod(clon,cyclic))
        return npsel.concatenate((a,clonmod),axis=axis)
    # process array(s)
    if len(arr) == 1:
        return _addcyclic_lon(arr[-1])
    else:
        return list(map(_addcyclic,arr[:-1])) + [_addcyclic_lon(arr[-1])]

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 invalid projection points on output
    if llcrnrlon > 1.e20 or urcrnrlon > 1.e20:
        raise ValueError('width and/or height too large for this projection, try smaller values')
    else:
        return corners

def _choosecornersllur(llcrnrx, llcrnry, urcrnrx, urcrnry,**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(urcrnrx, urcrnry, inverse=True)
    llcrnrlon, llcrnrlat = p(llcrnrx, llcrnry, inverse=True)
    corners = llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat
    # test for invalid projection points on output
    if llcrnrlon > 1.e20 or urcrnrlon > 1.e20:
        raise ValueError('width and/or height too large for this projection, try smaller values')
    else:
        return corners

def maskoceans(lonsin,latsin,datain,inlands=True,resolution='l',grid=5):
    """
    mask data (``datain``), defined on a grid with latitudes ``latsin``
    longitudes ``lonsin`` so that points over water will not be plotted.

    .. tabularcolumns:: |l|L|

    ==============   ====================================================
    Arguments        Description
    ==============   ====================================================
    lonsin, latsin   rank-2 arrays containing longitudes and latitudes of
                     grid.
    datain           rank-2 input array on grid defined by ``lonsin`` and
                     ``latsin``.
    inlands          if False, masked only ocean points and not inland
                     lakes (Default True).
    resolution       gshhs coastline resolution used to define land/sea
                     mask (default 'l', available 'c','l','i','h' or 'f')
    grid             land/sea mask grid spacing in minutes (Default 5;
                     10, 2.5 and 1.25 are also available).
    ==============   ====================================================

    returns a masked array the same shape as datain with "wet" points masked.
    """
    # read in land/sea mask.
    lsmask_lons, lsmask_lats, lsmask =\
    _readlsmask(lakes=inlands,resolution=resolution,grid=grid)
    # nearest-neighbor interpolation to output grid.
    lsmasko = interp(lsmask,lsmask_lons,lsmask_lats,lonsin,latsin,masked=True,order=0)
    # mask input data.
    mask = lsmasko == 0
    return ma.masked_array(datain,mask=mask)

def _readlsmask(lakes=True,resolution='l',grid=5):
    # read in land/sea mask.
    if grid == 10:
        nlons = 2160
    elif grid == 5:
        nlons = 4320
    elif grid == 2.5:
        nlons = 8640
    elif grid == 1.25:
        nlons = 17280
    else:
        raise ValueError('grid for land/sea mask must be 10,5,2.5 or 1.25')
    nlats = nlons//2
    import gzip
    lsmaskf =\
    gzip.open(os.path.join(basemap_datadir,'lsmask_%smin_%s.bin' %\
        (grid,resolution)), 'rb')
    lsmask =\
    np.reshape(np.frombuffer(lsmaskf.read(),dtype=np.uint8),(nlats,nlons))
    if lakes:
        lsmask =\
        np.where(lsmask==2,np.array(0,dtype=np.uint8),lsmask)
    lsmaskf.close()
    delta = 360./nlons
    lsmask_lons = np.linspace(-180+0.5*delta,180-0.5*delta,nlons).astype(np.float32)
    lsmask_lats = np.linspace(-90+0.5*delta,90-0.5*delta,nlats).astype(np.float32)
    return lsmask_lons, lsmask_lats, lsmask

class _tup(tuple):
    # tuple with an added remove method.
    # used for objects returned by drawparallels and drawmeridians.
    def remove(self):
        for item in self:
            for x in item:
                x.remove()
class _dict(dict):
    # override __delitem__ to first call remove method on values.
    def __delitem__(self,key):
        self[key].remove()
        super(_dict, self).__delitem__(key)


def _setlonlab(fmt, lon, labelstyle):
    """Set longitude label string (called by :meth:`Basemap.drawmeridians`)."""

    try:
        # `fmt` is a function that returns a formatted string.
        lonlab = fmt(lon)
    except:
        # `fmt` is a format string.
        degchar = b"\xc2\xb0".decode("utf-8")
        if lon > 180:
            if mpl.rcParams["text.usetex"]:
                if labelstyle == "+/-":
                    lonlabstr = r"${\/-%s\/^{\circ}}$" % fmt
                else:
                    lonlabstr = r"${%s\/^{\circ}\/W}$" % fmt
            else:
                if labelstyle == "+/-":
                    lonlabstr = r"-%s%s" % (fmt, degchar)
                else:
                    lonlabstr = r"%s%sW" % (fmt, degchar)
            lonlab = lonlabstr % np.fabs(lon - 360)
        elif lon < 180 and lon != 0:
            if mpl.rcParams["text.usetex"]:
                if labelstyle == "+/-":
                    lonlabstr = r"${\/+%s\/^{\circ}}$" % fmt
                else:
                    lonlabstr = r"${%s\/^{\circ}\/E}$" % fmt
            else:
                if labelstyle == "+/-":
                    lonlabstr = r"+%s%s" % (fmt, degchar)
                else:
                    lonlabstr = r"%s%sE" % (fmt, degchar)
            lonlab = lonlabstr % lon
        else:
            if mpl.rcParams["text.usetex"]:
                lonlabstr = r"${%s\/^{\circ}}$" % fmt
            else:
                lonlabstr = r"%s%s" % (fmt, degchar)
            lonlab = lonlabstr % lon
    return lonlab


def _setlatlab(fmt, lat, labelstyle):
    """Set latitude label string (called by :meth:`Basemap.drawparallels`)."""

    try:
        # `fmt` is a function that returns a formatted string.
        latlab = fmt(lat)
    except:
        # `fmt` is a format string.
        degchar = b"\xc2\xb0".decode("utf-8")
        if lat < 0:
            if mpl.rcParams["text.usetex"]:
                if labelstyle == "+/-":
                    latlabstr = r"${\/-%s\/^{\circ}}$" % fmt
                else:
                    latlabstr = r"${%s\/^{\circ}\/S}$" % fmt
            else:
                if labelstyle == "+/-":
                    latlabstr = r"-%s%s" % (fmt, degchar)
                else:
                    latlabstr = r"%s%sS" % (fmt, degchar)
            latlab = latlabstr % np.fabs(lat)
        elif lat > 0:
            if mpl.rcParams["text.usetex"]:
                if labelstyle == "+/-":
                    latlabstr = r"${\/+%s\/^{\circ}}$" % fmt
                else:
                    latlabstr = r"${%s\/^{\circ}\/N}$" % fmt
            else:
                if labelstyle == "+/-":
                    latlabstr = r"+%s%s" % (fmt, degchar)
                else:
                    latlabstr = r"%s%sN" % (fmt, degchar)
            latlab = latlabstr % lat
        else:
            if mpl.rcParams["text.usetex"]:
                latlabstr = r"${%s\/^{\circ}}$" % fmt
            else:
                latlabstr = r"%s%s" % (fmt, degchar)
            latlab = latlabstr % lat
    return latlab