from __future__ import (absolute_import, division, print_function)
"""
Module for plotting data on maps with matplotlib.
Contains the :class:`Basemap` class (which does most of the
heavy lifting), and the following functions:
:func:`interp`: bilinear interpolation between rectilinear grids.
:func:`maskoceans`: mask 'wet' points of an input array.
:func:`shiftgrid`: shifts global lat/lon grids east or west.
:func:`addcyclic`: Add cyclic (wraparound) point in longitude.
"""
import os
import sys
import math
import warnings
import functools
try:
from urllib2 import urlopen
from urllib import urlretrieve
except ImportError:
from urllib.request import urlopen
from urllib.request import urlretrieve
import numpy as np
import numpy.ma as ma
import matplotlib as mpl
from matplotlib.artist import Artist
from matplotlib.collections import LineCollection
from matplotlib.collections import PolyCollection
from matplotlib.image import imread
from matplotlib.lines import Line2D
from matplotlib.patches import Circle
from matplotlib.patches import Ellipse
from matplotlib.patches import FancyArrowPatch
from matplotlib.patches import Polygon
from matplotlib.transforms import Bbox
from mpl_toolkits.axes_grid1 import make_axes_locatable
import pyproj
import _geoslib
from . proj import Proj
__version__ = "1.4.1"
# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
# check to see if environment variable BASEMAPDATA set to a directory,
# and if so look for the data there.
if 'BASEMAPDATA' in os.environ:
basemap_datadir = os.environ['BASEMAPDATA']
if not os.path.isdir(basemap_datadir):
raise RuntimeError('Path in environment BASEMAPDATA not a directory')
else:
from mpl_toolkits import basemap_data
basemap_datadir = os.path.abspath(list(basemap_data.__path__)[0])
# module variable that sets the default value for the 'latlon' kwarg.
# can be set to True by user so plotting functions can take lons,lats
# in degrees by default, instead of x,y (map projection coords in meters).
latlon_default = False
# supported map projections.
_projnames = {'cyl' : 'Cylindrical Equidistant',
'merc' : 'Mercator',
'tmerc' : 'Transverse Mercator',
'omerc' : 'Oblique Mercator',
'mill' : 'Miller Cylindrical',
'gall' : 'Gall Stereographic Cylindrical',
'cea' : 'Cylindrical Equal Area',
'lcc' : 'Lambert Conformal',
'laea' : 'Lambert Azimuthal Equal Area',
'nplaea' : 'North-Polar Lambert Azimuthal',
'splaea' : 'South-Polar Lambert Azimuthal',
'eqdc' : 'Equidistant Conic',
'aeqd' : 'Azimuthal Equidistant',
'npaeqd' : 'North-Polar Azimuthal Equidistant',
'spaeqd' : 'South-Polar Azimuthal Equidistant',
'aea' : 'Albers Equal Area',
'stere' : 'Stereographic',
'npstere' : 'North-Polar Stereographic',
'spstere' : 'South-Polar Stereographic',
'cass' : 'Cassini-Soldner',
'poly' : 'Polyconic',
'ortho' : 'Orthographic',
'geos' : 'Geostationary',
'nsper' : 'Near-Sided Perspective',
'sinu' : 'Sinusoidal',
'moll' : 'Mollweide',
'hammer' : 'Hammer',
'robin' : 'Robinson',
'kav7' : 'Kavrayskiy VII',
'eck4' : 'Eckert IV',
'vandg' : 'van der Grinten',
'mbtfpq' : 'McBryde-Thomas Flat-Polar Quartic',
'gnom' : 'Gnomonic',
'rotpole' : 'Rotated Pole',
}
supported_projections = []
for _items in _projnames.items():
supported_projections.append(" %-17s%-40s\n" % (_items))
supported_projections = ''.join(supported_projections)
_cylproj = ['cyl','merc','mill','gall','cea']
_pseudocyl = ['moll','robin','eck4','kav7','sinu','mbtfpq','vandg','hammer']
_dg2rad = math.radians(1.)
_rad2dg = math.degrees(1.)
# projection specific parameters.
projection_params = {'cyl' : 'corners only (no width/height)',
'merc' : 'corners plus lat_ts (no width/height)',
'tmerc' : 'lon_0,lat_0,k_0',
'omerc' : 'lon_0,lat_0,lat_1,lat_2,lon_1,lon_2,no_rot,k_0',
'mill' : 'corners only (no width/height)',
'gall' : 'corners only (no width/height)',
'cea' : 'corners only plus lat_ts (no width/height)',
'lcc' : 'lon_0,lat_0,lat_1,lat_2,k_0',
'laea' : 'lon_0,lat_0',
'nplaea' : 'bounding_lat,lon_0,lat_0,no corners or width/height',
'splaea' : 'bounding_lat,lon_0,lat_0,no corners or width/height',
'eqdc' : 'lon_0,lat_0,lat_1,lat_2',
'aeqd' : 'lon_0,lat_0',
'npaeqd' : 'bounding_lat,lon_0,lat_0,no corners or width/height',
'spaeqd' : 'bounding_lat,lon_0,lat_0,no corners or width/height',
'aea' : 'lon_0,lat_0,lat_1',
'stere' : 'lon_0,lat_0,lat_ts,k_0',
'npstere' : 'bounding_lat,lon_0,lat_0,no corners or width/height',
'spstere' : 'bounding_lat,lon_0,lat_0,no corners or width/height',
'cass' : 'lon_0,lat_0',
'poly' : 'lon_0,lat_0',
'ortho' : 'lon_0,lat_0,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
'geos' : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
'nsper' : 'lon_0,satellite_height,llcrnrx,llcrnry,urcrnrx,urcrnry,no width/height',
'sinu' : 'lon_0,lat_0,no corners or width/height',
'moll' : 'lon_0,lat_0,no corners or width/height',
'hammer' : 'lon_0,lat_0,no corners or width/height',
'robin' : 'lon_0,lat_0,no corners or width/height',
'eck4' : 'lon_0,lat_0,no corners or width/height',
'kav7' : 'lon_0,lat_0,no corners or width/height',
'vandg' : 'lon_0,lat_0,no corners or width/height',
'mbtfpq' : 'lon_0,lat_0,no corners or width/height',
'gnom' : 'lon_0,lat_0',
'rotpole' : 'lon_0,o_lat_p,o_lon_p,corner lat/lon or corner x,y (no width/height)'
}
# create dictionary that maps epsg codes to Basemap kwargs.
epsgf = open(os.path.join(basemap_datadir, 'epsg'))
epsg_dict={}
for line in epsgf:
if line.startswith("#"):
continue
l = line.split()
code = l[0].strip("<>")
parms = ' '.join(l[1:-1])
_kw_args={}
for s in l[1:-1]:
try:
k,v = s.split('=')
except:
pass
k = k.strip("+")
if k=='proj':
if v == 'longlat': v = 'cyl'
if v not in _projnames:
continue
k='projection'
if k=='k':
k='k_0'
if k in ['projection','lat_1','lat_2','lon_0','lat_0',\
'a','b','k_0','lat_ts','ellps','datum']:
if k not in ['projection','ellps','datum']:
v = float(v)
_kw_args[k]=v
if 'projection' in _kw_args:
if 'a' in _kw_args:
if 'b' in _kw_args:
_kw_args['rsphere']=(_kw_args['a'],_kw_args['b'])
del _kw_args['b']
else:
_kw_args['rsphere']=_kw_args['a']
del _kw_args['a']
if 'datum' in _kw_args:
if _kw_args['datum'] == 'NAD83':
_kw_args['ellps'] = 'GRS80'
elif _kw_args['datum'] == 'NAD27':
_kw_args['ellps'] = 'clrk66'
elif _kw_args['datum'] == 'WGS84':
_kw_args['ellps'] = 'WGS84'
del _kw_args['datum']
# supported epsg projections.
# omerc not supported yet, since we can't handle
# alpha,gamma and lonc keywords.
if _kw_args['projection'] != 'omerc':
epsg_dict[code]=_kw_args
epsgf.close()
# The __init__ docstring is pulled out here because it is so long;
# Having it in the usual place makes it hard to get from the
# __init__ argument list to the code that uses the arguments.
_Basemap_init_doc = """
Sets up a basemap with specified map projection.
and creates the coastline data structures in map projection
coordinates.
Calling a Basemap class instance with the arguments lon, lat will
convert lon/lat (in degrees) to x/y map projection coordinates
(in meters). The inverse transformation is done if the optional keyword
``inverse`` is set to True.
The desired projection is set with the projection keyword. Default is ``cyl``.
Supported values for the projection keyword are:
============== ====================================================
Value Description
============== ====================================================
%(supported_projections)s
============== ====================================================
For most map projections, the map projection region can either be
specified by setting these keywords:
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
llcrnrlon longitude of lower left hand corner of the desired map
domain (degrees).
llcrnrlat latitude of lower left hand corner of the desired map
domain (degrees).
urcrnrlon longitude of upper right hand corner of the desired map
domain (degrees).
urcrnrlat latitude of upper right hand corner of the desired map
domain (degrees).
============== ====================================================
or these
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
width width of desired map domain in projection coordinates
(meters).
height height of desired map domain in projection coordinates
(meters).
lon_0 center of desired map domain (in degrees).
lat_0 center of desired map domain (in degrees).
============== ====================================================
For ``sinu``, ``moll``, ``hammer``, ``npstere``, ``spstere``, ``nplaea``, ``splaea``,
``npaeqd``, ``spaeqd``, ``robin``, ``eck4``, ``kav7``, or ``mbtfpq``, the values of
llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, width and height are ignored
(because either they are computed internally, or entire globe is
always plotted).
For the cylindrical projections (``cyl``, ``merc``, ``mill``, ``cea`` and ``gall``),
the default is to use
llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other
projections except ``ortho``, ``geos`` and ``nsper``, either the lat/lon values of the
corners or width and height must be specified by the user.
For ``ortho``, ``geos`` and ``nsper``, the lat/lon values of the corners may be specified,
or the x/y values of the corners (llcrnrx,llcrnry,urcrnrx,urcrnry) in the
coordinate system of the global projection (with x=0,y=0 at the center
of the global projection). If the corners are not specified,
the entire globe is plotted.
For ``rotpole``, the lat/lon values of the corners on the unrotated sphere
may be provided as llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat, or the lat/lon
values of the corners on the rotated sphere can be given as
llcrnrx,llcrnry,urcrnrx,urcrnry.
Other keyword arguments:
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
resolution resolution of boundary database to use. Can be ``c``
(crude), ``l`` (low), ``i`` (intermediate), ``h``
(high), ``f`` (full) or None.
If None, no boundary data will be read in (and
class methods such as drawcoastlines will raise an
if invoked).
Resolution drops off by roughly 80%% between datasets.
Higher res datasets are much slower to draw.
Default ``c``. Coastline data is from the GSHHS
(http://www.soest.hawaii.edu/wessel/gshhs/gshhs.html).
State, country and river datasets from the Generic
Mapping Tools (http://gmt.soest.hawaii.edu).
area_thresh coastline or lake with an area smaller than
area_thresh in km^2 will not be plotted.
Default 10000,1000,100,10,1 for resolution
``c``, ``l``, ``i``, ``h``, ``f``.
rsphere radius of the sphere used to define map projection
(default 6370997 meters, close to the arithmetic mean
radius of the earth). If given as a sequence, the
first two elements are interpreted as the radii
of the major and minor axes of an ellipsoid.
Note: sometimes an ellipsoid is specified by the
major axis and an inverse flattening parameter (if).
The minor axis (b) can be computed from the major
axis (a) and the inverse flattening parameter using
the formula if = a/(a-b).
ellps string describing ellipsoid ('GRS80' or 'WGS84',
for example). If both rsphere and ellps are given,
rsphere is ignored. Default None. See pyproj.pj_ellps
for allowed values.
suppress_ticks suppress automatic drawing of axis ticks and labels
in map projection coordinates. Default True,
so parallels and meridians can be labelled instead.
If parallel or meridian labelling is requested
(using drawparallels and drawmeridians methods),
automatic tick labelling will be supressed even if
suppress_ticks=False. suppress_ticks=False
is useful if you want to use your own custom tick
formatter, or if you want to let matplotlib label
the axes in meters using map projection
coordinates.
fix_aspect fix aspect ratio of plot to match aspect ratio
of map projection region (default True).
anchor determines how map is placed in axes rectangle
(passed to axes.set_aspect). Default is ``C``,
which means map is centered.
Allowed values are
``C``, ``SW``, ``S``, ``SE``, ``E``, ``NE``,
``N``, ``NW``, and ``W``.
celestial use astronomical conventions for longitude (i.e.
negative longitudes to the east of 0). Default False.
Implies resolution=None.
ax set default axes instance
(default None - matplotlib.pyplot.gca() may be used
to get the current axes instance).
If you do not want matplotlib.pyplot to be imported,
you can either set this to a pre-defined axes
instance, or use the ``ax`` keyword in each Basemap
method call that does drawing. In the first case,
all Basemap method calls will draw to the same axes
instance. In the second case, you can draw to
different axes with the same Basemap instance.
You can also use the ``ax`` keyword in individual
method calls to selectively override the default
axes instance.
============== ====================================================
The following keywords are map projection parameters which all default to
None. Not all parameters are used by all projections, some are ignored.
The module variable ``projection_params`` is a dictionary which
lists which parameters apply to which projections.
.. tabularcolumns:: |l|L|
================ ====================================================
Keyword Description
================ ====================================================
lat_ts latitude of true scale. Optional for stereographic,
cylindrical equal area and mercator projections.
default is lat_0 for stereographic projection.
default is 0 for mercator and cylindrical equal area
projections.
lat_1 first standard parallel for lambert conformal,
albers equal area and equidistant conic.
Latitude of one of the two points on the projection
centerline for oblique mercator. If lat_1 is not given, but
lat_0 is, lat_1 is set to lat_0 for lambert
conformal, albers equal area and equidistant conic.
lat_2 second standard parallel for lambert conformal,
albers equal area and equidistant conic.
Latitude of one of the two points on the projection
centerline for oblique mercator. If lat_2 is not
given it is set to lat_1 for lambert conformal,
albers equal area and equidistant conic.
lon_1 Longitude of one of the two points on the projection
centerline for oblique mercator.
lon_2 Longitude of one of the two points on the projection
centerline for oblique mercator.
k_0 Scale factor at natural origin (used
by 'tmerc', 'omerc', 'stere' and 'lcc').
no_rot only used by oblique mercator.
If set to True, the map projection coordinates will
not be rotated to true North. Default is False
(projection coordinates are automatically rotated).
lat_0 central latitude (y-axis origin) - used by all
projections.
lon_0 central meridian (x-axis origin) - used by all
projections.
o_lat_p latitude of rotated pole (only used by 'rotpole')
o_lon_p longitude of rotated pole (only used by 'rotpole')
boundinglat bounding latitude for pole-centered projections
(npstere,spstere,nplaea,splaea,npaeqd,spaeqd).
These projections are square regions centered
on the north or south pole.
The longitude lon_0 is at 6-o'clock, and the
latitude circle boundinglat is tangent to the edge
of the map at lon_0.
round cut off pole-centered projection at boundinglat
(so plot is a circle instead of a square). Only
relevant for npstere,spstere,nplaea,splaea,npaeqd
or spaeqd projections. Default False.
satellite_height height of satellite (in m) above equator -
only relevant for geostationary
and near-sided perspective (``geos`` or ``nsper``)
projections. Default 35,786 km.
================ ====================================================
Useful instance variables:
.. tabularcolumns:: |l|L|
================ ====================================================
Variable Name Description
================ ====================================================
projection map projection. Print the module variable
``supported_projections`` to see a list of allowed
values.
epsg EPSG code defining projection (see
http://spatialreference.org for a list of
EPSG codes and their definitions).
aspect map aspect ratio
(size of y dimension / size of x dimension).
llcrnrlon longitude of lower left hand corner of the
selected map domain.
llcrnrlat latitude of lower left hand corner of the
selected map domain.
urcrnrlon longitude of upper right hand corner of the
selected map domain.
urcrnrlat latitude of upper right hand corner of the
selected map domain.
llcrnrx x value of lower left hand corner of the
selected map domain in map projection coordinates.
llcrnry y value of lower left hand corner of the
selected map domain in map projection coordinates.
urcrnrx x value of upper right hand corner of the
selected map domain in map projection coordinates.
urcrnry y value of upper right hand corner of the
selected map domain in map projection coordinates.
rmajor equatorial radius of ellipsoid used (in meters).
rminor polar radius of ellipsoid used (in meters).
resolution resolution of boundary dataset being used (``c``
for crude, ``l`` for low, etc.).
If None, no boundary dataset is associated with the
Basemap instance.
proj4string the string describing the map projection that is
used by PROJ.4.
================ ====================================================
**Converting from Geographic (lon/lat) to Map Projection (x/y) Coordinates**
Calling a Basemap class instance with the arguments lon, lat will
convert lon/lat (in degrees) to x/y map projection
coordinates (in meters). If optional keyword ``inverse`` is
True (default is False), the inverse transformation from x/y
to lon/lat is performed.
For cylindrical equidistant projection (``cyl``), this
does nothing (i.e. x,y == lon,lat).
For non-cylindrical projections, the inverse transformation
always returns longitudes between -180 and 180 degrees. For
cylindrical projections (self.projection == ``cyl``, ``mill``,
``cea``, ``gall`` or ``merc``)
the inverse transformation will return longitudes between
self.llcrnrlon and self.llcrnrlat.
Input arguments lon, lat can be either scalar floats, sequences
or numpy arrays.
**Example Usage:**
>>> from mpl_toolkits.basemap import Basemap
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # read in topo data (on a regular lat/lon grid)
>>> etopo = np.loadtxt('etopo20data.gz')
>>> lons = np.loadtxt('etopo20lons.gz')
>>> lats = np.loadtxt('etopo20lats.gz')
>>> # create Basemap instance for Robinson projection.
>>> m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
>>> # compute map projection coordinates for lat/lon grid.
>>> x, y = m(*np.meshgrid(lons,lats))
>>> # make filled contour plot.
>>> cs = m.contourf(x,y,etopo,30,cmap=plt.cm.jet)
>>> m.drawcoastlines() # draw coastlines
>>> m.drawmapboundary() # draw a line around the map region
>>> m.drawparallels(np.arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels
>>> m.drawmeridians(np.arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
>>> plt.title('Robinson Projection') # add a title
>>> plt.show()
[this example (simpletest.py) plus many others can be found in the
examples directory of source distribution. The "OO" version of this
example (which does not use matplotlib.pyplot) is called "simpletest_oo.py".]
""" % locals()
# unsupported projection error message.
_unsupported_projection = ["'%s' is an unsupported projection.\n"]
_unsupported_projection.append("The supported projections are:\n")
_unsupported_projection.append(supported_projections)
_unsupported_projection = ''.join(_unsupported_projection)
def _validated_ll(param, name, minval, maxval):
item = np.squeeze(param)
param = float(param if item.shape else item)
if param > maxval or param < minval:
raise ValueError('%s must be between %f and %f degrees' %
(name, minval, maxval))
return param
def _validated_or_none(param, name, minval, maxval):
if param is None:
return None
return _validated_ll(param, name, minval, maxval)
def _insert_validated(d, param, name, minval, maxval):
if param is not None:
d[name] = _validated_ll(param, name, minval, maxval)
def _transform(plotfunc):
# shift data and longitudes to map projection region, then compute
# transformation to map projection coordinates.
@functools.wraps(plotfunc)
def with_transform(self,x,y,data,*args,**kwargs):
# input coordinates are latitude/longitude, not map projection coords.
if kwargs.pop('latlon', latlon_default):
# shift data to map projection region for
# cylindrical and pseudo-cylindrical projections.
if self.projection in _cylproj or self.projection in _pseudocyl:
x, data = self.shiftdata(x, data,
fix_wrap_around=plotfunc.__name__ not in ["scatter"])
# convert lat/lon coords to map projection coords.
x, y = self(x,y)
return plotfunc(self,x,y,data,*args,**kwargs)
return with_transform
def _transform1d(plotfunc):
# shift data and longitudes to map projection region, then compute
# transformation to map projection coordinates.
@functools.wraps(plotfunc)
def with_transform(self,x,y,*args,**kwargs):
x = np.asarray(x)
# input coordinates are latitude/longitude, not map projection coords.
if kwargs.pop('latlon', latlon_default):
# shift data to map projection region for
# cylindrical and pseudo-cylindrical projections.
if self.projection in _cylproj or self.projection in _pseudocyl:
if x.ndim == 1:
x = self.shiftdata(x, fix_wrap_around=plotfunc.__name__ not in ["scatter"])
elif x.ndim == 0:
if x > 180:
x = x - 360.
# convert lat/lon coords to map projection coords.
x, y = self(x,y)
return plotfunc(self,x,y,*args,**kwargs)
return with_transform
def _transformuv(plotfunc):
# shift data and longitudes to map projection region, then compute
# transformation to map projection coordinates. Works when call
# signature has two data arrays instead of one.
@functools.wraps(plotfunc)
def with_transform(self,x,y,u,v,*args,**kwargs):
# input coordinates are latitude/longitude, not map projection coords.
if kwargs.pop('latlon', latlon_default):
# shift data to map projection region for
# cylindrical and pseudo-cylindrical projections.
if self.projection in _cylproj or self.projection in _pseudocyl:
x1, u = self.shiftdata(x, u)
x, v = self.shiftdata(x, v)
# convert lat/lon coords to map projection coords.
x, y = self(x,y)
return plotfunc(self,x,y,u,v,*args,**kwargs)
return with_transform
[docs]class Basemap(object):
[docs] def __init__(self, llcrnrlon=None, llcrnrlat=None,
urcrnrlon=None, urcrnrlat=None,
llcrnrx=None, llcrnry=None,
urcrnrx=None, urcrnry=None,
width=None, height=None,
projection='cyl', resolution='c',
area_thresh=None, rsphere=6370997.0,
ellps=None, lat_ts=None,
lat_1=None, lat_2=None,
lat_0=None, lon_0=None,
lon_1=None, lon_2=None,
o_lon_p=None, o_lat_p=None,
k_0=None,
no_rot=False,
suppress_ticks=True,
satellite_height=35786000,
boundinglat=None,
fix_aspect=True,
anchor='C',
celestial=False,
round=False,
epsg=None,
ax=None):
# docstring is added after __init__ method definition
# set epsg code if given, set to 4326 for projection='cyl':
if epsg is not None:
self.epsg = epsg
elif projection == 'cyl':
self.epsg = 4326
# replace kwarg values with those implied by epsg code,
# if given.
if hasattr(self,'epsg'):
if str(self.epsg) not in epsg_dict:
raise ValueError('%s is not a supported EPSG code' %
self.epsg)
epsg_params = epsg_dict[str(self.epsg)]
for k in epsg_params:
if k == 'projection':
projection = epsg_params[k]
elif k == 'rsphere':
rsphere = epsg_params[k]
elif k == 'ellps':
ellps = epsg_params[k]
elif k == 'lat_1':
lat_1 = epsg_params[k]
elif k == 'lat_2':
lat_2 = epsg_params[k]
elif k == 'lon_0':
lon_0 = epsg_params[k]
elif k == 'lat_0':
lat_0 = epsg_params[k]
elif k == 'lat_ts':
lat_ts = epsg_params[k]
elif k == 'k_0':
k_0 = epsg_params[k]
# fix aspect to ratio to match aspect ratio of map projection
# region
self.fix_aspect = fix_aspect
# where to put plot in figure (default is 'C' or center)
self.anchor = anchor
# geographic or celestial coords?
self.celestial = celestial
# map projection.
self.projection = projection
# bounding lat (for pole-centered plots)
self.boundinglat = boundinglat
# is a round pole-centered plot desired?
self.round = round
# full disk projection?
self._fulldisk = False # default value
# set up projection parameter dict.
projparams = {}
projparams['proj'] = projection
# if ellps keyword specified, it over-rides rsphere.
if ellps is not None:
try:
elldict = pyproj.pj_ellps[ellps]
except KeyError:
raise ValueError(
'illegal ellps definition, allowed values are %s' %
pyproj.pj_ellps.keys())
projparams['a'] = elldict['a']
if 'b' in elldict:
projparams['b'] = elldict['b']
else:
projparams['b'] = projparams['a']*(1.0-(1.0/elldict['rf']))
else:
try:
if rsphere[0] > rsphere[1]:
projparams['a'] = rsphere[0]
projparams['b'] = rsphere[1]
else:
projparams['a'] = rsphere[1]
projparams['b'] = rsphere[0]
except:
if projection == 'tmerc':
# use bR_a instead of R because of obscure bug
# in proj4 for tmerc projection.
projparams['bR_a'] = rsphere
else:
projparams['R'] = rsphere
# set units to meters.
projparams['units']='m'
# check for sane values of lon_0, lat_0, lat_ts, lat_1, lat_2
lat_0 = _validated_or_none(lat_0, 'lat_0', -90, 90)
lat_1 = _validated_or_none(lat_1, 'lat_1', -90, 90)
lat_2 = _validated_or_none(lat_2, 'lat_2', -90, 90)
lat_ts = _validated_or_none(lat_ts, 'lat_ts', -90, 90)
lon_0 = _validated_or_none(lon_0, 'lon_0', -360, 720)
lon_1 = _validated_or_none(lon_1, 'lon_1', -360, 720)
lon_2 = _validated_or_none(lon_2, 'lon_2', -360, 720)
llcrnrlon = _validated_or_none(llcrnrlon, 'llcrnrlon', -360, 720)
urcrnrlon = _validated_or_none(urcrnrlon, 'urcrnrlon', -360, 720)
llcrnrlat = _validated_or_none(llcrnrlat, 'llcrnrlat', -90, 90)
urcrnrlat = _validated_or_none(urcrnrlat, 'urcrnrlat', -90, 90)
_insert_validated(projparams, lat_0, 'lat_0', -90, 90)
_insert_validated(projparams, lat_1, 'lat_1', -90, 90)
_insert_validated(projparams, lat_2, 'lat_2', -90, 90)
_insert_validated(projparams, lat_ts, 'lat_ts', -90, 90)
_insert_validated(projparams, lon_0, 'lon_0', -360, 720)
_insert_validated(projparams, lon_1, 'lon_1', -360, 720)
_insert_validated(projparams, lon_2, 'lon_2', -360, 720)
if projection in ['geos','nsper']:
projparams['h'] = satellite_height
# check for sane values of projection corners.
using_corners = (None not in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat])
if using_corners:
self.llcrnrlon = _validated_ll(llcrnrlon, 'llcrnrlon', -360, 720)
self.urcrnrlon = _validated_ll(urcrnrlon, 'urcrnrlon', -360, 720)
self.llcrnrlat = _validated_ll(llcrnrlat, 'llcrnrlat', -90, 90)
self.urcrnrlat = _validated_ll(urcrnrlat, 'urcrnrlat', -90, 90)
# for each of the supported projections,
# compute lat/lon of domain corners
# and set values in projparams dict as needed.
if projection in ['lcc', 'eqdc', 'aea']:
if projection == 'lcc' and k_0 is not None:
projparams['k_0']=k_0
# if lat_0 is given, but not lat_1,
# set lat_1=lat_0
if lat_1 is None and lat_0 is not None:
lat_1 = lat_0
projparams['lat_1'] = lat_1
if lat_1 is None or lon_0 is None:
raise ValueError('must specify lat_1 or lat_0 and lon_0 for %s basemap (lat_2 is optional)' % _projnames[projection])
if lat_2 is None:
projparams['lat_2'] = lat_1
if not using_corners:
using_cornersxy = (None not in [llcrnrx,llcrnry,urcrnrx,urcrnry])
if using_cornersxy:
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecornersllur(llcrnrx,llcrnry,urcrnrx,urcrnry,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
else:
if width is None or height is None:
raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
if lon_0 is None or lat_0 is None:
raise ValueError('must specify lon_0 and lat_0 when using width, height to specify projection region')
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'stere':
if k_0 is not None:
projparams['k_0']=k_0
if lat_0 is None or lon_0 is None:
raise ValueError('must specify lat_0 and lon_0 for Stereographic basemap (lat_ts is optional)')
if not using_corners:
if width is None or height is None:
raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
if lon_0 is None or lat_0 is None:
raise ValueError('must specify lon_0 and lat_0 when using width, height to specify projection region')
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection in ['spstere', 'npstere',
'splaea', 'nplaea',
'spaeqd', 'npaeqd']:
if (projection == 'splaea' and boundinglat >= 0) or\
(projection == 'nplaea' and boundinglat <= 0):
raise ValueError('boundinglat cannot extend into opposite hemisphere')
if boundinglat is None or lon_0 is None:
raise ValueError('must specify boundinglat and lon_0 for %s basemap' % _projnames[projection])
if projection[0] == 's':
sgn = -1
else:
sgn = 1
rootproj = projection[2:]
projparams['proj'] = rootproj
if rootproj == 'stere':
projparams['lat_ts'] = sgn * 90.
projparams['lat_0'] = sgn * 90.
self.llcrnrlon = lon_0 - sgn*45.
self.urcrnrlon = lon_0 + sgn*135.
proj = pyproj.Proj(projparams)
x,y = proj(lon_0,boundinglat)
lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
self.urcrnrlat = self.llcrnrlat
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[projection])
elif projection == 'laea':
if lat_0 is None or lon_0 is None:
raise ValueError('must specify lat_0 and lon_0 for Lambert Azimuthal basemap')
if not using_corners:
if width is None or height is None:
raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
if lon_0 is None or lat_0 is None:
raise ValueError('must specify lon_0 and lat_0 when using width, height to specify projection region')
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection in ['tmerc','gnom','cass','poly'] :
if projection == 'tmerc' and k_0 is not None:
projparams['k_0']=k_0
if projection == 'gnom' and 'R' not in projparams:
raise ValueError('gnomonic projection only works for perfect spheres - not ellipsoids')
if lat_0 is None or lon_0 is None:
raise ValueError('must specify lat_0 and lon_0 for Transverse Mercator, Gnomonic, Cassini-Soldnerr and Polyconic basemap')
if not using_corners:
if width is None or height is None:
raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
if lon_0 is None or lat_0 is None:
raise ValueError('must specify lon_0 and lat_0 when using width, height to specify projection region')
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'ortho':
if 'R' not in projparams:
raise ValueError('orthographic projection only works for perfect spheres - not ellipsoids')
if lat_0 is None or lon_0 is None:
raise ValueError('must specify lat_0 and lon_0 for Orthographic basemap')
if (lat_0 == 90 or lat_0 == -90) and\
None in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
# for ortho plot centered on pole, set boundinglat to equator.
# (so meridian labels can be drawn in this special case).
self.boundinglat = 0
self.round = True
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
urcrnrlat = 90.
self._fulldisk = True
else:
self._fulldisk = False
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
# FIXME: won't work for points exactly on equator??
if np.abs(lat_0) < 1.e-2: lat_0 = 1.e-2
projparams['lat_0'] = lat_0
elif projection == 'geos':
if lat_0 is not None and lat_0 != 0:
raise ValueError('lat_0 must be zero for Geostationary basemap')
if lon_0 is None:
raise ValueError('must specify lon_0 for Geostationary basemap')
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
urcrnrlat = 90.
self._fulldisk = True
else:
self._fulldisk = False
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'nsper':
if 'R' not in projparams:
raise ValueError('near-sided perspective projection only works for perfect spheres - not ellipsoids')
if lat_0 is None or lon_0 is None:
raise ValueError('must specify lon_0 and lat_0 for near-sided perspective Basemap')
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
urcrnrlat = 90.
self._fulldisk = True
else:
self._fulldisk = False
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection in _pseudocyl:
if lon_0 is None:
raise ValueError('must specify lon_0 for %s projection' % _projnames[self.projection])
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
llcrnrlon = lon_0-180.
llcrnrlat = -90.
urcrnrlon = lon_0+180
urcrnrlat = 90.
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'omerc':
if k_0 is not None:
projparams['k_0']=k_0
if lat_1 is None or lon_1 is None or lat_2 is None or lon_2 is None:
raise ValueError('must specify lat_1,lon_1 and lat_2,lon_2 for Oblique Mercator basemap')
projparams['lat_1'] = lat_1
projparams['lon_1'] = lon_1
projparams['lat_2'] = lat_2
projparams['lon_2'] = lon_2
projparams['lat_0'] = lat_0
if no_rot:
projparams['no_rot']=''
#if not using_corners:
# raise ValueError, 'cannot specify map region with width and height keywords for this projection, please specify lat/lon values of corners'
if not using_corners:
if width is None or height is None:
raise ValueError('must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat) in degrees or width and height in meters')
if lon_0 is None or lat_0 is None:
raise ValueError('must specify lon_0 and lat_0 when using width, height to specify projection region')
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'aeqd':
if lat_0 is None or lon_0 is None:
raise ValueError('must specify lat_0 and lon_0 for Azimuthal Equidistant basemap')
if not using_corners:
if width is None or height is None:
self._fulldisk = True
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
urcrnrlat = 90.
else:
self._fulldisk = False
if lon_0 is None or lat_0 is None:
raise ValueError('must specify lon_0 and lat_0 when using width, height to specify projection region')
if not self._fulldisk:
llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection in _cylproj:
if projection == 'merc' or projection == 'cea':
if lat_ts is None:
lat_ts = 0.
projparams['lat_ts'] = lat_ts
if not using_corners:
llcrnrlat = -90.
urcrnrlat = 90.
if lon_0 is not None:
llcrnrlon = lon_0-180.
urcrnrlon = lon_0+180.
else:
llcrnrlon = -180.
urcrnrlon = 180
if projection == 'merc':
# clip plot region to be within -89.99S to 89.99N
# (mercator is singular at poles)
if llcrnrlat < -89.99: llcrnrlat = -89.99
if llcrnrlat > 89.99: llcrnrlat = 89.99
if urcrnrlat < -89.99: urcrnrlat = -89.99
if urcrnrlat > 89.99: urcrnrlat = 89.99
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
if lon_0 is not None:
projparams['lon_0'] = lon_0
else:
projparams['lon_0']=0.5*(llcrnrlon+urcrnrlon)
elif projection == 'rotpole':
if lon_0 is None or o_lon_p is None or o_lat_p is None:
raise ValueError('must specify lon_0,o_lat_p,o_lon_p for rotated pole Basemap')
if width is not None or height is not None:
sys.stdout.write('warning: width and height keywords ignored for %s projection' % _projnames[self.projection])
projparams['lon_0']=lon_0
projparams['o_lon_p']=o_lon_p
projparams['o_lat_p']=o_lat_p
projparams['o_proj']='longlat'
projparams['proj']='ob_tran'
if not using_corners and None in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
raise ValueError('must specify lat/lon values of corners in degrees')
if None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
p = pyproj.Proj(projparams)
llcrnrx = _dg2rad*llcrnrx; llcrnry = _dg2rad*llcrnry
urcrnrx = _dg2rad*urcrnrx; urcrnry = _dg2rad*urcrnry
llcrnrlon, llcrnrlat = p(llcrnrx,llcrnry,inverse=True)
urcrnrlon, urcrnrlat = p(urcrnrx,urcrnry,inverse=True)
self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
else:
raise ValueError(_unsupported_projection % projection)
# initialize proj4
proj = Proj(projparams,self.llcrnrlon,self.llcrnrlat,self.urcrnrlon,self.urcrnrlat)
# make sure axis ticks are suppressed.
self.noticks = suppress_ticks
# map boundary not yet drawn.
self._mapboundarydrawn = False
# make Proj instance a Basemap instance variable.
self.projtran = proj
# copy some Proj attributes.
atts = ['rmajor','rminor','esq','flattening','ellipsoid','projparams']
for att in atts:
self.__dict__[att] = proj.__dict__[att]
# these only exist for geostationary projection.
if hasattr(proj,'_width'):
self.__dict__['_width'] = proj.__dict__['_width']
if hasattr(proj,'_height'):
self.__dict__['_height'] = proj.__dict__['_height']
# spatial reference string (useful for georeferencing output
# images with gdal_translate).
if hasattr(self,'_proj4'):
#self.srs = proj._proj4.srs
self.srs = proj._proj4.pjinitstring
else:
pjargs = []
for key,value in self.projparams.items():
# 'cyl' projection translates to 'eqc' in PROJ.4
if projection == 'cyl' and key == 'proj':
value = 'eqc'
# ignore x_0 and y_0 settings for 'cyl' projection
# (they are not consistent with what PROJ.4 uses)
elif projection == 'cyl' and key in ['x_0','y_0']:
continue
pjargs.append('+'+key+"="+str(value)+' ')
self.srs = ''.join(pjargs)
self.proj4string = self.srs
# set instance variables defining map region.
self.xmin = proj.xmin
self.xmax = proj.xmax
self.ymin = proj.ymin
self.ymax = proj.ymax
if projection == 'cyl':
self.aspect = (self.urcrnrlat-self.llcrnrlat)/(self.urcrnrlon-self.llcrnrlon)
else:
self.aspect = (proj.ymax-proj.ymin)/(proj.xmax-proj.xmin)
if projection in ['geos','ortho','nsper'] and \
None not in [llcrnrx,llcrnry,urcrnrx,urcrnry]:
self.llcrnrx = llcrnrx+0.5*proj.xmax
self.llcrnry = llcrnry+0.5*proj.ymax
self.urcrnrx = urcrnrx+0.5*proj.xmax
self.urcrnry = urcrnry+0.5*proj.ymax
self._fulldisk = False
else:
self.llcrnrx = proj.llcrnrx
self.llcrnry = proj.llcrnry
self.urcrnrx = proj.urcrnrx
self.urcrnry = proj.urcrnry
if self.projection == 'rotpole':
lon0,lat0 = self(0.5*(self.llcrnrx + self.urcrnrx),\
0.5*(self.llcrnry + self.urcrnry),\
inverse=True)
self.projparams['lat_0']=lat0
# if ax == None, pyplot.gca may be used.
self.ax = ax
self.lsmask = None
# This will record hashs of Axes instances.
self._initialized_axes = set()
# set defaults for area_thresh.
self.resolution = resolution
# celestial=True implies resolution=None (no coastlines).
if self.celestial:
self.resolution=None
if area_thresh is None and self.resolution is not None:
if resolution == 'c':
area_thresh = 10000.
elif resolution == 'l':
area_thresh = 1000.
elif resolution == 'i':
area_thresh = 100.
elif resolution == 'h':
area_thresh = 10.
elif resolution == 'f':
area_thresh = 1.
else:
raise ValueError("boundary resolution must be one of 'c','l','i','h' or 'f'")
self.area_thresh = area_thresh
# define map boundary polygon (in lat/lon coordinates)
blons, blats, self._boundarypolyll, self._boundarypolyxy = self._getmapboundary()
self.boundarylats = blats
self.boundarylons = blons
# set min/max lats for projection domain.
if self.projection in _cylproj:
self.latmin = self.llcrnrlat
self.latmax = self.urcrnrlat
self.lonmin = self.llcrnrlon
self.lonmax = self.urcrnrlon
elif self.projection in ['ortho','geos','nsper'] + _pseudocyl:
self.latmin = -90.
self.latmax = 90.
self.lonmin = self.llcrnrlon
self.lonmax = self.urcrnrlon
else:
lons, lats = self.makegrid(1001,1001)
lats = ma.masked_where(lats > 1.e20,lats)
lons = ma.masked_where(lons > 1.e20,lons)
self.latmin = lats.min()
self.latmax = lats.max()
self.lonmin = lons.min()
self.lonmax = lons.max()
NPole = _geoslib.Point(self(0.,90.))
SPole = _geoslib.Point(self(0.,-90.))
if lat_0 is None:
lon_0, lat_0 =\
self(0.5*(self.xmin+self.xmax),
0.5*(self.ymin+self.ymax),inverse=True)
Dateline = _geoslib.Point(self(180.,lat_0))
Greenwich = _geoslib.Point(self(0.,lat_0))
hasNP = NPole.within(self._boundarypolyxy)
hasSP = SPole.within(self._boundarypolyxy)
hasPole = hasNP or hasSP
hasDateline = Dateline.within(self._boundarypolyxy)
hasGreenwich = Greenwich.within(self._boundarypolyxy)
# projection crosses dateline (and not Greenwich or pole).
if not hasPole and hasDateline and not hasGreenwich:
if self.lonmin < 0 and self.lonmax > 0.:
lons = np.where(lons < 0, lons+360, lons)
self.lonmin = lons.min()
self.lonmax = lons.max()
# read in coastline polygons, only keeping those that
# intersect map boundary polygon.
if self.resolution is not None:
self.coastsegs, self.coastpolygontypes =\
self._readboundarydata('gshhs',as_polygons=True)
# reformat for use in matplotlib.patches.Polygon.
self.coastpolygons = []
for seg in self.coastsegs:
x, y = list(zip(*seg))
self.coastpolygons.append((x,y))
# replace coastsegs with line segments (instead of polygons)
self.coastsegs, types =\
self._readboundarydata('gshhs',as_polygons=False)
self.coastsegs = [sg for sg in self.coastsegs if len(sg) > 0]
# create geos Polygon structures for land areas.
# currently only used in is_land method.
self.landpolygons=[]
self.lakepolygons=[]
if self.resolution is not None and len(self.coastpolygons) > 0:
#self.islandinlakepolygons=[]
#self.lakeinislandinlakepolygons=[]
x, y = list(zip(*self.coastpolygons))
for x,y,typ in zip(x,y,self.coastpolygontypes):
b = np.asarray([x,y]).T
if typ == 1: self.landpolygons.append(_geoslib.Polygon(b))
if typ == 2: self.lakepolygons.append(_geoslib.Polygon(b))
#if typ == 3: self.islandinlakepolygons.append(_geoslib.Polygon(b))
#if typ == 4: self.lakeinislandinlakepolygons.append(_geoslib.Polygon(b))
# set __init__'s docstring
__init__.__doc__ = _Basemap_init_doc
def __call__(self,x,y,inverse=False):
"""
Calling a Basemap class instance with the arguments lon, lat will
convert lon/lat (in degrees) to x/y map projection
coordinates (in meters). If optional keyword ``inverse`` is
True (default is False), the inverse transformation from x/y
to lon/lat is performed.
For cylindrical equidistant projection (``cyl``), this
does nothing (i.e. x,y == lon,lat).
For non-cylindrical projections, the inverse transformation
always returns longitudes between -180 and 180 degrees. For
cylindrical projections (self.projection == ``cyl``,
``cea``, ``mill``, ``gall`` or ``merc``)
the inverse transformation will return longitudes between
self.llcrnrlon and self.llcrnrlat.
Input arguments lon, lat can be either scalar floats,
sequences, or numpy arrays.
"""
if self.celestial:
# don't assume center of map is at greenwich
# (only relevant for cyl or pseudo-cyl projections)
if self.projection in _pseudocyl or self.projection in _cylproj:
lon_0=self.projparams['lon_0']
else:
lon_0 = 0.
if self.celestial and not inverse:
try:
x = 2.*lon_0-x
except TypeError:
x = [2*lon_0-xx for xx in x]
if self.projection == 'rotpole' and inverse:
try:
x = _dg2rad*x
except TypeError:
x = [_dg2rad*xx for xx in x]
try:
y = _dg2rad*y
except TypeError:
y = [_dg2rad*yy for yy in y]
xout,yout = self.projtran(x,y,inverse=inverse)
if self.celestial and inverse:
try:
xout = -2.*lon_0-xout
except:
xout = [-2.*lon_0-xx for xx in xout]
if self.projection == 'rotpole' and not inverse:
try:
xout = _rad2dg*xout
xout = np.where(xout < 0., xout+360, xout)
except TypeError:
xout = [_rad2dg*xx for xx in xout]
xout = [xx+360. if xx < 0 else xx for xx in xout]
try:
yout = _rad2dg*yout
except TypeError:
yout = [_rad2dg*yy for yy in yout]
return xout,yout
[docs] def makegrid(self,nx,ny,returnxy=False):
"""
return arrays of shape (ny,nx) containing lon,lat coordinates of
an equally spaced native projection grid.
If ``returnxy = True``, the x,y values of the grid are returned also.
"""
return self.projtran.makegrid(nx,ny,returnxy=returnxy)
def _readboundarydata(self,name,as_polygons=False):
"""
read boundary data, clip to map projection region.
"""
# only gshhs coastlines can be polygons.
if name != 'gshhs': as_polygons=False
try:
bdatfile = open(os.path.join(basemap_datadir,name+'_'+self.resolution+'.dat'),'rb')
bdatmetafile = open(os.path.join(basemap_datadir,name+'meta_'+self.resolution+'.dat'),'r')
except:
raise IOError(
"Unable to open boundary dataset file. Only the 'crude', 'low' "
"and 'intermediate' resolution datasets are installed by default. "
"If you are requesting a 'high' or 'full' resolution dataset, "
"you need to install the `basemap-data-hires` package")
polygons = []
polygon_types = []
# coastlines are polygons, other boundaries are line segments.
if name == 'gshhs':
Shape = _geoslib.Polygon
else:
Shape = _geoslib.LineString
# see if map projection region polygon contains a pole.
NPole = _geoslib.Point(self(0.,90.))
SPole = _geoslib.Point(self(0.,-90.))
boundarypolyxy = self._boundarypolyxy
boundarypolyll = self._boundarypolyll
hasNP = NPole.within(boundarypolyxy)
hasSP = SPole.within(boundarypolyxy)
containsPole = hasNP or hasSP
# these projections cannot cross pole.
if containsPole and\
self.projection in _cylproj + _pseudocyl + ['geos']:
raise ValueError('%s projection cannot cross pole'%(self.projection))
# make sure some projections have has containsPole=True
# we will compute the intersections in stereographic
# coordinates, then transform back. This is
# because these projections are only defined on a hemisphere, and
# some boundary features (like Eurasia) would be undefined otherwise.
tostere =\
['omerc','ortho','gnom','nsper','nplaea','npaeqd','splaea','spaeqd']
if self.projection in tostere and name == 'gshhs':
containsPole = True
lon_0=self.projparams['lon_0']
lat_0=self.projparams['lat_0']
re = self.projparams['R']
# center of stereographic projection restricted to be
# nearest one of 6 points on the sphere (every 90 deg lat/lon).
lon0 = 90.*(np.around(lon_0/90.))
lat0 = 90.*(np.around(lat_0/90.))
if np.abs(int(lat0)) == 90: lon0=0.
maptran = pyproj.Proj(proj='stere',lon_0=lon0,lat_0=lat0,R=re)
# boundary polygon for ortho/gnom/nsper projection
# in stereographic coordinates.
b = self._boundarypolyll.boundary
blons = b[:,0]; blats = b[:,1]
b[:,0], b[:,1] = maptran(blons, blats)
boundarypolyxy = _geoslib.Polygon(b)
for line in bdatmetafile:
linesplit = line.split()
area = float(linesplit[1])
south = float(linesplit[3])
north = float(linesplit[4])
crossdatelineE=False; crossdatelineW=False
if name == 'gshhs':
id = linesplit[7]
if id.endswith('E'):
crossdatelineE = True
elif id.endswith('W'):
crossdatelineW = True
# make sure south/north limits of dateline crossing polygons
# (Eurasia) are the same, since they will be merged into one.
# (this avoids having one filtered out and not the other).
if crossdatelineE:
south_save=south
north_save=north
if crossdatelineW:
south=south_save
north=north_save
if area < 0.: area = 1.e30
useit = self.latmax>=south and self.latmin<=north and area>self.area_thresh
if useit:
typ = int(linesplit[0])
npts = int(linesplit[2])
offsetbytes = int(linesplit[5])
bytecount = int(linesplit[6])
bdatfile.seek(offsetbytes,0)
# read in binary string convert into an npts by 2
# numpy array (first column is lons, second is lats).
polystring = bdatfile.read(bytecount)
# binary data is little endian.
b = np.array(np.frombuffer(polystring,dtype='<f4'),'f8')
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[3:,:]
b2 = b.copy()
poly = Shape(b)
# if map boundary polygon is a valid one in lat/lon
# coordinates (i.e. it does not contain either pole),
# the intersections of the boundary geometries
# and the map projection region can be computed before
# transforming the boundary geometry to map projection
# coordinates (this saves time, especially for small map
# regions and high-resolution boundary geometries).
if not containsPole:
# close Antarctica.
if name == 'gshhs' and south < -89:
lons2 = b[:,0]
lats = b[:,1]
lons1 = lons2 - 360.
lons3 = lons2 + 360.
lons = lons1.tolist()+lons2.tolist()+lons3.tolist()
lats = lats.tolist()+lats.tolist()+lats.tolist()
lonstart,latstart = lons[0], lats[0]
lonend,latend = lons[-1], lats[-1]
lons.insert(0,lonstart)
lats.insert(0,-90.)
lons.append(lonend)
lats.append(-90.)
b = np.empty((len(lons),2),np.float64)
b[:,0] = lons; b[:,1] = lats
poly = Shape(b)
if not poly.is_valid(): poly=poly.fix()
# if polygon instersects map projection
# region, process it.
if poly.intersects(boundarypolyll):
if name != 'gshhs' or as_polygons:
geoms = poly.intersection(boundarypolyll)
else:
# convert polygons to line segments
poly = _geoslib.LineString(poly.boundary)
geoms = poly.intersection(boundarypolyll)
# iterate over geometries in intersection.
for psub in geoms:
b = psub.boundary
blons = b[:,0]; blats = b[:,1]
bx, by = self(blons, blats)
polygons.append(list(zip(bx,by)))
polygon_types.append(typ)
else:
# create duplicate polygons shifted by -360 and +360
# (so as to properly treat polygons that cross
# Greenwich meridian).
b2[:,0] = b[:,0]-360
poly1 = Shape(b2)
b2[:,0] = b[:,0]+360
poly2 = Shape(b2)
polys = [poly1,poly,poly2]
for poly in polys:
# try to fix "non-noded intersection" errors.
if not poly.is_valid(): poly=poly.fix()
# if polygon instersects map projection
# region, process it.
if poly.intersects(boundarypolyll):
if name != 'gshhs' or as_polygons:
geoms = poly.intersection(boundarypolyll)
else:
# convert polygons to line segments
# note: use fix method here or Eurasia
# line segments sometimes disappear.
poly = _geoslib.LineString(poly.fix().boundary)
geoms = poly.intersection(boundarypolyll)
# iterate over geometries in intersection.
for psub in geoms:
b = psub.boundary
blons = b[:,0]; blats = b[:,1]
# transformation from lat/lon to
# map projection coordinates.
bx, by = self(blons, blats)
if not as_polygons or len(bx) > 4:
polygons.append(list(zip(bx,by)))
polygon_types.append(typ)
# if map boundary polygon is not valid in lat/lon
# coordinates, compute intersection between map
# projection region and boundary geometries in map
# projection coordinates.
else:
# transform coordinates from lat/lon
# to map projection coordinates.
# special case for ortho/gnom/nsper, compute coastline polygon
# vertices in stereographic coords.
if name == 'gshhs' and as_polygons and self.projection in tostere:
b[:,0], b[:,1] = maptran(b[:,0], b[:,1])
else:
b[:,0], b[:,1] = self(b[:,0], b[:,1])
goodmask = np.logical_and(b[:,0]<1.e20,b[:,1]<1.e20)
# if less than two points are valid in
# map proj coords, skip this geometry.
if np.sum(goodmask) <= 1: continue
if name != 'gshhs' or (name == 'gshhs' and not as_polygons):
# if not a polygon,
# just remove parts of geometry that are undefined
# in this map projection.
bx = np.compress(goodmask, b[:,0])
by = np.compress(goodmask, b[:,1])
# split coastline segments that jump across entire plot.
xd = (bx[1:]-bx[0:-1])**2
yd = (by[1:]-by[0:-1])**2
dist = np.sqrt(xd+yd)
split = dist > 0.1*(self.xmax-self.xmin)
if np.sum(split) and self.projection not in _cylproj:
ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist()
iprev = 0
ind.append(len(xd))
for i in ind:
# don't add empty lists.
if len(list(range(iprev,i))):
polygons.append(list(zip(bx[iprev:i],by[iprev:i])))
iprev = i
else:
polygons.append(list(zip(bx,by)))
polygon_types.append(typ)
continue
# create a GEOS geometry object.
if name == 'gshhs' and not as_polygons:
# convert polygons to line segments
poly = _geoslib.LineString(poly.boundary)
else:
# this is a workaround to avoid
# GEOS_ERROR: CGAlgorithmsDD::orientationIndex encountered NaN/Inf numbers
b[np.isposinf(b)] = 1e20
b[np.isneginf(b)] = -1e20
poly = Shape(b)
# this is a workaround to avoid
# "GEOS_ERROR: TopologyException:
# found non-noded intersection between ..."
if not poly.is_valid(): poly=poly.fix()
# if geometry instersects map projection
# region, and doesn't have any invalid points, process it.
if goodmask.all() and poly.intersects(boundarypolyxy):
# if geometry intersection calculation fails,
# just move on.
try:
geoms = poly.intersection(boundarypolyxy)
except:
continue
# iterate over geometries in intersection.
for psub in geoms:
b = psub.boundary
# if projection in ['ortho','gnom','nsper'],
# transform polygon from stereographic
# to ortho/gnom/nsper coordinates.
if self.projection in tostere:
# if coastline polygon covers more than 99%
# of map region for fulldisk projection,
# it's probably bogus, so skip it.
#areafrac = psub.area()/boundarypolyxy.area()
#if self.projection == ['ortho','nsper']:
# if name == 'gshhs' and\
# self._fulldisk and\
# areafrac > 0.99: continue
# inverse transform from stereographic
# to lat/lon.
b[:,0], b[:,1] = maptran(b[:,0], b[:,1], inverse=True)
# orthographic/gnomonic/nsper.
b[:,0], b[:,1]= self(b[:,0], b[:,1])
if not as_polygons or len(b) > 4:
polygons.append(list(zip(b[:,0],b[:,1])))
polygon_types.append(typ)
bdatfile.close()
bdatmetafile.close()
return polygons, polygon_types
def _getmapboundary(self):
"""
create map boundary polygon (in lat/lon and x/y coordinates)
"""
nx = 100; ny = 100
maptran = self
if self.projection in ['ortho','geos','nsper']:
# circular region.
thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
rminor = self._height
rmajor = self._width
x = rmajor*np.cos(thetas) + rmajor
y = rminor*np.sin(thetas) + rminor
b = np.empty((len(x),2),np.float64)
b[:,0]=x; b[:,1]=y
boundaryxy = _geoslib.Polygon(b)
# compute proj instance for full disk, if necessary.
if not self._fulldisk:
projparms = self.projparams.copy()
del projparms['x_0']
del projparms['y_0']
if self.projection == 'ortho':
llcrnrx = -self.rmajor
llcrnry = -self.rmajor
urcrnrx = -llcrnrx
urcrnry = -llcrnry
else:
llcrnrx = -self._width
llcrnry = -self._height
urcrnrx = -llcrnrx
urcrnry = -llcrnry
projparms['x_0']=-llcrnrx
projparms['y_0']=-llcrnry
maptran = pyproj.Proj(projparms)
elif self.projection == 'aeqd' and self._fulldisk:
# circular region.
thetas = np.linspace(0.,2.*np.pi,2*nx*ny)[:-1]
rminor = self._height
rmajor = self._width
x = rmajor*np.cos(thetas) + rmajor
y = rminor*np.sin(thetas) + rminor
b = np.empty((len(x),2),np.float64)
b[:,0]=x; b[:,1]=y
boundaryxy = _geoslib.Polygon(b)
elif self.projection in _pseudocyl:
nx = 10*nx; ny = 10*ny
# quasi-elliptical region.
lon_0 = self.projparams['lon_0']
# left side
lats1 = np.linspace(-89.9999,89.9999,ny).tolist()
lons1 = len(lats1)*[lon_0-179.9]
# top.
lons2 = np.linspace(lon_0-179.9,lon_0+179.9,nx).tolist()
lats2 = len(lons2)*[89.9999]
# right side
lats3 = np.linspace(89.9999,-89.9999,ny).tolist()
lons3 = len(lats3)*[lon_0+179.9]
# bottom.
lons4 = np.linspace(lon_0+179.9,lon_0-179.9,nx).tolist()
lats4 = len(lons4)*[-89.9999]
lons = np.array(lons1+lons2+lons3+lons4,np.float64)
lats = np.array(lats1+lats2+lats3+lats4,np.float64)
x, y = maptran(lons,lats)
b = np.empty((len(x),2),np.float64)
b[:,0]=x; b[:,1]=y
boundaryxy = _geoslib.Polygon(b)
else: # all other projections are rectangular.
nx = 100*nx; ny = 100*ny
# left side (x = xmin, ymin <= y <= ymax)
yy = np.linspace(self.ymin, self.ymax, ny)[:-1]
x = len(yy)*[self.xmin]; y = yy.tolist()
# top (y = ymax, xmin <= x <= xmax)
xx = np.linspace(self.xmin, self.xmax, nx)[:-1]
x = x + xx.tolist()
y = y + len(xx)*[self.ymax]
# right side (x = xmax, ymin <= y <= ymax)
yy = np.linspace(self.ymax, self.ymin, ny)[:-1]
x = x + len(yy)*[self.xmax]; y = y + yy.tolist()
# bottom (y = ymin, xmin <= x <= xmax)
xx = np.linspace(self.xmax, self.xmin, nx)[:-1]
x = x + xx.tolist()
y = y + len(xx)*[self.ymin]
x = np.array(x,np.float64)
y = np.array(y,np.float64)
b = np.empty((4,2),np.float64)
b[:,0]=[self.xmin,self.xmin,self.xmax,self.xmax]
b[:,1]=[self.ymin,self.ymax,self.ymax,self.ymin]
boundaryxy = _geoslib.Polygon(b)
if self.projection in _cylproj:
# make sure map boundary doesn't quite include pole.
if self.urcrnrlat > 89.9999:
urcrnrlat = 89.9999
else:
urcrnrlat = self.urcrnrlat
if self.llcrnrlat < -89.9999:
llcrnrlat = -89.9999
else:
llcrnrlat = self.llcrnrlat
lons = [self.llcrnrlon, self.llcrnrlon, self.urcrnrlon, self.urcrnrlon]
lats = [llcrnrlat, urcrnrlat, urcrnrlat, llcrnrlat]
self.boundarylonmin = min(lons)
self.boundarylonmax = max(lons)
x, y = self(lons, lats)
b = np.empty((len(x),2),np.float64)
b[:,0]=x; b[:,1]=y
boundaryxy = _geoslib.Polygon(b)
else:
if self.projection not in _pseudocyl:
lons, lats = maptran(x,y,inverse=True)
# fix lons so there are no jumps.
n = 1
lonprev = lons[0]
for lon,lat in zip(lons[1:],lats[1:]):
if np.abs(lon-lonprev) > 90.:
if lonprev < 0:
lon = lon - 360.
else:
lon = lon + 360
lons[n] = lon
lonprev = lon
n = n + 1
self.boundarylonmin = lons.min()
self.boundarylonmax = lons.max()
# for circular full disk projections where boundary is
# a latitude circle, set boundarylonmax and boundarylonmin
# to cover entire world (so parallels will be drawn).
if self._fulldisk and \
np.abs(self.boundarylonmax-self.boundarylonmin) < 1.:
self.boundarylonmin = -180.
self.boundarylonmax = 180.
b = np.empty((len(lons),2),np.float64)
b[:,0] = lons; b[:,1] = lats
boundaryll = _geoslib.Polygon(b)
return lons, lats, boundaryll, boundaryxy
[docs] def drawmapboundary(self,color='k',linewidth=1.0,fill_color=None,\
zorder=None,ax=None):
"""
draw boundary around map projection region, optionally
filling interior of region.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
linewidth line width for boundary (default 1.)
color color of boundary line (default black)
fill_color fill the map region background with this
color (default is to fill with axis
background color). If set to the string
'none', no filling is done.
zorder sets the zorder for filling map background
(default 0).
ax axes instance to use
(default None, use default axes instance).
============== ====================================================
returns matplotlib.collections.PatchCollection representing map boundary.
"""
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# if no fill_color given, use axes background color.
# if fill_color is string 'none', really don't fill.
if fill_color is None:
mpl_version = tuple(map(int, mpl.__version__.split(".")[:2]))
if mpl_version >= (2, 0):
fill_color = ax.get_facecolor()
else:
fill_color = ax.get_axis_bgcolor()
elif fill_color == 'none' or fill_color == 'None':
fill_color = None
limb = None
if self.projection in ['ortho','geos','nsper'] or (self.projection=='aeqd' and\
self._fulldisk):
limb = Ellipse((self._width,self._height),2.*self._width,2.*self._height)
if self.projection in ['ortho','geos','nsper','aeqd'] and self._fulldisk:
# elliptical region.
ax.set_frame_on(False)
elif self.projection in _pseudocyl: # elliptical region.
ax.set_frame_on(False)
nx = 100; ny = 100
if self.projection == 'vandg':
nx = 10*nx; ny = 10*ny
# quasi-elliptical region.
lon_0 = self.projparams['lon_0']
# left side
lats1 = np.linspace(-89.9999,89.99999,ny).tolist()
lons1 = len(lats1)*[lon_0-179.9]
# top.
lons2 = np.linspace(lon_0-179.9999,lon_0+179.9999,nx).tolist()
lats2 = len(lons2)*[89.9999]
# right side
lats3 = np.linspace(89.9999,-89.9999,ny).tolist()
lons3 = len(lats3)*[lon_0+179.9999]
# bottom.
lons4 = np.linspace(lon_0+179.9999,lon_0-179.9999,nx).tolist()
lats4 = len(lons4)*[-89.9999]
lons = np.array(lons1+lons2+lons3+lons4,np.float64)
lats = np.array(lats1+lats2+lats3+lats4,np.float64)
x, y = self(lons,lats)
xy = list(zip(x,y))
limb = Polygon(xy)
elif self.round:
ax.set_frame_on(False)
limb = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
radius=0.5*(self.xmax-self.xmin),fc='none')
else: # all other projections are rectangular.
ax.set_frame_on(True)
for spine in ax.spines.values():
spine.set_linewidth(linewidth)
spine.set_edgecolor(color)
if zorder is not None:
spine.set_zorder(zorder)
if self.projection not in ['geos','ortho','nsper']:
limb = ax.patch
if limb is not None:
if limb is not ax.patch:
ax.add_patch(limb)
self._mapboundarydrawn = limb
if fill_color is None:
limb.set_fill(False)
else:
limb.set_facecolor(fill_color)
limb.set_zorder(0)
limb.set_edgecolor(color)
limb.set_linewidth(linewidth)
if zorder is not None:
limb.set_zorder(zorder)
limb.set_clip_on(True)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
return limb
[docs] def fillcontinents(self,color='0.8',lake_color=None,ax=None,zorder=None,alpha=None):
"""
Fill continents.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
color color to fill continents (default gray).
lake_color color to fill inland lakes (default axes background).
ax axes instance (overrides default axes instance).
zorder sets the zorder for the continent polygons (if not
specified, uses default zorder for a Polygon patch).
Set to zero if you want to paint over the filled
continents).
alpha sets alpha transparency for continent polygons
============== ====================================================
After filling continents, lakes are re-filled with
axis background color.
returns a list of matplotlib.patches.Polygon objects.
"""
if self.resolution is None:
raise AttributeError('there are no boundary datasets associated with this Basemap instance')
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# get axis background color.
mpl_version = tuple(map(int, mpl.__version__.split(".")[:2]))
if mpl_version >= (2, 0):
axisbgc = ax.get_facecolor()
else:
axisbgc = ax.get_axis_bgcolor()
npoly = 0
polys = []
for x,y in self.coastpolygons:
xa = np.array(x,np.float32)
ya = np.array(y,np.float32)
# check to see if all four corners of domain in polygon (if so,
# don't draw since it will just fill in the whole map).
# ** turn this off for now since it prevents continents that
# fill the whole map from being filled **
#delx = 10; dely = 10
#if self.projection in ['cyl']:
# delx = 0.1
# dely = 0.1
#test1 = np.fabs(xa-self.urcrnrx) < delx
#test2 = np.fabs(xa-self.llcrnrx) < delx
#test3 = np.fabs(ya-self.urcrnry) < dely
#test4 = np.fabs(ya-self.llcrnry) < dely
#hasp1 = np.sum(test1*test3)
#hasp2 = np.sum(test2*test3)
#hasp4 = np.sum(test2*test4)
#hasp3 = np.sum(test1*test4)
#if not hasp1 or not hasp2 or not hasp3 or not hasp4:
if 1:
xy = list(zip(xa.tolist(),ya.tolist()))
if self.coastpolygontypes[npoly] not in [2,4]:
poly = Polygon(xy,facecolor=color,edgecolor=color,linewidth=0)
else: # lakes filled with background color by default
if lake_color is None:
poly = Polygon(xy,facecolor=axisbgc,edgecolor=axisbgc,linewidth=0)
else:
poly = Polygon(xy,facecolor=lake_color,edgecolor=lake_color,linewidth=0)
if zorder is not None:
poly.set_zorder(zorder)
if alpha is not None:
poly.set_alpha(alpha)
ax.add_patch(poly)
polys.append(poly)
npoly = npoly + 1
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip continent polygons to map limbs
polys,c = self._cliplimb(ax,polys)
return polys
def _cliplimb(self,ax,coll):
if not self._mapboundarydrawn:
return coll, None
c = self._mapboundarydrawn
if c not in ax.patches:
p = ax.add_patch(c)
#p.set_clip_on(False)
try:
coll.set_clip_path(c)
except:
for item in coll:
item.set_clip_path(c)
return coll,c
[docs] def drawcoastlines(self,linewidth=1.,linestyle='solid',color='k',antialiased=1,ax=None,zorder=None):
"""
Draw coastlines.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
linewidth coastline width (default 1.)
linestyle coastline linestyle (default solid)
color coastline color (default black)
antialiased antialiasing switch for coastlines (default True).
ax axes instance (overrides default axes instance)
zorder sets the zorder for the coastlines (if not specified,
uses default zorder for
matplotlib.patches.LineCollections).
============== ====================================================
returns a matplotlib.patches.LineCollection object.
"""
if self.resolution is None:
raise AttributeError('there are no boundary datasets associated with this Basemap instance')
# get current axes instance (if none specified).
ax = ax or self._check_ax()
coastlines = LineCollection(self.coastsegs,antialiaseds=(antialiased,))
coastlines.set_color(color)
coastlines.set_linestyle(linestyle)
coastlines.set_linewidth(linewidth)
coastlines.set_label('_nolabel_')
if zorder is not None:
coastlines.set_zorder(zorder)
ax.add_collection(coastlines)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
coastlines,c = self._cliplimb(ax,coastlines)
return coastlines
[docs] def drawcountries(self,linewidth=0.5,linestyle='solid',color='k',antialiased=1,ax=None,zorder=None):
"""
Draw country boundaries.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
linewidth country boundary line width (default 0.5)
linestyle coastline linestyle (default solid)
color country boundary line color (default black)
antialiased antialiasing switch for country boundaries (default
True).
ax axes instance (overrides default axes instance)
zorder sets the zorder for the country boundaries (if not
specified uses default zorder for
matplotlib.patches.LineCollections).
============== ====================================================
returns a matplotlib.patches.LineCollection object.
"""
if self.resolution is None:
raise AttributeError('there are no boundary datasets associated with this Basemap instance')
# read in country line segments, only keeping those that
# intersect map boundary polygon.
if not hasattr(self,'cntrysegs'):
self.cntrysegs, types = self._readboundarydata('countries')
# get current axes instance (if none specified).
ax = ax or self._check_ax()
countries = LineCollection(self.cntrysegs,antialiaseds=(antialiased,))
countries.set_color(color)
countries.set_linestyle(linestyle)
countries.set_linewidth(linewidth)
countries.set_label('_nolabel_')
if zorder is not None:
countries.set_zorder(zorder)
ax.add_collection(countries)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip countries to map limbs
countries,c = self._cliplimb(ax,countries)
return countries
[docs] def drawstates(self,linewidth=0.5,linestyle='solid',color='k',antialiased=1,ax=None,zorder=None):
"""
Draw state boundaries in Americas.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
linewidth state boundary line width (default 0.5)
linestyle coastline linestyle (default solid)
color state boundary line color (default black)
antialiased antialiasing switch for state boundaries
(default True).
ax axes instance (overrides default axes instance)
zorder sets the zorder for the state boundaries (if not
specified, uses default zorder for
matplotlib.patches.LineCollections).
============== ====================================================
returns a matplotlib.patches.LineCollection object.
"""
if self.resolution is None:
raise AttributeError('there are no boundary datasets associated with this Basemap instance')
# read in state line segments, only keeping those that
# intersect map boundary polygon.
if not hasattr(self,'statesegs'):
self.statesegs, types = self._readboundarydata('states')
# get current axes instance (if none specified).
ax = ax or self._check_ax()
states = LineCollection(self.statesegs,antialiaseds=(antialiased,))
states.set_color(color)
states.set_linestyle(linestyle)
states.set_linewidth(linewidth)
states.set_label('_nolabel_')
if zorder is not None:
states.set_zorder(zorder)
ax.add_collection(states)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip states to map limbs
states,c = self._cliplimb(ax,states)
return states
[docs] def drawcounties(self,linewidth=0.1,linestyle='solid',color='k',antialiased=1,
facecolor='none',ax=None,zorder=None,drawbounds=False):
"""
Draw county boundaries in US. The county boundary shapefile
originates with the NOAA Coastal Geospatial Data Project
(http://coastalgeospatial.noaa.gov/data_gis.html).
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
linewidth county boundary line width (default 0.1)
linestyle coastline linestyle (default solid)
color county boundary line color (default black)
facecolor fill color of county (default is no fill)
antialiased antialiasing switch for county boundaries
(default True).
ax axes instance (overrides default axes instance)
zorder sets the zorder for the county boundaries (if not
specified, uses default zorder for
matplotlib.patches.LineCollections).
============== ====================================================
returns a matplotlib.patches.LineCollection object.
"""
ax = ax or self._check_ax()
gis_file = os.path.join(basemap_datadir,'UScounties')
county_info = self.readshapefile(gis_file,'counties',\
default_encoding='latin-1',drawbounds=drawbounds)
counties = [coords for coords in self.counties]
counties = PolyCollection(counties, antialiaseds=(antialiased,))
counties.set_linestyle(linestyle)
counties.set_linewidth(linewidth)
counties.set_edgecolor(color)
counties.set_facecolor(facecolor)
counties.set_label('counties')
if zorder:
counties.set_zorder(zorder)
ax.add_collection(counties)
return counties
[docs] def drawrivers(self,linewidth=0.5,linestyle='solid',color='k',antialiased=1,ax=None,zorder=None):
"""
Draw major rivers.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
linewidth river boundary line width (default 0.5)
linestyle coastline linestyle (default solid)
color river boundary line color (default black)
antialiased antialiasing switch for river boundaries (default
True).
ax axes instance (overrides default axes instance)
zorder sets the zorder for the rivers (if not
specified uses default zorder for
matplotlib.patches.LineCollections).
============== ====================================================
returns a matplotlib.patches.LineCollection object.
"""
if self.resolution is None:
raise AttributeError('there are no boundary datasets associated with this Basemap instance')
# read in river line segments, only keeping those that
# intersect map boundary polygon.
if not hasattr(self,'riversegs'):
self.riversegs, types = self._readboundarydata('rivers')
# get current axes instance (if none specified).
ax = ax or self._check_ax()
rivers = LineCollection(self.riversegs,antialiaseds=(antialiased,))
rivers.set_color(color)
rivers.set_linestyle(linestyle)
rivers.set_linewidth(linewidth)
rivers.set_label('_nolabel_')
if zorder is not None:
rivers.set_zorder(zorder)
ax.add_collection(rivers)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip rivers to map limbs
rivers,c = self._cliplimb(ax,rivers)
return rivers
[docs] 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
[docs] def readshapefile(self,shapefile,name,drawbounds=True,zorder=None,
linewidth=0.5,color='k',antialiased=1,ax=None,
default_encoding='utf-8',encoding_errors='strict'):
"""
Read in shape file, optionally draw boundaries on map.
.. note::
- Assumes shapes are 2D
- only works for Point, MultiPoint, Polyline and Polygon shapes.
- vertices/points must be in geographic (lat/lon) coordinates.
Mandatory Arguments:
.. tabularcolumns:: |l|L|
================ ====================================================
Argument Description
================ ====================================================
shapefile path to shapefile components. Example:
shapefile='/home/jeff/esri/world_borders' assumes
that world_borders.shp, world_borders.shx and
world_borders.dbf live in /home/jeff/esri.
name name for Basemap attribute to hold the shapefile
vertices or points in map projection
coordinates. Class attribute name+'_info' is a list
of dictionaries, one for each shape, containing
attributes of each shape from dbf file, For
example, if name='counties', self.counties
will be a list of x,y vertices for each shape in
map projection coordinates and self.counties_info
will be a list of dictionaries with shape
attributes. Rings in individual Polygon
shapes are split out into separate polygons, and
additional keys 'RINGNUM' and 'SHAPENUM' are added
to the shape attribute dictionary.
================ ====================================================
The following optional keyword arguments are only relevant for Polyline
and Polygon shape types, for Point and MultiPoint shapes they are
ignored.
.. tabularcolumns:: |l|L|
================ =====================================================
Keyword Description
================ =====================================================
drawbounds draw boundaries of shapes (default True).
zorder shape boundary zorder (if not specified,
default for 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)
default_encoding encoding used to parse properties from .dbf files
(default utf-8)
encoding_errors encoding error handling (default strict), other
possible values: ignore, replace and backslashreplace
================ =====================================================
A tuple (num_shapes, type, min, max) containing shape file info
is returned.
num_shapes is the number of shapes, type is the type code (one of
the SHPT* constants defined in the shapelib module, see
http://shapelib.maptools.org/shp_api.html) and min and
max are 4-element lists with the minimum and maximum values of the
vertices. If ``drawbounds=True`` a
matplotlib.patches.LineCollection object is appended to the tuple.
"""
import shapefile as shp
from shapefile import Reader
shp.default_encoding = default_encoding
if not os.path.exists('%s.shp'%shapefile):
raise IOError('cannot locate %s.shp'%shapefile)
if not os.path.exists('%s.shx'%shapefile):
raise IOError('cannot locate %s.shx'%shapefile)
if not os.path.exists('%s.dbf'%shapefile):
raise IOError('cannot locate %s.dbf'%shapefile)
# open shapefile, read vertices for each object, convert
# to map projection coordinates (only works for 2D shape types).
try:
shf = Reader(shapefile, encoding=default_encoding,
encodingErrors=encoding_errors)
except:
raise IOError('error reading shapefile %s.shp' % shapefile)
fields = shf.fields
coords = []; attributes = []
msg = " ".join([
"shapefile must have lat/lon vertices - it looks like this one",
"has vertices in map projection coordinates. You can convert the",
"shapefile to geographic coordinates using the shpproj utility",
"from the shapelib tools (http://shapelib.maptools.org/shapelib-tools.html)"])
shptype = shf.shapes()[0].shapeType
bbox = shf.bbox.tolist()
info = (shf.numRecords,shptype,bbox[0:2]+[0.,0.],bbox[2:]+[0.,0.])
npoly = 0
for shprec in shf.shapeRecords():
shp = shprec.shape; rec = shprec.record
npoly = npoly + 1
if shptype != shp.shapeType:
raise ValueError('readshapefile can only handle a single shape type per file')
if shptype not in [1,3,5,8]:
raise ValueError('readshapefile can only handle 2D shape types')
verts = shp.points
if shptype in [1,8]: # a Point or MultiPoint shape.
lons, lats = list(zip(*verts))
if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
raise ValueError(msg)
# if latitude is slightly greater than 90, truncate to 90
lats = [max(min(lat, 90.0), -90.0) for lat in lats]
if len(verts) > 1: # MultiPoint
x,y = self(lons, lats)
coords.append(list(zip(x,y)))
else: # single Point
x,y = self(lons[0], lats[0])
coords.append((x,y))
attdict={}
for r,key in zip(rec,fields[1:]):
attdict[key[0]]=r
attributes.append(attdict)
else: # a Polyline or Polygon shape.
parts = shp.parts.tolist()
ringnum = 0
for indx1,indx2 in zip(parts,parts[1:]+[len(verts)]):
ringnum = ringnum + 1
lons, lats = list(zip(*verts[indx1:indx2]))
if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
raise ValueError(msg)
# if latitude is slightly greater than 90, truncate to 90
lats = [max(min(lat, 90.0), -90.0) for lat in lats]
x, y = self(lons, lats)
coords.append(list(zip(x,y)))
attdict={}
for r,key in zip(rec,fields[1:]):
attdict[key[0]]=r
# add information about ring number to dictionary.
attdict['RINGNUM'] = ringnum
attdict['SHAPENUM'] = npoly
attributes.append(attdict)
# draw shape boundaries for polylines, polygons using LineCollection.
if shptype not in [1,8] and drawbounds:
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# make LineCollections for each polygon.
lines = LineCollection(coords,antialiaseds=(antialiased,))
lines.set_color(color)
lines.set_linewidth(linewidth)
lines.set_label('_nolabel_')
if zorder is not None:
lines.set_zorder(zorder)
ax.add_collection(lines)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip boundaries to map limbs
lines,c = self._cliplimb(ax,lines)
info = info + (lines,)
self.__dict__[name]=coords
self.__dict__[name+'_info']=attributes
return info
[docs] def drawparallels(self,circles,color='k',textcolor='k',linewidth=1.,zorder=None, \
dashes=[1,1],labels=[0,0,0,0],labelstyle=None, \
fmt='%g',xoffset=None,yoffset=None,ax=None,latmax=None,
**text_kwargs):
r"""
Draw and label parallels (latitude lines) for values (in degrees)
given in the sequence ``circles``.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
color color to draw parallels (default black).
textcolor color to draw labels (default black).
linewidth line width for parallels (default 1.)
zorder sets the zorder for parallels (if not specified,
uses default zorder for matplotlib.lines.Line2D
objects).
dashes dash pattern for parallels (default [1,1], i.e.
1 pixel on, 1 pixel off).
labels list of 4 values (default [0,0,0,0]) that control
whether parallels are labelled where they intersect
the left, right, top or bottom of the plot. For
example labels=[1,0,0,1] will cause parallels
to be labelled where they intersect the left and
and bottom of the plot, but not the right and top.
labelstyle if set to "+/-", north and south latitudes are
labelled with "+" and "-", otherwise they are
labelled with "N" and "S".
fmt a format string to format the parallel labels
(default '%g') **or** a function that takes a
latitude value in degrees as it's only argument
and returns a formatted string.
xoffset label offset from edge of map in x-direction
(default is 0.01 times width of map in map
projection coordinates).
yoffset label offset from edge of map in y-direction
(default is 0.01 times height of map in map
projection coordinates).
ax axes instance (overrides default axes instance)
latmax absolute value of latitude to which meridians are drawn
(default is 80).
\**text_kwargs additional keyword arguments controlling text
for labels that are passed on to
the text method of the axes instance (see
matplotlib.pyplot.text documentation).
============== ====================================================
returns a dictionary whose keys are the parallel values, and
whose values are tuples containing lists of the
matplotlib.lines.Line2D and matplotlib.text.Text instances
associated with each parallel. Deleting an item from the
dictionary removes the corresponding parallel from the plot.
"""
text_kwargs['color']=textcolor # pass textcolor kwarg on to ax.text
# if celestial=True, don't use "N" and "S" labels.
if labelstyle is None and self.celestial:
labelstyle="+/-"
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# don't draw meridians past latmax, always draw parallel at latmax.
if latmax is None: latmax = 80.
# offset for labels.
if yoffset is None:
yoffset = (self.urcrnry-self.llcrnry)/100.
if self.aspect > 1:
yoffset = self.aspect*yoffset
else:
yoffset = yoffset/self.aspect
if xoffset is None:
xoffset = (self.urcrnrx-self.llcrnrx)/100.
if self.projection in _cylproj + _pseudocyl:
lons = np.linspace(self.llcrnrlon, self.urcrnrlon, 10001)
elif self.projection in ['tmerc']:
lon_0 = self.projparams['lon_0']
# tmerc only defined within +/- 90 degrees of lon_0
lons = np.linspace(lon_0-90,lon_0+90,100001)
else:
lonmin = self.boundarylonmin; lonmax = self.boundarylonmax
lons = np.linspace(lonmin, lonmax, 10001)
# make sure latmax degree parallel is drawn if projection not merc or cyl or miller
try:
circlesl = list(circles)
except:
circlesl = circles
if self.projection not in _cylproj + _pseudocyl:
if max(circlesl) > 0 and latmax not in circlesl:
circlesl.append(latmax)
if min(circlesl) < 0 and -latmax not in circlesl:
circlesl.append(-latmax)
xdelta = 0.01*(self.xmax-self.xmin)
ydelta = 0.01*(self.ymax-self.ymin)
linecolls = {}
for circ in circlesl:
lats = circ*np.ones(len(lons),np.float32)
x,y = self(lons,lats)
# remove points outside domain.
# leave a little slop around edges (3*xdelta)
# don't really know why, but this appears to be needed to
# or lines sometimes don't reach edge of plot.
testx = np.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta)
x = np.compress(testx, x)
y = np.compress(testx, y)
testy = np.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta)
x = np.compress(testy, x)
y = np.compress(testy, y)
lines = []
if len(x) > 1 and len(y) > 1:
# split into separate line segments if necessary.
# (not necessary for cylindrical or pseudocylindricl projections)
xd = (x[1:]-x[0:-1])**2
yd = (y[1:]-y[0:-1])**2
dist = np.sqrt(xd+yd)
if self.projection not in ['cyl','rotpole']:
split = dist > self.rmajor/10.
else:
split = dist > 1.
if np.sum(split) and self.projection not in _cylproj:
ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist()
xl = []
yl = []
iprev = 0
ind.append(len(xd))
for i in ind:
xl.append(x[iprev:i])
yl.append(y[iprev:i])
iprev = i
else:
xl = [x]
yl = [y]
# draw each line segment.
for x,y in zip(xl,yl):
# skip if only a point.
if len(x) > 1 and len(y) > 1:
l = Line2D(x,y,linewidth=linewidth)
l.set_color(color)
l.set_dashes(dashes)
l.set_label('_nolabel_')
if zorder is not None:
l.set_zorder(zorder)
ax.add_line(l)
lines.append(l)
linecolls[circ] = (lines,[])
# draw labels for parallels
# parallels not labelled for fulldisk orthographic or geostationary
if self.projection in ['ortho','geos','nsper','vandg','aeqd'] and max(labels):
if self.projection == 'vandg' or self._fulldisk:
sys.stdout.write('Warning: Cannot label parallels on %s basemap' % _projnames[self.projection])
labels = [0,0,0,0]
# search along edges of map to see if parallels intersect.
# if so, find x,y location of intersection and draw a label there.
dx = (self.xmax-self.xmin)/1000.
dy = (self.ymax-self.ymin)/1000.
if self.projection in _pseudocyl:
lon_0 = self.projparams['lon_0']
for dolab,side in zip(labels,['l','r','t','b']):
if not dolab: continue
# for cylindrical projections, don't draw parallels on top or bottom.
if self.projection in _cylproj + _pseudocyl and side in ['t','b']: continue
if side in ['l','r']:
nmax = int((self.ymax-self.ymin)/dy+1)
yy = np.linspace(self.llcrnry,self.urcrnry,nmax)
if side == 'l':
if self.projection in _pseudocyl:
lats = np.linspace(-89.99,89,99,nmax)
if self.celestial:
lons = (self.projparams['lon_0']+180.)*np.ones(len(lats),lats.dtype)
else:
lons = (self.projparams['lon_0']-180.)*np.ones(len(lats),lats.dtype)
xx, yy = self(lons, lats)
else:
xx = self.llcrnrx*np.ones(yy.shape,yy.dtype)
lons,lats = self(xx,yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
if self.projection in _pseudocyl:
lats = np.linspace(-89.99,89,99,nmax)
if self.celestial:
lons = (self.projparams['lon_0']-180.)*np.ones(len(lats),lats.dtype)
else:
lons = (self.projparams['lon_0']+180.)*np.ones(len(lats),lats.dtype)
xx, yy = self(lons, lats)
else:
xx = self.urcrnrx*np.ones(yy.shape,yy.dtype)
lons,lats = self(xx,yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError('inverse transformation undefined - please adjust the map projection region')
# adjust so 0 <= lons < 360
lons = [(lon+360) % 360 for lon in lons]
else:
nmax = int((self.xmax-self.xmin)/dx+1)
xx = np.linspace(self.llcrnrx,self.urcrnrx,nmax)
if side == 'b':
lons,lats = self(xx,self.llcrnry*np.ones(xx.shape,np.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
lons,lats = self(xx,self.urcrnry*np.ones(xx.shape,np.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError('inverse transformation undefined - please adjust the map projection region')
# adjust so 0 <= lons < 360
lons = [(lon+360) % 360 for lon in lons]
for lat in circles:
# don't label parallels for round polar plots
if self.round: continue
# find index of parallel (there may be two, so
# search from left and right).
nl = _searchlist(lats,lat)
nr = _searchlist(lats[::-1],lat)
if nr != -1: nr = len(lons)-nr-1
latlab = _setlatlab(fmt,lat,labelstyle)
# parallels can intersect each map edge twice.
for i,n in enumerate([nl,nr]):
# don't bother if close to the first label.
if i and abs(nr-nl) < 100: continue
if n >= 0:
t = None
if side == 'l':
if self.projection in _pseudocyl:
if self.celestial:
xlab,ylab = self(lon_0+179.9,lat)
else:
xlab,ylab = self(lon_0-179.9,lat)
else:
xlab = self.llcrnrx
xlab = xlab-xoffset
if self.projection in _pseudocyl:
if lat>0:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='bottom',**text_kwargs)
elif lat<0:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='top',**text_kwargs)
else:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='center',**text_kwargs)
else:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='right',verticalalignment='center',**text_kwargs)
elif side == 'r':
if self.projection in _pseudocyl:
if self.celestial:
xlab,ylab = self(lon_0-179.9,lat)
else:
xlab,ylab = self(lon_0+179.9,lat)
else:
xlab = self.urcrnrx
xlab = xlab+xoffset
if self.projection in _pseudocyl:
if lat>0:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='bottom',**text_kwargs)
elif lat<0:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='top',**text_kwargs)
else:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='center',**text_kwargs)
else:
t=ax.text(xlab,yy[n],latlab,horizontalalignment='left',verticalalignment='center',**text_kwargs)
elif side == 'b':
t = ax.text(xx[n],self.llcrnry-yoffset,latlab,horizontalalignment='center',verticalalignment='top',**text_kwargs)
else:
t = ax.text(xx[n],self.urcrnry+yoffset,latlab,horizontalalignment='center',verticalalignment='bottom',**text_kwargs)
if t is not None: linecolls[lat][1].append(t)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
keys = list(linecolls.keys()); vals = list(linecolls.values())
for k,v in zip(keys,vals):
if v == ([], []):
del linecolls[k]
# add a remove method to each tuple.
else:
linecolls[k] = _tup(linecolls[k])
# override __delitem__ in dict to call remove() on values.
pardict = _dict(linecolls)
# clip parallels for round polar plots (and delete labels).
for lines, _ in pardict.values():
self._cliplimb(ax, lines)
return pardict
[docs] def drawmeridians(self,meridians,color='k',textcolor='k',linewidth=1., zorder=None,\
dashes=[1,1],labels=[0,0,0,0],labelstyle=None,\
fmt='%g',xoffset=None,yoffset=None,ax=None,latmax=None,
**text_kwargs):
r"""
Draw and label meridians (longitude lines) for values (in degrees)
given in the sequence ``meridians``.
.. tabularcolumns:: |l|L|
============== ====================================================
Keyword Description
============== ====================================================
color color to draw meridians (default black).
textcolor color to draw labels (default black).
linewidth line width for meridians (default 1.)
zorder sets the zorder for meridians (if not specified,
uses default zorder for matplotlib.lines.Line2D
objects).
dashes dash pattern for meridians (default [1,1], i.e.
1 pixel on, 1 pixel off).
labels list of 4 values (default [0,0,0,0]) that control
whether meridians are labelled where they intersect
the left, right, top or bottom of the plot. For
example labels=[1,0,0,1] will cause meridians
to be labelled where they intersect the left and
and bottom of the plot, but not the right and top.
labelstyle if set to "+/-", east and west longitudes are
labelled with "+" and "-", otherwise they are
labelled with "E" and "W".
fmt a format string to format the meridian labels
(default '%g') **or** a function that takes a
longitude value in degrees as it's only argument
and returns a formatted string.
xoffset label offset from edge of map in x-direction
(default is 0.01 times width of map in map
projection coordinates).
yoffset label offset from edge of map in y-direction
(default is 0.01 times height of map in map
projection coordinates).
ax axes instance (overrides default axes instance)
latmax absolute value of latitude to which meridians are drawn
(default is 80).
\**text_kwargs additional keyword arguments controlling text
for labels that are passed on to
the text method of the axes instance (see
matplotlib.pyplot.text documentation).
============== ====================================================
returns a dictionary whose keys are the meridian values, and
whose values are tuples containing lists of the
matplotlib.lines.Line2D and matplotlib.text.Text instances
associated with each meridian. Deleting an item from the
dictionary removes the correpsonding meridian from the plot.
"""
text_kwargs['color']=textcolor # pass textcolor kwarg on to ax.text
# for cylindrical projections, try to handle wraparound (i.e. if
# projection is defined in -180 to 0 and user asks for meridians from
# 180 to 360 to be drawn, it should work)
if self.projection in _cylproj or self.projection in _pseudocyl:
def addlon(meridians,madd):
minside = (madd >= self.llcrnrlon and madd <= self.urcrnrlon)
if minside and madd not in meridians: meridians.append(madd)
return meridians
merids = list(meridians)
meridians = []
for m in merids:
meridians = addlon(meridians,m)
meridians = addlon(meridians,m+360)
meridians = addlon(meridians,m-360)
meridians.sort()
# if celestial=True, don't use "E" and "W" labels.
if labelstyle is None and self.celestial:
labelstyle="+/-"
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# don't draw meridians past latmax, always draw parallel at latmax.
if latmax is None: latmax = 80. # unused w/ cyl, merc or miller proj.
# offset for labels.
if yoffset is None:
yoffset = (self.urcrnry-self.llcrnry)/100.
if self.aspect > 1:
yoffset = self.aspect*yoffset
else:
yoffset = yoffset/self.aspect
if xoffset is None:
xoffset = (self.urcrnrx-self.llcrnrx)/100.
lats = np.linspace(self.latmin,self.latmax,10001)
if self.projection not in _cylproj + _pseudocyl:
testlat = np.logical_and(lats>-latmax,lats<latmax)
lats = np.compress(testlat,lats)
xdelta = 0.01*(self.xmax-self.xmin)
ydelta = 0.01*(self.ymax-self.ymin)
linecolls = {}
for merid in meridians:
lons = merid*np.ones(len(lats),np.float32)
x,y = self(lons,lats)
# remove points outside domain.
# leave a little slop around edges (3*xdelta)
# don't really know why, but this appears to be needed to
# or lines sometimes don't reach edge of plot.
testx = np.logical_and(x>=self.xmin-3*xdelta,x<=self.xmax+3*xdelta)
x = np.compress(testx, x)
y = np.compress(testx, y)
testy = np.logical_and(y>=self.ymin-3*ydelta,y<=self.ymax+3*ydelta)
x = np.compress(testy, x)
y = np.compress(testy, y)
lines = []
if len(x) > 1 and len(y) > 1:
# split into separate line segments if necessary.
# (not necessary for mercator or cylindrical or miller).
xd = (x[1:]-x[0:-1])**2
yd = (y[1:]-y[0:-1])**2
dist = np.sqrt(xd+yd)
if self.projection not in ['cyl','rotpole']:
split = dist > self.rmajor/10.
else:
split = dist > 1.
if np.sum(split) and self.projection not in _cylproj:
ind = (np.compress(split,np.squeeze(split*np.indices(xd.shape)))+1).tolist()
xl = []
yl = []
iprev = 0
ind.append(len(xd))
for i in ind:
xl.append(x[iprev:i])
yl.append(y[iprev:i])
iprev = i
else:
xl = [x]
yl = [y]
# draw each line segment.
for x,y in zip(xl,yl):
# skip if only a point.
if len(x) > 1 and len(y) > 1:
l = Line2D(x,y,linewidth=linewidth)
l.set_color(color)
l.set_dashes(dashes)
l.set_label('_nolabel_')
if zorder is not None:
l.set_zorder(zorder)
ax.add_line(l)
lines.append(l)
linecolls[merid] = (lines,[])
# draw labels for meridians.
# meridians not labelled for sinusoidal, hammer, mollweide,
# VanDerGrinten or full-disk orthographic/geostationary.
if self.projection in ['sinu','moll','hammer','vandg'] and max(labels):
sys.stdout.write('Warning: Cannot label meridians on %s basemap' % _projnames[self.projection])
labels = [0,0,0,0]
if self.projection in ['ortho','geos','nsper','aeqd'] and max(labels):
if self._fulldisk and self.boundinglat is None:
sys.stdout.write(" ".join([
"Warning: Cannot label meridians on full-disk Geostationary,"
"Orthographic or Azimuthal equidistant basemap"]))
labels = [0,0,0,0]
# search along edges of map to see if parallels intersect.
# if so, find x,y location of intersection and draw a label there.
dx = (self.xmax-self.xmin)/1000.
dy = (self.ymax-self.ymin)/1000.
if self.projection in _pseudocyl:
lon_0 = self.projparams['lon_0']
xmin,ymin = self(lon_0-179.9,-90)
xmax,ymax = self(lon_0+179.9,90)
for dolab,side in zip(labels,['l','r','t','b']):
if not dolab or self.round: continue
# for cylindrical projections, don't draw meridians on left or right.
if self.projection in _cylproj + _pseudocyl and side in ['l','r']: continue
if side in ['l','r']:
nmax = int((self.ymax-self.ymin)/dy+1)
yy = np.linspace(self.llcrnry,self.urcrnry,nmax)
if side == 'l':
lons,lats = self(self.llcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
lons,lats = self(self.urcrnrx*np.ones(yy.shape,np.float32),yy,inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError('inverse transformation undefined - please adjust the map projection region')
# adjust so 0 <= lons < 360
lons = [(lon+360) % 360 for lon in lons]
else:
nmax = int((self.xmax-self.xmin)/dx+1)
if self.projection in _pseudocyl:
xx = np.linspace(xmin,xmax,nmax)
else:
xx = np.linspace(self.llcrnrx,self.urcrnrx,nmax)
if side == 'b':
lons,lats = self(xx,self.llcrnry*np.ones(xx.shape,np.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
else:
lons,lats = self(xx,self.urcrnry*np.ones(xx.shape,np.float32),inverse=True)
lons = lons.tolist(); lats = lats.tolist()
if max(lons) > 1.e20 or max(lats) > 1.e20:
raise ValueError('inverse transformation undefined - please adjust the map projection region')
# adjust so 0 <= lons < 360
lons = [(lon+360) % 360 for lon in lons]
for lon in meridians:
# adjust so 0 <= lon < 360
lon2 = (lon+360) % 360
# find index of meridian (there may be two, so
# search from left and right).
nl = _searchlist(lons,lon2)
nr = _searchlist(lons[::-1],lon2)
if nr != -1: nr = len(lons)-nr-1
lonlab = _setlonlab(fmt,lon2,labelstyle)
# meridians can intersect each map edge twice.
for i,n in enumerate([nl,nr]):
lat = lats[n]/100.
# no meridians > latmax for projections other than merc,cyl,miller.
if self.projection not in _cylproj and lat > latmax: continue
# don't bother if close to the first label.
if i and abs(nr-nl) < 100: continue
if n >= 0:
t = None
if side == 'l':
t = ax.text(self.llcrnrx-xoffset,yy[n],lonlab,horizontalalignment='right',verticalalignment='center',**text_kwargs)
elif side == 'r':
t = ax.text(self.urcrnrx+xoffset,yy[n],lonlab,horizontalalignment='left',verticalalignment='center',**text_kwargs)
elif side == 'b':
t = ax.text(xx[n],self.llcrnry-yoffset,lonlab,horizontalalignment='center',verticalalignment='top',**text_kwargs)
else:
t = ax.text(xx[n],self.urcrnry+yoffset,lonlab,horizontalalignment='center',verticalalignment='bottom',**text_kwargs)
if t is not None: linecolls[lon][1].append(t)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# remove empty values from linecolls dictionary
keys = list(linecolls.keys()); vals = list(linecolls.values())
for k,v in zip(keys,vals):
if v == ([], []):
del linecolls[k]
else:
# add a remove method to each tuple.
linecolls[k] = _tup(linecolls[k])
# override __delitem__ in dict to call remove() on values.
meridict = _dict(linecolls)
# for round polar plots, clip meridian lines and label them.
if self.round:
# label desired?
label = False
for lab in labels:
if lab: label = True
for merid in meridict:
if not label: continue
# label
lonlab = _setlonlab(fmt,merid,labelstyle)
x,y = self(merid,self.boundinglat)
r = np.sqrt((x-0.5*(self.xmin+self.xmax))**2+
(y-0.5*(self.ymin+self.ymax))**2)
r = r + np.sqrt(xoffset**2+yoffset**2)
if self.projection.startswith('np'):
pole = 1
elif self.projection.startswith('sp'):
pole = -1
elif self.projection == 'ortho' and self.round:
pole = 1
if pole == 1:
theta = (np.pi/180.)*(merid-self.projparams['lon_0']-90)
if self.projection == 'ortho' and\
self.projparams['lat_0'] == -90:
theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90)
x = r*np.cos(theta)+0.5*(self.xmin+self.xmax)
y = r*np.sin(theta)+0.5*(self.ymin+self.ymax)
if x > 0.5*(self.xmin+self.xmax)+xoffset:
horizalign = 'left'
elif x < 0.5*(self.xmin+self.xmax)-xoffset:
horizalign = 'right'
else:
horizalign = 'center'
if y > 0.5*(self.ymin+self.ymax)+yoffset:
vertalign = 'bottom'
elif y < 0.5*(self.ymin+self.ymax)-yoffset:
vertalign = 'top'
else:
vertalign = 'center'
# labels [l,r,t,b]
if labels[0] and not labels[1] and x >= 0.5*(self.xmin+self.xmax)+xoffset: continue
if labels[1] and not labels[0] and x <= 0.5*(self.xmin+self.xmax)-xoffset: continue
if labels[2] and not labels[3] and y <= 0.5*(self.ymin+self.ymax)-yoffset: continue
if labels[3] and not labels[2]and y >= 0.5*(self.ymin+self.ymax)+yoffset: continue
elif pole == -1:
theta = (np.pi/180.)*(-merid+self.projparams['lon_0']+90)
x = r*np.cos(theta)+0.5*(self.xmin+self.xmax)
y = r*np.sin(theta)+0.5*(self.ymin+self.ymax)
if x > 0.5*(self.xmin+self.xmax)-xoffset:
horizalign = 'right'
elif x < 0.5*(self.xmin+self.xmax)+xoffset:
horizalign = 'left'
else:
horizalign = 'center'
if y > 0.5*(self.ymin+self.ymax)-yoffset:
vertalign = 'top'
elif y < 0.5*(self.ymin+self.ymax)+yoffset:
vertalign = 'bottom'
else:
vertalign = 'center'
# labels [l,r,t,b]
if labels[0] and not labels[1] and x <= 0.5*(self.xmin+self.xmax)+xoffset: continue
if labels[1] and not labels[0] and x >= 0.5*(self.xmin+self.xmax)-xoffset: continue
if labels[2] and not labels[3] and y >= 0.5*(self.ymin+self.ymax)-yoffset: continue
if labels[3] and not labels[2] and y <= 0.5*(self.ymin+self.ymax)+yoffset: continue
t=ax.text(x,y,lonlab,horizontalalignment=horizalign,verticalalignment=vertalign,**text_kwargs)
meridict[merid][1].append(t)
for lines, _ in meridict.values():
self._cliplimb(ax, lines)
return meridict
[docs] def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs):
r"""
Draw a polygon centered at ``lon_0,lat_0``. The polygon
approximates a circle on the surface of the earth with radius
``radius_deg`` degrees latitude along longitude ``lon_0``,
made up of ``npts`` vertices.
The polygon represents a Tissot's indicatrix
(http://en.wikipedia.org/wiki/Tissot's_Indicatrix),
which when drawn on a map shows the distortion
inherent in the map projection.
.. note::
Cannot handle situations in which the polygon intersects
the edge of the map projection domain, and then re-enters the domain.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \**kwargs passed on to matplotlib.patches.Polygon.
returns a matplotlib.patches.Polygon object."""
ax = kwargs.pop('ax', None) or self._check_ax()
g = pyproj.Geod(a=self.rmajor,b=self.rminor)
az12,az21,dist = g.inv(lon_0,lat_0,lon_0,lat_0+radius_deg)
seg = [self(lon_0,lat_0+radius_deg)]
delaz = 360./npts
az = az12
for n in range(npts):
az = az+delaz
lon, lat, az21 = g.fwd(lon_0, lat_0, az, dist)
x,y = self(lon,lat)
# add segment if it is in the map projection region.
if x < 1.e20 and y < 1.e20:
seg.append((x,y))
poly = Polygon(seg,**kwargs)
ax.add_patch(poly)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip polygons to map limbs
poly,c = self._cliplimb(ax,poly)
return poly
[docs] 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
[docs] def drawgreatcircle(self,lon1,lat1,lon2,lat2,del_s=100.,**kwargs):
r"""
Draw a great circle on the map from the longitude-latitude
pair ``lon1,lat1`` to ``lon2,lat2``
.. tabularcolumns:: |l|L|
============== =======================================================
Keyword Description
============== =======================================================
del_s points on great circle computed every del_s kilometers
(default 100).
\**kwargs other keyword arguments are passed on to :meth:`plot`
method of Basemap instance.
============== =======================================================
Returns a list with a single ``matplotlib.lines.Line2D`` object like a
call to ``pyplot.plot()``.
"""
# use great circle formula for a perfect sphere.
gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
az12,az21,dist = gc.inv(lon1,lat1,lon2,lat2)
npoints = int((dist+0.5*1000.*del_s)/(1000.*del_s))
lonlats = gc.npts(lon1,lat1,lon2,lat2,npoints)
lons = [lon1]; lats = [lat1]
for lon, lat in lonlats:
lons.append(lon)
lats.append(lat)
lons.append(lon2); lats.append(lat2)
x, y = self(lons, lats)
# Correct wrap around effect of great circles
# get points
_p = self.plot(x,y,**kwargs)
p = _p[0].get_path()
# since we know the difference between any two points, we can use this to find wrap arounds on the plot
max_dist = 1000*del_s*2
# calculate distances and compare with max allowable distance
dists = np.abs(np.diff(p.vertices[:,0]))
cuts = np.where( dists > max_dist )[0]
# if there are any cut points, cut them and begin again at the next point
for i,k in enumerate(cuts):
# vertex to cut at
cut_point = cuts[i]
# create new vertices with a nan inbetween and set those as the path's vertices
verts = np.concatenate(
[p.vertices[:cut_point, :],
[[np.nan, np.nan]],
p.vertices[cut_point+1:, :]]
)
p.codes = None
p.vertices = verts
return _p
[docs] 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
[docs] def set_axes_limits(self,ax=None):
"""
Final step in Basemap method wrappers of Axes plotting methods:
Set axis limits, fix aspect ratio for map domain using current
or specified axes instance. This is done only once per axes
instance.
In interactive mode, this method always calls draw_if_interactive
before returning.
"""
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# If we have already set the axes limits, and if the user
# has not defeated this by turning autoscaling back on,
# then all we need to do is plot if interactive.
if (hash(ax) in self._initialized_axes
and not ax.get_autoscalex_on()
and not ax.get_autoscaley_on()):
if mpl.is_interactive():
import matplotlib.pyplot as plt
plt.draw_if_interactive()
return
self._initialized_axes.add(hash(ax))
# Take control of axis scaling:
ax.set_autoscale_on(False)
# update data limits for map domain.
corners = ((self.llcrnrx, self.llcrnry), (self.urcrnrx, self.urcrnry))
ax.update_datalim(corners)
ax.set_xlim((self.llcrnrx, self.urcrnrx))
ax.set_ylim((self.llcrnry, self.urcrnry))
# if map boundary not yet drawn for elliptical maps, draw it with default values.
if not self._mapboundarydrawn or self._mapboundarydrawn not in ax.patches:
# elliptical map, draw boundary manually.
if ((self.projection in ['ortho', 'geos', 'nsper', 'aeqd'] and
self._fulldisk) or self.round or
self.projection in _pseudocyl):
# first draw boundary, no fill
limb1 = self.drawmapboundary(fill_color='none', ax=ax)
# draw another filled patch, with no boundary.
limb2 = self.drawmapboundary(fill_color='none', linewidth=0, ax=ax)
self._mapboundarydrawn = limb2
# for elliptical map, always turn off axis_frame.
if ((self.projection in ['ortho', 'geos', 'nsper', 'aeqd'] and
self._fulldisk) or self.round or
self.projection in _pseudocyl):
# turn off axes frame.
ax.set_frame_on(False)
# make sure aspect ratio of map preserved.
# plot is re-centered in bounding rectangle.
# (anchor instance var determines where plot is placed)
if self.fix_aspect:
ax.set_aspect('equal',anchor=self.anchor)
else:
ax.set_aspect('auto',anchor=self.anchor)
# make sure axis ticks are turned off.
if self.noticks:
ax.set_xticks([])
ax.set_yticks([])
# force draw if in interactive mode.
if mpl.is_interactive():
import matplotlib.pyplot as plt
plt.draw_if_interactive()
def _save_use_hold(self, ax, kwargs):
h = kwargs.pop('hold', None)
if hasattr(ax, '_hold'):
self._tmp_hold = ax._hold
if h is not None:
ax._hold = h
def _restore_hold(self, ax):
if hasattr(ax, '_hold'):
ax._hold = self._tmp_hold
[docs] @_transform1d
def scatter(self, *args, **kwargs):
r"""
Plot points with markers on the map
(see matplotlib.pyplot.scatter documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axes instance.
Other \**kwargs passed on to matplotlib.pyplot.scatter.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
ret = ax.scatter(*args, **kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt:
plt.sci(ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret,c = self._cliplimb(ax,ret)
return ret
[docs] @_transform1d
def plot(self, *args, **kwargs):
r"""
Draw lines and/or markers on the map
(see matplotlib.pyplot.plot documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \**kwargs passed on to matplotlib.pyplot.plot.
"""
ax = kwargs.pop('ax', None) or self._check_ax()
self._save_use_hold(ax, kwargs)
try:
ret = ax.plot(*args, **kwargs)
finally:
self._restore_hold(ax)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret,c = self._cliplimb(ax,ret)
return ret
[docs] def imshow(self, *args, **kwargs):
r"""
Display an image over the map
(see matplotlib.pyplot.imshow documentation).
``extent`` and ``origin`` keywords set automatically so image
will be drawn over map region.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \**kwargs passed on to matplotlib.pyplot.plot.
returns an matplotlib.image.AxesImage instance.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
kwargs['extent']=(self.llcrnrx,self.urcrnrx,self.llcrnry,self.urcrnry)
# use origin='lower', unless overridden.
if 'origin' not in kwargs:
kwargs['origin']='lower'
self._save_use_hold(ax, kwargs)
try:
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=DeprecationWarning)
ret = ax.imshow(*args, **kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt:
plt.sci(ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip image to map limbs
ret,c = self._cliplimb(ax,ret)
return ret
[docs] @_transform
def pcolor(self,x,y,data,**kwargs):
r"""
Make a pseudo-color plot over the map
(see matplotlib.pyplot.pcolor documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
If x or y are outside projection limb (i.e. they have values > 1.e20)
they will be convert to masked arrays with those values masked.
As a result, those values will not be plotted.
If ``tri`` is set to ``True``, an unstructured grid is assumed
(x,y,data must be 1-d) and matplotlib.pyplot.tripcolor is used.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \**kwargs passed on to matplotlib.pyplot.pcolor (or tripcolor if
``tri=True``).
Note: (taken from matplotlib.pyplot.pcolor documentation)
Ideally the dimensions of x and y should be one greater than those of data;
if the dimensions are the same, then the last row and column of data will be ignored.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
if kwargs.pop('tri', False):
import matplotlib.tri as tri
# for unstructured grids, toss out points outside
# projection limb (don't use those points in triangulation).
if ma.isMA(data):
data = data.filled(fill_value=1.e30)
masked=True
else:
masked=False
mask = np.logical_or(x<1.e20,y<1.e20)
x = np.compress(mask,x)
y = np.compress(mask,y)
data = np.compress(mask,data)
if masked:
triang = tri.Triangulation(x, y)
z = data[triang.triangles]
mask = (z > 1.e20).sum(axis=-1)
triang.set_mask(mask)
ret = ax.tripcolor(triang,data,**kwargs)
else:
ret = ax.tripcolor(x,y,data,**kwargs)
else:
# make x,y masked arrays
# (masked where data is outside of projection limb)
x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
ret = ax.pcolor(x,y,data,**kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt:
plt.sci(ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret,c = self._cliplimb(ax,ret)
if self.round:
# for some reason, frame gets turned on.
ax.set_frame_on(False)
return ret
[docs] @_transform
def pcolormesh(self,x,y,data,**kwargs):
r"""
Make a pseudo-color plot over the map
(see matplotlib.pyplot.pcolormesh documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \**kwargs passed on to matplotlib.pyplot.pcolormesh.
Note: (taken from matplotlib.pyplot.pcolor documentation)
Ideally the dimensions of x and y should be one greater than those of data;
if the dimensions are the same, then the last row and column of data will be ignored.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
# fix for invalid grid points
if ((np.any(x > 1e20) or np.any(y > 1e20)) and
x.ndim == 2 and y.ndim == 2):
if x.shape != y.shape:
raise ValueError('pcolormesh: x and y need same dimension')
nx,ny = x.shape
if nx < data.shape[0] or ny < data.shape[1]:
raise ValueError('pcolormesh: data dimension needs to be at least that of x and y.')
mask = (
(x[:-1,:-1] > 1e20) |
(x[1:,:-1] > 1e20) |
(x[:-1,1:] > 1e20) |
(x[1:,1:] > 1e20) |
(y[:-1,:-1] > 1e20) |
(y[1:,:-1] > 1e20) |
(y[:-1,1:] > 1e20) |
(y[1:,1:] > 1e20)
)
# we do not want to overwrite original array
data = data[:nx-1,:ny-1].copy()
data[mask] = np.nan
self._save_use_hold(ax, kwargs)
try:
ret = ax.pcolormesh(x,y,data,**kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt:
plt.sci(ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret,c = self._cliplimb(ax,ret)
if self.round:
# for some reason, frame gets turned on.
ax.set_frame_on(False)
return ret
[docs] def hexbin(self,x,y,**kwargs):
r"""
Make a hexagonal binning plot of x versus y, where x, y are 1-D
sequences of the same length, N. If C is None (the default), this is a
histogram of the number of occurences of the observations at
(x[i],y[i]).
If C is specified, it specifies values at the coordinate (x[i],y[i]).
These values are accumulated for each hexagonal bin and then reduced
according to reduce_C_function, which defaults to the numpy mean function
(np.mean). (If C is specified, it must also be a 1-D sequence of the
same length as x and y.)
x, y and/or C may be masked arrays, in which case only unmasked points
will be plotted.
(see matplotlib.pyplot.hexbin documentation).
Extra keyword ``ax`` can be used to override the default axis instance.
Other \**kwargs passed on to matplotlib.pyplot.hexbin
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
# make x,y masked arrays
# (masked where data is outside of projection limb)
x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
ret = ax.hexbin(x,y,**kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt:
plt.sci(ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret,c = self._cliplimb(ax,ret)
return ret
[docs] @_transform
def contour(self,x,y,data,*args,**kwargs):
r"""
Make a contour plot over the map
(see matplotlib.pyplot.contour documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
If ``tri`` is set to ``True``, an unstructured grid is assumed
(x,y,data must be 1-d) and matplotlib.pyplot.tricontour is used.
Other \*args and \**kwargs passed on to matplotlib.pyplot.contour
(or tricontour if ``tri=True``).
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
if kwargs.pop('tri', False):
import matplotlib.tri as tri
# for unstructured grids, toss out points outside
# projection limb (don't use those points in triangulation).
if ma.isMA(data):
data = data.filled(fill_value=1.e30)
masked=True
else:
masked=False
mask = np.logical_or(x<self.xmin,y<self.xmin) +\
np.logical_or(x>self.xmax,y>self.xmax)
x = np.compress(mask,x)
y = np.compress(mask,y)
data = np.compress(mask,data)
if masked:
triang = tri.Triangulation(x, y)
z = data[triang.triangles]
mask = (z > 1.e20).sum(axis=-1)
triang.set_mask(mask)
CS = ax.tricontour(triang,data,*args,**kwargs)
else:
CS = ax.tricontour(x,y,data,*args,**kwargs)
else:
# make sure x is monotonically increasing - if not,
# print warning suggesting that the data be shifted in longitude
# with the shiftgrid function.
# only do this check for global projections.
if self.projection in _cylproj + _pseudocyl:
xx = x[x.shape[0]//2,:]
condition = (xx >= self.xmin) & (xx <= self.xmax)
xl = xx.compress(condition).tolist()
xs = xl[:]
xs.sort()
if xl != xs:
sys.stdout.write(" ".join([
"WARNING: x coordinate not montonically increasing",
"- contour plot may not be what you expect. If it",
"looks odd, you can either adjust the map projection",
"region to be consistent with your data, or (if your",
"data is on a global lat/lon grid) use the shiftdata",
"method to adjust the data to be consistent with the",
"map projection region (see examples/shiftdata.py)"]))
# mask for points more than one grid length outside projection limb.
xx = ma.masked_where(x > 1.e20, x)
yy = ma.masked_where(y > 1.e20, y)
epsx = np.abs(xx[:,1:]-xx[:,0:-1]).max()
epsy = np.abs(yy[1:,:]-yy[0:-1,:]).max()
xymask = \
np.logical_or(np.greater(x,self.xmax+epsx),np.greater(y,self.ymax+epsy))
xymask = xymask + \
np.logical_or(np.less(x,self.xmin-epsx),np.less(y,self.ymin-epsy))
data = ma.asarray(data)
# combine with data mask.
mask = np.logical_or(ma.getmaskarray(data),xymask)
data = ma.masked_array(data,mask=mask)
CS = ax.contour(x,y,data,*args,**kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt and CS.get_array() is not None:
plt.sci(CS)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
if isinstance(CS, Artist):
# Since MPL 3.8, `QuadContourSet` objects are `Artist` objects too.
CS, c = self._cliplimb(ax, CS)
else:
CS.collections, c = self._cliplimb(ax, CS.collections)
return CS
[docs] @_transform
def contourf(self,x,y,data,*args,**kwargs):
r"""
Make a filled contour plot over the map
(see matplotlib.pyplot.contourf documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
If x or y are outside projection limb (i.e. they have values > 1.e20),
the corresponing data elements will be masked.
Extra keyword 'ax' can be used to override the default axis instance.
If ``tri`` is set to ``True``, an unstructured grid is assumed
(x,y,data must be 1-d) and matplotlib.pyplot.tricontourf is used.
Other \*args and \**kwargs passed on to matplotlib.pyplot.contourf
(or tricontourf if ``tri=True``).
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
if kwargs.get('tri', False):
import matplotlib.tri as tri
# for unstructured grids, toss out points outside
# projection limb (don't use those points in triangulation).
if ma.isMA(data):
data = data.filled(fill_value=1.e30)
masked=True
else:
masked=False
mask = np.logical_or(x<1.e20,y<1.e20)
x = np.compress(mask,x)
y = np.compress(mask,y)
data = np.compress(mask,data)
if masked:
triang = tri.Triangulation(x, y)
z = data[triang.triangles]
mask = (z > 1.e20).sum(axis=-1)
triang.set_mask(mask)
CS = ax.tricontourf(triang,data,*args,**kwargs)
else:
CS = ax.tricontourf(x,y,data,*args,**kwargs)
else:
# make sure x is monotonically increasing - if not,
# print warning suggesting that the data be shifted in longitude
# with the shiftgrid function.
# only do this check for global projections.
if self.projection in _cylproj + _pseudocyl:
xx = x[x.shape[0]//2,:]
condition = (xx >= self.xmin) & (xx <= self.xmax)
xl = xx.compress(condition).tolist()
xs = xl[:]
xs.sort()
if xl != xs:
sys.stdout.write(" ".join([
"WARNING: x coordinate not montonically increasing",
"- contour plot may not be what you expect. If it",
"looks odd, you can either adjust the map projection",
"region to be consistent with your data, or (if your",
"data is on a global lat/lon grid) use the shiftgrid",
"function to adjust the data to be consistent with the",
"map projection region (see examples/contour_demo.py)"]))
# mask for points more than one grid length outside projection limb.
xx = ma.masked_where(x > 1.e20, x)
yy = ma.masked_where(y > 1.e20, y)
if self.projection != 'omerc':
epsx = np.abs(xx[:,1:]-xx[:,0:-1]).max()
epsy = np.abs(yy[1:,:]-yy[0:-1,:]).max()
else: # doesn't work for omerc (FIXME)
epsx = 0.; epsy = 0
xymask = \
np.logical_or(np.greater(x,self.xmax+epsx),np.greater(y,self.ymax+epsy))
xymask = xymask + \
np.logical_or(np.less(x,self.xmin-epsx),np.less(y,self.ymin-epsy))
data = ma.asarray(data)
# combine with data mask.
mask = np.logical_or(ma.getmaskarray(data),xymask)
data = ma.masked_array(data,mask=mask)
CS = ax.contourf(x,y,data,*args,**kwargs)
finally:
self._restore_hold(ax)
# reset current active image (only if pyplot is imported).
if plt and CS.get_array() is not None:
plt.sci(CS)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
if isinstance(CS, Artist):
# Since MPL 3.8, `QuadContourSet` objects are `Artist` objects too.
CS, c = self._cliplimb(ax, CS)
else:
CS.collections, c = self._cliplimb(ax, CS.collections)
return CS
[docs] @_transformuv
def quiver(self, x, y, u, v, *args, **kwargs):
r"""
Make a vector plot (u, v) with arrows on the map.
Arguments may be 1-D or 2-D arrays or sequences
(see matplotlib.pyplot.quiver documentation for details).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \*args and \**kwargs passed on to matplotlib.pyplot.quiver.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
ret = ax.quiver(x,y,u,v,*args,**kwargs)
finally:
self._restore_hold(ax)
if plt is not None and ret.get_array() is not None:
plt.sci(ret)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret,c = self._cliplimb(ax,ret)
return ret
[docs] @_transformuv
def streamplot(self, x, y, u, v, *args, **kwargs):
r"""
Draws streamlines of a vector flow.
(see matplotlib.pyplot.streamplot documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \*args and \**kwargs passed on to matplotlib.pyplot.streamplot.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
self._save_use_hold(ax, kwargs)
try:
ret = ax.streamplot(x,y,u,v,*args,**kwargs)
finally:
self._restore_hold(ax)
if plt is not None and ret.lines.get_array() is not None:
plt.sci(ret.lines)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
ret.lines,c = self._cliplimb(ax,ret.lines)
ret.arrows,c = self._cliplimb(ax,ret.arrows)
# streamplot arrows not returned in matplotlib 1.1.1, so clip all
# FancyArrow patches attached to axes instance.
if c is not None:
for p in ax.patches:
if isinstance(p,FancyArrowPatch): p.set_clip_path(c)
return ret
[docs] @_transformuv
def barbs(self, x, y, u, v, *args, **kwargs):
r"""
Make a wind barb plot (u, v) with on the map.
(see matplotlib.pyplot.barbs documentation).
If ``latlon`` keyword is set to True, x,y are intrepreted as
longitude and latitude in degrees. Data and longitudes are
automatically shifted to match map projection region for cylindrical
and pseudocylindrical projections, and x,y are transformed to map
projection coordinates. If ``latlon`` is False (default), x and y
are assumed to be map projection coordinates.
Extra keyword ``ax`` can be used to override the default axis instance.
Other \*args and \**kwargs passed on to matplotlib.pyplot.barbs
Returns two matplotlib.axes.Barbs instances, one for the Northern
Hemisphere and one for the Southern Hemisphere.
"""
ax, plt = self._ax_plt_from_kw(kwargs)
lons, lats = self(x, y, inverse=True)
unh = ma.masked_where(lats <= 0, u)
vnh = ma.masked_where(lats <= 0, v)
ush = ma.masked_where(lats > 0, u)
vsh = ma.masked_where(lats > 0, v)
self._save_use_hold(ax, kwargs)
try:
retnh = ax.barbs(x,y,unh,vnh,*args,**kwargs)
kwargs['flip_barb']=True
retsh = ax.barbs(x,y,ush,vsh,*args,**kwargs)
finally:
self._restore_hold(ax)
# Because there are two collections returned in general,
# we can't set the current image...
#if plt is not None and ret.get_array() is not None:
# plt.sci(retnh)
# set axes limits to fit map region.
self.set_axes_limits(ax=ax)
# clip to map limbs
retnh,c = self._cliplimb(ax,retnh)
retsh,c = self._cliplimb(ax,retsh)
return retnh,retsh
[docs] def drawlsmask(self,land_color="0.8",ocean_color="w",lsmask=None,
lsmask_lons=None,lsmask_lats=None,lakes=True,resolution='l',grid=5,**kwargs):
r"""
Draw land-sea mask image.
.. note::
The land-sea mask image cannot be overlaid on top
of other images, due to limitations in matplotlib image handling
(you can't specify the zorder of an image).
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
land_color desired land color (color name or rgba tuple).
Default gray ("0.8").
ocean_color desired water color (color name or rgba tuple).
Default white.
lsmask An array of 0's for ocean pixels, 1's for
land pixels and 2's for lake/pond pixels.
Default is None
(default 5-minute resolution land-sea mask is used).
lakes Plot lakes and ponds (Default True)
lsmask_lons 1d array of longitudes for lsmask (ignored
if lsmask is None). Longitudes must be ordered
from -180 W eastward.
lsmask_lats 1d array of latitudes for lsmask (ignored
if lsmask is None). Latitudes must be ordered
from -90 S northward.
resolution gshhs coastline resolution used to define land/sea
mask (default 'l', available 'c','l','i','h' or 'f')
grid land/sea mask grid spacing in minutes (Default 5;
10, 2.5 and 1.25 are also available).
\**kwargs extra keyword arguments passed on to
:meth:`imshow`
============== ====================================================
If any of the lsmask, lsmask_lons or lsmask_lats keywords are not
set, the built in GSHHS land-sea mask datasets are used.
Extra keyword ``ax`` can be used to override the default axis instance.
returns a matplotlib.image.AxesImage instance.
"""
# convert land and water colors to integer rgba tuples with
# values between 0 and 255.
from matplotlib.colors import ColorConverter
c = ColorConverter()
# if conversion fails, assume it's because the color
# given is already an rgba tuple with values between 0 and 255.
try:
cl = c.to_rgba(land_color)
rgba_land = tuple([int(255*x) for x in cl])
except:
rgba_land = land_color
try:
co = c.to_rgba(ocean_color)
rgba_ocean = tuple([int(255*x) for x in co])
except:
rgba_ocean = ocean_color
# look for axes instance (as keyword, an instance variable
# or from plt.gca().
ax = kwargs.pop('ax', None) or self._check_ax()
# Clear saved lsmask if new lsmask is passed
if lsmask is not None or lsmask_lons is not None \
or lsmask_lats is not None:
# Make sure passed lsmask is not the same as cached mask
if lsmask is not self.lsmask:
self.lsmask = None
# if lsmask,lsmask_lons,lsmask_lats keywords not given,
# read default land-sea mask in from file.
if lsmask is None or lsmask_lons is None or lsmask_lats is None:
# if lsmask instance variable already set, data already
# read in.
if self.lsmask is None:
# read in land/sea mask.
lsmask_lons, lsmask_lats, lsmask =\
_readlsmask(lakes=lakes,resolution=resolution,grid=grid)
# instance variable lsmask is set on first invocation,
# it contains the land-sea mask interpolated to the native
# projection grid. Further calls to drawlsmask will not
# redo the interpolation (unless a new land-sea mask is passed
# in via the lsmask, lsmask_lons, lsmask_lats keywords).
# is it a cylindrical projection whose limits lie
# outside the limits of the image?
cylproj = self.projection in _cylproj and \
(self.urcrnrlon > lsmask_lons[-1] or \
self.llcrnrlon < lsmask_lons[0])
if cylproj:
# stack grids side-by-side (in longitiudinal direction), so
# any range of longitudes may be plotted on a world map.
# in versions of NumPy later than 1.10.0, concatenate will
# not stack these arrays as expected. If axis 1 is outside
# the dimensions of the array, concatenate will now raise
# an IndexError. Using hstack instead.
lsmask_lons = \
np.hstack((lsmask_lons,lsmask_lons[1:] + 360))
lsmask = \
np.hstack((lsmask,lsmask[:,1:]))
else:
if lakes: lsmask = np.where(lsmask==2,np.array(0,np.uint8),lsmask)
# transform mask to nx x ny regularly spaced native projection grid
# nx and ny chosen to have roughly the same horizontal
# resolution as mask.
if self.lsmask is None:
nlons = len(lsmask_lons)
nlats = len(lsmask_lats)
if self.projection == 'cyl':
dx = lsmask_lons[1]-lsmask_lons[0]
else:
dx = (np.pi/180.)*(lsmask_lons[1]-lsmask_lons[0])*self.rmajor
nx = int((self.xmax-self.xmin)/dx)+1; ny = int((self.ymax-self.ymin)/dx)+1
# interpolate rgba values from proj='cyl' (geographic coords)
# to a rectangular map projection grid.
mask,x,y = self.transform_scalar(lsmask,lsmask_lons,\
lsmask_lats,nx,ny,returnxy=True,order=0,masked=255)
lsmask_lats.dtype
# for these projections, points outside the projection
# limb have to be set to transparent manually.
if self.projection in _pseudocyl:
lons, lats = self(x, y, inverse=True)
lon_0 = self.projparams['lon_0']
lats = lats[:,nx//2]
lons1 = (lon_0+180.)*np.ones(lons.shape[0],np.float64)
lons2 = (lon_0-180.)*np.ones(lons.shape[0],np.float64)
xmax,ytmp = self(lons1,lats)
xmin,ytmp = self(lons2,lats)
for j in range(lats.shape[0]):
xx = x[j,:]
mask[j,:]=np.where(np.logical_or(xx<xmin[j],xx>xmax[j]),\
255,mask[j,:])
self.lsmask = mask
ny, nx = self.lsmask.shape
rgba = np.ones((ny,nx,4),np.uint8)
rgba_land = np.array(rgba_land,np.uint8)
rgba_ocean = np.array(rgba_ocean,np.uint8)
for k in range(4):
rgba[:,:,k] = np.where(self.lsmask,rgba_land[k],rgba_ocean[k])
# make points outside projection limb transparent.
rgba[:,:,3] = np.where(self.lsmask==255,0,rgba[:,:,3])
# plot mask as rgba image.
im = self.imshow(rgba,interpolation='nearest',ax=ax,**kwargs)
# clip to map limbs.
im,c = self._cliplimb(ax,im)
return im
[docs] def bluemarble(self, ax=None, scale=None, **kwargs):
r"""Display blue marble image as map background.
The original image is available at::
https://visibleearth.nasa.gov/
The default image size is 5400x2700, which can be a bit slow
memory-consuming. The ``scale`` keyword allows to downsample
the image (e.g. ``scale=0.5`` downsamples to 2700x1350).
Parameters
----------
ax : matplotlib.image.AxesImage, optional
if given, alternative axes instance where the image is drawn
scale : float, optional
if given, rescale the image by the given factor
\**kwargs : dict, optional
keyword arguments passed on to :meth:`Basemap.imshow`.
Returns
-------
ax : matplotlib.image.AxesImage
axes instance
"""
image = "bluemarble"
if ax is not None:
return self.warpimage(image=image, ax=ax, scale=scale, **kwargs)
return self.warpimage(image=image, scale=scale, **kwargs)
[docs] def shadedrelief(self, ax=None, scale=None, **kwargs):
r"""Display shaded relief image as map background.
The original image is available at::
https://www.shadedrelief.com/
The default image size is 5400x2700, which can be a bit slow
memory-consuming. The ``scale`` keyword allows to downsample
the image (e.g. ``scale=0.5`` downsamples to 2700x1350).
Parameters
----------
ax : matplotlib.image.AxesImage, optional
if given, alternative axes instance where the image is drawn
scale : float, optional
if given, rescale the image by the given factor
\**kwargs : dict, optional
keyword arguments passed on to :meth:`Basemap.imshow`.
Returns
-------
ax : matplotlib.image.AxesImage
axes instance
"""
image = "shadedrelief"
if ax is not None:
return self.warpimage(image=image, ax=ax, scale=scale, **kwargs)
return self.warpimage(image=image, scale=scale, **kwargs)
[docs] def etopo(self, ax=None, scale=None, **kwargs):
r"""Display ETOPO relief image as map background.
The original image is available at::
http://www.ngdc.noaa.gov/mgg/global/global.html
The default image size is 5400x2700, which can be a bit slow
memory-consuming. The ``scale`` keyword allows to downsample
the image (e.g. ``scale=0.5`` downsamples to 2700x1350).
Parameters
----------
ax : matplotlib.image.AxesImage, optional
if given, alternative axes instance where the image is drawn
scale : float, optional
if given, rescale the image by the given factor
\**kwargs : dict, optional
keyword arguments passed on to :meth:`Basemap.imshow`.
Returns
-------
ax : matplotlib.image.AxesImage
axes instance
"""
image = "etopo"
if ax is not None:
return self.warpimage(image=image, ax=ax, scale=scale, **kwargs)
return self.warpimage(image=image, scale=scale, **kwargs)
[docs] def warpimage(self,image="bluemarble",scale=None,**kwargs):
r"""
Display an image (filename given by ``image`` keyword) as a map background.
If image is a URL (starts with 'http'), it is downloaded to a temp
file using urllib.urlretrieve.
Default (if ``image`` not specified) is to display
'blue marble next generation' image from http://visibleearth.nasa.gov/.
Specified image must have pixels covering the whole globe in a regular
lat/lon grid, starting and -180W and the South Pole.
Works with the global images from
http://earthobservatory.nasa.gov/Features/BlueMarble/BlueMarble_monthlies.php.
The ``scale`` keyword can be used to downsample (rescale) the image.
Values less than 1.0 will speed things up at the expense of image
resolution.
Extra keyword ``ax`` can be used to override the default axis instance.
\**kwargs passed on to :meth:`imshow`.
returns a matplotlib.image.AxesImage instance.
"""
# fix PIL import on some versions of OSX and scipy
try:
from PIL import Image
except ImportError:
try:
import Image
except ImportError:
raise ImportError("warpimage method requires PIL "
"(http://pillow.readthedocs.io)")
import matplotlib.image as mpimg
def pil_to_array(*args, **kwargs):
"""Call :func:`~mpimg.pil_to_array` ignoring deprecation warnings."""
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=DeprecationWarning)
return mpimg.pil_to_array(*args, **kwargs)
if self.celestial:
raise ValueError("warpimage does not work in celestial coordinates")
ax = kwargs.pop('ax', None) or self._check_ax()
# default image file is blue marble next generation
# from NASA (http://visibleearth.nasa.gov).
if image == "bluemarble":
file = os.path.join(basemap_datadir,'bmng.jpg')
# display shaded relief image (from
# http://www.shadedreliefdata.com)
elif image == "shadedrelief":
file = os.path.join(basemap_datadir,'shadedrelief.jpg')
# display etopo image (from
# http://www.ngdc.noaa.gov/mgg/image/globalimages.html)
elif image == "etopo":
file = os.path.join(basemap_datadir,'etopo1.jpg')
else:
file = image
# if image is same as previous invocation, used cached data.
# if not, regenerate rgba data.
if not hasattr(self,'_bm_file') or self._bm_file != file:
newfile = True
else:
newfile = False
if file.startswith('http'):
self._bm_file, headers = urlretrieve(file)
else:
self._bm_file = file
# bmproj is True if map projection region is same as
# image region.
bmproj = self.projection == 'cyl' and \
self.llcrnrlon == -180 and self.urcrnrlon == 180 and \
self.llcrnrlat == -90 and self.urcrnrlat == 90
# read in jpeg image to rgba array of normalized floats.
if not hasattr(self,'_bm_rgba') or newfile:
pilImage = Image.open(self._bm_file)
if scale is not None:
w, h = pilImage.size
width = int(np.round(w*scale))
height = int(np.round(h*scale))
pilImage = pilImage.resize((width,height),Image.LANCZOS)
# orientation of arrays returned by pil_to_array changed in
# matplotlib 1.2 (https://github.com/matplotlib/matplotlib/pull/616)
self._bm_rgba = pil_to_array(pilImage)[::-1,:]
# define lat/lon grid that image spans.
nlons = self._bm_rgba.shape[1]; nlats = self._bm_rgba.shape[0]
delta = 360./float(nlons)
self._bm_lons = np.arange(-180.+0.5*delta,180.,delta)
self._bm_lats = np.arange(-90.+0.5*delta,90.,delta)
# is it a cylindrical projection whose limits lie
# outside the limits of the image?
cylproj = self.projection in _cylproj and \
(self.urcrnrlon > self._bm_lons[-1] or \
self.llcrnrlon < self._bm_lons[0])
# if pil_to_array returns a 2D array, it's a grayscale image.
# create an RGB image, with R==G==B.
if self._bm_rgba.ndim == 2:
tmp = np.empty(self._bm_rgba.shape+(3,),np.uint8)
for k in range(3):
tmp[:,:,k] = self._bm_rgba
self._bm_rgba = tmp
if cylproj and not bmproj:
# stack grids side-by-side (in longitiudinal direction), so
# any range of longitudes may be plotted on a world map.
self._bm_lons = \
np.concatenate((self._bm_lons,self._bm_lons+360),0)
self._bm_rgba = \
np.concatenate((self._bm_rgba,self._bm_rgba),1)
# convert to normalized floats.
self._bm_rgba = self._bm_rgba.astype(np.float32)/255.
if not bmproj: # interpolation necessary.
if newfile or not hasattr(self,'_bm_rgba_warped'):
# transform to nx x ny regularly spaced native
# projection grid.
# nx and ny chosen to have roughly the
# same horizontal res as original image.
if self.projection != 'cyl':
dx = 2.*np.pi*self.rmajor/float(nlons)
nx = int((self.xmax-self.xmin)/dx)+1
ny = int((self.ymax-self.ymin)/dx)+1
else:
dx = 360./float(nlons)
nx = int((self.urcrnrlon-self.llcrnrlon)/dx)+1
ny = int((self.urcrnrlat-self.llcrnrlat)/dx)+1
self._bm_rgba_warped = np.ones((ny,nx,4),np.float64)
# interpolate rgba values from geographic coords (proj='cyl')
# to map projection coords.
# if masked=True, values outside of
# projection limb will be masked.
for k in range(self._bm_rgba.shape[2]):
self._bm_rgba_warped[:,:,k],x,y = \
self.transform_scalar(self._bm_rgba[:,:,k],\
self._bm_lons,self._bm_lats,nx,ny,returnxy=True)
# for ortho,geos mask pixels outside projection limb.
if self.projection in ['geos','ortho','nsper'] or \
(self.projection == 'aeqd' and self._fulldisk):
lonsr,latsr = self(x,y,inverse=True)
mask = ma.zeros((ny,nx,4),np.int8)
mask[:,:,0] = np.logical_or(lonsr>1.e20,latsr>1.e30)
for k in range(1,4):
mask[:,:,k] = mask[:,:,0]
self._bm_rgba_warped = \
ma.masked_array(self._bm_rgba_warped,mask=mask)
# make points outside projection limb transparent.
self._bm_rgba_warped = self._bm_rgba_warped.filled(0.)
# treat pseudo-cyl projections such as mollweide, robinson and sinusoidal.
elif self.projection in _pseudocyl and \
self.projection != 'hammer':
lonsr,latsr = self(x,y,inverse=True)
mask = ma.zeros((ny,nx,4),np.int8)
lon_0 = self.projparams['lon_0']
lonright = lon_0 + 180. - 1E-10
lonleft = lon_0 - 180. + 1E-10
x1 = np.array(ny*[0.5*(self.xmax + self.xmin)],np.float64)
y1 = np.linspace(self.ymin, self.ymax, ny)
lons1, lats1 = self(x1,y1,inverse=True)
lats1 = np.where(lats1 < -89.999999, -89.999999, lats1)
lats1 = np.where(lats1 > 89.999999, 89.999999, lats1)
for j,lat in enumerate(lats1):
xmax,ymax = self(lonright,lat)
xmin,ymin = self(lonleft,lat)
mask[j,:,0] = np.logical_or(x[j,:]>xmax,x[j,:]<xmin)
for k in range(1,4):
mask[:,:,k] = mask[:,:,0]
self._bm_rgba_warped = \
ma.masked_array(self._bm_rgba_warped,mask=mask)
# make points outside projection limb transparent.
# FIXME: Probably not needed anymore
self._bm_rgba_warped = self._bm_rgba_warped.filled(0.)
# plot warped rgba image.
im = self.imshow(self._bm_rgba_warped,ax=ax,**kwargs)
# for hammer projection, use clip path defined by
# projection limb (patch created in drawmapboundary).
# FIXME: Is this now redundant?
if self.projection == 'hammer':
if not self._mapboundarydrawn:
self.drawmapboundary(color='none',linewidth=None)
im.set_clip_path(self._mapboundarydrawn)
else:
# bmproj True, no interpolation necessary.
im = self.imshow(self._bm_rgba,ax=ax,**kwargs)
# clip to map limbs
im,c = self._cliplimb(ax,im)
return im
[docs] def arcgisimage(self, server="http://server.arcgisonline.com/ArcGIS",
service="World_Imagery", xpixels=400, ypixels=None,
dpi=96, cachedir=None, verbose=False, **kwargs):
r"""Display background image using ArcGIS Server REST API.
In order to use this method, the :class:`Basemap` instance
must be created using the ``epsg`` keyword to define the
map projection, unless the "cyl" projection is used (in
which case the EPSG code 4326 is assumed).
Parameters
----------
server : str, optional
base URL of the web map server
service : str, optional
service (image type) hosted by the server
xpixels : int, optional
requested number of image pixels in the `x`-direction
ypixels : int, optional
requested number of image pixels in the `y`-direction;
if not given, it is inferred from ``xpixels`` and the
aspect ratio of the map projection region
dpi : int, optional
device resolution of the exported image
cachedir : str, optional
if given, directory to use as cache folder for the images
retrieved from the server
verbose : bool, optional
if True, print debugging information
\**kwargs : dict, optional
keyword-only arguments; currently, only ``ax`` is supported
to override the default :class:`matplotlib.axes.Axes`
instance
Returns
-------
aximg : matplotlib.image.AxesImage
image axes instance
"""
# Fix PIL import on some versions of OSX and scipy.
try:
from PIL import Image
except ImportError:
try:
import Image
except ImportError:
raise ImportError("arcgisimage method requires PIL "
"(http://pillow.readthedocs.io)")
if not hasattr(self, "epsg"):
raise ValueError("the Basemap instance must be created using "
"an EPSG code (http://spatialreference.org) "
"in order to use the wmsmap method")
ax = kwargs.pop("ax", None) or self._check_ax()
# Find the `(x, y)` values at the corner points.
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=FutureWarning)
p = pyproj.Proj(init="epsg:%s" % self.epsg, preserve_units=True)
xmin, ymin = p(self.llcrnrlon, self.llcrnrlat)
xmax, ymax = p(self.urcrnrlon, self.urcrnrlat)
if self.projection in _cylproj:
dateline = _geoslib.Point(self(180., 0.5 * (self.llcrnrlat + self.urcrnrlat)))
if dateline.within(self._boundarypolyxy):
raise ValueError("arcgisimage cannot handle images that cross "
"the dateline for cylindrical projections")
# If ypixels is not given, compute it with xpixels and aspect ratio.
if ypixels is None:
ypixels = int(self.aspect * xpixels)
# Construct a URL using the ArcGIS Server REST API.
basemap_url = "".join([
"%s/rest/services/%s/MapServer/export?",
"bbox=%s,%s,%s,%s&",
"bboxSR=%s&",
"imageSR=%s&",
"size=%s,%s&",
"dpi=%s&",
"format=png32&",
"transparent=true&",
"f=image",
]) % (server, service, xmin, ymin, xmax, ymax, self.epsg, self.epsg, xpixels, ypixels, dpi)
# Print URL in verbose mode.
if verbose: # pragma: no cover
print(basemap_url)
# Try to return fast if cache is enabled.
if cachedir is not None:
# Generate a filename for the cached file.
filename = "%s-bbox-%s-%s-%s-%s-bboxsr%s-imagesr%s-size-%s-%s-dpi%s.png" % \
(service, xmin, ymin, xmax, ymax, self.epsg, self.epsg, xpixels, ypixels, dpi)
# Return fast if the image is already in the cache.
cache_path = os.path.join(cachedir, filename)
if os.path.isfile(cache_path):
if verbose: # pragma: no cover
print("Image already in cache")
img = Image.open(cache_path)
return self.imshow(img, ax=ax, origin="upper")
# Retrieve image from the remote server.
import contextlib
conn = urlopen(basemap_url)
with contextlib.closing(conn):
img = Image.open(conn)
# Save to cache if requested.
if cachedir is not None:
# Check if the cache directory exists, if not create it.
if not os.path.exists(cachedir):
os.makedirs(cachedir)
img.save(cache_path)
return self.imshow(img, ax=ax, origin="upper")
[docs] def wmsimage(self,server,\
xpixels=400,ypixels=None,\
format='png',alpha=None,verbose=False,**kwargs):
r"""
Retrieve an image using from a WMS server using the
Open Geospatial Consortium (OGC) standard interface
and display on the map. Requires OWSLib
(http://pypi.python.org/pypi/OWSLib).
In order to use this method, the Basemap instance must be
created using the ``epsg`` keyword to define the map projection, unless
the ``cyl`` projection is used (in which case the epsg code 4326 is
assumed).
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
server WMS server URL.
xpixels requested number of image pixels in x-direction
(default 400).
ypixels requested number of image pixels in y-direction.
Default (None) is to infer the number from
from xpixels and the aspect ratio of the
map projection region.
format desired image format (default 'png')
alpha The alpha blending value,
between 0 (transparent) and 1 (opaque) (default None)
verbose if True, print WMS server info (default
False).
\**kwargs extra keyword arguments passed on to
OWSLib.wms.WebMapService.getmap.
============== ====================================================
Extra keyword ``ax`` can be used to override the default axis instance.
returns a matplotlib.image.AxesImage instance.
"""
try:
from owslib.wms import WebMapService
except ImportError:
raise ImportError('OWSLib required to use wmsimage method')
import io
ax = kwargs.pop('ax', None) or self._check_ax()
if not hasattr(self,'epsg'):
raise ValueError("the Basemap instance must be created using "
"an EPSG code (http://spatialreference.org) "
"in order to use the wmsmap method")
if 'layers' not in kwargs:
raise ValueError('no layers specified')
# find the x,y values at the corner points.
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=FutureWarning)
p = pyproj.Proj(init="epsg:%s" % self.epsg, preserve_units=True)
xmin,ymin = p(self.llcrnrlon,self.llcrnrlat)
xmax,ymax = p(self.urcrnrlon,self.urcrnrlat)
if self.projection in _cylproj:
Dateline =\
_geoslib.Point(self(180.,0.5*(self.llcrnrlat+self.urcrnrlat)))
hasDateline = Dateline.within(self._boundarypolyxy)
if hasDateline:
raise ValueError("wmsimage cannot handle images that cross "
"the dateline for cylindrical projections")
if self.projection == 'cyl':
xmin = (180./np.pi)*xmin; xmax = (180./np.pi)*xmax
ymin = (180./np.pi)*ymin; ymax = (180./np.pi)*ymax
# ypixels not given, find by scaling xpixels by the map aspect ratio.
if ypixels is None:
ypixels = int(self.aspect*xpixels)
if verbose: print(server)
wms = WebMapService(server)
if verbose:
print('id: %s, version: %s' %
(wms.identification.type,wms.identification.version))
print('title: %s, abstract: %s' %
(wms.identification.title,wms.identification.abstract))
print('available layers:')
print(list(wms.contents))
print('projection options:')
print(wms[kwargs['layers'][0]].crsOptions)
# remove keys from kwargs that are over-ridden
for k in ['format','bbox','service','size','srs']:
if 'format' in kwargs: del kwargs['format']
img = wms.getmap(service='wms',bbox=(xmin,ymin,xmax,ymax),
size=(xpixels,ypixels),format='image/%s'%format,
srs='EPSG:%s' % self.epsg, **kwargs)
# return AxesImage instance.
# this works for png and jpeg.
return self.imshow(imread(io.BytesIO(img.read()),
format=format),origin='upper',alpha=alpha,ax=ax)
[docs] def drawmapscale(self,lon,lat,lon0,lat0,length,barstyle='simple',\
units='km',fontsize=9,yoffset=None,labelstyle='simple',\
fontcolor='k',fillcolor1='w',fillcolor2='k',ax=None,\
format='%d',zorder=None,linecolor=None,linewidth=None):
"""
Draw a map scale at ``lon,lat`` of length ``length``
representing distance in the map
projection coordinates at ``lon0,lat0``.
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
units the units of the length argument (Default km).
barstyle ``simple`` or ``fancy`` (roughly corresponding
to the styles provided by Generic Mapping Tools).
Default ``simple``.
fontsize for map scale annotations, default 9.
fontcolor for map scale annotations, default black.
labelstyle ``simple`` (default) or ``fancy``. For
``fancy`` the map scale factor (ratio betwee
the actual distance and map projection distance
at lon0,lat0) and the value of lon0,lat0 are also
displayed on the top of the scale bar. For
``simple``, just the units are display on top
and the distance below the scale bar.
If equal to False, plot an empty label.
format a string formatter to format numeric values
yoffset yoffset controls how tall the scale bar is,
and how far the annotations are offset from the
scale bar. Default is 0.02 times the height of
the map (0.02*(self.ymax-self.ymin)).
fillcolor1(2) colors of the alternating filled regions
(default white and black). Only relevant for
'fancy' barstyle.
zorder sets the zorder for the map scale.
linecolor sets the color of the scale, by default, fontcolor
is used
linewidth linewidth for scale and ticks
============== ====================================================
Extra keyword ``ax`` can be used to override the default axis instance.
"""
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# not valid for cylindrical projection
if self.projection == 'cyl':
raise ValueError("cannot draw map scale for projection='cyl'")
# convert length to meters
lenlab = length
if units == 'km':
length = length*1000
elif units == 'mi':
length = length*1609.344
elif units == 'nmi':
length = length*1852
elif units == 'ft':
length = length*0.3048
elif units != 'm':
raise KeyError("units must be 'm' (meters), 'km' (kilometers), "
"'mi' (miles), 'nmi' (nautical miles), or 'ft' (feet)")
# reference point and center of scale.
x0,y0 = self(lon0,lat0)
xc,yc = self(lon,lat)
# make sure lon_0 between -180 and 180
lon_0 = ((lon0+360) % 360) - 360
degchar = b"\xc2\xb0".decode("utf-8")
if lat0>0:
if lon>0:
lonlatstr = '%g%sN, %g%sE' % (lat0, degchar, lon_0, degchar)
elif lon<0:
lonlatstr = '%g%sN, %g%sW' % (lat0, degchar, lon_0, degchar)
else:
lonlatstr = '%g%s, %g%sW' % (lat0, degchar, lon_0, degchar)
else:
if lon>0:
lonlatstr = '%g%sS, %g%sE' % (lat0, degchar, lon_0, degchar)
elif lon<0:
lonlatstr = '%g%sS, %g%sW' % (lat0, degchar, lon_0, degchar)
else:
lonlatstr = '%g%sS, %g%s' % (lat0, degchar, lon_0, degchar)
# left edge of scale
lon1,lat1 = self(x0-length/2,y0,inverse=True)
x1,y1 = self(lon1,lat1)
# right edge of scale
lon4,lat4 = self(x0+length/2,y0,inverse=True)
x4,y4 = self(lon4,lat4)
x1 = x1-x0+xc; y1 = y1-y0+yc
x4 = x4-x0+xc; y4 = y4-y0+yc
if x1 > 1.e20 or x4 > 1.e20 or y1 > 1.e20 or y4 > 1.e20:
raise ValueError("scale bar positioned outside projection limb")
# scale factor for true distance
gc = pyproj.Geod(a=self.rmajor,b=self.rminor)
az12,az21,dist = gc.inv(lon1,lat1,lon4,lat4)
scalefact = dist/length
# label to put on top of scale bar.
if labelstyle=='simple':
labelstr = units
elif labelstyle == 'fancy':
labelstr = units+" (scale factor %4.2f at %s)"%(scalefact,lonlatstr)
elif labelstyle == False:
labelstr = ''
else:
raise KeyError("labelstyle must be 'simple' or 'fancy'")
# default y offset is 2 percent of map height.
if yoffset is None: yoffset = 0.02*(self.ymax-self.ymin)
rets = [] # will hold all plot objects generated.
# set linecolor
if linecolor is None:
linecolor = fontcolor
# 'fancy' style
if barstyle == 'fancy':
#we need 5 sets of x coordinates (in map units)
#quarter scale
lon2,lat2 = self(x0-length/4,y0,inverse=True)
x2,y2 = self(lon2,lat2)
x2 = x2-x0+xc; y2 = y2-y0+yc
#three quarter scale
lon3,lat3 = self(x0+length/4,y0,inverse=True)
x3,y3 = self(lon3,lat3)
x3 = x3-x0+xc; y3 = y3-y0+yc
#plot top line
ytop = yc+yoffset/2
ybottom = yc-yoffset/2
ytick = ybottom - yoffset/2
ytext = ytick - yoffset/2
rets.append(self.plot([x1,x4],[ytop,ytop],color=linecolor, linewidth=linewidth)[0])
#plot bottom line
rets.append(self.plot([x1,x4],[ybottom,ybottom],color=linecolor, linewidth=linewidth)[0])
#plot left edge
rets.append(self.plot([x1,x1],[ybottom,ytop],color=linecolor, linewidth=linewidth)[0])
#plot right edge
rets.append(self.plot([x4,x4],[ybottom,ytop],color=linecolor, linewidth=linewidth)[0])
#make a filled black box from left edge to 1/4 way across
rets.append(ax.fill([x1,x2,x2,x1,x1],[ytop,ytop,ybottom,ybottom,ytop],\
ec=fontcolor,fc=fillcolor1)[0])
#make a filled white box from 1/4 way across to 1/2 way across
rets.append(ax.fill([x2,xc,xc,x2,x2],[ytop,ytop,ybottom,ybottom,ytop],\
ec=fontcolor,fc=fillcolor2)[0])
#make a filled white box from 1/2 way across to 3/4 way across
rets.append(ax.fill([xc,x3,x3,xc,xc],[ytop,ytop,ybottom,ybottom,ytop],\
ec=fontcolor,fc=fillcolor1)[0])
#make a filled white box from 3/4 way across to end
rets.append(ax.fill([x3,x4,x4,x3,x3],[ytop,ytop,ybottom,ybottom,ytop],\
ec=fontcolor,fc=fillcolor2)[0])
#plot 3 tick marks at left edge, center, and right edge
rets.append(self.plot([x1,x1],[ytick,ybottom],color=linecolor, linewidth=linewidth)[0])
rets.append(self.plot([xc,xc],[ytick,ybottom],color=linecolor, linewidth=linewidth)[0])
rets.append(self.plot([x4,x4],[ytick,ybottom],color=linecolor, linewidth=linewidth)[0])
#label 3 tick marks
rets.append(ax.text(x1,ytext,format % (0),\
horizontalalignment='center',\
verticalalignment='top',\
fontsize=fontsize,color=fontcolor))
rets.append(ax.text(xc,ytext,format % (0.5*lenlab),\
horizontalalignment='center',\
verticalalignment='top',\
fontsize=fontsize,color=fontcolor))
rets.append(ax.text(x4,ytext,format % (lenlab),\
horizontalalignment='center',\
verticalalignment='top',\
fontsize=fontsize,color=fontcolor))
#put units, scale factor on top
rets.append(ax.text(xc,ytop+yoffset/2,labelstr,\
horizontalalignment='center',\
verticalalignment='bottom',\
fontsize=fontsize,color=fontcolor))
# 'simple' style
elif barstyle == 'simple':
rets.append(self.plot([x1,x4],[yc,yc],color=linecolor, linewidth=linewidth)[0])
rets.append(self.plot([x1,x1],[yc-yoffset,yc+yoffset],color=linecolor, linewidth=linewidth)[0])
rets.append(self.plot([x4,x4],[yc-yoffset,yc+yoffset],color=linecolor, linewidth=linewidth)[0])
rets.append(ax.text(xc,yc-yoffset,format % lenlab,\
verticalalignment='top',horizontalalignment='center',\
fontsize=fontsize,color=fontcolor))
#put units, scale factor on top
rets.append(ax.text(xc,yc+yoffset,labelstr,\
horizontalalignment='center',\
verticalalignment='bottom',\
fontsize=fontsize,color=fontcolor))
else:
raise KeyError("barstyle must be 'simple' or 'fancy'")
if zorder is not None:
for ret in rets:
try:
ret.set_zorder(zorder)
except:
pass
return rets
[docs] def colorbar(self,mappable=None,location='right',size="5%",pad='2%',fig=None,ax=None,**kwargs):
r"""
Add colorbar to axes associated with a map.
The colorbar axes instance is created using the axes_grid toolkit.
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
mappable the Image, ContourSet, etc. to which the colorbar
applies. Default None, matplotlib.pyplot.gci() is
used to retrieve the current image mappable.
location where to put colorbar ('top','bottom','left','right')
Default 'right'.
size width of colorbar axes (string 'N%', where N is
an integer describing the fractional width of the parent
axes). Default '5%'.
pad Padding between parent axes and colorbar axes in
same units as size. Default '2%'.
fig Figure instance the map axes instance is associated
with. Default None, and matplotlib.pyplot.gcf() is used
to retrieve the current active figure instance.
ax The axes instance which the colorbar will be
associated with. Default None, searches for self.ax,
and if None uses matplotlib.pyplot.gca().
\**kwargs extra keyword arguments passed on to
colorbar method of the figure instance.
============== ====================================================
Returns a matplotlib colorbar instance.
"""
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# get current figure instance (if none specified).
if fig is None or mappable is None:
import matplotlib.pyplot as plt
if fig is None:
fig = plt.gcf()
# get current mappable if none specified.
if mappable is None:
mappable = plt.gci()
# create colorbar axes uses axes_grid toolkit.
divider = make_axes_locatable(ax)
if location in ['left','right']:
orientation = 'vertical'
elif location in ['top','bottom']:
orientation = 'horizontal'
else:
raise ValueError('location must be top,bottom,left or right')
cax = divider.append_axes(location, size=size, pad=pad)
# create colorbar.
cb = fig.colorbar(mappable,orientation=orientation,cax=cax,**kwargs)
fig.sca(ax) # reset parent axes as current axes.
return cb
[docs] def nightshade(self,date,color="k",delta=0.25,alpha=0.5,ax=None,zorder=2):
"""
Shade the regions of the map that are in darkness at the time
specifed by ``date``. ``date`` is a datetime instance,
assumed to be UTC.
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
color color to shade night regions (default black).
delta day/night terminator is computed with a
a resolution of ``delta`` degrees (default 0.25).
alpha alpha transparency for shading (default 0.5, so
map background shows through).
zorder zorder for shading (default 2).
============== ====================================================
Extra keyword ``ax`` can be used to override the default axis instance.
returns a matplotlib.contour.ContourSet instance.
"""
from .solar import daynight_grid
# make sure date is UTC, or naive with repect to time zones
if date.utcoffset():
raise ValueError('datetime instance must be UTC, not {0}'.format(date.tzname()))
# get current axes instance (if none specified).
ax = ax or self._check_ax()
# create grid of day=0, night=1
lons,lats,daynight = daynight_grid(date,delta,self.lonmin,self.lonmax)
x,y = self(lons,lats)
# contour the day-night grid, coloring the night area
# with the specified color and transparency.
CS = self.contourf(x,y,daynight,1,colors=[color],alpha=alpha,ax=ax)
if isinstance(CS, Artist):
# Since MPL 3.8, `QuadContourSet` objects are `Artist` objects too.
CS.set_zorder(zorder)
# clip to map limbs
CS, c = self._cliplimb(ax, CS)
else:
# set zorder on ContourSet collections show night shading
# is on top.
for c in CS.collections:
c.set_zorder(zorder)
# clip to map limbs
CS.collections, c = self._cliplimb(ax, CS.collections)
return CS
def _check_ax(self):
"""
Returns the axis on which to draw.
Returns self.ax, or if self.ax=None returns plt.gca().
"""
if self.ax is None:
try:
ax = plt.gca()
except:
import matplotlib.pyplot as plt
ax = plt.gca()
# associate an axes instance with this Basemap instance
# the first time this method is called.
#self.ax = ax
else:
ax = self.ax
return ax
def _ax_plt_from_kw(self, kw):
"""
Return (ax, plt), where ax is the current axes, and plt is
None or a reference to the pyplot module.
plt will be None if ax was popped from kw or taken from self.ax;
otherwise, pyplot was used and is returned.
"""
plt = None
_ax = kw.pop('ax', None)
if _ax is None:
_ax = self.ax
if _ax is None:
import matplotlib.pyplot as plt
_ax = plt.gca()
return _ax, plt
[docs] def shiftdata(self,lonsin,datain=None,lon_0=None,fix_wrap_around=True):
"""
Shift longitudes (and optionally data) so that they match map projection region.
Only valid for cylindrical/pseudo-cylindrical global projections and data
on regular lat/lon grids. longitudes and data can be 1-d or 2-d, if 2-d
it is assumed longitudes are 2nd (rightmost) dimension.
.. tabularcolumns:: |l|L|
================ ======================================================
Arguments Description
================ ======================================================
lonsin original 1-d or 2-d longitudes.
================ ======================================================
.. tabularcolumns:: |l|L|
================ ======================================================
Keywords Description
================ ======================================================
datain original 1-d or 2-d data. Default None.
lon_0 center of map projection region. Defaut None,
given by current map projection.
fix_wrap_around if True reindex (if required) longitudes (and data) to
avoid jumps caused by remapping of longitudes of
points from outside of the [lon_0-180, lon_0+180]
interval back into the interval.
If False do not reindex longitudes and data, but do
make sure that longitudes are in the
[lon_0-180, lon_0+180] range.
================ ======================================================
if datain given, returns ``dataout,lonsout`` (data and longitudes shifted to fit in interval
[lon_0-180,lon_0+180]), otherwise just returns longitudes. If
transformed longitudes lie outside map projection region, data is
masked and longitudes are set to 1.e30.
"""
if lon_0 is None and 'lon_0' not in self.projparams:
raise ValueError('lon_0 keyword must be provided')
lonsin = np.asarray(lonsin)
if lonsin.ndim not in [1,2]:
raise ValueError('1-d or 2-d longitudes required')
if datain is not None:
# if it's a masked array, leave it alone.
if not ma.isMA(datain): datain = np.asarray(datain)
if datain.ndim not in [1,2]:
raise ValueError('1-d or 2-d data required')
if lon_0 is None:
lon_0 = self.projparams['lon_0']
# 2-d data.
if lonsin.ndim == 2:
nlats = lonsin.shape[0]
nlons = lonsin.shape[1]
lonsin1 = lonsin[0,:]
lonsin1 = np.where(lonsin1 > lon_0+180, lonsin1-360 ,lonsin1)
lonsin1 = np.where(lonsin1 < lon_0-180, lonsin1+360 ,lonsin1)
if nlons > 1:
londiff = np.abs(lonsin1[0:-1]-lonsin1[1:])
londiff_sort = np.sort(londiff)
thresh = 360.-londiff_sort[-2] if nlons > 2 else 360.-londiff_sort[-1]
itemindex = nlons-np.where(londiff>=thresh)[0]
else:
lonsin[0, :] = lonsin1
itemindex = 0
# if no shift necessary, itemindex will be
# empty, so don't do anything
if fix_wrap_around and np.all(itemindex != 0) and itemindex.size:
# check to see if cyclic (wraparound) point included
# if so, remove it.
if np.abs(lonsin1[0]-lonsin1[-1]) < 1.e-4:
hascyclic = True
lonsin_save = lonsin.copy()
lonsin = lonsin[:,1:]
if datain is not None:
datain_save = datain.copy()
datain = datain[:,1:]
else:
hascyclic = False
lonsin = np.where(lonsin > lon_0+180, lonsin-360 ,lonsin)
lonsin = np.where(lonsin < lon_0-180, lonsin+360 ,lonsin)
lonsin = np.roll(lonsin,itemindex-1,axis=1)
if datain is not None:
# np.roll works on ndarrays and on masked arrays
datain = np.roll(datain,itemindex-1,axis=1)
# add cyclic point back at beginning.
if hascyclic:
lonsin_save[:,1:] = lonsin
lonsin_save[:,0] = lonsin[:,-1]-360.
lonsin = lonsin_save
if datain is not None:
datain_save[:,1:] = datain
datain_save[:,0] = datain[:,-1]
datain = datain_save
# 1-d data.
elif lonsin.ndim == 1:
nlons = len(lonsin)
lonsin = np.where(lonsin > lon_0+180, lonsin-360 ,lonsin)
lonsin = np.where(lonsin < lon_0-180, lonsin+360 ,lonsin)
if nlons > 1:
londiff = np.abs(lonsin[0:-1]-lonsin[1:])
londiff_sort = np.sort(londiff)
thresh = 360.-londiff_sort[-2] if nlons > 2 else 360.0 - londiff_sort[-1]
itemindex = len(lonsin)-np.where(londiff>=thresh)[0]
else:
itemindex = 0
if fix_wrap_around and np.all(itemindex != 0) and itemindex.size:
# check to see if cyclic (wraparound) point included
# if so, remove it.
if np.abs(lonsin[0]-lonsin[-1]) < 1.e-4:
hascyclic = True
lonsin_save = lonsin.copy()
lonsin = lonsin[1:]
if datain is not None:
datain_save = datain.copy()
datain = datain[1:]
else:
hascyclic = False
lonsin = np.roll(lonsin,itemindex-1)
if datain is not None:
datain = np.roll(datain,itemindex-1)
# add cyclic point back at beginning.
if hascyclic:
lonsin_save[1:] = lonsin
lonsin_save[0] = lonsin[-1]-360.
lonsin = lonsin_save
if datain is not None:
datain_save[1:] = datain
datain_save[0] = datain[-1]
datain = datain_save
# mask points outside
# map region so they don't wrap back in the domain.
mask = np.logical_or(lonsin<lon_0-180,lonsin>lon_0+180)
lonsin = np.where(mask,1.e30,lonsin)
if datain is not None and mask.any():
datain = ma.masked_where(mask, datain)
if datain is not None:
return lonsin, datain
else:
return lonsin
### End of Basemap class
def _searchlist(a,x):
"""
like bisect, but works for lists that are not sorted,
and are not in increasing order.
returns -1 if x does not fall between any two elements"""
# make sure x is a float (and not an array scalar)
x = float(x)
itemprev = a[0]
nslot = -1
eps = 180.
for n,item in enumerate(a[1:]):
if item < itemprev:
if itemprev-item>eps:
if ((x>itemprev and x<=360.) or (x<item and x>=0.)):
nslot = n+1
break
elif x <= itemprev and x > item and itemprev:
nslot = n+1
break
else:
if item-itemprev>eps:
if ((x<itemprev and x>=0.) or (x>item and x<=360.)):
nslot = n+1
break
elif x >= itemprev and x < item:
nslot = n+1
break
itemprev = item
return nslot
[docs]def interp(datain,xin,yin,xout,yout,checkbounds=False,masked=False,order=1):
"""
Interpolate data (``datain``) on a rectilinear grid (with x = ``xin``
y = ``yin``) to a grid with x = ``xout``, y= ``yout``.
.. tabularcolumns:: |l|L|
============== ====================================================
Arguments Description
============== ====================================================
datain a rank-2 array with 1st dimension corresponding to
y, 2nd dimension x.
xin, yin rank-1 arrays containing x and y of
datain grid in increasing order.
xout, yout rank-2 arrays containing x and y of desired output grid.
============== ====================================================
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
checkbounds If True, values of xout and yout are checked to see
that they lie within the range specified by xin
and xin.
If False, and xout,yout are outside xin,yin,
interpolated values will be clipped to values on
boundary of input grid (xin,yin)
Default is False.
masked If True, points outside the range of xin and yin
are masked (in a masked array).
If masked is set to a number, then
points outside the range of xin and yin will be
set to that number. Default False.
order 0 for nearest-neighbor interpolation, 1 for
bilinear interpolation, 3 for cublic spline
(default 1). order=3 requires scipy.ndimage.
============== ====================================================
.. note::
If datain is a masked array and order=1 (bilinear interpolation) is
used, elements of dataout will be masked if any of the four surrounding
points in datain are masked. To avoid this, do the interpolation in two
passes, first with order=1 (producing dataout1), then with order=0
(producing dataout2). Then replace all the masked values in dataout1
with the corresponding elements in dataout2 (using numpy.where).
This effectively uses nearest neighbor interpolation if any of the
four surrounding points in datain are masked, and bilinear interpolation
otherwise.
Returns ``dataout``, the interpolated data on the grid ``xout, yout``.
"""
# xin and yin must be monotonically increasing.
if xin[-1]-xin[0] < 0 or yin[-1]-yin[0] < 0:
raise ValueError('xin and yin must be increasing!')
if xout.shape != yout.shape:
raise ValueError('xout and yout must have same shape!')
# check that xout,yout are
# within region defined by xin,yin.
if checkbounds:
if xout.min() < xin.min() or \
xout.max() > xin.max() or \
yout.min() < yin.min() or \
yout.max() > yin.max():
raise ValueError('yout or xout outside range of yin or xin')
# compute grid coordinates of output grid.
delx = xin[1:]-xin[0:-1]
dely = yin[1:]-yin[0:-1]
if max(delx)-min(delx) < 1.e-4 and max(dely)-min(dely) < 1.e-4:
# regular input grid.
xcoords = (len(xin)-1)*(xout-xin[0])/(xin[-1]-xin[0])
ycoords = (len(yin)-1)*(yout-yin[0])/(yin[-1]-yin[0])
else:
# irregular (but still rectilinear) input grid.
xoutflat = xout.flatten(); youtflat = yout.flatten()
ix = (np.searchsorted(xin,xoutflat)-1).tolist()
iy = (np.searchsorted(yin,youtflat)-1).tolist()
xoutflat = xoutflat.tolist(); xin = xin.tolist()
youtflat = youtflat.tolist(); yin = yin.tolist()
xcoords = []; ycoords = []
for n,i in enumerate(ix):
if i < 0:
xcoords.append(-1) # outside of range on xin (lower end)
elif i >= len(xin)-1:
xcoords.append(len(xin)) # outside range on upper end.
else:
xcoords.append(float(i)+(xoutflat[n]-xin[i])/(xin[i+1]-xin[i]))
for m,j in enumerate(iy):
if j < 0:
ycoords.append(-1) # outside of range of yin (on lower end)
elif j >= len(yin)-1:
ycoords.append(len(yin)) # outside range on upper end
else:
ycoords.append(float(j)+(youtflat[m]-yin[j])/(yin[j+1]-yin[j]))
xcoords = np.reshape(xcoords,xout.shape)
ycoords = np.reshape(ycoords,yout.shape)
# data outside range xin,yin will be clipped to
# values on boundary.
if masked:
xmask = np.logical_or(np.less(xcoords,0),np.greater(xcoords,len(xin)-1))
ymask = np.logical_or(np.less(ycoords,0),np.greater(ycoords,len(yin)-1))
xymask = np.logical_or(xmask,ymask)
xcoords = np.clip(xcoords,0,len(xin)-1)
ycoords = np.clip(ycoords,0,len(yin)-1)
# interpolate to output grid using bilinear interpolation.
if order == 1:
xi = xcoords.astype(np.int32)
yi = ycoords.astype(np.int32)
xip1 = xi+1
yip1 = yi+1
xip1 = np.clip(xip1,0,len(xin)-1)
yip1 = np.clip(yip1,0,len(yin)-1)
delx = xcoords-xi.astype(np.float32)
dely = ycoords-yi.astype(np.float32)
dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
delx*dely*datain[yip1,xip1] + \
(1.-delx)*dely*datain[yip1,xi] + \
delx*(1.-dely)*datain[yi,xip1]
elif order == 0:
xcoordsi = np.around(xcoords).astype(np.int32)
ycoordsi = np.around(ycoords).astype(np.int32)
dataout = datain[ycoordsi,xcoordsi]
elif order == 3:
try:
from scipy.ndimage import map_coordinates
except ImportError:
raise ValueError('scipy.ndimage must be installed if order=3')
coords = [ycoords,xcoords]
dataout = map_coordinates(datain,coords,order=3,mode='nearest')
else:
raise ValueError('order keyword must be 0, 1 or 3')
if masked:
newmask = ma.mask_or(ma.getmask(dataout), xymask)
dataout = ma.masked_array(dataout, mask=newmask)
if not isinstance(masked, bool):
dataout = dataout.filled(masked)
return dataout
[docs]def shiftgrid(lon0,datain,lonsin,start=True,cyclic=360.0):
"""
Shift global lat/lon grid east or west.
.. tabularcolumns:: |l|L|
============== ====================================================
Arguments Description
============== ====================================================
lon0 starting longitude for shifted grid
(ending longitude if start=False). lon0 must be on
input grid (within the range of lonsin).
datain original data with longitude the right-most
dimension.
lonsin original longitudes.
============== ====================================================
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
start if True, lon0 represents the starting longitude
of the new grid. if False, lon0 is the ending
longitude. Default True.
cyclic width of periodic domain (default 360)
============== ====================================================
returns ``dataout,lonsout`` (data and longitudes on shifted grid).
"""
if np.fabs(lonsin[-1]-lonsin[0]-cyclic) > 1.e-4:
# Use all data instead of raise ValueError, 'cyclic point not included'
start_idx = 0
else:
# If cyclic, remove the duplicate point
start_idx = 1
if lon0 < lonsin[0] or lon0 > lonsin[-1]:
raise ValueError('lon0 outside of range of lonsin')
i0 = np.argmin(np.fabs(lonsin-lon0))
i0_shift = len(lonsin)-i0
if ma.isMA(datain):
dataout = ma.zeros(datain.shape,datain.dtype)
else:
dataout = np.zeros(datain.shape,datain.dtype)
if ma.isMA(lonsin):
lonsout = ma.zeros(lonsin.shape,lonsin.dtype)
else:
lonsout = np.zeros(lonsin.shape,lonsin.dtype)
if start:
lonsout[0:i0_shift] = lonsin[i0:]
else:
lonsout[0:i0_shift] = lonsin[i0:]-cyclic
dataout[...,0:i0_shift] = datain[...,i0:]
if start:
lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]+cyclic
else:
lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]
dataout[...,i0_shift:] = datain[...,start_idx:i0+start_idx]
return dataout,lonsout
[docs]def addcyclic(*arr,**kwargs):
"""
Adds cyclic (wraparound) points in longitude to one or several arrays,
the last array being longitudes in degrees. e.g.
``data1out, data2out, lonsout = addcyclic(data1,data2,lons)``
============== ====================================================
Keywords Description
============== ====================================================
axis the dimension representing longitude (default -1,
or right-most)
cyclic width of periodic domain (default 360)
============== ====================================================
"""
# get (default) keyword arguments
axis = kwargs.get('axis',-1)
cyclic = kwargs.get('cyclic',360)
# define functions
def _addcyclic(a):
"""addcyclic function for a single data array"""
npsel = np.ma if np.ma.is_masked(a) else np
slicer = [slice(None)] * np.ndim(a)
try:
slicer[axis] = slice(0, 1)
except IndexError:
raise ValueError('The specified axis does not correspond to an '
'array dimension.')
return npsel.concatenate((a,a[tuple(slicer)]),axis=axis)
def _addcyclic_lon(a):
"""addcyclic function for a single longitude array"""
# select the right numpy functions
npsel = np.ma if np.ma.is_masked(a) else np
# get cyclic longitudes
clon = (np.take(a,[0],axis=axis)
+ cyclic * np.sign(np.diff(np.take(a,[0,-1],axis=axis),axis=axis)))
# ensure the values do not exceed cyclic
clonmod = npsel.where(clon<=cyclic,clon,np.mod(clon,cyclic))
return npsel.concatenate((a,clonmod),axis=axis)
# process array(s)
if len(arr) == 1:
return _addcyclic_lon(arr[-1])
else:
return list(map(_addcyclic,arr[:-1])) + [_addcyclic_lon(arr[-1])]
def _choosecorners(width,height,**kwargs):
"""
private function to determine lat/lon values of projection region corners,
given width and height of projection region in meters.
"""
p = pyproj.Proj(kwargs)
urcrnrlon, urcrnrlat = p(0.5*width,0.5*height, inverse=True)
llcrnrlon, llcrnrlat = p(-0.5*width,-0.5*height, inverse=True)
corners = llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat
# test for invalid projection points on output
if llcrnrlon > 1.e20 or urcrnrlon > 1.e20:
raise ValueError('width and/or height too large for this projection, try smaller values')
else:
return corners
def _choosecornersllur(llcrnrx, llcrnry, urcrnrx, urcrnry,**kwargs):
"""
private function to determine lat/lon values of projection region corners,
given width and height of projection region in meters.
"""
p = pyproj.Proj(kwargs)
urcrnrlon, urcrnrlat = p(urcrnrx, urcrnry, inverse=True)
llcrnrlon, llcrnrlat = p(llcrnrx, llcrnry, inverse=True)
corners = llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat
# test for invalid projection points on output
if llcrnrlon > 1.e20 or urcrnrlon > 1.e20:
raise ValueError('width and/or height too large for this projection, try smaller values')
else:
return corners
[docs]def maskoceans(lonsin,latsin,datain,inlands=True,resolution='l',grid=5):
"""
mask data (``datain``), defined on a grid with latitudes ``latsin``
longitudes ``lonsin`` so that points over water will not be plotted.
.. tabularcolumns:: |l|L|
============== ====================================================
Arguments Description
============== ====================================================
lonsin, latsin rank-2 arrays containing longitudes and latitudes of
grid.
datain rank-2 input array on grid defined by ``lonsin`` and
``latsin``.
inlands if False, masked only ocean points and not inland
lakes (Default True).
resolution gshhs coastline resolution used to define land/sea
mask (default 'l', available 'c','l','i','h' or 'f')
grid land/sea mask grid spacing in minutes (Default 5;
10, 2.5 and 1.25 are also available).
============== ====================================================
returns a masked array the same shape as datain with "wet" points masked.
"""
# read in land/sea mask.
lsmask_lons, lsmask_lats, lsmask =\
_readlsmask(lakes=inlands,resolution=resolution,grid=grid)
# nearest-neighbor interpolation to output grid.
lsmasko = interp(lsmask,lsmask_lons,lsmask_lats,lonsin,latsin,masked=True,order=0)
# mask input data.
mask = lsmasko == 0
return ma.masked_array(datain,mask=mask)
def _readlsmask(lakes=True,resolution='l',grid=5):
# read in land/sea mask.
if grid == 10:
nlons = 2160
elif grid == 5:
nlons = 4320
elif grid == 2.5:
nlons = 8640
elif grid == 1.25:
nlons = 17280
else:
raise ValueError('grid for land/sea mask must be 10,5,2.5 or 1.25')
nlats = nlons//2
import gzip
lsmaskf =\
gzip.open(os.path.join(basemap_datadir,'lsmask_%smin_%s.bin' %\
(grid,resolution)), 'rb')
lsmask =\
np.reshape(np.frombuffer(lsmaskf.read(),dtype=np.uint8),(nlats,nlons))
if lakes:
lsmask =\
np.where(lsmask==2,np.array(0,dtype=np.uint8),lsmask)
lsmaskf.close()
delta = 360./nlons
lsmask_lons = np.linspace(-180+0.5*delta,180-0.5*delta,nlons).astype(np.float32)
lsmask_lats = np.linspace(-90+0.5*delta,90-0.5*delta,nlats).astype(np.float32)
return lsmask_lons, lsmask_lats, lsmask
class _tup(tuple):
# tuple with an added remove method.
# used for objects returned by drawparallels and drawmeridians.
def remove(self):
for item in self:
for x in item:
x.remove()
class _dict(dict):
# override __delitem__ to first call remove method on values.
def __delitem__(self,key):
self[key].remove()
super(_dict, self).__delitem__(key)
def _setlonlab(fmt, lon, labelstyle):
"""Set longitude label string (called by :meth:`Basemap.drawmeridians`)."""
try:
# `fmt` is a function that returns a formatted string.
lonlab = fmt(lon)
except:
# `fmt` is a format string.
degchar = b"\xc2\xb0".decode("utf-8")
if lon > 180:
if mpl.rcParams["text.usetex"]:
if labelstyle == "+/-":
lonlabstr = r"${\/-%s\/^{\circ}}$" % fmt
else:
lonlabstr = r"${%s\/^{\circ}\/W}$" % fmt
else:
if labelstyle == "+/-":
lonlabstr = r"-%s%s" % (fmt, degchar)
else:
lonlabstr = r"%s%sW" % (fmt, degchar)
lonlab = lonlabstr % np.fabs(lon - 360)
elif lon < 180 and lon != 0:
if mpl.rcParams["text.usetex"]:
if labelstyle == "+/-":
lonlabstr = r"${\/+%s\/^{\circ}}$" % fmt
else:
lonlabstr = r"${%s\/^{\circ}\/E}$" % fmt
else:
if labelstyle == "+/-":
lonlabstr = r"+%s%s" % (fmt, degchar)
else:
lonlabstr = r"%s%sE" % (fmt, degchar)
lonlab = lonlabstr % lon
else:
if mpl.rcParams["text.usetex"]:
lonlabstr = r"${%s\/^{\circ}}$" % fmt
else:
lonlabstr = r"%s%s" % (fmt, degchar)
lonlab = lonlabstr % lon
return lonlab
def _setlatlab(fmt, lat, labelstyle):
"""Set latitude label string (called by :meth:`Basemap.drawparallels`)."""
try:
# `fmt` is a function that returns a formatted string.
latlab = fmt(lat)
except:
# `fmt` is a format string.
degchar = b"\xc2\xb0".decode("utf-8")
if lat < 0:
if mpl.rcParams["text.usetex"]:
if labelstyle == "+/-":
latlabstr = r"${\/-%s\/^{\circ}}$" % fmt
else:
latlabstr = r"${%s\/^{\circ}\/S}$" % fmt
else:
if labelstyle == "+/-":
latlabstr = r"-%s%s" % (fmt, degchar)
else:
latlabstr = r"%s%sS" % (fmt, degchar)
latlab = latlabstr % np.fabs(lat)
elif lat > 0:
if mpl.rcParams["text.usetex"]:
if labelstyle == "+/-":
latlabstr = r"${\/+%s\/^{\circ}}$" % fmt
else:
latlabstr = r"${%s\/^{\circ}\/N}$" % fmt
else:
if labelstyle == "+/-":
latlabstr = r"+%s%s" % (fmt, degchar)
else:
latlabstr = r"%s%sN" % (fmt, degchar)
latlab = latlabstr % lat
else:
if mpl.rcParams["text.usetex"]:
latlabstr = r"${%s\/^{\circ}}$" % fmt
else:
latlabstr = r"%s%s" % (fmt, degchar)
latlab = latlabstr % lat
return latlab