8000 Numpification and some reformatting of pyproj and proj · matplotlib/basemap@2044601 · GitHub
[go: up one dir, main page]
Skip to content

Commit 2044601

Browse files
efiringefiring
committed
Numpification and some reformatting of pyproj and proj
svn path=/trunk/toolkits/basemap/; revision=4071
1 parent 3b88110 commit 2044601

File tree

2 files changed

+425
-348
lines changed

2 files changed

+425
-348
lines changed

‎lib/matplotlib/toolkits/basemap/proj.py

Lines changed: 126 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -1,41 +1,49 @@
1-
import matplotlib.numerix as NX
1+
import numpy as npy
22
import pyproj
33
import math
44

55
__version__ = '1.2.2'
66
_dg2rad = math.radians(1.)
77
_rad2dg = math.degrees(1.)
88

9+
_upper_right_out_of_bounds = (
10+
'the upper right corner of the plot is not in the map projection region')
11+
12+
_lower_left_out_of_bounds = (
13+
'the lower left corner of the plot is not in the map projection region')
14+
15+
916
class Proj(object):
1017
"""
11-
peforms cartographic transformations (converts from longitude,latitude
12-
to native map projection x,y coordinates and vice versa) using proj
13-
(http://proj.maptools.org/)
14-
Uses a pyrex generated C-interface to libproj.
15-
16-
__init__ method sets up projection information.
17-
__call__ method compute transformations.
18-
See docstrings for __init__ and __call__ for details.
18+
peforms cartographic transformations (converts from longitude,latitude
19+
to native map projection x,y coordinates and vice versa) using proj
20+
(http://proj.maptools.org/)
21+
Uses a pyrex generated C-interface to libproj.
1922
20-
Contact: Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
23+
__init__ method sets up projection information.
24+
__call__ method compute transformations.
25+
See docstrings for __init__ and __call__ for details.
26+
27+
Contact: Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
2128
"""
2229

