Menu
▾
▴
SF.net SVN: matplotlib: [4102] trunk/toolkits/basemap/lib/matplotlib/ toolkits/basemap/basemap.py
|
From: <ef...@us...> - 2007-11-04 06:09:24
|
Revision: 4102
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4102&view=rev
Author: efiring
Date: 2007-11-03 23:09:22 -0700 (Sat, 03 Nov 2007)
Log Message:
-----------
Basic numpification, reformatting, some refactoring
Modified Paths:
--------------
trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/basemap.py
Modified: trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/basemap.py
===================================================================
--- trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-04 01:02:08 UTC (rev 4101)
+++ trunk/toolkits/basemap/lib/matplotlib/toolkits/basemap/basemap.py 2007-11-04 06:09:22 UTC (rev 4102)
@@ -10,94 +10,28 @@
from matplotlib.lines import Line2D
import pyproj, sys, os, math, dbflib
from proj import Proj
-import matplotlib.numerix as NX
-from matplotlib.numerix import ma
+from matplotlib.numerix import npyma as ma
import numpy as npy
from numpy import linspace
from matplotlib.numerix.mlab import squeeze
-from matplotlib.cbook import popd, is_scalar
+from matplotlib.cbook import is_scalar, dedent
+
from shapelib import ShapeFile
+import time
+
+
# basemap data files now installed in lib/matplotlib/toolkits/basemap/data
basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
__version__ = '0.9.7'
-# test to see numerix set to use numpy (if not, raise an error)
-if NX.which[0] != 'numpy':
- raise ImportError("numerix must be set to numpy")
-class Basemap(object):
- """
- Set up a basemap with one of 19 supported map projections
- (cylindrical equidistant, mercator, polyconic, oblique mercator,
- transverse mercator, miller cylindrical, lambert conformal conic,
- azimuthal equidistant, equidistant conic, lambert azimuthal equal area,
- albers equal area conic, gnomonic, orthographic, sinusoidal, mollweide,
- geostationary, robinson, cassini-soldner or stereographic).
- Doesn't actually draw anything, but sets up the map projection class and
- creates the coastline, lake river and political boundary data
- structures in native map projection coordinates.
- Uses a pyrex interface to C-code from proj.4 (http://proj.maptools.org).
-
- Useful instance variables:
-
- projection - map projection ('cyl','merc','mill','lcc','eqdc','aea',
- 'laea', 'nplaea', 'splaea', 'tmerc', 'omerc', 'cass', 'gnom', 'poly',
- 'sinu', 'moll', 'ortho', 'robin', 'aeqd', 'npaeqd', 'spaeqd', 'stere',
- 'geos', 'npstere' or 'spstere')
- (projections prefixed with 'np' or 'sp' are special case polar-centric
- versions of the parent projection)
- aspect - map aspect ratio (size of y dimension / size of x dimension).
- llcrnrlon - longitude of lower left hand corner of the desired map domain.
- llcrnrlon - latitude of lower left hand corner of the desired map domain.
- urcrnrlon - longitude of upper right hand corner of the desired map domain.
- urcrnrlon - latitude of upper right hand corner of the desired map domain.
- llcrnrx,llcrnry,urcrnrx,urcrnry - corners of map domain in projection coordinates.
- rmajor,rminor - equatorial and polar radii of ellipsoid used (in meters).
- resolution - resolution of boundary dataset being used ('c' for crude,
- 'l' for low, etc.). If None, no boundary dataset is associated with the
- Basemap instance.
- srs - a string representing the 'spatial reference system' for the map
- projection as defined by PROJ.4.
-
- Example Usage:
-
->>> from matplotlib.toolkits.basemap import Basemap
->>> from pylab import load, meshgrid, title, arange, show
->>> # read in topo data (on a regular lat/lon grid)
->>> etopo = load('etopo20data.gz')
->>> lons = load('etopo20lons.gz')
->>> lats = load('etopo20lats.gz')
->>> # create Basemap instance for Robinson projection.
->>> m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
->>> # compute native map projection coordinates for lat/lon grid.
->>> x, y = m(*meshgrid(lons,lats))
->>> # make filled contour plot.
->>> cs = m.contourf(x,y,etopo,30,cmap=cm.jet)
->>> m.drawcoastlines() # draw coastlines
->>> m.drawmapboundary() # draw a line around the map region
->>> m.drawparallels(arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels
->>> m.drawmeridians(arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
->>> title('Robinson Projection') # add a title
->>> show()
-
- [this example (simpletest.py) plus many others can be found in the
- examples directory of source distribution. The "OO" version of this
- example (which does not use pylab) is called "simpletest_oo.py".]
-
- Contact: Jeff Whitaker <jef...@no...>
- """
-
- def __init__(self,llcrnrlon=None,llcrnrlat=None,
- urcrnrlon=None,urcrnrlat=None,\
- width=None,height=None,\
- projection='cyl',resolution='c',area_thresh=None,rsphere=6370997.0,\
- lat_ts=None,lat_1=None,lat_2=None,lat_0=None,lon_0=None,\
- lon_1=None,lon_2=None,suppress_ticks=True,\
- satellite_height=None,boundinglat=None,anchor='C',ax=None):
- """
+# 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 = """
create a Basemap instance.
arguments:
@@ -133,12 +67,12 @@
lon_0 - center of desired map domain (in degrees).
lat_0 - center of desired map domain (in degrees).
- For 'sinu', 'moll', 'npstere', 'spstere', 'nplaea', 'splaea', 'nplaea',
- 'splaea', 'npaeqd', 'spaeqd' or 'robin', the values of
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,width and height are ignored (because
- either they are computed internally, or entire globe is always plotted). For the
- cylindrical projections ('cyl','merc' and 'mill'), the default is to use
- llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other
+ For 'sinu', 'moll', 'npstere', 'spstere', 'nplaea', 'splaea', 'nplaea',
+ 'splaea', 'npaeqd', 'spaeqd' or 'robin', the values of
+ llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,width and height are ignored (because
+ either they are computed internally, or entire globe is always plotted). For the
+ cylindrical projections ('cyl','merc' and 'mill'), the default is to use
+ llcrnrlon=-180,llcrnrlat=-90, urcrnrlon=180 and urcrnrlat=90). For all other
projections except 'ortho' and 'geos', either the lat/lon values of the
corners or width and height must be specified by the user.
For 'ortho' and 'geos', the lat/lon values of the corners may be specified,
@@ -191,7 +125,7 @@
The following parameters are map projection parameters which all default to
None. Not all parameters are used by all projections, some are ignored.
- lat_ts - latitude of natural origin (used for mercator, and
+ lat_ts - latitude of natural origin (used for mercator, and
optionally for stereographic projection).
lat_1 - first standard parallel for lambert conformal, albers
equal area projection and equidistant conic projections. Latitude of one
@@ -201,13 +135,13 @@
lat_2 - second standard parallel for lambert conformal, albers
equal area projection and equidistant conic projections. Latitude of one
of the two points on the projection centerline for oblique mercator.
- If lat_2 is not given, it is set to lat_1 for
+ If lat_2 is not given, it is set to lat_1 for
lambert conformal, albers equal area and equidistant conic.
lon_1 - longitude of one of the two points on the projection centerline
for oblique mercator.
lon_2 - longitude of one of the two points on the projection centerline
for oblique mercator.
- lat_0 - central latitude (y-axis origin) - used by all projections,
+ lat_0 - central latitude (y-axis origin) - used by all projections,
lon_0 - central meridian (x-axis origin) - used by all projections,
boundinglat - bounding latitude for pole-centered projections (npstere,spstere,
nplaea,splaea,npaeqd,spaeqd). These projections are square regions centered
@@ -215,10 +149,147 @@
latitude circle boundinglat is tangent to the edge of the map at lon_0.
satellite_height - height of satellite (in m) above equator -
only relevant for geostationary projections ('geos').
-
+
"""
+_unsupported_projection = """
+ unsupported projection, use 'cyl' - cylindrical equidistant, 'merc' -
+ mercator, 'lcc' - lambert conformal conic, 'stere' - stereographic,
+ 'npstere' - stereographic, special case centered on north pole.
+ 'spstere' - stereographic, special case centered on south pole,
+ 'aea' - albers equal area conic, 'tmerc' - transverse mercator,
+ 'aeqd' - azimuthal equidistant, 'mill' - miller cylindrical,
+ 'npaeqd' - azimuthal equidistant, special case centered on north pole,
+ 'spaeqd' - azimuthal equidistant, special case centered on south pole,
+ 'eqdc' - equidistant conic, 'laea' - lambert azimuthal equal area,
+ 'nplaea' - lambert azimuthal, special case centered on north pole,
+ 'splaea' - lambert azimuthal, special case centered on south pole,
+ 'cass' - cassini-soldner (transverse cylindrical equidistant),
+ 'poly' - polyconic, 'omerc' - oblique mercator, 'ortho' - orthographic,
+ 'geos' - geostationary, 'sinu' - sinusoidal, 'moll' - mollweide,
+ 'robin' - robinson, or 'gnom' - gnomonic. You tried '%s'
+ """
+
+# This allows substitution of longer names into error messages.
+projnames = {'cyl' : 'Cylindrical Equidistant',
+ 'merc' : 'Mercator',
+ 'tmerc' : 'Transverse Mercator',
+ 'omerc' : 'Oblique Mercator',
+ 'mill' : 'Miller Cylindrical',
+ 'llc' : 'Lambert Conformal',
+ 'laea' : 'Lambert Azimuthal Equal Area',
+ 'nplaea' : 'North-Polar Lambert Azimuthal',
+ 'splaea' : 'South-Polar Lambert Azimuthal',
+ 'eqdc' : 'Equidistant Conic',
+ 'eaqd' : 'Azimuthal Equidistant',
+ 'npaeqd' : 'North-Polar Azimuthal Equidistant',
+ 'spaeqd' : 'South-Polar Azimuthal Equidistant',
+ 'aea' : 'Albers Equal Area',
+ 'stere' : 'Stereographic',
+ 'npstere' : 'Nouth-Polar Stereographic',
+ 'spstere' : 'South-Polar Stereographic',
+ 'cass' : 'Cassini-Soldner',
+ 'poly' : 'Polyconic',
+ 'ortho' : 'Orthographic',
+ 'geos' : 'Geostationary',
+ 'sinu' : 'Sinusoidal',
+ 'moll' : 'Mollweide',
+ 'robin' : 'Robinson',
+ 'gnom' : 'Gnomonic',
+ }
+
+def _validated_ll(param, name, minval, maxval):
+ param = float(param)
+ if param > maxval or param < minval:
+ raise ValueError('%s must be between %f and %f degrees' %
+ (name, minval, maxval))
+ return param
+
+def _insert_validated(d, param, name, minval, maxval):
+ if param is not None:
+ d[name] = _validated_ll(param, name, minval, maxval)
+
+
+class Basemap(object):
+ """
+ Set up a basemap with one of 19 supported map projections
+ (cylindrical equidistant, mercator, polyconic, oblique mercator,
+ transverse mercator, miller cylindrical, lambert conformal conic,
+ azimuthal equidistant, equidistant conic, lambert azimuthal equal area,
+ albers equal area conic, gnomonic, orthographic, sinusoidal, mollweide,
+ geostationary, robinson, cassini-soldner or stereographic).
+ Doesn't actually draw anything, but sets up the map projection class and
+ creates the coastline, lake river and political boundary data
+ structures in native map projection coordinates.
+ Uses a pyrex interface to C-code from proj.4 (http://proj.maptools.org).
+
+ Useful instance variables:
+
+ projection - map projection ('cyl','merc','mill','lcc','eqdc','aea',
+ 'laea', 'nplaea', 'splaea', 'tmerc', 'omerc', 'cass', 'gnom', 'poly',
+ 'sinu', 'moll', 'ortho', 'robin', 'aeqd', 'npaeqd', 'spaeqd', 'stere',
+ 'geos', 'npstere' or 'spstere')
+ (projections prefixed with 'np' or 'sp' are special case polar-centric
+ versions of the parent projection)
+ aspect - map aspect ratio (size of y dimension / size of x dimension).
+ llcrnrlon - longitude of lower left hand corner of the desired map domain.
+ llcrnrlon - latitude of lower left hand corner of the desired map domain.
+ urcrnrlon - longitude of upper right hand corner of the desired map domain.
+ urcrnrlon - latitude of upper right hand corner of the desired map domain.
+ llcrnrx,llcrnry,urcrnrx,urcrnry - corners of map domain in projection coordinates.
+ rmajor,rminor - equatorial and polar radii of ellipsoid used (in meters).
+ resolution - resolution of boundary dataset being used ('c' for crude,
+ 'l' for low, etc.). If None, no boundary dataset is associated with the
+ Basemap instance.
+ srs - a string representing the 'spatial reference system' for the map
+ projection as defined by PROJ.4.
+
+ Example Usage:
+
+ >>> from matplotlib.toolkits.basemap import Basemap
+ >>> from pylab import load, meshgrid, title, arange, show
+ >>> # read in topo data (on a regular lat/lon grid)
+ >>> etopo = load('etopo20data.gz')
+ >>> lons = load('etopo20lons.gz')
+ >>> lats = load('etopo20lats.gz')
+ >>> # create Basemap instance for Robinson projection.
+ >>> m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
+ >>> # compute native map projection coordinates for lat/lon grid.
+ >>> x, y = m(*meshgrid(lons,lats))
+ >>> # make filled contour plot.
+ >>> cs = m.contourf(x,y,etopo,30,cmap=cm.jet)
+ >>> m.drawcoastlines() # draw coastlines
+ >>> m.drawmapboundary() # draw a line around the map region
+ >>> m.drawparallels(arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels
+ >>> m.drawmeridians(arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
+ >>> title('Robinson Projection') # add a title
+ >>> show()
+
+ [this example (simpletest.py) plus many others can be found in the
+ examples directory of source distribution. The "OO" version of this
+ example (which does not use pylab) is called "simpletest_oo.py".]
+
+ Contact: Jeff Whitaker <jef...@no...>
+ """
+
+
+ def __init__(self, llcrnrlon=None, llcrnrlat=None,
+ urcrnrlon=None, urcrnrlat=None,
+ width=None, height=None,
+ projection='cyl', resolution='c',
+ area_thresh=None, rsphere=6370997.0,
+ lat_ts=None,
+ lat_1=None, lat_2=None,
+ lat_0=None, lon_0=None,
+ lon_1=None, lon_2=None,
+ suppress_ticks=True,
+ satellite_height=None,
+ boundinglat=None,
+ anchor='C',
+ ax=None):
+ # docstring is added after definition
+ #print "starting: ", time.clock()
# where to put plot in figure (default is 'C' or center)
self.anchor = anchor
# map projection.
@@ -228,241 +299,103 @@
projparams = {}
projparams['proj'] = projection
try:
- if rsphere[0] > rsphere[1]:
- projparams['a'] = rsphere[0]
- projparams['b'] = rsphere[1]
- else:
- projparams['a'] = rsphere[1]
- projparams['b'] = rsphere[0]
- except:
+ projparams['a'] = max(rsphere)
+ projparams['b'] = min(rsphere)
+ except TypeError:
projparams['a'] = rsphere
projparams['b'] = rsphere
# set units to meters.
projparams['units']='m'
# check for sane values of lon_0, lat_0, lat_ts, lat_1, lat_2
- if lat_0 is not None:
- if lat_0 > 90. or lat_0 < -90.:
- raise ValueError, 'lat_0 must be between -90 and +90 degrees'
- else:
- projparams['lat_0'] = lat_0
- if lon_0 is not None:
- if lon_0 < -720. or lon_0 > 720.:
- raise ValueError, 'lon_0 must be between -720 and +720 degrees'
- else:
- projparams['lon_0'] = lon_0
- if lon_1 is not None:
- if lon_1 < -720. or lon_1 > 720.:
- raise ValueError, 'lon_1 must be between -720 and +720 degrees'
- else:
- projparams['lon_1'] = lon_1
- if lon_2 is not None:
- if lon_2 < -720. or lon_2 > 720.:
- raise ValueError, 'lon_2 must be between -720 and +720 degrees'
- else:
- projparams['lon_2'] = lon_2
- if lat_1 is not None:
- if lat_1 > 90. or lat_1 < -90.:
- raise ValueError, 'lat_1 must be between -90 and +90 degrees'
- else:
- projparams['lat_1'] = lat_1
- if lat_2 is not None:
- if lat_2 > 90. or lat_2 < -90.:
- raise ValueError, 'lat_2 must be between -90 and +90 degrees'
- else:
- projparams['lat_2'] = lat_2
- if lat_ts is not None:
- if lat_ts > 90. or lat_ts < -90.:
- raise ValueError, 'lat_ts must be between -90 and +90 degrees'
- else:
- projparams['lat_ts'] = lat_ts
+ _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 satellite_height is not None:
projparams['h'] = satellite_height
- if None not in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
- # make sure lat/lon limits are converted to floats.
- self.llcrnrlon = float(llcrnrlon)
- self.llcrnrlat = float(llcrnrlat)
- self.urcrnrlon = float(urcrnrlon)
- self.urcrnrlat = float(urcrnrlat)
- # check values of urcrnrlon,urcrnrlat and llcrnrlon,llcrnrlat
- if self.urcrnrlat > 90.0 or self.llcrnrlat > 90.0:
- raise ValueError, 'urcrnrlat and llcrnrlat must be less than 90'
- if self.urcrnrlat < -90.0 or self.llcrnrlat < -90.0:
- raise ValueError, 'urcrnrlat and llcrnrlat must be greater than -90'
- if self.urcrnrlon > 720. or self.llcrnrlon > 720.:
- raise ValueError, 'urcrnrlon and llcrnrlon must be less than 720'
- if self.urcrnrlon < -360. or self.llcrnrlon < -360.:
- raise ValueError, 'urcrnrlon and llcrnrlon must be greater than -360'
- # for each of the supported projections, compute lat/lon of domain corners
+ 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 == 'lcc':
+ if projection in ['lcc', 'eqdc', 'aea']:
# if lat_0 is given, but not lat_1,
# set lat_1=lat_0
if lat_1 is None and lat_0 is not None:
lat_1 = lat_0
projparams['lat_1'] = lat_1
if lat_1 is None or lon_0 is None:
- raise ValueError, 'must specify lat_1 or lat_0 and lon_0 for Lambert Conformal basemap (lat_2 is optional)'
+ raise ValueError('must specify lat_1 or lat_0 and lon_0 for %(projection)s basemap (lat_2 is optional)' % projnames)
if lat_2 is None:
projparams['lat_2'] = lat_1
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
if width is None or height is None:
raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
-
- elif projection == 'eqdc':
- # if lat_0 is given, but not lat_1,
- # set lat_1=lat_0
- if lat_1 is None and lat_0 is not None:
- lat_1 = lat_0
- projparams['lat_1'] = lat_1
- if lat_1 is None or lon_0 is None:
- raise ValueError, 'must specify lat_1 or lat_0 and lon_0 for Equidistant Conic basemap (lat_2 is optional)'
- if lat_2 is None:
- projparams['lat_2'] = lat_1
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
- if width is None or height is None:
- raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
- elif projection == 'aea':
- # if lat_0 is given, but not lat_1,
- # set lat_1=lat_0
- if lat_1 is None and lat_0 is not None:
- lat_1 = lat_0
- projparams['lat_1'] = lat_1
- if lat_1 is None or lon_0 is None:
- raise ValueError, 'must specify lat_1 or lat_0 and lon_0 for Albers Equal Area basemap (lat_2 is optional)'
- if lat_2 is None:
- projparams['lat_2'] = lat_1
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
- if width is None or height is None:
- raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
+ 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
+
+ # skipping over the following for now; it can be beautified and
+ # consolidated later
elif projection == 'stere':
if lat_0 is None or lon_0 is None:
raise ValueError, 'must specify lat_0 and lon_0 for Stereographic basemap (lat_ts is optional)'
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
if width is None or height is None:
raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
- elif projection == 'spstere':
+ if 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 boundinglat is None or lon_0 is None:
- raise ValueError, 'must specify boundinglat and lon_0 for South-Polar Stereographic basemap'
- projparams['lat_ts'] = -90.
- projparams['lat_0'] = -90.
- projparams['proj'] = 'stere'
- self.llcrnrlon = lon_0+45.
- self.urcrnrlon = lon_0-135.
+ raise ValueError('must specify boundinglat and lon_0 for %(projection) basemap' % projnames)
+ 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:
print 'warning: width and height keywords ignored for %s projection' % self.projection
- elif projection == 'npstere':
- if boundinglat is None or lon_0 is None:
- raise ValueError, 'must specify boundinglat and lon_0 for North-Polar Stereographic basemap'
- projparams['lat_ts'] = 90.
- projparams['lat_0'] = 90.
- projparams['proj'] = 'stere'
- self.llcrnrlon = lon_0-45.
- self.urcrnrlon = lon_0+135.
- proj = pyproj.Proj(projparams)
- x,y = proj(lon_0,boundinglat)
- lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
- self.urcrnrlat = self.llcrnrlat
- if width is not None or height is not None:
- print 'warning: width and height keywords ignored for %s projection' % self.projection
- elif projection == 'splaea':
- if boundinglat is None or lon_0 is None:
- raise ValueError, 'must specify boundinglat and lon_0 for South-Polar Lambert Azimuthal basemap'
- projparams['lat_0'] = -90.
- projparams['proj'] = 'laea'
- self.llcrnrlon = lon_0+45.
- self.urcrnrlon = lon_0-135.
- proj = pyproj.Proj(projparams)
- x,y = proj(lon_0,boundinglat)
- lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
- self.urcrnrlat = self.llcrnrlat
- if width is not None or height is not None:
- print 'warning: width and height keywords ignored for %s projection' % self.projection
- elif projection == 'nplaea':
- if boundinglat is None or lon_0 is None:
- raise ValueError, 'must specify boundinglat and lon_0 for North-Polar Lambert Azimuthal basemap'
- projparams['lat_0'] = 90.
- projparams['proj'] = 'laea'
- self.llcrnrlon = lon_0-45.
- self.urcrnrlon = lon_0+135.
- proj = pyproj.Proj(projparams)
- x,y = proj(lon_0,boundinglat)
- lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
- self.urcrnrlat = self.llcrnrlat
- if width is not None or height is not None:
- print 'warning: width and height keywords ignored for %s projection' % self.projection
- elif projection == 'spaeqd':
- if boundinglat is None or lon_0 is None:
- raise ValueError, 'must specify boundinglat and lon_0 for South-Polar Azimuthal Equidistant basemap'
- projparams['lat_0'] = -90.
- projparams['proj'] = 'aeqd'
- self.llcrnrlon = lon_0+45.
- self.urcrnrlon = lon_0-135.
- proj = pyproj.Proj(projparams)
- x,y = proj(lon_0,boundinglat)
- lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
- self.urcrnrlat = self.llcrnrlat
- if width is not None or height is not None:
- print 'warning: width and height keywords ignored for %s projection' % self.projection
- elif projection == 'npaeqd':
- if boundinglat is None or lon_0 is None:
- raise ValueError, 'must specify boundinglat and lon_0 for North-Polar Azimuthal Equidistant basemap'
- projparams['lat_0'] = 90.
- projparams['proj'] = 'aeqd'
- self.llcrnrlon = lon_0-45.
- self.urcrnrlon = lon_0+135.
- proj = pyproj.Proj(projparams)
- x,y = proj(lon_0,boundinglat)
- lon,self.llcrnrlat = proj(math.sqrt(2.)*y,0.,inverse=True)
- self.urcrnrlat = self.llcrnrlat
- if width is not None or height is not None:
- print 'warning: width and height keywords ignored for %s projection' % self.projection
elif projection == 'laea':
if lat_0 is None or lon_0 is None:
raise ValueError, 'must specify lat_0 and lon_0 for Lambert Azimuthal basemap'
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
if width is None or height is None:
raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
+ if lon_0 is None or lat_0 is None:
+ raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
+ llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
+ self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
+ self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'merc':
if lat_ts is None:
raise ValueError, 'must specify lat_ts for Mercator basemap'
# clip plot region to be within -89.99S to 89.99N
# (mercator is singular at poles)
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
@@ -478,16 +411,15 @@
elif projection in ['tmerc','gnom','cass','poly'] :
if lat_0 is None or lon_0 is None:
raise ValueError, 'must specify lat_0 and lon_0 for Transverse Mercator, Gnomonic, Cassini-Soldnerr Polyconic basemap'
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
if width is None or height is None:
raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
-
+ 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 projparams['a'] != projparams['b']:
raise ValueError, 'orthographic projection only works for perfect spheres - not ellipsoids'
@@ -495,7 +427,7 @@
raise ValueError, 'must specify lat_0 and lon_0 for Orthographic basemap'
if width is not None or height is not None:
print 'warning: width and height keywords ignored for %s projection' % self.projection
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
@@ -510,7 +442,7 @@
raise ValueError, 'must specify lon_0 and satellite_height for Geostationary basemap'
if width is not None or height is not None:
print 'warning: width and height keywords ignored for %s projection' % self.projection
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
@@ -538,29 +470,27 @@
projparams['lon_1'] = lon_1
projparams['lat_2'] = lat_2
projparams['lon_2'] = lon_2
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
if width is None or height is None:
raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
+ if lon_0 is None or lat_0 is None:
+ raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
+ llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
+ self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
+ self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'aeqd':
if lat_0 is None or lon_0 is None:
raise ValueError, 'must specify lat_0 and lon_0 for Azimuthal Equidistant basemap'
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
if width is None or height is None:
raise ValueError, 'must either specify lat/lon values of corners (llcrnrlon,llcrnrlat,ucrnrlon,urcrnrlat) in degrees or width and height in meters'
- else:
- if lon_0 is None or lat_0 is None:
- raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
- llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
- self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
- self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
+ if lon_0 is None or lat_0 is None:
+ raise ValueError, 'must specify lon_0 and lat_0 when using width, height to specify projection region'
+ llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat = _choosecorners(width,height,**projparams)
+ self.llcrnrlon = llcrnrlon; self.llcrnrlat = llcrnrlat
+ self.urcrnrlon = urcrnrlon; self.urcrnrlat = urcrnrlat
elif projection == 'mill':
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
@@ -570,7 +500,7 @@
if width is not None or height is not None:
print 'warning: width and height keywords ignored for %s projection' % self.projection
elif projection == 'cyl':
- if None in [llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat]:
+ if not using_corners:
llcrnrlon = -180.
llcrnrlat = -90.
urcrnrlon = 180
@@ -580,26 +510,12 @@
if width is not None or height is not None:
print 'warning: width and height keywords ignored for %s projection' % self.projection
else:
- raise ValueError, """
- unsupported projection, use 'cyl' - cylindrical equidistant, 'merc' -
- mercator, 'lcc' - lambert conformal conic, 'stere' - stereographic,
- 'npstere' - stereographic, special case centered on north pole.
- 'spstere' - stereographic, special case centered on south pole,
- 'aea' - albers equal area conic, 'tmerc' - transverse mercator,
- 'aeqd' - azimuthal equidistant, 'mill' - miller cylindrical,
- 'npaeqd' - azimuthal equidistant, special case centered on north pole,
- 'spaeqd' - azimuthal equidistant, special case centered on south pole,
- 'eqdc' - equidistant conic, 'laea' - lambert azimuthal equal area,
- 'nplaea' - lambert azimuthal, special case centered on north pole,
- 'splaea' - lambert azimuthal, special case centered on south pole,
- 'cass' - cassini-soldner (transverse cylindrical equidistant),
- 'poly' - polyconic, 'omerc' - oblique mercator, 'ortho' - orthographic,
- 'geos' - geostationary, 'sinu' - sinusoidal, 'moll' - mollweide,
- 'robin' - robinson, or 'gnom' - gnomonic. You tried '%s'""" % projection
+ 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
@@ -665,6 +581,8 @@
# if no boundary data needed, we are done.
if self.resolution is None:
return
+
+ ################## starting boundary processing ###################
if area_thresh is None:
if resolution == 'c':
area_thresh = 10000.
@@ -681,8 +599,10 @@
msg = """
Unable to open boundary dataset file. Only the 'crude', 'low'
and 'intermediate' resolution datasets are installed by default. If you
-are requesting a 'high' resolution dataset, you need to download
+are requesting a 'high' resolution dataset, you need to download
and install those files manually (see the basemap README for details)."""
+
+ #print "reading data file:", time.clock()
try:
bdatfile = open(os.path.join(basemap_datadir,'gshhs_'+resolution+'.txt'))
except:
@@ -705,6 +625,7 @@
coastlats.append(lat)
coastsegtype.append(typ)
coastsegind.append(len(coastlons))
+ #print "read coasts", time.clock()
# read in country boundary data.
cntrylons = []; cntrylats = []; cntrysegind = []
@@ -770,7 +691,7 @@
# so valid longitudes can range from -360 to 720.
# This means a lot of redundant processing is done when
# creating the class instance, but it a lot easier to figure
- # out what to do when the projection domain straddles the
+ # out what to do when the projection domain straddles the
# Greenwich meridian.
coastlons2 = [lon+360. for lon in coastlons]
cntrylons2 = [lon+360. for lon in cntrylons]
@@ -780,35 +701,50 @@
cntrylons3 = [lon-360. for lon in cntrylons]
statelons3 = [lon-360. for lon in statelons]
riverlons3 = [lon-360. for lon in riverlons]
-
+ #print "starting to make coast polygons", time.clock()
# transform coastline polygons to native map coordinates.
- xc,yc = proj(NX.array(coastlons),NX.array(coastlats))
+ xc,yc = proj(npy.array(coastlons),npy.array(coastlats))
xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(coastlons2),NX.array(coastlats))
- xc3,yc3 = proj(NX.array(coastlons3),NX.array(coastlats))
+ xc2,yc2 = proj(npy.array(coastlons2),npy.array(coastlats))
+ xc3,yc3 = proj(npy.array(coastlons3),npy.array(coastlats))
xc2 = xc2.tolist(); yc2 = yc2.tolist()
xc3 = xc3.tolist(); yc3 = yc3.tolist()
# set up segments in form needed for LineCollection,
# ignoring 'inf' values that are off the map.
- segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:])]
- segmentsll = [zip(coastlons[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:])]
+ segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in
+ zip(coastsegind[:-1],coastsegind[1:])]
+ segmentsll = [zip(coastlons[i0:i1],coastlats[i0:i1]) for i0,i1 in
+ zip(coastsegind[:-1],coastsegind[1:])]
segtypes = [i for i in coastsegtype[:-1]]
- segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segmentsll2 = [zip(coastlons2[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
- segtypes2 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1]) if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
+ segments2 = [zip(xc2[i0:i1],yc2[i0:i1]) for i0,i1 in
+ zip(coastsegind[:-1],coastsegind[1:])
+ if max(xc2[i0:i1]) < 1.e20
+ and max(yc2[i0:i1]) < 1.e20]
+ segmentsll2 = [zip(coastlons2[i0:i1],coastlats[i0:i1]) for i0,i1 in
+ zip(coastsegind[:-1],coastsegind[1:])
+ if max(xc2[i0:i1]) < 1.e20
+ and max(yc2[i0:i1]) < 1.e20]
+ segtypes2 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1])
+ if max(xc2[i0:i1]) < 1.e20 and max(yc2[i0:i1]) < 1.e20]
+
segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
segmentsll3 = [zip(coastlons3[i0:i1],coastlats[i0:i1]) for i0,i1 in zip(coastsegind[:-1],coastsegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
segtypes3 = [i for i0,i1,i in zip(coastsegind[:-1],coastsegind[1:],coastsegtype[:-1]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
- self.coastsegs = segments+segments2+segments3
- self.coastsegsll = segmentsll+segmentsll2+segmentsll3
- self.coastsegtypes = segtypes+segtypes2+segtypes3
+ self.coastsegs = segments +segments2+segments3
+ self.coastsegsll = segmentsll +segmentsll2+segmentsll3
+ self.coastsegtypes = segtypes +segtypes2+segtypes3
+ #print len(coastsegind)
+ #print len(segments), len(segments2), len(segments3)
+ #print len(self.coastsegs), len(self.coastsegsll), len(self.coastsegtypes)
+ #print "made segments", time.clock()
+
# same as above for country segments.
- xc,yc = proj(NX.array(cntrylons),NX.array(cntrylats))
+ xc,yc = proj(npy.array(cntrylons),npy.array(cntrylats))
xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(cntrylons2),NX.array(cntrylats))
- xc3,yc3 = proj(NX.array(cntrylons3),NX.array(cntrylats))
+ xc2,yc2 = proj(npy.array(cntrylons2),npy.array(cntrylats))
+ xc3,yc3 = proj(npy.array(cntrylons3),npy.array(cntrylats))
xc2 = xc2.tolist(); yc2 = yc2.tolist()
xc3 = xc3.tolist(); yc3 = yc3.tolist()
segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(cntrysegind[:-1],cntrysegind[1:])]
@@ -817,10 +753,10 @@
self.cntrysegs = segments+segments2+segments3
# same as above for state segments.
- xc,yc = proj(NX.array(statelons),NX.array(statelats))
+ xc,yc = proj(npy.array(statelons),npy.array(statelats))
xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(statelons2),NX.array(statelats))
- xc3,yc3 = proj(NX.array(statelons3),NX.array(statelats))
+ xc2,yc2 = proj(npy.array(statelons2),npy.array(statelats))
+ xc3,yc3 = proj(npy.array(statelons3),npy.array(statelats))
xc2 = xc2.tolist(); yc2 = yc2.tolist()
xc3 = xc3.tolist(); yc3 = yc3.tolist()
segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(statesegind[:-1],statesegind[1:])]
@@ -829,10 +765,10 @@
self.statesegs = segments+segments2+segments3
# same as above for river segments.
- xc,yc = proj(NX.array(riverlons),NX.array(riverlats))
+ xc,yc = proj(npy.array(riverlons),npy.array(riverlats))
xc = xc.tolist(); yc = yc.tolist()
- xc2,yc2 = proj(NX.array(riverlons2),NX.array(riverlats))
- xc3,yc3 = proj(NX.array(riverlons3),NX.array(riverlats))
+ xc2,yc2 = proj(npy.array(riverlons2),npy.array(riverlats))
+ xc3,yc3 = proj(npy.array(riverlons3),npy.array(riverlats))
xc2 = xc2.tolist(); yc2 = yc2.tolist()
xc3 = xc3.tolist(); yc3 = yc3.tolist()
segments = [zip(xc[i0:i1],yc[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:])]
@@ -840,6 +776,7 @@
segments3 = [zip(xc3[i0:i1],yc3[i0:i1]) for i0,i1 in zip(riversegind[:-1],riversegind[1:]) if max(xc3[i0:i1]) < 1.e20 and max(yc3[i0:i1]) < 1.e20]
self.riversegs = segments+segments2+segments3
+ #print "Making final set of polygons", time.clock()
# store coast polygons for filling.
self.coastpolygons = []
coastpolygonsll = []
@@ -939,6 +876,7 @@
self.coastpolygons = polygons
coastpolygonsll = polygonsll
self.coastpolygontypes = polygontypes
+ #print "made coastpolygons", time.clock()
states = []; rivers = []; countries = []
for seg in self.cntrysegs:
if self._insidemap_seg(seg):
@@ -953,12 +891,13 @@
self.riversegs = rivers
self.cntryegs = countries
+ #print "starting projection limb checks", time.clock()
# split up segments that go outside projection limb
coastsegs = []
coastsegtypes = []
for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
+ xx = npy.array([x for x,y in seg],npy.float32)
+ yy = npy.array([y for x,y in seg],npy.float32)
i1,i2 = self._splitseg(xx,yy)
if i1 and i2:
for i,j in zip(i1,i2):
@@ -970,10 +909,11 @@
coastsegs.append(segtype)
self.coastsegs = coastsegs
self.coastsegtypes = coastsegtypes
+ #print "finished p l checks", time.clock()
states = []
for seg in self.statesegs:
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
+ xx = npy.array([x for x,y in seg],npy.float32)
+ yy = npy.array([y for x,y in seg],npy.float32)
i1,i2 = self._splitseg(xx,yy)
if i1 and i2:
for i,j in zip(i1,i2):
@@ -984,8 +924,8 @@
self.statesegs = states
countries = []
for seg in self.cntrysegs:
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
+ xx = npy.array([x for x,y in seg],npy.float32)
+ yy = npy.array([y for x,y in seg],npy.float32)
i1,i2 = self._splitseg(xx,yy)
if i1 and i2:
for i,j in zip(i1,i2):
@@ -996,8 +936,8 @@
self.cntrysegs = countries
rivers = []
for seg in self.riversegs:
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
+ xx = npy.array([x for x,y in seg],npy.float32)
+ yy = npy.array([y for x,y in seg],npy.float32)
i1,i2 = self._splitseg(xx,yy)
if i1 and i2:
for i,j in zip(i1,i2):
@@ -1007,18 +947,19 @@
rivers.append(seg)
self.riversegs = rivers
+ "Starting remaining coast processing", time.clock()
# split coastline segments that jump across entire plot.
coastsegs = []
coastsegtypes = []
for seg,segtype in zip(self.coastsegs,self.coastsegtypes):
- xx = NX.array([x for x,y in seg],NX.Float32)
- yy = NX.array([y for x,y in seg],NX.Float32)
+ xx = npy.array([x for x,y in seg],npy.float32)
+ yy = npy.array([y for x,y in seg],npy.float32)
xd = (xx[1:]-xx[0:-1])**2
yd = (yy[1:]-yy[0:-1])**2
- dist = NX.sqrt(xd+yd)
+ dist = npy.sqrt(xd+yd)
split = dist > 5000000.
- if NX.sum(split) and self.projection not in ['merc','cyl','mill']:
- ind = (NX.compress(split,squeeze(split*NX.indices(xd.shape)))+1).tolist()
+ if npy.sum(split) and self.projection not in ['merc','cyl','mill']:
+ ind = (npy.compress(split,squeeze(split*npy.indices(xd.shape)))+1).tolist()
iprev = 0
ind.append(len(xd))
for i in ind:
@@ -1055,11 +996,11 @@
y = poly[1]
lons = polyll[0]
lats = polyll[1]
- mask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
+ mask = npy.logical_or(npy.greater(x,1.e20),npy.greater(y,1.e20))
# replace values in polygons that are over the horizon.
xsave = False
ysave = False
- if NX.sum(mask):
+ if npy.sum(mask):
i1,i2 = self._splitseg(x,y,mask=mask)
# loop over segments of polygon that are outside projection limb.
for i,j in zip(i1,i2):
@@ -1147,12 +1088,12 @@
lons.reverse()
lats.reverse()
xx,yy = self(lons,lats)
- xx = NX.array(xx); yy = NX.array(yy)
- xdist = NX.fabs(xx[1:]-xx[0:-1])
+ xx = npy.array(xx); yy = npy.array(yy)
+ xdist = npy.fabs(xx[1:]-xx[0:-1])
if max(xdist) > 1000000:
- nmin = NX.argmax(xdist)+1
- xnew = NX.zeros(len(xx),NX.Float64)
- ynew = NX.zeros(len(xx),NX.Float64)
+ nmin = npy.argmax(xdist)+1
+ xnew = npy.zeros(len(xx),npy.float64)
+ ynew = npy.zeros(len(xx),npy.float64)
lonsnew = len(xx)*[0]
latsnew = len(xx)*[0]
xnew[0:len(xx)-nmin] = xx[nmin:]
@@ -1169,12 +1110,12 @@
x.reverse()
y.reverse()
# close polygon (add lines along left and right edges down to S pole)
- for phi in NX.arange(-89.999,lats[0],0.1):
+ for phi in npy.arange(-89.999,lats[0],0.1):
xx,yy = self(lon_0-179.99,phi)
xn.append(xx); yn.append(yy)
xn = xn+x
yn = yn+y
- for phi in NX.arange(lats[-1],-89.999,-0.1):
+ for phi in npy.arange(lats[-1],-89.999,-0.1):
xx,yy = self(lon_0+179.99,phi)
xn.append(xx); yn.append(yy)
# move points outside map to edge of map
@@ -1190,11 +1131,17 @@
xn.append(x); yn.append(y)
coastpolygons.append((xn,yn))
self.coastpolygons = coastpolygons
+ #print "finished init", time.clock()
+ __init__.__doc__ = _Basemap_init_doc
+ #### End of the world's longest __init__
+
+
+
def _splitseg(self,xx,yy,mask=None):
"""split segment up around missing values (outside projection limb)"""
if mask is None:
- mask = (NX.logical_or(NX.greater_equal(xx,1.e20),NX.greater_equal(yy,1.e20))).tolist()
+ mask = (npy.logical_or(npy.greater_equal(xx,1.e20),npy.greater_equal(yy,1.e20))).tolist()
i1=[]; i2=[]
mprev = 1
for i,m in enumerate(mask):
@@ -1239,37 +1186,38 @@
def __call__(self,x,y,inverse=False):
"""
- Calling a Basemap class instance with the arguments lon, lat will
- convert lon/lat (in degrees) to x/y native map projection
- coordinates (in meters). If optional keyword 'inverse' is
- True (default is False), the inverse transformation from x/y
- to lon/lat is performed.
+ Calling a Basemap class instance with the arguments lon, lat will
+ convert lon/lat (in degrees) to x/y native map projection
+ coordinates (in meters). If optional keyword 'inverse' is
+ True (default is False), the inverse transformation from x/y
+ to lon/lat is performed.
- For cylindrical equidistant projection ('cyl'), this
- does nothing (i.e. x,y == lon,lat).
+ For cylindrical equidistant projection ('cyl'), this
+ does nothing (i.e. x,y == lon,lat).
- For non-cylindrical projections, the inverse transformation
- always returns longitudes between -180 and 180 degrees. For
- cylindrical projections (self.projection == 'cyl','mill' or 'merc')
- the inverse transformation will return longitudes between
- self.llcrnrlon and self.llcrnrlat.
+ For non-cylindrical projections, the inverse transformation
+ always returns longitudes between -180 and 180 degrees. For
+ cylindrical projections (self.projection == 'cyl','mill' or 'merc')
+ the inverse transformation will return longitudes between
+ self.llcrnrlon and self.llcrnrlat.
- input arguments lon, lat can be either scalar floats or N arrays.
+ input arguments lon, lat can be either scalar floats or N arrays.
"""
return self.projtran(x,y,inverse=inverse)
def makegrid(self,nx,ny,returnxy=False):
"""
- return arrays of shape (ny,nx) containing lon,lat coordinates of
- an equally spaced native projection grid.
- if returnxy = True, the x,y values of the grid are returned also.
+ return arrays of shape (ny,nx) containing lon,lat coordinates of
+ an equally spaced native projection grid.
+ if returnxy = True, the x,y values of the grid are returned also.
"""
return self.projtran.makegrid(nx,ny,returnxy=returnxy)
def drawmapboundary(self,color='k',linewidth=1.0,ax=None):
"""
- draw boundary around map projection region. If ax=None (default),
- default axis instance is used, otherwise specified axis instance is used.
+ draw boundary around map projection region. If ax=None (default),
+ default axis instance is used, otherwise specified axis
+ instance is used.
"""
# get current axes instance (if none specified).
if ax is None and self.ax is None:
@@ -1301,21 +1249,21 @@
ellps.set_clip_on(False)
elif self.projection in ['moll','robin','sinu']: # elliptical region.
# left side
- lats = NX.arange(-89.9,89.9+dtheta,dtheta).tolist()
+ lats = npy.arange(-89.9,89.9+dtheta,dtheta).tolist()
lons = len(lats)*[self.projparams['lon_0']-179.9]
x,y = self(lons,lats)
# top.
- lons = NX.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
+ lons = npy.arange(self.projparams['lon_0']-179.9,self.projparams['lon_0']+179+dtheta,dtheta).tolist()
lats = len(lons)*[89.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
# right side
- lats = NX.arange(89.9,-89.9-dtheta,-dtheta).tolist()
+ lats = npy.arange(89.9,-89.9-dtheta,-dtheta).tolist()
lons = len(lats)*[self.projparams['lon_0']+179.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
# bottom.
- lons = NX.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
+ lons = npy.arange(self.projparams['lon_0']+179.9,self.projparams['lon_0']-180-dtheta,-dtheta).tolist()
lats = len(lons)*[-89.9]
xx,yy = self(lons,lats)
x = x+xx; y = y+yy
@@ -1334,18 +1282,18 @@
def fillcontinents(self,color='0.8',ax=None,zorder=None):
"""
- Fill continents.
+ Fill continents.
- color - color to fill continents (default gray).
- 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).
+ color - color to fill continents (default gray).
+ 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).
- After filling continents, lakes are re-filled with axis background color.
+ After filling continents, lakes are re-filled with axis background color.
"""
if self.resolution is None:
- raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
+ raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
# get current axes instance (if none specified).
if ax is None and self.ax is None:
try:
@@ -1359,18 +1307,18 @@
axisbgc = ax.get_axis_bgcolor()
np = 0
for x,y in self.coastpolygons:
- xa = NX.array(x,NX.Flo...
[truncated message content] |
