diff --git a/.requirements-2.6.txt b/.requirements-2.6.txt new file mode 100644 index 000000000..9a23970d3 --- /dev/null +++ b/.requirements-2.6.txt @@ -0,0 +1 @@ +unittest2 diff --git a/.travis.yml b/.travis.yml index 4e5183250..08f68ef45 100644 --- a/.travis.yml +++ b/.travis.yml @@ -32,6 +32,10 @@ install: pip install --upgrade setuptools pip install $NUMPY pip install $MPL + pip install -r requirements.txt + + - if [[ $TRAVIS_PYTHON_VERSION == '2.6' ]]; then pip install -r .requirements-2.6.txt; fi + - | cd geos-3.3.3 export GEOS_DIR=$HOME/.local/ diff --git a/Changelog b/Changelog index 8372d02cc..88965d52e 100644 --- a/Changelog +++ b/Changelog @@ -1,5 +1,7 @@ version 1.0.8 (not yet released) -------------------------------- +* removes inline pyshp and pyproj code and uses makes these packages required + external dependencies (handling issue 230). * add 'textcolor' kwarg to drawmeridians and drawparallels (issue 145). * don't assume grid is regular when adding cyclic point in addcyclic. New kwargs 'axis' and 'cyclic' added. More than one array diff --git a/LICENSE_pyshp b/LICENSE_pyshp deleted file mode 100644 index 56e4c61a1..000000000 --- a/LICENSE_pyshp +++ /dev/null @@ -1,3 +0,0 @@ -pyshp is released under the MIT license, which can be found at: - -http://www.opensource.org/licenses/mit-license.php diff --git a/README.md b/README.md index 5ea5b3b66..579e97b53 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,10 @@ using matplotlib. * numpy +* [pyproj](https://github.com/jswhit/pyproj) + +* [pyshp](https://github.com/GeospatialPython/pyshp) + * The GEOS (Geometry Engine - Open Source) library (version 3.1.1 or higher). Source code is included in the geos-3.3.3 directory. @@ -31,9 +35,6 @@ source code for the GEOS library is included in the 'geos-3.3.3' directory under the terms given in LICENSE_geos. -shapefile.py from pyshp.googlecode.com is included under the terms given -in LICENSE_pyshp. - the land-sea mask, coastline, lake, river and political boundary data are extracted from datasets provided with the [Generic Mapping Tools (GMT)](http://gmt.soest.hawaii.edu) and are included under the terms given in LICENSE_data. diff --git a/examples/allskymap.py b/examples/allskymap.py index 2aaefe1c3..cb356da3c 100644 --- a/examples/allskymap.py +++ b/examples/allskymap.py @@ -23,8 +23,9 @@ from numpy import * import matplotlib.pyplot as pl from matplotlib.pyplot import * -from mpl_toolkits.basemap import Basemap, pyproj -from mpl_toolkits.basemap.pyproj import Geod +from mpl_toolkits.basemap import Basemap +import pyproj +from pyproj import Geod __all__ = ['AllSkyMap'] diff --git a/examples/testwmsimage.py b/examples/testwmsimage.py index b9ad0d1f3..d8bd7e8be 100644 --- a/examples/testwmsimage.py +++ b/examples/testwmsimage.py @@ -3,7 +3,8 @@ from a WMS server and display it on a map (using the wmsimage convenience method) """ -from mpl_toolkits.basemap import Basemap, pyproj +from mpl_toolkits.basemap import Basemap +import pyproj from datetime import datetime import numpy as np import matplotlib.pyplot as plt diff --git a/examples/utmtest.py b/examples/utmtest.py index 8b8766fc9..bbf0ea09a 100644 --- a/examples/utmtest.py +++ b/examples/utmtest.py @@ -1,5 +1,5 @@ from __future__ import print_function -from mpl_toolkits.basemap import pyproj +import pyproj from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt # proj4 definition of UTM zone 17. diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py index 68fea3de1..0177a1634 100644 --- a/lib/mpl_toolkits/basemap/__init__.py +++ b/lib/mpl_toolkits/basemap/__init__.py @@ -27,7 +27,7 @@ from matplotlib.patches import Ellipse, Circle, Polygon, FancyArrowPatch from matplotlib.lines import Line2D from matplotlib.transforms import Bbox -from . import pyproj +import pyproj from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib.image import imread import sys, os, math @@ -2128,7 +2128,7 @@ def readshapefile(self,shapefile,name,drawbounds=True,zorder=None, matplotlib.patches.LineCollection object is appended to the tuple. """ import shapefile as shp - from .shapefile import Reader + from shapefile import Reader shp.default_encoding = default_encoding if not os.path.exists('%s.shp'%shapefile): raise IOError('cannot locate %s.shp'%shapefile) diff --git a/lib/mpl_toolkits/basemap/proj.py b/lib/mpl_toolkits/basemap/proj.py index 16510d11d..eaaf5f689 100644 --- a/lib/mpl_toolkits/basemap/proj.py +++ b/lib/mpl_toolkits/basemap/proj.py @@ -1,5 +1,5 @@ import numpy as np -from . import pyproj +import pyproj import math from matplotlib.cbook import dedent diff --git a/lib/mpl_toolkits/basemap/pyproj.py b/lib/mpl_toolkits/basemap/pyproj.py deleted file mode 100644 index cd0950582..000000000 --- a/lib/mpl_toolkits/basemap/pyproj.py +++ /dev/null @@ -1,868 +0,0 @@ -""" -Cython wrapper to provide python interfaces to -PROJ.4 (http://trac.osgeo.org/proj/) functions. - -Performs cartographic transformations and geodetic computations. - -The Proj class can convert from geographic (longitude,latitude) -to native map projection (x,y) coordinates and vice versa, or -from one map projection coordinate system directly to another. -The module variable pj_list is a dictionary containing all the -available projections and their descriptions. - -The Geod class can perform forward and inverse geodetic, or -Great Circle, computations. The forward computation involves -determining latitude, longitude and back azimuth of a terminus -point given the latitude and longitude of an initial point, plus -azimuth and distance. The inverse computation involves -determining the forward and back azimuths and distance given the -latitudes and longitudes of an initial and terminus point. - -Input coordinates can be given as python arrays, lists/tuples, -scalars or numpy/Numeric/numarray arrays. Optimized for objects -that support the Python buffer protocol (regular python and -numpy array objects). - -Download: http://code.google.com/p/pyproj/downloads/list - -Requirements: python 2.4 or higher. - -Example scripts are in 'test' subdirectory of source distribution. -The 'test()' function will run the examples in the docstrings. - -Contact: Jeffrey Whitaker >> from mpl_toolkits.basemap import Proj - >>> p = Proj(proj='utm',zone=10,ellps='WGS84') # use kwargs - >>> x,y = p(-120.108, 34.36116666) - >>> 'x=%9.3f y=%11.3f' % (x,y) - 'x=765975.641 y=3805993.134' - >>> 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True) - 'lon=-120.108 lat=34.361' - >>> # do 3 cities at a time in a tuple (Fresno, LA, SF) - >>> lons = (-119.72,-118.40,-122.38) - >>> lats = (36.77, 33.93, 37.62 ) - >>> x,y = p(lons, lats) - >>> 'x: %9.3f %9.3f %9.3f' % x - 'x: 792763.863 925321.537 554714.301' - >>> 'y: %9.3f %9.3f %9.3f' % y - 'y: 4074377.617 3763936.941 4163835.303' - >>> lons, lats = p(x, y, inverse=True) # inverse transform - >>> 'lons: %8.3f %8.3f %8.3f' % lons - 'lons: -119.720 -118.400 -122.380' - >>> 'lats: %8.3f %8.3f %8.3f' % lats - 'lats: 36.770 33.930 37.620' - >>> p2 = Proj('+proj=utm +zone=10 +ellps=WGS84') # use proj4 string - >>> x,y = p2(-120.108, 34.36116666) - >>> 'x=%9.3f y=%11.3f' % (x,y) - 'x=765975.641 y=3805993.134' - >>> p = Proj(init="epsg:32667") - >>> 'x=%12.3f y=%12.3f (meters)' % p(-114.057222, 51.045) - 'x=-1783486.760 y= 6193833.196 (meters)' - >>> p = Proj("+init=epsg:32667",preserve_units=True) - >>> 'x=%12.3f y=%12.3f (feet)' % p(-114.057222, 51.045) - 'x=-5851322.810 y=20320934.409 (feet)' - """ - # if projparams is None, use kwargs. - if projparams is None: - if len(kwargs) == 0: - raise RuntimeError('no projection control parameters specified') - else: - projstring = _dict2string(kwargs) - elif type(projparams) == str: - # if projparams is a string, interpret as a proj4 init string. - projstring = projparams - else: # projparams a dict - projstring = _dict2string(projparams) - # make sure units are meters if preserve_units is False. - if not projstring.count('+units=') and not preserve_units: - projstring = '+units=m '+projstring - else: - kvpairs = [] - for kvpair in projstring.split(): - if kvpair.startswith('+units') and not preserve_units: - k,v = kvpair.split('=') - kvpairs.append(k+'=m ') - else: - kvpairs.append(kvpair+' ') - projstring = ''.join(kvpairs) - # look for EPSG, replace with epsg (EPSG only works - # on case-insensitive filesystems). - projstring = projstring.replace('EPSG','epsg') - return _proj.Proj.__new__(self, projstring) - - def __call__(self, *args, **kw): - #,lon,lat,inverse=False,radians=False,errcheck=False): - """ - Calling a Proj 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. If optional keyword 'radians' is True - (default is False) the units of lon/lat are radians instead of - degrees. If optional keyword 'errcheck' is True (default is - False) an exception is raised if the transformation is invalid. - If errcheck=False and the transformation is invalid, no - exception is raised and 1.e30 is returned. - - Instead of calling with lon, lat, a single ndarray of - shape n,2 may be used, and one of the same shape will - be returned; this is more efficient. - - Inputs should be doubles (they will be cast to doubles if they - are not, causing a slight performance hit). - - Works with numpy and regular python array objects, python - sequences and scalars, but is fastest for array objects. - """ - inverse = kw.get('inverse', False) - radians = kw.get('radians', False) - errcheck = kw.get('errcheck', False) - #if len(args) == 1: - # latlon = np.array(args[0], copy=True, - # order='C', dtype=float, ndmin=2) - # if inverse: - # _proj.Proj._invn(self, latlon, radians=radians, errcheck=errcheck) - # else: - # _proj.Proj._fwdn(self, latlon, radians=radians, errcheck=errcheck) - # return latlon - lon, lat = args - # process inputs, making copies that support buffer API. - inx, xisfloat, xislist, xistuple = _copytobuffer(lon) - iny, yisfloat, yislist, yistuple = _copytobuffer(lat) - # call proj4 functions. inx and iny modified in place. - if inverse: - _proj.Proj._inv(self, inx, iny, radians=radians, errcheck=errcheck) - else: - _proj.Proj._fwd(self, inx, iny, radians=radians, errcheck=errcheck) - # if inputs were lists, tuples or floats, convert back. - outx = _convertback(xisfloat,xislist,xistuple,inx) - outy = _convertback(yisfloat,yislist,xistuple,iny) - return outx, outy - - def is_latlong(self): - """returns True if projection in geographic (lon/lat) coordinates""" - return _proj.Proj.is_latlong(self) - - def is_geocent(self): - """returns True if projection in geocentric (x/y) coordinates""" - return _proj.Proj.is_geocent(self) - -def transform(p1, p2, x, y, z=None, radians=False): - """ - x2, y2, z2 = transform(p1, p2, x1, y1, z1, radians=False) - - Transform points between two coordinate systems defined by the - Proj instances p1 and p2. - - The points x1,y1,z1 in the coordinate system defined by p1 are - transformed to x2,y2,z2 in the coordinate system defined by p2. - - z1 is optional, if it is not set it is assumed to be zero (and - only x2 and y2 are returned). - - In addition to converting between cartographic and geographic - projection coordinates, this function can take care of datum - shifts (which cannot be done using the __call__ method of the - Proj instances). It also allows for one of the coordinate - systems to be geographic (proj = 'latlong'). - - If optional keyword 'radians' is True (default is False) and p1 - is defined in geographic coordinate (pj.is_latlong() is True), - x1,y1 is interpreted as radians instead of the default degrees. - Similarly, if p2 is defined in geographic coordinates and - radians=True, x2, y2 are returned in radians instead of degrees. - if p1.is_latlong() and p2.is_latlong() both are False, the - radians keyword has no effect. - - x,y and z can be numpy or regular python arrays, python - lists/tuples or scalars. Arrays are fastest. For projections in - geocentric coordinates, values of x and y are given in meters. - z is always meters. - - Example usage: - - >>> # projection 1: UTM zone 15, grs80 ellipse, NAD83 datum - >>> # (defined by epsg code 26915) - >>> p1 = Proj(init='epsg:26915') - >>> # projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum - >>> p2 = Proj(init='epsg:26715') - >>> # find x,y of Jefferson City, MO. - >>> x1, y1 = p1(-92.199881,38.56694) - >>> # transform this point to projection 2 coordinates. - >>> x2, y2 = transform(p1,p2,x1,y1) - >>> '%9.3f %11.3f' % (x1,y1) - '569704.566 4269024.671' - >>> '%9.3f %11.3f' % (x2,y2) - '569722.342 4268814.027' - >>> '%8.3f %5.3f' % p2(x2,y2,inverse=True) - ' -92.200 38.567' - >>> # process 3 points at a time in a tuple - >>> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri - >>> lons = (-92.22,-94.72,-90.37) - >>> x1, y1 = p1(lons,lats) - >>> x2, y2 = transform(p1,p2,x1,y1) - >>> xy = x1+y1 - >>> '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy - '567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005' - >>> xy = x2+y2 - >>> '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy - '567721.149 351747.558 728569.133 4297989.112 4353489.644 4292106.305' - >>> lons, lats = p2(x2,y2,inverse=True) - >>> xy = lons+lats - >>> '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy - ' -92.220 -94.720 -90.370 38.830 39.320 38.750' - >>> # test datum shifting, installation of extra datum grid files. - >>> p1 = Proj(proj='latlong',datum='WGS84') - >>> x1 = -111.5; y1 = 45.25919444444 - >>> p2 = Proj(proj="utm",zone=10,datum='NAD27') - >>> x2, y2 = transform(p1, p2, x1, y1) - >>> "%12.3f %12.3f" % (x2,y2) - ' 1402285.991 5076292.423' - """ - # process inputs, making copies that support buffer API. - inx, xisfloat, xislist, xistuple = _copytobuffer(x) - iny, yisfloat, yislist, yistuple = _copytobuffer(y) - if z is not None: - inz, zisfloat, zislist, zistuple = _copytobuffer(z) - else: - inz = None - # call pj_transform. inx,iny,inz buffers modified in place. - _proj._transform(p1,p2,inx,iny,inz,radians) - # if inputs were lists, tuples or floats, convert back. - outx = _convertback(xisfloat,xislist,xistuple,inx) - outy = _convertback(yisfloat,yislist,xistuple,iny) - if inz is not None: - outz = _convertback(zisfloat,zislist,zistuple,inz) - return outx, outy, outz - else: - return outx, outy - -def _copytobuffer_return_scalar(x): - try: - # inx,isfloat,islist,istuple - return array('d',(float(x),)),True,False,False - except: - raise TypeError('input must be an array, list, tuple or scalar') - -def _copytobuffer(x): - """ - return a copy of x as an object that supports the python Buffer - API (python array if input is float, list or tuple, numpy array - if input is a numpy array). returns copyofx, isfloat, islist, - istuple (islist is True if input is a list, istuple is true if - input is a tuple, isfloat is true if input is a float). - """ - # make sure x supports Buffer API and contains doubles. - isfloat = False; islist = False; istuple = False - # first, if it's a numpy array scalar convert to float - # (array scalars don't support buffer API) - if hasattr(x,'shape'): - if x.shape == (): - return _copytobuffer_return_scalar(x) - else: - try: - # typecast numpy arrays to double. - # (this makes a copy - which is crucial - # since buffer is modified in place) - x.dtype.char - inx = x.copy(order="C").astype('d') - - # inx,isfloat,islist,istuple - return inx,False,False,False - except: - try: # perhaps they are Numeric/numarrays? - # sorry, not tested yet. - # i don't know Numeric/numarrays has `shape'. - x.typecode() - inx = x.astype('d') - # inx,isfloat,islist,istuple - return inx,False,False,False - except: - raise TypeError('input must be an array, list, tuple or scalar') - else: - # perhaps they are regular python arrays? - if hasattr(x, 'typecode'): - #x.typecode - inx = array('d',x) - # try to convert to python array - # a list. - elif type(x) == list: - inx = array('d',x) - islist = True - # a tuple. - elif type(x) == tuple: - inx = array('d',x) - istuple = True - # a scalar? - else: - return _copytobuffer_return_scalar(x) - return inx,isfloat,islist,istuple - -def _convertback(isfloat,islist,istuple,inx): - # if inputs were lists, tuples or floats, convert back to original type. - if isfloat: - return inx[0] - elif islist: - return inx.tolist() - elif istuple: - return tuple(inx) - else: - return inx - -def _dict2string(projparams): - # convert a dict to a proj4 string. - pjargs = [] - for key,value in projparams.items(): - pjargs.append('+'+key+"="+str(value)+' ') - return ''.join(pjargs) - -class Geod(_proj.Geod): - """ - performs forward and inverse geodetic, or Great Circle, - computations. The forward computation (using the 'fwd' method) - involves determining latitude, longitude and back azimuth of a - computations. The forward computation (using the 'fwd' method) - involves determining latitude, longitude and back azimuth of a - terminus point given the latitude and longitude of an initial - point, plus azimuth and distance. The inverse computation (using - the 'inv' method) involves determining the forward and back - azimuths and distance given the latitudes and longitudes of an - initial and terminus point. - """ - def __new__(self, initstring=None, **kwargs): - """ - initialize a Geod class instance. - - Geodetic parameters for specifying the ellipsoid - can be given in a dictionary 'initparams', as keyword arguments, - or as as proj4 geod initialization string. - Following is a list of the ellipsoids that may be defined using the - 'ellps' keyword (these are stored in the model variable pj_ellps):: - - MERIT a=6378137.0 rf=298.257 MERIT 1983 - SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85 - GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980) - IAU76 a=6378140.0 rf=298.257 IAU 1976 - airy a=6377563.396 b=6356256.910 Airy 1830 - APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965 - airy a=6377563.396 b=6356256.910 Airy 1830 - APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965 - NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965 - mod_airy a=6377340.189 b=6356034.446 Modified Airy - andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.) - aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer. 1969 - GRS67 a=6378160.0 rf=298.247167427 GRS 67(IUGG 1967) - bessel a=6377397.155 rf=299.1528128 Bessel 1841 - bess_nam a=6377483.865 rf=299.1528128 Bessel 1841 (Namibia) - clrk66 a=6378206.4 b=6356583.8 Clarke 1866 - clrk80 a=6378249.145 rf=293.4663 Clarke 1880 mod. - CPM a=6375738.7 rf=334.29 Comm. des Poids et Mesures 1799 - delmbr a=6376428. rf=311.5 Delambre 1810 (Belgium) - engelis a=6378136.05 rf=298.2566 Engelis 1985 - evrst30 a=6377276.345 rf=300.8017 Everest 1830 - evrst48 a=6377304.063 rf=300.8017 Everest 1948 - evrst56 a=6377301.243 rf=300.8017 Everest 1956 - evrst69 a=6377295.664 rf=300.8017 Everest 1969 - evrstSS a=6377298.556 rf=300.8017 Everest (Sabah & Sarawak) - fschr60 a=6378166. rf=298.3 Fischer (Mercury Datum) 1960 - fschr60m a=6378155. rf=298.3 Modified Fischer 1960 - fschr68 a=6378150. rf=298.3 Fischer 1968 - helmert a=6378200. rf=298.3 Helmert 1906 - hough a=6378270.0 rf=297. Hough - helmert a=6378200. rf=298.3 Helmert 1906 - hough a=6378270.0 rf=297. Hough - intl a=6378388.0 rf=297. International 1909 (Hayford) - krass a=6378245.0 rf=298.3 Krassovsky, 1942 - kaula a=6378163. rf=298.24 Kaula 1961 - lerch a=6378139. rf=298.257 Lerch 1979 - mprts a=6397300. rf=191. Maupertius 1738 - new_intl a=6378157.5 b=6356772.2 New International 1967 - plessis a=6376523. b=6355863. Plessis 1817 (France) - SEasia a=6378155.0 b=6356773.3205 Southeast Asia - walbeck a=6376896.0 b=6355834.8467 Walbeck - WGS60 a=6378165.0 rf=298.3 WGS 60 - WGS66 a=6378145.0 rf=298.25 WGS 66 - WGS72 a=6378135.0 rf=298.26 WGS 72 - WGS84 a=6378137.0 rf=298.257223563 WGS 84 - sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997) - - The parameters of the ellipsoid may also be set directly using - the 'a' (semi-major or equatorial axis radius) keyword, and - any one of the following keywords: 'b' (semi-minor, - or polar axis radius), 'e' (eccentricity), 'es' (eccentricity - squared), 'f' (flattening), or 'rf' (reciprocal flattening). - - See the proj documentation (http://trac.osgeo.org/proj/) for more - - See the proj documentation (http://trac.osgeo.org/proj/) for more - information about specifying ellipsoid parameters (specifically, - the chapter 'Specifying the Earth's figure' in the main Proj - users manual). - - Example usage: - - >>> from mpl_toolkits.basemap.pyproj import Geod - >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid. - >>> # specify the lat/lons of some cities. - >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.) - >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.) - >>> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.) - >>> london_lat = 51.+(32./60.); london_lon = -(5./60.) - >>> # compute forward and back azimuths, plus distance - >>> # between Boston and Portland. - >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat) - >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist) - '-66.531 75.654 4164192.708' - >>> # compute latitude, longitude and back azimuth of Portland, - >>> # given Boston lat/lon, forward azimuth and distance to Portland. - >>> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist) - >>> "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz) - '45.517 -123.683 75.654' - >>> # compute the azimuths, distances from New York to several - >>> # cities (pass a list) - >>> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat] - >>> lons2 = [boston_lon, portland_lon, london_lon] - >>> lats2 = [boston_lat, portland_lat, london_lat] - >>> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2) - >>> for faz,baz,d in list(zip(az12,az21,dist)): "%7.3f %7.3f %9.3f" % (faz,baz,d) - ' 54.663 -123.448 288303.720' - '-65.463 79.342 4013037.318' - ' 51.254 -71.576 5579916.651' - >>> g2 = Geod('+ellps=clrk66') # use proj4 style initialization string - >>> az12,az21,dist = g2.inv(boston_lon,boston_lat,portland_lon,portland_lat) - >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist) - '-66.531 75.654 4164192.708' - """ - # if initparams is a proj-type init string, - # convert to dict. - ellpsd = {} - if initstring is not None: - for kvpair in initstring.split(): - k,v = kvpair.split('=') - k = k.lstrip('+') - if k in ['a','b','rf','f','es','e']: - v = float(v) - ellpsd[k] = v - # merge this dict with kwargs dict. - kwargs = dict(list(kwargs.items()) + list(ellpsd.items())) - self.sphere = False - if 'ellps' in kwargs: - # ellipse name given, look up in pj_ellps dict - ellps_dict = pj_ellps[kwargs['ellps']] - a = ellps_dict['a'] - if ellps_dict['description']=='Normal Sphere': - self.sphere = True - if 'b' in ellps_dict: - b = ellps_dict['b'] - es = 1. - (b * b) / (a * a) - f = (a - b)/a - elif 'rf' in ellps_dict: - f = 1./ellps_dict['rf'] - b = a*(1. - f) - es = 1. - (b * b) / (a * a) - else: - # a (semi-major axis) and one of - # b the semi-minor axis - # rf the reciprocal flattening - # f flattening - # es eccentricity squared - # must be given. - a = kwargs['a'] - if 'b' in kwargs: - b = kwargs['b'] - es = 1. - (b * b) / (a * a) - f = (a - b)/a - elif 'rf' in kwargs: - f = 1./kwargs['rf'] - b = a*(1. - f) - es = 1. - (b * b) / (a * a) - elif 'f' in kwargs: - f = kwargs['f'] - b = a*(1. - f) - es = 1. - (b/a)**2 - elif 'es' in kwargs: - es = kwargs['es'] - b = math.sqrt(a**2 - es*a**2) - f = (a - b)/a - elif 'e' in kwargs: - es = kwargs['e']**2 - b = math.sqrt(a**2 - es*a**2) - f = (a - b)/a - else: - b = a - f = 0. - es = 0. - #msg='ellipse name or a, plus one of f,es,b must be given' - #raise ValueError(msg) - if math.fabs(f) < 1.e-8: self.sphere = True - self.a = a - self.b = b - self.f = f - self.es = es - return _proj.Geod.__new__(self, a, f) - - def fwd(self, lons, lats, az, dist, radians=False): - """ - forward transformation - Returns longitudes, latitudes and back - azimuths of terminus points given longitudes (lons) and - latitudes (lats) of initial points, plus forward azimuths (az) - and distances (dist). - latitudes (lats) of initial points, plus forward azimuths (az) - and distances (dist). - - Works with numpy and regular python array objects, python - sequences and scalars. - - if radians=True, lons/lats and azimuths are radians instead of - degrees. Distances are in meters. - """ - # process inputs, making copies that support buffer API. - inx, xisfloat, xislist, xistuple = _copytobuffer(lons) - iny, yisfloat, yislist, yistuple = _copytobuffer(lats) - inz, zisfloat, zislist, zistuple = _copytobuffer(az) - ind, disfloat, dislist, distuple = _copytobuffer(dist) - _proj.Geod._fwd(self, inx, iny, inz, ind, radians=radians) - # if inputs were lists, tuples or floats, convert back. - outx = _convertback(xisfloat,xislist,xistuple,inx) - outy = _convertback(yisfloat,yislist,xistuple,iny) - outz = _convertback(zisfloat,zislist,zistuple,inz) - return outx, outy, outz - - def inv(self,lons1,lats1,lons2,lats2,radians=False): - """ - inverse transformation - Returns forward and back azimuths, plus - distances between initial points (specified by lons1, lats1) and - terminus points (specified by lons2, lats2). - - Works with numpy and regular python array objects, python - sequences and scalars. - - if radians=True, lons/lats and azimuths are radians instead of - degrees. Distances are in meters. - """ - # process inputs, making copies that support buffer API. - inx, xisfloat, xislist, xistuple = _copytobuffer(lons1) - iny, yisfloat, yislist, yistuple = _copytobuffer(lats1) - inz, zisfloat, zislist, zistuple = _copytobuffer(lons2) - ind, disfloat, dislist, distuple = _copytobuffer(lats2) - _proj.Geod._inv(self,inx,iny,inz,ind,radians=radians) - # if inputs were lists, tuples or floats, convert back. - outx = _convertback(xisfloat,xislist,xistuple,inx) - outy = _convertback(yisfloat,yislist,xistuple,iny) - outz = _convertback(zisfloat,zislist,zistuple,inz) - return outx, outy, outz - - def npts(self, lon1, lat1, lon2, lat2, npts, radians=False): - """ - Given a single initial point and terminus point (specified by - python floats lon1,lat1 and lon2,lat2), returns a list of - longitude/latitude pairs describing npts equally spaced - intermediate points along the geodesic between the initial and - terminus points. - - if radians=True, lons/lats are radians instead of degrees. - - Example usage: - - >>> from mpl_toolkits.basemap.pyproj import Geod - >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid. - >>> # specify the lat/lons of Boston and Portland. - >>> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid. - >>> # specify the lat/lons of Boston and Portland. - >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.) - >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.) - >>> # find ten equally spaced points between Boston and Portland. - >>> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10) - >>> for lon,lat in lonlats: '%6.3f %7.3f' % (lat, lon) - '43.528 -75.414' - '44.637 -79.883' - '45.565 -84.512' - '46.299 -89.279' - '46.830 -94.156' - '47.149 -99.112' - '47.251 -104.106' - '47.136 -109.100' - '46.805 -114.051' - '46.262 -118.924' - >>> # test with radians=True (inputs/outputs in radians, not degrees) - >>> import math - >>> dg2rad = math.radians(1.) - >>> rad2dg = math.degrees(1.) - >>> lonlats = g.npts(dg2rad*boston_lon,dg2rad*boston_lat,dg2rad*portland_lon,dg2rad*portland_lat,10,radians=True) - >>> for lon,lat in lonlats: '%6.3f %7.3f' % (rad2dg*lat, rad2dg*lon) - '43.528 -75.414' - '44.637 -79.883' - '45.565 -84.512' - '46.299 -89.279' - '46.830 -94.156' - '47.149 -99.112' - '47.251 -104.106' - '47.136 -109.100' - '46.805 -114.051' - '46.262 -118.924' - """ - lons, lats = _proj.Geod._npts(self, lon1, lat1, lon2, lat2, npts, radians=radians) - return list(zip(lons, lats)) - -def test(): - """run the examples in the docstrings using the doctest module""" - import doctest, pyproj - doctest.testmod(pyproj,verbose=True) - -if __name__ == "__main__": test() diff --git a/lib/mpl_toolkits/basemap/shapefile.py b/lib/mpl_toolkits/basemap/shapefile.py deleted file mode 100644 index bf2de5214..000000000 --- a/lib/mpl_toolkits/basemap/shapefile.py +++ /dev/null @@ -1,1167 +0,0 @@ -""" -shapefile.py -Provides read and write support for ESRI Shapefiles. -author: jlawheadgeospatialpython.com -date: 20130727 -version: 1.2.0 -Compatible with Python versions 2.4-3.x -""" - -__version__ = "1.2.0" - -from struct import pack, unpack, calcsize, error -import os -import sys -import time -import array -import tempfile - -# -# Constants for shape types -default_encoding = 'utf-8' -NULL = 0 -POINT = 1 -POLYLINE = 3 -POLYGON = 5 -MULTIPOINT = 8 -POINTZ = 11 -POLYLINEZ = 13 -POLYGONZ = 15 -MULTIPOINTZ = 18 -POINTM = 21 -POLYLINEM = 23 -POLYGONM = 25 -MULTIPOINTM = 28 -MULTIPATCH = 31 - -PYTHON3 = sys.version_info[0] == 3 - -if PYTHON3: - xrange = range - -def b(v): - if PYTHON3: - if isinstance(v, str): - # For python 3 encode str to bytes. - return v.encode(default_encoding) - elif isinstance(v, bytes): - # Already bytes. - return v - else: - # Error. - raise Exception('Unknown input type') - else: - # For python 2 assume str passed in and return str. - return v - -def u(v): - if PYTHON3: - if isinstance(v, bytes): - # For python 3 decode bytes to str. - return v.decode(default_encoding, errors="replace") - elif isinstance(v, str): - # Already str. - return v - else: - # Error. - raise Exception('Unknown input type') - else: - # For python 2 assume str passed in and return str. - return v - -def is_string(v): - if PYTHON3: - return isinstance(v, str) - else: - return isinstance(v, basestring) - -class _Array(array.array): - """Converts python tuples to lits of the appropritate type. - Used to unpack different shapefile header parts.""" - def __repr__(self): - return str(self.tolist()) - -def signed_area(coords): - """Return the signed area enclosed by a ring using the linear time - algorithm at http://www.cgafaq.info/wiki/Polygon_Area. A value >= 0 - indicates a counter-clockwise oriented ring. - """ - xs, ys = map(list, zip(*coords)) - xs.append(xs[1]) - ys.append(ys[1]) - return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(coords)))/2.0 - -class _Shape: - def __init__(self, shapeType=None): - """Stores the geometry of the different shape types - specified in the Shapefile spec. Shape types are - usually point, polyline, or polygons. Every shape type - except the "Null" type contains points at some level for - example verticies in a polygon. If a shape type has - multiple shapes containing points within a single - geometry record then those shapes are called parts. Parts - are designated by their starting index in geometry record's - list of shapes.""" - self.shapeType = shapeType - self.points = [] - - @property - def __geo_interface__(self): - if self.shapeType in [POINT, POINTM, POINTZ]: - return { - 'type': 'Point', - 'coordinates': tuple(self.points[0]) - } - elif self.shapeType in [MULTIPOINT, MULTIPOINTM, MULTIPOINTZ]: - return { - 'type': 'MultiPoint', - 'coordinates': tuple([tuple(p) for p in self.points]) - } - elif self.shapeType in [POLYLINE, POLYLINEM, POLYLINEZ]: - if len(self.parts) == 1: - return { - 'type': 'LineString', - 'coordinates': tuple([tuple(p) for p in self.points]) - } - else: - ps = None - coordinates = [] - for part in self.parts: - if ps == None: - ps = part - continue - else: - coordinates.append(tuple([tuple(p) for p in self.points[ps:part]])) - ps = part - else: - coordinates.append(tuple([tuple(p) for p in self.points[part:]])) - return { - 'type': 'MultiLineString', - 'coordinates': tuple(coordinates) - } - elif self.shapeType in [POLYGON, POLYGONM, POLYGONZ]: - if len(self.parts) == 1: - return { - 'type': 'Polygon', - 'coordinates': (tuple([tuple(p) for p in self.points]),) - } - else: - ps = None - coordinates = [] - for part in self.parts: - if ps == None: - ps = part - continue - else: - coordinates.append(tuple([tuple(p) for p in self.points[ps:part]])) - ps = part - else: - coordinates.append(tuple([tuple(p) for p in self.points[part:]])) - polys = [] - poly = [coordinates[0]] - for coord in coordinates[1:]: - if signed_area(coord) < 0: - polys.append(poly) - poly = [coord] - else: - poly.append(coord) - polys.append(poly) - if len(polys) == 1: - return { - 'type': 'Polygon', - 'coordinates': tuple(polys[0]) - } - elif len(polys) > 1: - return { - 'type': 'MultiPolygon', - 'coordinates': polys - } - -class _ShapeRecord: - """A shape object of any type.""" - def __init__(self, shape=None, record=None): - self.shape = shape - self.record = record - -class ShapefileException(Exception): - """An exception to handle shapefile specific problems.""" - pass - -class Reader: - """Reads the three files of a shapefile as a unit or - separately. If one of the three files (.shp, .shx, - .dbf) is missing no exception is thrown until you try - to call a method that depends on that particular file. - The .shx index file is used if available for efficiency - but is not required to read the geometry from the .shp - file. The "shapefile" argument in the constructor is the - name of the file you want to open. - - You can instantiate a Reader without specifying a shapefile - and then specify one later with the load() method. - - Only the shapefile headers are read upon loading. Content - within each file is only accessed when required and as - efficiently as possible. Shapefiles are usually not large - but they can be. - """ - def __init__(self, *args, **kwargs): - self.shp = None - self.shx = None - self.dbf = None - self.shapeName = "Not specified" - self._offsets = [] - self.shpLength = None - self.numRecords = None - self.fields = [] - self.__dbfHdrLength = 0 - # See if a shapefile name was passed as an argument - if len(args) > 0: - if is_string(args[0]): - self.load(args[0]) - return - if "shp" in kwargs.keys(): - if hasattr(kwargs["shp"], "read"): - self.shp = kwargs["shp"] - if hasattr(self.shp, "seek"): - self.shp.seek(0) - if "shx" in kwargs.keys(): - if hasattr(kwargs["shx"], "read"): - self.shx = kwargs["shx"] - if hasattr(self.shx, "seek"): - self.shx.seek(0) - if "dbf" in kwargs.keys(): - if hasattr(kwargs["dbf"], "read"): - self.dbf = kwargs["dbf"] - if hasattr(self.dbf, "seek"): - self.dbf.seek(0) - if self.shp or self.dbf: - self.load() - else: - raise ShapefileException("Shapefile Reader requires a shapefile or file-like object.") - - def load(self, shapefile=None): - """Opens a shapefile from a filename or file-like - object. Normally this method would be called by the - constructor with the file object or file name as an - argument.""" - if shapefile: - (shapeName, ext) = os.path.splitext(shapefile) - self.shapeName = shapeName - try: - self.shp = open("%s.shp" % shapeName, "rb") - except IOError: - raise ShapefileException("Unable to open %s.shp" % shapeName) - try: - self.shx = open("%s.shx" % shapeName, "rb") - except IOError: - raise ShapefileException("Unable to open %s.shx" % shapeName) - try: - self.dbf = open("%s.dbf" % shapeName, "rb") - except IOError: - raise ShapefileException("Unable to open %s.dbf" % shapeName) - if self.shp: - self.__shpHeader() - if self.dbf: - self.__dbfHeader() - - def __getFileObj(self, f): - """Checks to see if the requested shapefile file object is - available. If not a ShapefileException is raised.""" - if not f: - raise ShapefileException("Shapefile Reader requires a shapefile or file-like object.") - if self.shp and self.shpLength is None: - self.load() - if self.dbf and len(self.fields) == 0: - self.load() - return f - - def __restrictIndex(self, i): - """Provides list-like handling of a record index with a clearer - error message if the index is out of bounds.""" - if self.numRecords: - rmax = self.numRecords - 1 - if abs(i) > rmax: - raise IndexError("Shape or Record index out of range.") - if i < 0: i = range(self.numRecords)[i] - return i - - def __shpHeader(self): - """Reads the header information from a .shp or .shx file.""" - if not self.shp: - raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no shp file found") - shp = self.shp - # File length (16-bit word * 2 = bytes) - shp.seek(24) - self.shpLength = unpack(">i", shp.read(4))[0] * 2 - # Shape type - shp.seek(32) - self.shapeType= unpack("2i", f.read(8)) - # Determine the start of the next record - next = f.tell() + (2 * recLength) - shapeType = unpack(" -10e38: - record.m.append(m) - else: - record.m.append(None) - # Read a single point - if shapeType in (1,11,21): - record.points = [_Array('d', unpack("<2d", f.read(16)))] - # Read a single Z value - if shapeType == 11: - record.z = unpack("i", shx.read(4))[0] * 2) - 100 - numRecords = shxRecordLength // 8 - # Jump to the first record. - shx.seek(100) - for r in range(numRecords): - # Offsets are 16-bit words just like the file length - self._offsets.append(unpack(">i", shx.read(4))[0] * 2) - shx.seek(shx.tell() + 4) - if not i == None: - return self._offsets[i] - - def shape(self, i=0): - """Returns a shape object for a shape in the the geometry - record file.""" - shp = self.__getFileObj(self.shp) - i = self.__restrictIndex(i) - offset = self.__shapeIndex(i) - if not offset: - # Shx index not available so iterate the full list. - for j,k in enumerate(self.iterShapes()): - if j == i: - return k - shp.seek(offset) - return self.__shape() - - def shapes(self): - """Returns all shapes in a shapefile.""" - shp = self.__getFileObj(self.shp) - # Found shapefiles which report incorrect - # shp file length in the header. Can't trust - # that so we seek to the end of the file - # and figure it out. - shp.seek(0,2) - self.shpLength = shp.tell() - shp.seek(100) - shapes = [] - while shp.tell() < self.shpLength: - shapes.append(self.__shape()) - return shapes - - def iterShapes(self): - """Serves up shapes in a shapefile as an iterator. Useful - for handling large shapefiles.""" - shp = self.__getFileObj(self.shp) - shp.seek(0,2) - self.shpLength = shp.tell() - shp.seek(100) - while shp.tell() < self.shpLength: - yield self.__shape() - - def __dbfHeaderLength(self): - """Retrieves the header length of a dbf file header.""" - if not self.__dbfHdrLength: - if not self.dbf: - raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no dbf file found)") - dbf = self.dbf - (self.numRecords, self.__dbfHdrLength) = \ - unpack("6i", 9994,0,0,0,0,0)) - # File length (Bytes / 2 = 16-bit words) - if headerType == 'shp': - f.write(pack(">i", self.__shpFileLength())) - elif headerType == 'shx': - f.write(pack('>i', ((100 + (len(self._shapes) * 8)) // 2))) - # Version, Shape type - f.write(pack("<2i", 1000, self.shapeType)) - # The shapefile's bounding box (lower left, upper right) - if self.shapeType != 0: - try: - f.write(pack("<4d", *self.bbox())) - except error: - raise ShapefileException("Failed to write shapefile bounding box. Floats required.") - else: - f.write(pack("<4d", 0,0,0,0)) - # Elevation - z = self.zbox() - # Measure - m = self.mbox() - try: - f.write(pack("<4d", z[0], z[1], m[0], m[1])) - except error: - raise ShapefileException("Failed to write shapefile elevation and measure values. Floats required.") - - def __dbfHeader(self): - """Writes the dbf header and field descriptors.""" - f = self.__getFileObj(self.dbf) - f.seek(0) - version = 3 - year, month, day = time.localtime()[:3] - year -= 1900 - # Remove deletion flag placeholder from fields - for field in self.fields: - if field[0].startswith("Deletion"): - self.fields.remove(field) - numRecs = len(self.records) - numFields = len(self.fields) - headerLength = numFields * 32 + 33 - recordLength = sum([int(field[2]) for field in self.fields]) + 1 - header = pack('2i", recNum, 0)) - recNum += 1 - start = f.tell() - # Shape Type - if self.shapeType != 31: - s.shapeType = self.shapeType - f.write(pack("i", length)) - f.seek(finish) - - def __shxRecords(self): - """Writes the shx records.""" - f = self.__getFileObj(self.shx) - f.seek(100) - for i in range(len(self._shapes)): - f.write(pack(">i", self._offsets[i] // 2)) - f.write(pack(">i", self._lengths[i])) - - def __dbfRecords(self): - """Writes the dbf records.""" - f = self.__getFileObj(self.dbf) - for record in self.records: - if not self.fields[0][0].startswith("Deletion"): - f.write(b(' ')) # deletion flag - for (fieldName, fieldType, size, dec), value in zip(self.fields, record): - fieldType = fieldType.upper() - size = int(size) - if fieldType.upper() == "N": - value = str(value).rjust(size) - elif fieldType == 'L': - value = str(value)[0].upper() - else: - value = str(value)[:size].ljust(size) - assert len(value) == size - value = b(value) - f.write(value) - - def null(self): - """Creates a null shape.""" - self._shapes.append(_Shape(NULL)) - - def point(self, x, y, z=0, m=0): - """Creates a point shape.""" - pointShape = _Shape(self.shapeType) - pointShape.points.append([x, y, z, m]) - self._shapes.append(pointShape) - - def line(self, parts=[], shapeType=POLYLINE): - """Creates a line shape. This method is just a convienience method - which wraps 'poly()'. - """ - self.poly(parts, shapeType, []) - - def poly(self, parts=[], shapeType=POLYGON, partTypes=[]): - """Creates a shape that has multiple collections of points (parts) - including lines, polygons, and even multipoint shapes. If no shape type - is specified it defaults to 'polygon'. If no part types are specified - (which they normally won't be) then all parts default to the shape type. - """ - polyShape = _Shape(shapeType) - polyShape.parts = [] - polyShape.points = [] - # Make sure polygons are closed - if shapeType in (5,15,25,31): - for part in parts: - if part[0] != part[-1]: - part.append(part[0]) - for part in parts: - polyShape.parts.append(len(polyShape.points)) - for point in part: - # Ensure point is list - if not isinstance(point, list): - point = list(point) - # Make sure point has z and m values - while len(point) < 4: - point.append(0) - polyShape.points.append(point) - if polyShape.shapeType == 31: - if not partTypes: - for part in parts: - partTypes.append(polyShape.shapeType) - polyShape.partTypes = partTypes - self._shapes.append(polyShape) - - def field(self, name, fieldType="C", size="50", decimal=0): - """Adds a dbf field descriptor to the shapefile.""" - self.fields.append((name, fieldType, size, decimal)) - - def record(self, *recordList, **recordDict): - """Creates a dbf attribute record. You can submit either a sequence of - field values or keyword arguments of field names and values. Before - adding records you must add fields for the record values using the - fields() method. If the record values exceed the number of fields the - extra ones won't be added. In the case of using keyword arguments to specify - field/value pairs only fields matching the already registered fields - will be added.""" - record = [] - fieldCount = len(self.fields) - # Compensate for deletion flag - if self.fields[0][0].startswith("Deletion"): fieldCount -= 1 - if recordList: - [record.append(recordList[i]) for i in range(fieldCount)] - elif recordDict: - for field in self.fields: - if field[0] in recordDict: - val = recordDict[field[0]] - if val is None: - record.append("") - else: - record.append(val) - if record: - self.records.append(record) - - def shape(self, i): - return self._shapes[i] - - def shapes(self): - """Return the current list of shapes.""" - return self._shapes - - def saveShp(self, target): - """Save an shp file.""" - if not hasattr(target, "write"): - target = os.path.splitext(target)[0] + '.shp' - if not self.shapeType: - self.shapeType = self._shapes[0].shapeType - self.shp = self.__getFileObj(target) - self.__shapefileHeader(self.shp, headerType='shp') - self.__shpRecords() - - def saveShx(self, target): - """Save an shx file.""" - if not hasattr(target, "write"): - target = os.path.splitext(target)[0] + '.shx' - if not self.shapeType: - self.shapeType = self._shapes[0].shapeType - self.shx = self.__getFileObj(target) - self.__shapefileHeader(self.shx, headerType='shx') - self.__shxRecords() - - def saveDbf(self, target): - """Save a dbf file.""" - if not hasattr(target, "write"): - target = os.path.splitext(target)[0] + '.dbf' - self.dbf = self.__getFileObj(target) - self.__dbfHeader() - self.__dbfRecords() - - def save(self, target=None, shp=None, shx=None, dbf=None): - """Save the shapefile data to three files or - three file-like objects. SHP and DBF files can also - be written exclusively using saveShp, saveShx, and saveDbf respectively. - If target is specified but not shp,shx, or dbf then the target path and - file name are used. If no options or specified, a unique base file name - is generated to save the files and the base file name is returned as a - string. - """ - # Create a unique file name if one is not defined - if shp: - self.saveShp(shp) - if shx: - self.saveShx(shx) - if dbf: - self.saveDbf(dbf) - elif not shp and not shx and not dbf: - generated = False - if not target: - temp = tempfile.NamedTemporaryFile(prefix="shapefile_",dir=os.getcwd()) - target = temp.name - generated = True - self.saveShp(target) - self.shp.close() - self.saveShx(target) - self.shx.close() - self.saveDbf(target) - self.dbf.close() - if generated: - return target -class Editor(Writer): - def __init__(self, shapefile=None, shapeType=POINT, autoBalance=1): - self.autoBalance = autoBalance - if not shapefile: - Writer.__init__(self, shapeType) - elif is_string(shapefile): - base = os.path.splitext(shapefile)[0] - if os.path.isfile("%s.shp" % base): - r = Reader(base) - Writer.__init__(self, r.shapeType) - self._shapes = r.shapes() - self.fields = r.fields - self.records = r.records() - - def select(self, expr): - """Select one or more shapes (to be implemented)""" - # TODO: Implement expressions to select shapes. - pass - - def delete(self, shape=None, part=None, point=None): - """Deletes the specified part of any shape by specifying a shape - number, part number, or point number.""" - # shape, part, point - if shape and part and point: - del self._shapes[shape][part][point] - # shape, part - elif shape and part and not point: - del self._shapes[shape][part] - # shape - elif shape and not part and not point: - del self._shapes[shape] - # point - elif not shape and not part and point: - for s in self._shapes: - if s.shapeType == 1: - del self._shapes[point] - else: - for part in s.parts: - del s[part][point] - # part, point - elif not shape and part and point: - for s in self._shapes: - del s[part][point] - # part - elif not shape and part and not point: - for s in self._shapes: - del s[part] - - def point(self, x=None, y=None, z=None, m=None, shape=None, part=None, point=None, addr=None): - """Creates/updates a point shape. The arguments allows - you to update a specific point by shape, part, point of any - shape type.""" - # shape, part, point - if shape and part and point: - try: self._shapes[shape] - except IndexError: self._shapes.append([]) - try: self._shapes[shape][part] - except IndexError: self._shapes[shape].append([]) - try: self._shapes[shape][part][point] - except IndexError: self._shapes[shape][part].append([]) - p = self._shapes[shape][part][point] - if x: p[0] = x - if y: p[1] = y - if z: p[2] = z - if m: p[3] = m - self._shapes[shape][part][point] = p - # shape, part - elif shape and part and not point: - try: self._shapes[shape] - except IndexError: self._shapes.append([]) - try: self._shapes[shape][part] - except IndexError: self._shapes[shape].append([]) - points = self._shapes[shape][part] - for i in range(len(points)): - p = points[i] - if x: p[0] = x - if y: p[1] = y - if z: p[2] = z - if m: p[3] = m - self._shapes[shape][part][i] = p - # shape - elif shape and not part and not point: - try: self._shapes[shape] - except IndexError: self._shapes.append([]) - - # point - # part - if addr: - shape, part, point = addr - self._shapes[shape][part][point] = [x, y, z, m] - else: - Writer.point(self, x, y, z, m) - if self.autoBalance: - self.balance() - - def validate(self): - """An optional method to try and validate the shapefile - as much as possible before writing it (not implemented).""" - #TODO: Implement validation method - pass - - def balance(self): - """Adds a corresponding empty attribute or null geometry record depending - on which type of record was created to make sure all three files - are in synch.""" - if len(self.records) > len(self._shapes): - self.null() - elif len(self.records) < len(self._shapes): - self.record() - - def __fieldNorm(self, fieldName): - """Normalizes a dbf field name to fit within the spec and the - expectations of certain ESRI software.""" - if len(fieldName) > 11: fieldName = fieldName[:11] - fieldName = fieldName.upper() - fieldName.replace(' ', '_') - -# Begin Testing -def test(): - import doctest - doctest.NORMALIZE_WHITESPACE = 1 - doctest.testfile("README.txt", verbose=1) - -if __name__ == "__main__": - """ - Doctests are contained in the file 'README.txt'. This library was originally developed - using Python 2.3. Python 2.4 and above have some excellent improvements in the built-in - testing libraries but for now unit testing is done using what's available in - 2.3. - """ - test() diff --git a/lib/mpl_toolkits/basemap/test.py b/lib/mpl_toolkits/basemap/test.py index 0a7b7a672..d03a10cdf 100644 --- a/lib/mpl_toolkits/basemap/test.py +++ b/lib/mpl_toolkits/basemap/test.py @@ -1,9 +1,20 @@ +import sys from mpl_toolkits.basemap import Basemap, shiftgrid import numpy as np +import pyproj # beginnings of a test suite. -from numpy.testing import TestCase,assert_almost_equal +from numpy.testing import TestCase, assert_almost_equal + +try: + from unittest import skipIf +except ImportError: + # for new features, fallback to unittest backport for Python 2.4 - 2.6 + from unittest2 import skipIf + +# For Python 3.x this will be true +PY3 = (sys.version_info[0] == 3) class TestRotateVector(TestCase): @@ -148,7 +159,8 @@ def test_less_than_n_by_3_points_should_work(self): lonsout = bm.shiftdata(lonsin[:, :2]) assert_almost_equal(lonsout_expected, lonsout) - +@skipIf(PY3 and pyproj.__version__ <= "1.9.4", + "Test skipped in Python 3.x with pyproj version 1.9.4 and below.") class TestProjectCoords(TestCase): def get_data(self): lons, lats = np.arange(-180, 180, 20), np.arange(-90, 90, 10) @@ -170,7 +182,6 @@ def test_convert(self): def test_results_should_be_same_for_c_and_f_order_arrays(self): - lons, lats, bmp = self.get_data() xx1, yy1 = bmp(lons.copy(order="C"), lats.copy(order="C")) diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 000000000..440c34361 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +pyproj>=1.9.3 +pyshp>=1.2.0 diff --git a/setup.py b/setup.py index 8b804cd89..790439314 100644 --- a/setup.py +++ b/setup.py @@ -1,8 +1,7 @@ import sys, glob, os, subprocess -major, minor1, minor2, s, tmp = sys.version_info -if major==2 and minor1<4 or major<2: - raise SystemExit("""matplotlib and the basemap toolkit require Python 2.4 or later.""") +if sys.version_info < (2, 6): + raise SystemExit("""matplotlib and the basemap toolkit require Python 2.6 or later.""") from distutils.dist import Distribution from distutils.util import convert_path @@ -81,57 +80,41 @@ def checkversion(GEOS_dir): geos_include_dirs=[os.path.join(GEOS_dir,'include'),inc_dirs] geos_library_dirs=[os.path.join(GEOS_dir,'lib'),os.path.join(GEOS_dir,'lib64')] -# proj4 and geos extensions. -deps = glob.glob('src/*.c') -deps.remove(os.path.join('src','_proj.c')) -deps.remove(os.path.join('src','_geoslib.c')) - packages = ['mpl_toolkits','mpl_toolkits.basemap'] namespace_packages = ['mpl_toolkits'] package_dirs = {'':'lib'} -extensions = [Extension("mpl_toolkits.basemap._proj",deps+['src/_proj.c'],include_dirs = ['src'],)] + # can't install _geoslib in mpl_toolkits.basemap namespace, # or Basemap objects won't be pickleable. -if sys.platform == 'win32': + # don't use runtime_library_dirs on windows (workaround # for a distutils bug - http://bugs.python.org/issue2437). - #extensions.append(Extension("mpl_toolkits.basemap._geoslib",['src/_geoslib.c'], - extensions.append(Extension("_geoslib",['src/_geoslib.c'], - library_dirs=geos_library_dirs, - include_dirs=geos_include_dirs, - libraries=['geos'])) +if sys.platform == 'win32': + runtime_lib_dirs = [] else: - #extensions.append(Extension("mpl_toolkits.basemap._geoslib",['src/_geoslib.c'], - extensions.append(Extension("_geoslib",['src/_geoslib.c'], - library_dirs=geos_library_dirs, - runtime_library_dirs=geos_library_dirs, - include_dirs=geos_include_dirs, - libraries=['geos_c'])) + runtime_lib_dirs = geos_library_dirs + +extensions = [ Extension("_geoslib",['src/_geoslib.c'], + library_dirs=geos_library_dirs, + runtime_library_dirs=runtime_lib_dirs, + include_dirs=geos_include_dirs, + libraries=['geos_c']) ] # Specify all the required mpl data -# create pyproj binary datum shift grid files. pathout =\ os.path.join('lib',os.path.join('mpl_toolkits',os.path.join('basemap','data'))) -if sys.argv[1] not in ['sdist','clean']: - cc = ccompiler.new_compiler() - sysconfig.get_config_vars() - sysconfig.customize_compiler(cc) - cc.set_include_dirs(['src']) - objects = cc.compile(['nad2bin.c', 'src/pj_malloc.c']) - execname = 'nad2bin' - cc.link_executable(objects, execname) - llafiles = glob.glob('datumgrid/*.lla') - cmd = os.path.join(os.getcwd(),execname) - for f in llafiles: - fout = os.path.basename(f.split('.lla')[0]) - fout = os.path.join(pathout,fout) - strg = '%s %s < %s' % (cmd, fout, f) - sys.stdout.write('executing %s\n' % strg) - subprocess.call(strg,shell=True) + datafiles = glob.glob(os.path.join(pathout,'*')) datafiles = [os.path.join('data',os.path.basename(f)) for f in datafiles] package_data = {'mpl_toolkits.basemap':datafiles} +requirements = [ + "numpy>=1.2.1", + "matplotlib>=1.0.0", + "pyproj >= 1.9.3", + "pyshp >= 1.2.0" +] + __version__ = "1.0.8" setup( name = "basemap", @@ -146,7 +129,7 @@ def checkversion(GEOS_dir): download_url = "https://downloads.sourceforge.net/project/matplotlib/matplotlib-toolkits/basemap-{0}/basemap-{0}.tar.gz".format(__version__), author = "Jeff Whitaker", author_email = "jeffrey.s.whitaker@noaa.gov", - install_requires = ["numpy>=1.2.1", "matplotlib>=1.0.0"], + install_requires = requirements, platforms = ["any"], license = "OSI Approved", keywords = ["python","plotting","plots","graphs","charts","GIS","mapping","map projections","maps"], @@ -161,7 +144,7 @@ def checkversion(GEOS_dir): packages = packages, namespace_packages = namespace_packages, package_dir = package_dirs, - ext_modules = extensions, + ext_modules = extensions, cmdclass = {'build_py': build_py}, package_data = package_data ) diff --git a/utils/readboundaries_shp.py b/utils/readboundaries_shp.py index 6beb5120f..4b153f027 100644 --- a/utils/readboundaries_shp.py +++ b/utils/readboundaries_shp.py @@ -1,5 +1,5 @@ import numpy as np -from mpl_toolkits.basemap.shapefile import Reader +from shapefile import Reader lsd = 5