23-
def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislatlon=True):
30+
def __init__(self,projparams,llcrnrlon,llcrnrlat,
31+
urcrnrlon,urcrnrlat,urcrnrislatlon=True):
2432
"""
25-
initialize a Proj class instance.
33+
initialize a Proj class instance.
2634
27-
Input 'projparams' is a dictionary containing proj map
28-
projection control parameter key/value pairs.
29-
See the proj documentation (http://www.remotesensing.org/proj/)
30-
for details.
35+
Input 'projparams' is a dictionary containing proj map
36+
projection control parameter key/value pairs.
37+
See the proj documentation (http://www.remotesensing.org/proj/)
38+
for details.
3139
32-
llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower left hand corner
33-
of projection region.
40+
llcrnrlon,llcrnrlat are lon and lat (in degrees) of lower
41+
left hand corner of projection region.
3442
35-
urcrnrlon,urcrnrlat are lon and lat (in degrees) of upper right hand corner
36-
of projection region if urcrnrislatlon=True (default). Otherwise,
37-
urcrnrlon,urcrnrlat are x,y in projection coordinates (units meters),
38-
assuming the lower left corner is x=0,y=0.
43+
urcrnrlon,urcrnrlat are lon and lat (in degrees) of upper
44+
right hand corner of projection region if urcrnrislatlon=True
45+
(default). Otherwise, urcrnrlon,urcrnrlat are x,y in projection
46+
coordinates (units meters), assuming the lower left corner is x=0,y=0.
3947
"""
4048
self.projparams = projparams
4149
self.projection = projparams['proj']
@@ -63,8 +71,9 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislat
6371
llcrnrx = llcrnrlon
6472
llcrnry = llcrnrlat
6573
elif self.projection == 'ortho':
66-
if llcrnrlon == -180 and llcrnrlat == -90 and urcrnrlon == 180 and urcrnrlat == 90:
67-
self._fulldisk = True
74+
if (llcrnrlon == -180 and llcrnrlat == -90 and
75+
urcrnrlon == 180 and urcrnrlat == 90):
76+
self._fulldisk = True
6877
self._proj4 = pyproj.Proj(projparams)
6978
llcrnrx = -self.rmajor
7079
llcrnry = -self.rmajor
@@ -75,27 +84,28 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislat
7584
self._proj4 = pyproj.Proj(projparams)
7685
llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
7786
if llcrnrx > 1.e20 or llcrnry > 1.e20:
78-
raise ValueError('the lower left corner of the plot is not in the map projection region')
87+
raise ValueError(_lower_left_out_of_bounds)
7988
elif self.projection == 'geos':
8089
self._proj4 = pyproj.Proj(projparams)
8190
# find major and minor axes of ellipse defining map proj region.
8291
delta = 0.01
83-
lats = NX.arange(0,90,delta)
92+
lats = npy.arange(0,90,delta)
8493
lon_0 = projparams['lon_0']
85-
lons = lon_0*NX.ones(len(lats),'d')
94+
lons = lon_0*npy.ones(len(lats),'d')
8695
x, y = self._proj4(lons, lats)
8796
yi = (y > 1.e20).tolist()
8897
ny = yi.index(1)-1
8998
height = y[ny]
90-
lons = NX.arange(lon_0,lon_0+90,delta)
91-
lats = NX.zeros(len(lons),'d')
99+
lons = npy.arange(lon_0,lon_0+90,delta)
100+
lats = npy.zeros(len(lons),'d')
92101
x, y = self(lons, lats)
93102
xi = (x > 1.e20).tolist()
94103
nx = xi.index(1)-1
95104
width = x[nx]
96105
self._height = height
97106
self._width = width
98-
if llcrnrlon == -180 and llcrnrlat == -90 and urcrnrlon == 180 and urcrnrlat == 90:
107+
if (llcrnrlon == -180 and llcrnrlat == -90 and
108+
urcrnrlon == 180 and urcrnrlat == 90):
99109
self._fulldisk = True
100110
llcrnrx = -width
101111
llcrnry = -height
@@ -105,7 +115,7 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislat
105115
self._fulldisk = False
106116
llcrnrx, llcrnry = self(llcrnrlon,llcrnrlat)
107117
if llcrnrx > 1.e20 or llcrnry > 1.e20:
108-
raise ValueError('the lower left corner of the plot is not in the map projection region')
118+
raise ValueError(_lower_left_out_of_bounds)
109119
elif self.projection in ['moll','robin','sinu']:
110120
self._proj4 = pyproj.Proj(projparams)
111121
xtmp,urcrnry = self(projparams['lon_0'],90.)
@@ -116,7 +126,7 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislat
116126
# note that for 'cyl' x,y == lon,lat
117127
self.projparams['x_0']=-llcrnrx
118128
self.projparams['y_0']=-llcrnry
119-
# reset with x_0, y_0.
129+
# reset with x_0, y_0.
120130
if self.projection != 'cyl':
121131
self._proj4 = pyproj.Proj(projparams)
122132
llcrnry = 0.
@@ -136,15 +146,15 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislat
136146
else:
137147
urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
138148
if urcrnrx > 1.e20 or urcrnry > 1.e20:
139-
raise ValueError('the upper right corner of the plot is not in the map projection region')
149+
raise ValueError(_upper_right_out_of_bounds)
140150
elif self.projection == 'geos':
141151
if self._fulldisk:
142152
urcrnrx = 2.*self._width
143153
urcrnry = 2.*self._height
144154
else:
145155
urcrnrx,urcrnry = self(urcrnrlon,urcrnrlat)
146156
if urcrnrx > 1.e20 or urcrnry > 1.e20:
147-
raise ValueError('the upper right corner of the plot is not in the map projection region')
157+
raise ValueError(_upper_right_out_of_bounds)
148158
elif self.projection in ['moll','robin','sinu']:
149159
xtmp,urcrnry = self(projparams['lon_0'],90.)
150160
urcrnrx,xtmp = self(projparams['lon_0']+180.,0)
@@ -172,24 +182,38 @@ def __init__(self,projparams,llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat,urcrnrislat
172182
self.ymax = llcrnry
173183
self.ymin = urcrnry
174184

175-
def __call__(self,x,y,inverse=False):
185+
def __call__(self, *args, **kw):
186+
# x,y,inverse=False):
176187
"""
177-
Calling a Proj class instance with the arguments lon, lat will
178-
convert lon/lat (in degrees) to x/y native map projection
179-
coordinates (in meters). If optional keyword 'inverse' is
180-
True (default is False), the inverse transformation from x/y
181-
to lon/lat is performed.
188+
Calling a Proj class instance with the arguments lon, lat will
189+
convert lon/lat (in degrees) to x/y native map projection
190+
coordinates (in meters). If optional keyword 'inverse' is
191+
True (default is False), the inverse transformation from x/y
192+
to lon/lat is performed.
182193
183-
For cylindrical equidistant projection ('cyl'), this
184-
does nothing (i.e. x,y == lon,lat).
194+
For cylindrical equidistant projection ('cyl'), this
195+
does nothing (i.e. x,y == lon,lat).
185196
186-
lon,lat can be either scalar floats or N arrays.
197+
lon,lat can be either scalar floats or N arrays.
187198
"""
199+
if len(args) == 1:
200+
xy = args[0]
201+
onearray = True
202+
else:
203+
x,y = args
204+
onearray = False
188205
if self.projection == 'cyl': # for cyl x,y == lon,lat
189-
return x,y
206+
if onearray:
207+
return xy
208+
else:
209+
return x,y
210+
inverse = kw.get('inverse', False)
211+
if onearray:
212+
outxy = self._proj4(xy, inverse=inverse)
213+
else:
214+
outx,outy = self._proj4(x, y, inverse=inverse)
190215
if inverse:
191-
outx,outy = self._proj4(x,y,inverse=True)
192-
if self.projection in ['merc','mill']:
216+
if self.projection in ['merc','mill']:
193217
if self.projection == 'merc':
194218
coslat = math.cos(math.radians(self.projparams['lat_ts']))
195219
sinlat = math.sin(math.radians(self.projparams['lat_ts']))
@@ -199,12 +223,14 @@ def __call__(self,x,y,inverse=False):
199223
# radius of curvature of the ellipse perpendicular to
200224
# the plane of the meridian.
201225
rcurv = self.rmajor*coslat/math.sqrt(1.-self.esq*sinlat**2)
202-
try: # x a scalar or an array
203-
outx = _rad2dg*(x/rcurv) + self.llcrnrlon
204-
except: # x a sequence
205-
outx = [_rad2dg*(xi/rcurv) + self.llcrnrlon for xi in x]
226+
if onearray:
227+
outxy[:,0] = _rad2dg*(xy[:,0]/rcurv) + self.llcrnrlon
228+
else:
229+
try: # x a scalar or an array
230+
outx = _rad2dg*(x/rcurv) + self.llcrnrlon
231+
except: # x a sequence
232+
outx = [_rad2dg*(xi/rcurv) + self.llcrnrlon for xi in x]
206233
else:
207-
outx,outy = self._proj4(x,y)
208234
if self.projection in ['merc','mill']:
209235
if self.projection == 'merc':
210236
coslat = math.cos(math.radians(self.projparams['lat_ts']))
@@ -215,28 +241,52 @@ def __call__(self,x,y,inverse=False):
215241
# radius of curvature of the ellipse perpendicular to
216242
# the plane of the meridian.
217243
rcurv = self.rmajor*coslat/math.sqrt(1.-self.esq*sinlat**2)
218-
try: # x is a scalar or an array
219-
outx = rcurv*_dg2rad*(x-self.llcrnrlon)
220-
except: # x is a sequence.
221-
outx = [rcurv*_dg2rad*(xi-self.llcrnrlon) for xi in x]
222-
return outx,outy
244+
if onearray:
245+
outxy[:,0] = rcurv*_dg2rad*(xy[:,0]-self.llcrnrlon)
246+
else:
247+
try: # x is a scalar or an array
248+
outx = rcurv*_dg2rad*(x-self.llcrnrlon)
249+
except: # x is a sequence.
250+
outx = [rcurv*_dg2rad*(xi-self.llcrnrlon) for xi in x]
251+
if onearray:
252+
return outxy
253+
else:
254+
return outx, outy
223255

224256
def makegrid(self,nx,ny,returnxy=False):
225257
"""
226-
return arrays of shape (ny,nx) containing lon,lat coordinates of
227-
an equally spaced native projection grid.
228-
if returnxy=True, the x,y values of the grid are returned also.
258+
return arrays of shape (ny,nx) containing lon,lat coordinates of
259+
an equally spaced native projection grid.
260+
if returnxy=True, the x,y values of the grid are returned also.
229261
"""
230262
dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
231263
dy = (self.urcrnry-self.llcrnry)/(ny-1)
232-
x = self.llcrnrx+dx*NX.indices((ny,nx),NX.Float32)[1,:,:]
233-
y = self.llcrnry+dy*NX.indices((ny,nx),NX.Float32)[0,:,:]
264+
x = self.llcrnrx+dx*npy.indices((ny,nx),npy.float32)[1,:,:]
265+
y = self.llcrnry+dy*npy.indices((ny,nx),npy.float32)[0,:,:]
234266
lons, lats = self(x, y, inverse=True)
235267
if returnxy:
236268
return lons, lats, x, y
237269
else:
238270
return lons, lats
239271

272+
def makegrid3d(self,nx,ny,returnxy=False):
273+
"""
274+
return array of shape (ny,nx, 2) containing lon,lat coordinates of
275+
an equally spaced native projection grid.
276+
if returnxy=True, the x,y values of the grid are returned also.
277+
"""
278+
dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
279+
dy = (self.urcrnry-self.llcrnry)/(ny-1)
280+
xy = npy.empty((ny,nx,2), npy.float64)
281+
xy[...,0] = self.llcrnrx+dx*npy.indices((ny,nx),npy.float32)[1,:,:]
282+
xy[...,1] = self.llcrnry+dy*npy.indices((ny,nx),npy.float32)[0,:,:]
283+
lonlat = self(xy, inverse=True)
284+
if returnxy:
285+
return lonlat, xy
286+
else:
287+
return lonlat
288+
289+
240290
if __name__ == "__main__":
241291

242292
params = {}
@@ -247,10 +297,10 @@ def makegrid(self,nx,ny,returnxy=False):
247297
params['lon_0'] = -107
248298
nx = 349; ny = 277; dx = 32463.41; dy = dx
249299
awips221 = Proj(params,-145.5,1.0,(nx-1)*dx,(ny-1)*dy,urcrnrislatlon=False)
250-
# AWIPS grid 221 parameters
251-
# (from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html)
300+
# AWIPS grid 221 parameters
301+
# (from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html)
252302
llcornerx, llcornery = awips221(-145.5,1.)
253-
# find 4 lon/lat corners of AWIPS grid 221.
303+
# find 4 lon/lat corners of AWIPS grid 221.
254304
llcornerx = 0.; llcornery = 0.
255305
lrcornerx = dx*(nx-1); lrcornery = 0.
256306
ulcornerx = 0.; ulcornery = dy*(ny-1)
@@ -270,13 +320,22 @@ def makegrid(self,nx,ny,returnxy=False):
270320
print ' -68.318 0.897'
271321
print ' -2.566 46.352'
272322
print ' 148.639 46.635'
273-
# compute lons and lats for the whole AWIPS grid 221 (377x249).
323+
# compute lons and lats for the whole AWIPS grid 221 (377x249).
274324
import time; t1 = time.clock()
275325
lons, lats = awips221.makegrid(nx,ny)
276326
t2 = time.clock()
277327
print 'compute lats/lons for all points on AWIPS 221 grid (%sx%s)' %(nx,ny)
278328
print 'max/min lons'
279-
print min(NX.ravel(lons)),max(NX.ravel(lons))
329+
print min(npy.ravel(lons)),max(npy.ravel(lons))
280330
print 'max/min lats'
281-
print min(NX.ravel(lats)),max(NX.ravel(lats))
331+
print min(npy.ravel(lats)),max(npy.ravel(lats))
332+
print 'took',t2-t1,'secs'
333+
print 'Same thing but with a single 3-D array'
334+
t1 = time.clock()
335+
lonlat, xy = awips221.makegrid3d(nx,ny, returnxy=True)
336+
t2 = time.clock()
282337
print 'took',t2-t1,'secs'
338+
339+
assert (lons==lonlat[...,0]).all(), "The longitudes are different"
340+
assert (lats==lonlat[...,1]).all(), "The latitudes are different"
341+

0 commit comments

Comments
 (0)
0