Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

5315 lines (4996 sloc) 253.356 kb
"""
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.
"""
from matplotlib import __version__ as _matplotlib_version
from matplotlib.cbook import is_scalar, dedent
# check to make sure matplotlib is not too old.
_mpl_required_version = '0.98'
if _matplotlib_version < _mpl_required_version:
msg = dedent("""
your matplotlib is too old - basemap requires version %s or
higher, you have version %s""" %
(_mpl_required_version,_matplotlib_version))
raise ImportError(msg)
from matplotlib import rcParams, is_interactive
from matplotlib.collections import LineCollection, PolyCollection
from matplotlib.patches import Ellipse, Circle, Polygon, FancyArrowPatch
from matplotlib.lines import Line2D
from matplotlib.transforms import Bbox
from mpl_toolkits.basemap import pyproj
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.image import imread
import sys, os, math
from .proj import Proj
import numpy as np
import numpy.ma as ma
import _geoslib
import functools
# 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:
basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
__version__ = '1.0.8'
# 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 False,
so parallels and meridians can be labelled instead.
If parallel or meridian labelling is requested
(using drawparallels and drawmeridians methods),
automatic tick labelling will be supressed even 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 don 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):
param = float(param)
if param > maxval or param < minval:
raise ValueError('%s must be between %f and %f degrees' %
(name, minval, maxval))
return param
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)
# 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)
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
_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,ucrnrlon,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,ucrnrlon,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):
msg='boundinglat cannot extend into opposite hemisphere'
raise ValueError(msg)
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,ucrnrlon,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,ucrnrlon,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:
msg='must specify lon_0 and lat_0 for near-sided perspective Basemap'
raise ValueError(msg)
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,ucrnrlon,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:
msg='must specify lon_0,o_lat_p,o_lon_p for rotated pole Basemap'
raise ValueError(msg)
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)
# 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.
"""
msg = dedent("""
Unable to open boundary dataset file. Only the 'crude', 'low',
'intermediate' and 'high' resolution datasets are installed by default.
If you are requesting a 'full' resolution dataset, you may need to
download and install those files separately
(see the basemap README for details).""")
# 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(msg)
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.fromstring(polystring,dtype='<f4'),'f8')
b.shape = (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[4:,:]
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:
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)
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:
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:
ax.set_frame_on(False)
# elliptical region.
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)
limb.set_clip_on(False)
if zorder is not None:
limb.set_zorder(zorder)
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,edgecolor=color,linewidth=linewidth)
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_clip_on(False)
if zorder is not None:
limb.set_zorder(zorder)
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')
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_clip_on(False)
if zorder is not None:
limb.set_zorder(zorder)
else: # all other projections are rectangular.
# use axesPatch for fill_color, spine for border line props.
for spine in ax.spines.values():
spine.set_linewidth(linewidth)
if self.projection not in ['geos','ortho','nsper']:
limb = ax.axesPatch
if fill_color is not None:
limb.set_facecolor(fill_color)
for spine in ax.spines.values():
spine.set_edgecolor(color)
ax.set_frame_on(True)
# FIXME? should zorder be set separately for edge and background?
if zorder is not None:
limb.set_zorder(zorder)
for spine in ax.spines.values():
spine.set_zorder(zorder)
else:
# use axesPatch for fill_color, spine for border line props.
for spine in ax.spines.values():
spine.set_edgecolor(color)
ax.set_frame_on(True)
# FIXME? should zorder be set separately for edge and background?
if zorder is not None:
ax.axesPatch.set_zorder(zorder)
for spine in ax.spines.values():
spine.set_zorder(zorder)
# for geos or ortho projections, also
# draw and fill map projection limb, clipped
# to rectangular region.
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.
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 for round polar plots.
if self.round: polys,c = self._clipcircle(ax,polys)
return polys
def _clipcircle(self,ax,coll):
c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
radius=0.5*(self.xmax-self.xmin),fc='none')
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)
# clip coastlines for round polar plots.
if self.round: coastlines,c = self._clipcircle(ax,coastlines)
ax.add_collection(coastlines)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
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)
# clip countries for round polar plots.
if self.round: countries,c = self._clipcircle(ax,countries)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
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)
# clip states for round polar plots.
if self.round: states,c = self._clipcircle(ax,states)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
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)
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)
# clip rivers for round polar plots.
if self.round: rivers,c = self._clipcircle(ax,rivers)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
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'):
"""
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 mathplotlib.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)
============== ====================================================
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)
except:
raise IOError('error reading shapefile %s.shp' % shapefile)
fields = shf.fields
coords = []; attributes = []
msg=dedent("""
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=(1,))
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)
# clip boundaries for round polar plots.
if self.round: lines,c = self._clipcircle(ax,lines)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
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):
"""
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).
if self.round:
c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
radius=0.5*(self.xmax-self.xmin),fc='none')
if c not in ax.patches:
p = ax.add_patch(c)
p.set_clip_on(False)
for par in pardict:
lines,labs = pardict[par]
for l in lines:
l.set_clip_path(c)
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):
"""
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(dedent(
"""'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:
c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
radius=0.5*(self.xmax-self.xmin),fc='none')
if c not in ax.patches:
p = ax.add_patch(c)
p.set_clip_on(False)
# label desired?
label = False
for lab in labels:
if lab: label = True
for merid in meridict:
lines,labs = meridict[merid]
# clip lines.
for l in lines:
l.set_clip_path(c)
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)
return meridict
def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs):
"""
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)
# clip polygons for round polar plots.
if self.round: poly,c = self._clipcircle(ax,poly)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
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):
"""
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 matplotlib.lines.Line2D object.
"""
# 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)[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 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')
# draw another filled patch, with no boundary.
limb2 = self.drawmapboundary(linewidth=0)
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 is_interactive():
import matplotlib.pyplot as plt
plt.draw_if_interactive()
@_transform1d
def scatter(self, *args, **kwargs):
"""
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)
# allow callers to override the hold state by passing hold=True|False
b = ax.ishold()
h = kwargs.pop('hold',None)
if h is not None:
ax.hold(h)
try:
ret = ax.scatter(*args, **kwargs)
except:
ax.hold(b)
raise
ax.hold(b)
# reset current active image (only if pyplot is imported).
if plt:
plt.sci(ret)
# clip for round polar plots.
if self.round: ret,c = self._clipcircle(ax,ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
return ret
@_transform1d
def plot(self, *args, **kwargs):
"""
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 projectio