matplotlib-devel Mailing List for matplotlib

On Wed, Sep 12, 2007 at 01:11:54PM -0500, John Hunter wrote:
> Then the front end (eg Line2D) could do something like
> def on_renderer_change(self, renderer):
> # on renderer change; path data already has the
> # separable/nonlinear part handled
> self.pathid = renderer.push_path(pathdata)
> Additionally, you would need to track when either the data or
> nonlinear mapping function are changed in order to remove the old
> path and push out a new one. Then at draw time, you would not need
> to push any data to the backend
> def on_draw(self, renderer):
> # on draw the line just needs to inform the backend to draw the
> # cached path with the current separable transformation
> renderer.draw_path(self.pathid, gc, septrans)
I am a bit tired, and I haven't been following the discussion too
closely, but I have the feeling this is the kind of pattern that Traits
makes both obvious and optimizes a lot. If you give me a concrete minimal
example I would be able to say more, and maybe try to see how it flows in
Traits.
I guess I am just trying to sell Traits, but it seems to me it solves the
problem very well, and has been well optimized (quicker than vanilla
Python for these kind of things).
If you are interested checking this line out, posting a minimal example
of what you are trying to achieve on the enthought-dev mailing-list would
get you higher quality answers.
Gaël

On 9/11/07, James Evans <jre...@ea...> wrote:
> When using the 'bar' plot command, unit information is lost when changing
> the units for a given axis. The attached script demonstrates this.
> This is not the case for line plots because of the use of Line2D's 'recache'
> method.
While the script clearly reveals a bug, but I don't think it is the
recache method (or its absence) that is the culprit. Rather it is the
conversion of the Axe.bar arguments (left, bottom, etc...) to arrays
before they can get stored array in their respective Rectangle
instances. I removed this conversion, and bar now should create the
Rectangle instances with unit info intact.
But there is an additional problem -- bar specifies the width and
bottom as optional scalar arguments, and to be safe, I suggest you
pass in dimensional values for these as well when using a unitized bar
chart. Otherwise they will be scalar and will default to the current
units, which is probably not what you want under a unit change (your
bar widths will be different under changes of units).
Here is a modified test script that appears to be working properly
under svn (r3844) with my recent change:
from pylab import figure, show, nx
cms = cm *nx.arange(0, 10, 2)
bottom=0*cm
width=0.8*cm
fig = figure()
ax1 = fig.add_subplot(2,2,1)
ax1.bar(cms, cms, bottom=bottom)
ax2 = fig.add_subplot(2,2,2)
ax2.bar(cms, cms, bottom=bottom, width=width, xunits=cm, yunits=inch)
ax3 = fig.add_subplot(2,2,3)
ax3.bar(cms, cms, bottom=bottom, width=width, xunits=inch, yunits=cm)
ax3.set_xlim(3, 6) # scalars are interpreted in current units
ax4 = fig.add_subplot(2,2,4)
ax4.bar(cms, cms, bottom=bottom, width=width, xunits=inch, yunits=inch)
#fig.savefig('simple_conversion_plot.png')
ax4.set_xlim(3*cm, 6*cm) # cm are converted to inches
show()

Jeff,
all tests pass.
Thanks again,
David
2007年9月12日, Jeff Whitaker <js...@fa...>:
>
> David Huard wrote:
> > Hi Jeff,
> >
> > Here is the output. Thanks for looking into this.
> >
> >
> > In [1]: from matplotlib.toolkits.basemap import pyproj
> >
> > In [2]: pyproj.test()
> > Trying:
> > from pyproj import Geod
> > Expecting nothing
> > ok
> > Trying:
> > g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
> > Expecting nothing
> > ok
> > Trying:
> > boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
> > Expecting nothing
> > ok
> > Trying:
> > portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
> > Expecting nothing
> > ok
> > Trying:
> > newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
> > Expecting nothing
> > ok
> > Trying:
> > london_lat = 51.+(32./60.); london_lon = -(5./60.)
> > Expecting nothing
> > ok
> > Trying:
> > az12,az21,dist =
> > g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> > Expecting nothing
> > **********************************************************************
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 401, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> > Failed example:
> > az12,az21,dist =
> > g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> > Exception raised:
> > Traceback (most recent call last):
> > File " doctest.py", line 1248, in __run
> > compileflags, 1) in test.globs
> > File "<doctest
> > matplotlib.toolkits.basemap.pyproj.Geod.__new__[6]>", line 1, in ?
> > az12,az21,dist =
> > g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 478, in inv
> > _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> > File "_geod.pyx", line 123, in _geod.Geod._inv
> > ValueError: undefined inverse geodesic (may be an antipodal point)
> > Trying:
> > print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> > Expecting:
> > -66.531 75.654 4164192.708
> > **********************************************************************
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 402, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> > Failed example:
> > print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> > Exception raised:
> > Traceback (most recent call last):
> > File "doctest.py", line 1248, in __run
> > compileflags, 1) in test.globs
> > File "<doctest matplotlib.toolkits.basemap.pyproj.Geod.__new_
> > _[7]>", line 1, in ?
> > print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> > NameError: name 'az12' is not defined
> > Trying:
> > endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> > Expecting nothing
> > **********************************************************************
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 406, in matplotlib.toolkits.basemap.pyproj.Geod.__new_ _
> > Failed example:
> > endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> > Exception raised:
> > Traceback (most recent call last):
> > File "doctest.py", line 1248, in __run
> > compileflags, 1) in test.globs
> > File "<doctest
> > matplotlib.toolkits.basemap.pyproj.Geod.__new__[8]>", line 1, in ?
> > endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12,
> dist)
> > NameError: name 'az12' is not defined
> > Trying:
> > print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> > Expecting:
> > 45.517 -123.683 75.654
> > **********************************************************************
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 407, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> > Failed example:
> > print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> > Exception raised:
> > Traceback (most recent call last):
> > File "doctest.py ", line 1248, in __run
> > compileflags, 1) in test.globs
> > File "<doctest
> > matplotlib.toolkits.basemap.pyproj.Geod.__new__[9]>", line 1, in ?
> > print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> > NameError: name 'endlat' is not defined
> > Trying:
> > lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
> > Expecting nothing
> > ok
> > Trying:
> > lons2 = [boston_lon, portland_lon, london_lon]
> > Expecting nothing
> > ok
> > Trying:
> > lats2 = [boston_lat, portland_lat, london_lat]
> > Expecting nothing
> > ok
> > Trying:
> > az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> > Expecting nothing
> > **********************************************************************
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 414, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> > Failed example:
> > az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> > Exception raised:
> > Traceback (most recent call last):
> > File "doctest.py", line 1248, in __run
> > compileflags, 1) in test.globs
> > File "<doctest matplotlib.toolkits.basemap.pyproj.Geod.__new_
> > _[13]>", line 1, in ?
> > az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 478, in inv
> > _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> > File "_geod.pyx", line 123, in _geod.Geod._inv
> > ValueError: undefined inverse geodesic (may be an antipodal point)
> > Trying:
> > for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" %
> > (faz,baz,d)
> > Expecting:
> > 54.663 -123.448 288303.720
> > -65.463 79.342 4013037.318
> > 51.254 -71.576 5579916.649
> > **********************************************************************
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> > line 415, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> > Failed example:
> > for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" %
> > (faz,baz,d)
> > Exception raised:
> > Traceback (most recent call last):
> > File "doctest.py", line 1248, in __run
> > compileflags, 1) in test.globs
> > File "<doctest
> > matplotlib.toolkits.basemap.pyproj.Geod.__new__[14]>", line 1, in ?
> > for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f
> > %9.3f" % (faz,baz,d)
> > NameError: name 'az12' is not defined
> > Trying:
> > from pyproj import Geod
> > Expecting nothing
> > ok
> > Trying:
> > g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
> > Expecting nothing
> > ok
> > Trying:
> > boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
> > Expecting nothing
> > ok
> > Trying:
> > portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
> > Expecting nothing
> > ok
> > Trying:
> > lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
> > Expecting nothing
> > ok
> > Trying:
> > for lon,lat in lonlats: print '%6.3f %7.3f' % (lat, lon)
> > Expecting:
> > 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
> > ok
> > Trying:
> > from pyproj import Proj
> > Expecting nothing
> > ok
> > Trying:
> > p = Proj(proj='utm',zone=10,ellps='WGS84')
> > Expecting nothing
> > ok
> > Trying:
> > x,y = p(-120.108, 34.36116666)
> > Expecting nothing
> > ok
> > Trying:
> > print 'x=%9.3f y=%11.3f' % (x,y)
> > Expecting:
> > x=765975.641 y=3805993.134
> > ok
> > Trying:
> > print 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
> > Expecting:
> > lon=-120.108 lat=34.361
> > ok
> > Trying:
> > lons = (-119.72,-118.40,-122.38)
> > Expecting nothing
> > ok
> > Trying:
> > lats = (36.77, 33.93, 37.62 )
> > Expecting nothing
> > ok
> > Trying:
> > x,y = p(lons, lats)
> > Expecting nothing
> > ok
> > Trying:
> > print 'x: %9.3f %9.3f %9.3f' % x
> > Expecting:
> > x: 792763.863 925321.537 554714.301
> > ok
> > Trying:
> > print 'y: %9.3f %9.3f %9.3f' % y
> > Expecting:
> > y: 4074377.617 3763936.941 4163835.303
> > ok
> > Trying:
> > lons, lats = p(x, y, inverse=True) # inverse transform
> > Expecting nothing
> > ok
> > Trying:
> > print 'lons: %8.3f %8.3f %8.3f' % lons
> > Expecting:
> > lons: -119.720 -118.400 -122.380
> > ok
> > Trying:
> > print 'lats: %8.3f %8.3f %8.3f' % lats
> > Expecting:
> > lats: 36.770 33.930 37.620
> > ok
> > Trying:
> > p1 = Proj(init='epsg:26915')
> > Expecting nothing
> > ok
> > Trying:
> > p2 = Proj(init='epsg:26715')
> > Expecting nothing
> > ok
> > Trying:
> > x1, y1 = p1(-92.199881,38.56694)
> > Expecting nothing
> > ok
> > Trying:
> > x2, y2 = transform(p1,p2,x1,y1)
> > Expecting nothing
> > ok
> > Trying:
> > print '%9.3f %11.3f' % (x1,y1)
> > Expecting:
> > 569704.566 4269024.671
> > ok
> > Trying:
> > print '%9.3f %11.3f' % (x2,y2)
> > Expecting:
> > 569706.333 4268817.680
> > ok
> > Trying:
> > print '%8.3f %5.3f' % p2(x2,y2,inverse=True)
> > Expecting:
> > -92.200 38.567
> > ok
> > Trying:
> > lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
> > Expecting nothing
> > ok
> > Trying:
> > lons = (-92.22,-94.72,-90.37)
> > Expecting nothing
> > ok
> > Trying:
> > x1, y1 = p1(lons,lats)
> > Expecting nothing
> > ok
> > Trying:
> > x2, y2 = transform(p1,p2,x1,y1)
> > Expecting nothing
> > ok
> > Trying:
> > xy = x1+y1
> > Expecting nothing
> > ok
> > Trying:
> > print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
> > Expecting:
> > 567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005
> > ok
> > Trying:
> > xy = x2+y2
> > Expecting nothing
> > ok
> > Trying:
> > print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
> > Expecting:
> > 567705.072 351727.113 728558.917 4297993.157 4353490.111 4292111.678
> > ok
> > Trying:
> > lons, lats = p2(x2,y2,inverse=True)
> > Expecting nothing
> > ok
> > Trying:
> > xy = lons+lats
> > Expecting nothing
> > ok
> > Trying:
> > print '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy
> > Expecting:
> > -92.220 -94.720 -90.370 38.830 39.320 38.750
> > ok
> > 11 items had no tests:
> > matplotlib.toolkits.basemap.pyproj
> > matplotlib.toolkits.basemap.pyproj.Geod
> > matplotlib.toolkits.basemap.pyproj.Geod.fwd
> > matplotlib.toolkits.basemap.pyproj.Geod.inv
> > matplotlib.toolkits.basemap.pyproj.Proj
> > matplotlib.toolkits.basemap.pyproj.Proj.__call__
> > matplotlib.toolkits.basemap.pyproj.Proj.is_geocent
> > matplotlib.toolkits.basemap.pyproj.Proj.is_latlong
> > matplotlib.toolkits.basemap.pyproj._convertback
> > matplotlib.toolkits.basemap.pyproj._copytobuffer
> > matplotlib.toolkits.basemap.pyproj.test
> > 3 items passed all tests:
> > 6 tests in matplotlib.toolkits.basemap.pyproj.Geod.npts
> > 13 tests in matplotlib.toolkits.basemap.pyproj.Proj.__new__
> > 18 tests in matplotlib.toolkits.basemap.pyproj.transform
> > **********************************************************************
> > 1 items had failures:
> > 6 of 15 in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> > 52 tests in 15 items.
> > 46 passed and 6 failed.
> > ***Test Failed*** 6 failures.
> >
> >
> > 2007年9月12日, Jeff Whitaker < js...@fa...
> > <mailto:js...@fa...>>:
> >
> > David Huard wrote:
> > > Hi, the pyproj package seems to cause a problem in the polarmap
> > > example of the basemap toolkit.
> > >
> > > Thanks,
> > >
> > > david
> > >
> > > huardda@allard:~/svnrepos/toolkits/basemap/examples$ python
> > polarmaps.py
> > > min/max etopo20 data:
> > > -9026.625 6228.8125
> > > plotting North Polar Lambert Azimuthal Equal Area basemap ...
> > > plotting North Polar Stereographic basemap ...
> > > plotting North Polar Azimuthal Equidistant basemap ...
> > > Traceback (most recent call last):
> > > File "polarmaps.py", line 51, in ?
> > >
> resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
> > > File
> > >
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/basemap.py",
> >
> > > line 1111, in __init__
> > >
> > >
> > az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[0]),
> math.radians(latsnew[0]),radians=True)
> > > File
> > >
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
> >
> > > line 478, in inv
> > > _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> > > File "_geod.pyx", line 123, in _geod.Geod._inv
> > > ValueError: undefined inverse geodesic (may be an antipodal point)
> > >
> > >
> > >
> > > This is from a fresh SVN version of both matplotlib and basemap.
> > > Linux, Ubuntu edgy, Xeon-64.
> >
> > David: Odd - I can't reproduce that on my mac. Can you try this
> >
> > >>> from matplotlib.toolkits.basemap import pyproj
> > >>> pyproj.test()
> >
> > and let me know if any of the tests fail?
> >
> > -Jeff
> >
> > --
> > Jeffrey S. Whitaker Phone : (303)497-6313
> > Meteorologist FAX : (303)497-6449
> > NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
> > <mailto:Jef...@no...>
> > 325 Broadway Office : Skaggs Research Cntr 1D-124
> > Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
> >
> >
> David: Should be fixed now in svn. Let me know if it isn't.
>
> -Jeff
>
> --
> Jeffrey S. Whitaker Phone : (303)497-6313
> Meteorologist FAX : (303)497-6449
> NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
> 325 Broadway Office : Skaggs Research Cntr 1D-124
> Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
>
>

On Wed, Sep 12, 2007 at 01:11:54PM -0500, John Hunter wrote:
> On 9/12/07, Michael Droettboom <md...@st...> wrote:
> 
> > So, I feel like I'm going in a bit of a circle here, and I might need a
> > reality check. I thought I'd better check in and see where you guys
> > (who've thought about this a lot longer than I have) see this going. A
> > statement of objectives of this part of the task would be helpful.
> > (e.g. what's the biggest problem with how transforms work now, and what
> > model would be a better fit). John, I know you've mentioned some to me
> > before, e.g. the LazyValue concept is quirky and relies on C and the PDF
> > stateful transforms model is close, but not quite what we need, etc. I
> > feel I have a better sense of the overall code structure now, but you
> > guys may have a better "gut" sense of what will fit best.
> 
> Here is a brief summary of what I see some of the problems to be with
> the existing approach to transformations, and what I would like to see
> improved in a refactoring. The three major objectives are clarity,
> extensibility and efficiency.
> 
> Clarity:
> 
> The existing transformation framework, written in C++ and
> making extensive use of deferred evaluation of binary operation
> trees and values by reference, is difficult for most developers to
> understand (and hence enhance). Additionally, since all the heavy
> lifting is done in C++, python developers who are not versed in C++
> have an additional barrier to making contributions.
Indeed!
> Extensibilty:
> 
> We would like to make it fairly easy for users to add additional
> non-linear transformations. The current framework requires adding a
> new function at the C++ layer, and hacking into axes.py to support
> additional functions. We would like the existing nonlinear
> transformations (log and polar) to be part of a general
> infrastructure where users could supply their own nonlinear
> functions which map (possibly nonseparable) (xhat, yhat) ->
> separable (x, y). There are two parts to this: one pretty easy and
> one pretty hard.
> 
> The easy part is supporting a transformation which has a separation
> callable that takes, eg an Nx2 array and returns and Nx2 array. For
> log, this will simply be log(XY), for polar, it will be
> r*cos(X[:,0]), r*sin(X[:,1]). Presumably we will want to take
> advantage of masked arrays to support invalid transformations, eg
> log of nonpositive data.
> 
> The harder part is to support axis, tick and label layout
> generically. Currently we do this by special casing log and polar,
> either with special tick locators and formatters (log) or special
> derived Axes (polar).
Another hard part is grids. More generally, a straight line in 
x,y becomes curved in x',y'. Ideally, a sequence of points plotted
on a straight line should lie directly on the transformed line. This
would make the caps on the polar_bar demo follow the arcs of the grid.
The extreme case is map projections, where for some projections, a 
straight line will not even be connected.
Another issue is zooming and panning. For amusement, try it with 
polar_demo.
> Efficiency:	
> 
> There are three parts to the efficiency question: the efficiency of
> the transformation itself, the efficiency with which transformation
> data structures are updated in the presence of viewlim changes
> (panning and zooming, window resizing) and the efficiency in getting
> transformed data to the backends. My guess is that the new design
> may be slower or not dramatically faster for the first two (which
> are not the bottleneck in most cases anyhow) but you might get
> sigificant savings on the 3rd.
Changing the internal representation of things like collections so that 
the transform can be done using numpy vectors will help a lot.
	- Paul

David Huard wrote:
> Hi Jeff,
>
> Here is the output. Thanks for looking into this.
>
>
> In [1]: from matplotlib.toolkits.basemap import pyproj
>
> In [2]: pyproj.test()
> Trying:
> from pyproj import Geod
> Expecting nothing
> ok
> Trying:
> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
> Expecting nothing
> ok
> Trying:
> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
> Expecting nothing
> ok
> Trying:
> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
> Expecting nothing
> ok
> Trying:
> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
> Expecting nothing
> ok
> Trying:
> london_lat = 51.+(32./60.); london_lon = -(5./60.)
> Expecting nothing
> ok
> Trying:
> az12,az21,dist = 
> g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> Expecting nothing
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 401, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> az12,az21,dist = 
> g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> Exception raised:
> Traceback (most recent call last):
> File " doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[6]>", line 1, in ?
> az12,az21,dist = 
> g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 478, in inv
> _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> File "_geod.pyx", line 123, in _geod.Geod._inv
> ValueError: undefined inverse geodesic (may be an antipodal point)
> Trying:
> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> Expecting:
> -66.531 75.654 4164192.708
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 402, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest matplotlib.toolkits.basemap.pyproj.Geod.__new_ 
> _[7]>", line 1, in ?
> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> NameError: name 'az12' is not defined
> Trying:
> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> Expecting nothing
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 406, in matplotlib.toolkits.basemap.pyproj.Geod.__new_ _
> Failed example:
> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[8]>", line 1, in ?
> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> NameError: name 'az12' is not defined
> Trying:
> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> Expecting:
> 45.517 -123.683 75.654
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 407, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py ", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[9]>", line 1, in ?
> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> NameError: name 'endlat' is not defined
> Trying:
> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
> Expecting nothing
> ok
> Trying:
> lons2 = [boston_lon, portland_lon, london_lon]
> Expecting nothing
> ok
> Trying:
> lats2 = [boston_lat, portland_lat, london_lat]
> Expecting nothing
> ok
> Trying:
> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> Expecting nothing
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 414, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest matplotlib.toolkits.basemap.pyproj.Geod.__new_ 
> _[13]>", line 1, in ?
> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 478, in inv
> _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> File "_geod.pyx", line 123, in _geod.Geod._inv
> ValueError: undefined inverse geodesic (may be an antipodal point)
> Trying:
> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" % 
> (faz,baz,d)
> Expecting:
> 54.663 -123.448 288303.720
> -65.463 79.342 4013037.318
> 51.254 -71.576 5579916.649
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 415, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" % 
> (faz,baz,d)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[14]>", line 1, in ?
> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f 
> %9.3f" % (faz,baz,d)
> NameError: name 'az12' is not defined
> Trying:
> from pyproj import Geod
> Expecting nothing
> ok
> Trying:
> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
> Expecting nothing
> ok
> Trying:
> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
> Expecting nothing
> ok
> Trying:
> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
> Expecting nothing
> ok
> Trying:
> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
> Expecting nothing
> ok
> Trying:
> for lon,lat in lonlats: print '%6.3f %7.3f' % (lat, lon)
> Expecting:
> 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
> ok
> Trying:
> from pyproj import Proj
> Expecting nothing
> ok
> Trying:
> p = Proj(proj='utm',zone=10,ellps='WGS84')
> Expecting nothing
> ok
> Trying:
> x,y = p(-120.108, 34.36116666)
> Expecting nothing
> ok
> Trying:
> print 'x=%9.3f y=%11.3f' % (x,y)
> Expecting:
> x=765975.641 y=3805993.134
> ok
> Trying:
> print 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
> Expecting:
> lon=-120.108 lat=34.361
> ok
> Trying:
> lons = (-119.72,-118.40,-122.38)
> Expecting nothing
> ok
> Trying:
> lats = (36.77, 33.93, 37.62 )
> Expecting nothing
> ok
> Trying:
> x,y = p(lons, lats)
> Expecting nothing
> ok
> Trying:
> print 'x: %9.3f %9.3f %9.3f' % x
> Expecting:
> x: 792763.863 925321.537 554714.301
> ok
> Trying:
> print 'y: %9.3f %9.3f %9.3f' % y
> Expecting:
> y: 4074377.617 3763936.941 4163835.303
> ok
> Trying:
> lons, lats = p(x, y, inverse=True) # inverse transform
> Expecting nothing
> ok
> Trying:
> print 'lons: %8.3f %8.3f %8.3f' % lons
> Expecting:
> lons: -119.720 -118.400 -122.380
> ok
> Trying:
> print 'lats: %8.3f %8.3f %8.3f' % lats
> Expecting:
> lats: 36.770 33.930 37.620
> ok
> Trying:
> p1 = Proj(init='epsg:26915')
> Expecting nothing
> ok
> Trying:
> p2 = Proj(init='epsg:26715')
> Expecting nothing
> ok
> Trying:
> x1, y1 = p1(-92.199881,38.56694)
> Expecting nothing
> ok
> Trying:
> x2, y2 = transform(p1,p2,x1,y1)
> Expecting nothing
> ok
> Trying:
> print '%9.3f %11.3f' % (x1,y1)
> Expecting:
> 569704.566 4269024.671
> ok
> Trying:
> print '%9.3f %11.3f' % (x2,y2)
> Expecting:
> 569706.333 4268817.680
> ok
> Trying:
> print '%8.3f %5.3f' % p2(x2,y2,inverse=True)
> Expecting:
> -92.200 38.567
> ok
> Trying:
> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
> Expecting nothing
> ok
> Trying:
> lons = (-92.22,-94.72,-90.37)
> Expecting nothing
> ok
> Trying:
> x1, y1 = p1(lons,lats)
> Expecting nothing
> ok
> Trying:
> x2, y2 = transform(p1,p2,x1,y1)
> Expecting nothing
> ok
> Trying:
> xy = x1+y1
> Expecting nothing
> ok
> Trying:
> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
> Expecting:
> 567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005
> ok
> Trying:
> xy = x2+y2
> Expecting nothing
> ok
> Trying:
> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
> Expecting:
> 567705.072 351727.113 728558.917 4297993.157 4353490.111 4292111.678
> ok
> Trying:
> lons, lats = p2(x2,y2,inverse=True)
> Expecting nothing
> ok
> Trying:
> xy = lons+lats
> Expecting nothing
> ok
> Trying:
> print '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy
> Expecting:
> -92.220 -94.720 -90.370 38.830 39.320 38.750
> ok
> 11 items had no tests:
> matplotlib.toolkits.basemap.pyproj
> matplotlib.toolkits.basemap.pyproj.Geod
> matplotlib.toolkits.basemap.pyproj.Geod.fwd
> matplotlib.toolkits.basemap.pyproj.Geod.inv
> matplotlib.toolkits.basemap.pyproj.Proj
> matplotlib.toolkits.basemap.pyproj.Proj.__call__
> matplotlib.toolkits.basemap.pyproj.Proj.is_geocent
> matplotlib.toolkits.basemap.pyproj.Proj.is_latlong
> matplotlib.toolkits.basemap.pyproj._convertback
> matplotlib.toolkits.basemap.pyproj._copytobuffer
> matplotlib.toolkits.basemap.pyproj.test
> 3 items passed all tests:
> 6 tests in matplotlib.toolkits.basemap.pyproj.Geod.npts
> 13 tests in matplotlib.toolkits.basemap.pyproj.Proj.__new__
> 18 tests in matplotlib.toolkits.basemap.pyproj.transform
> **********************************************************************
> 1 items had failures:
> 6 of 15 in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> 52 tests in 15 items.
> 46 passed and 6 failed.
> ***Test Failed*** 6 failures.
>
>
> 2007年9月12日, Jeff Whitaker < js...@fa... 
> <mailto:js...@fa...>>:
>
> David Huard wrote:
> > Hi, the pyproj package seems to cause a problem in the polarmap
> > example of the basemap toolkit.
> >
> > Thanks,
> >
> > david
> >
> > huardda@allard:~/svnrepos/toolkits/basemap/examples$ python
> polarmaps.py
> > min/max etopo20 data:
> > -9026.625 6228.8125
> > plotting North Polar Lambert Azimuthal Equal Area basemap ...
> > plotting North Polar Stereographic basemap ...
> > plotting North Polar Azimuthal Equidistant basemap ...
> > Traceback (most recent call last):
> > File "polarmaps.py", line 51, in ?
> > resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/basemap.py",
>
> > line 1111, in __init__
> >
> >
> az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[0]),math.radians(latsnew[0]),radians=True)
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
>
> > line 478, in inv
> > _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> > File "_geod.pyx", line 123, in _geod.Geod._inv
> > ValueError: undefined inverse geodesic (may be an antipodal point)
> >
> >
> >
> > This is from a fresh SVN version of both matplotlib and basemap.
> > Linux, Ubuntu edgy, Xeon-64.
>
> David: Odd - I can't reproduce that on my mac. Can you try this
>
> >>> from matplotlib.toolkits.basemap import pyproj
> >>> pyproj.test()
>
> and let me know if any of the tests fail?
>
> -Jeff
>
> --
> Jeffrey S. Whitaker Phone : (303)497-6313
> Meteorologist FAX : (303)497-6449
> NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
> <mailto:Jef...@no...>
> 325 Broadway Office : Skaggs Research Cntr 1D-124
> Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
>
>
David: Should be fixed now in svn. Let me know if it isn't.
-Jeff
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

On 9/12/07, Michael Droettboom <md...@st...> wrote:
> So, I feel like I'm going in a bit of a circle here, and I might need a
> reality check. I thought I'd better check in and see where you guys
> (who've thought about this a lot longer than I have) see this going. A
> statement of objectives of this part of the task would be helpful.
> (e.g. what's the biggest problem with how transforms work now, and what
> model would be a better fit). John, I know you've mentioned some to me
> before, e.g. the LazyValue concept is quirky and relies on C and the PDF
> stateful transforms model is close, but not quite what we need, etc. I
> feel I have a better sense of the overall code structure now, but you
> guys may have a better "gut" sense of what will fit best.
Here is a brief summary of what I see some of the problems to be with
the existing approach to transformations, and what I would like to see
improved in a refactoring. The three major objectives are clarity,
extensibility and efficiency.
Clarity:
 The existing transformation framework, written in C++ and
 making extensive use of deferred evaluation of binary operation
 trees and values by reference, is difficult for most developers to
 understand (and hence enhance). Additionally, since all the heavy
 lifting is done in C++, python developers who are not versed in C++
 have an additional barrier to making contributions.
Extensibilty:
 We would like to make it fairly easy for users to add additional
 non-linear transformations. The current framework requires adding a
 new function at the C++ layer, and hacking into axes.py to support
 additional functions. We would like the existing nonlinear
 transformations (log and polar) to be part of a general
 infrastructure where users could supply their own nonlinear
 functions which map (possibly nonseparable) (xhat, yhat) ->
 separable (x, y). There are two parts to this: one pretty easy and
 one pretty hard.
 The easy part is supporting a transformation which has a separation
 callable that takes, eg an Nx2 array and returns and Nx2 array. For
 log, this will simply be log(XY), for polar, it will be
 r*cos(X[:,0]), r*sin(X[:,1]). Presumably we will want to take
 advantage of masked arrays to support invalid transformations, eg
 log of nonpositive data.
 The harder part is to support axis, tick and label layout
 generically. Currently we do this by special casing log and polar,
 either with special tick locators and formatters (log) or special
 derived Axes (polar).
Efficiency:	
 There are three parts to the efficiency question: the efficiency of
 the transformation itself, the efficiency with which transformation
 data structures are updated in the presence of viewlim changes
 (panning and zooming, window resizing) and the efficiency in getting
 transformed data to the backends. My guess is that the new design
 may be slower or not dramatically faster for the first two (which
 are not the bottleneck in most cases anyhow) but you might get
 sigificant savings on the 3rd.
 What we would like to support is something like an operation which
 pushes the partially transformed data to the backend, and the
 backend then stores this data in a path or other data structure, and
 when the viewlimits are changed or the window is resized, the
 backend merely needs to get an updated affine to redraw the data. I
 say "partially transformed" because in the case of nonlinear or
 separable transformations, we will probably want to do the
 nonlinear/separation part first, and then push this to the backend
 which can build a path (eg an agg::path in agg) and on pan and zoom
 we would only need to let the backend know what the current affine
 is. There is more than one way to solve this problem: in mpl1 I
 used a path dictionary keyed off of a path id which was updated on
 renderer changes.
 Then the front end (eg Line2D) could do something like
 def on_renderer_change(self, renderer):
 # on renderer change; path data already has the
 # separable/nonlinear part handled
 self.pathid = renderer.push_path(pathdata)
 Additionally, you would need to track when either the data or
 nonlinear mapping function are changed in order to remove the old
 path and push out a new one. Then at draw time, you would not need
 to push any data to the backend
 def on_draw(self, renderer):
 # on draw the line just needs to inform the backend to draw the
 # cached path with the current separable transformation
 renderer.draw_path(self.pathid, gc, septrans)
 Ken I believe used some meta class magic to solve this problem. It
 may be that the ideal approach is different from either of these so
 feel free to experiment and I'll give it some more thought too. Ken
 will hopefully pipe in with his perspective too.

David Huard wrote:
> Hi Jeff,
>
> Here is the output. Thanks for looking into this.
>
>
> In [1]: from matplotlib.toolkits.basemap import pyproj
>
> In [2]: pyproj.test()
> Trying:
> from pyproj import Geod
> Expecting nothing
> ok
> Trying:
> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
> Expecting nothing
> ok
> Trying:
> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
> Expecting nothing
> ok
> Trying:
> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
> Expecting nothing
> ok
> Trying:
> newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
> Expecting nothing
> ok
> Trying:
> london_lat = 51.+(32./60.); london_lon = -(5./60.)
> Expecting nothing
> ok
> Trying:
> az12,az21,dist = 
> g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> Expecting nothing
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 401, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> az12,az21,dist = 
> g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> Exception raised:
> Traceback (most recent call last):
> File " doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[6]>", line 1, in ?
> az12,az21,dist = 
> g.inv(boston_lon,boston_lat,portland_lon,portland_lat)
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 478, in inv
> _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> File "_geod.pyx", line 123, in _geod.Geod._inv
> ValueError: undefined inverse geodesic (may be an antipodal point)
> Trying:
> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> Expecting:
> -66.531 75.654 4164192.708
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 402, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest matplotlib.toolkits.basemap.pyproj.Geod.__new_ 
> _[7]>", line 1, in ?
> print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
> NameError: name 'az12' is not defined
> Trying:
> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> Expecting nothing
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 406, in matplotlib.toolkits.basemap.pyproj.Geod.__new_ _
> Failed example:
> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[8]>", line 1, in ?
> endlon, endlat, backaz = g.fwd(boston_lon, boston_lat, az12, dist)
> NameError: name 'az12' is not defined
> Trying:
> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> Expecting:
> 45.517 -123.683 75.654
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 407, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py ", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[9]>", line 1, in ?
> print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
> NameError: name 'endlat' is not defined
> Trying:
> lons1 = 3*[newyork_lon]; lats1 = 3*[newyork_lat]
> Expecting nothing
> ok
> Trying:
> lons2 = [boston_lon, portland_lon, london_lon]
> Expecting nothing
> ok
> Trying:
> lats2 = [boston_lat, portland_lat, london_lat]
> Expecting nothing
> ok
> Trying:
> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> Expecting nothing
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 414, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest matplotlib.toolkits.basemap.pyproj.Geod.__new_ 
> _[13]>", line 1, in ?
> az12,az21,dist = g.inv(lons1,lats1,lons2,lats2)
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 478, in inv
> _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> File "_geod.pyx", line 123, in _geod.Geod._inv
> ValueError: undefined inverse geodesic (may be an antipodal point)
> Trying:
> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" % 
> (faz,baz,d)
> Expecting:
> 54.663 -123.448 288303.720
> -65.463 79.342 4013037.318
> 51.254 -71.576 5579916.649
> **********************************************************************
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 415, in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> Failed example:
> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f %9.3f" % 
> (faz,baz,d)
> Exception raised:
> Traceback (most recent call last):
> File "doctest.py", line 1248, in __run
> compileflags, 1) in test.globs
> File "<doctest 
> matplotlib.toolkits.basemap.pyproj.Geod.__new__[14]>", line 1, in ?
> for faz,baz,d in zip(az12,az21,dist): print "%7.3f %7.3f 
> %9.3f" % (faz,baz,d)
> NameError: name 'az12' is not defined
> Trying:
> from pyproj import Geod
> Expecting nothing
> ok
> Trying:
> g = Geod(ellps='clrk66') # Use Clarke 1966 ellipsoid.
> Expecting nothing
> ok
> Trying:
> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
> Expecting nothing
> ok
> Trying:
> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
> Expecting nothing
> ok
> Trying:
> lonlats = g.npts(boston_lon,boston_lat,portland_lon,portland_lat,10)
> Expecting nothing
> ok
> Trying:
> for lon,lat in lonlats: print '%6.3f %7.3f' % (lat, lon)
> Expecting:
> 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
> ok
> Trying:
> from pyproj import Proj
> Expecting nothing
> ok
> Trying:
> p = Proj(proj='utm',zone=10,ellps='WGS84')
> Expecting nothing
> ok
> Trying:
> x,y = p(-120.108, 34.36116666)
> Expecting nothing
> ok
> Trying:
> print 'x=%9.3f y=%11.3f' % (x,y)
> Expecting:
> x=765975.641 y=3805993.134
> ok
> Trying:
> print 'lon=%8.3f lat=%5.3f' % p(x,y,inverse=True)
> Expecting:
> lon=-120.108 lat=34.361
> ok
> Trying:
> lons = (-119.72,-118.40,-122.38)
> Expecting nothing
> ok
> Trying:
> lats = (36.77, 33.93, 37.62 )
> Expecting nothing
> ok
> Trying:
> x,y = p(lons, lats)
> Expecting nothing
> ok
> Trying:
> print 'x: %9.3f %9.3f %9.3f' % x
> Expecting:
> x: 792763.863 925321.537 554714.301
> ok
> Trying:
> print 'y: %9.3f %9.3f %9.3f' % y
> Expecting:
> y: 4074377.617 3763936.941 4163835.303
> ok
> Trying:
> lons, lats = p(x, y, inverse=True) # inverse transform
> Expecting nothing
> ok
> Trying:
> print 'lons: %8.3f %8.3f %8.3f' % lons
> Expecting:
> lons: -119.720 -118.400 -122.380
> ok
> Trying:
> print 'lats: %8.3f %8.3f %8.3f' % lats
> Expecting:
> lats: 36.770 33.930 37.620
> ok
> Trying:
> p1 = Proj(init='epsg:26915')
> Expecting nothing
> ok
> Trying:
> p2 = Proj(init='epsg:26715')
> Expecting nothing
> ok
> Trying:
> x1, y1 = p1(-92.199881,38.56694)
> Expecting nothing
> ok
> Trying:
> x2, y2 = transform(p1,p2,x1,y1)
> Expecting nothing
> ok
> Trying:
> print '%9.3f %11.3f' % (x1,y1)
> Expecting:
> 569704.566 4269024.671
> ok
> Trying:
> print '%9.3f %11.3f' % (x2,y2)
> Expecting:
> 569706.333 4268817.680
> ok
> Trying:
> print '%8.3f %5.3f' % p2(x2,y2,inverse=True)
> Expecting:
> -92.200 38.567
> ok
> Trying:
> lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
> Expecting nothing
> ok
> Trying:
> lons = (-92.22,-94.72,-90.37)
> Expecting nothing
> ok
> Trying:
> x1, y1 = p1(lons,lats)
> Expecting nothing
> ok
> Trying:
> x2, y2 = transform(p1,p2,x1,y1)
> Expecting nothing
> ok
> Trying:
> xy = x1+y1
> Expecting nothing
> ok
> Trying:
> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
> Expecting:
> 567703.344 351730.944 728553.093 4298200.739 4353698.725 4292319.005
> ok
> Trying:
> xy = x2+y2
> Expecting nothing
> ok
> Trying:
> print '%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy
> Expecting:
> 567705.072 351727.113 728558.917 4297993.157 4353490.111 4292111.678
> ok
> Trying:
> lons, lats = p2(x2,y2,inverse=True)
> Expecting nothing
> ok
> Trying:
> xy = lons+lats
> Expecting nothing
> ok
> Trying:
> print '%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy
> Expecting:
> -92.220 -94.720 -90.370 38.830 39.320 38.750
> ok
> 11 items had no tests:
> matplotlib.toolkits.basemap.pyproj
> matplotlib.toolkits.basemap.pyproj.Geod
> matplotlib.toolkits.basemap.pyproj.Geod.fwd
> matplotlib.toolkits.basemap.pyproj.Geod.inv
> matplotlib.toolkits.basemap.pyproj.Proj
> matplotlib.toolkits.basemap.pyproj.Proj.__call__
> matplotlib.toolkits.basemap.pyproj.Proj.is_geocent
> matplotlib.toolkits.basemap.pyproj.Proj.is_latlong
> matplotlib.toolkits.basemap.pyproj._convertback
> matplotlib.toolkits.basemap.pyproj._copytobuffer
> matplotlib.toolkits.basemap.pyproj.test
> 3 items passed all tests:
> 6 tests in matplotlib.toolkits.basemap.pyproj.Geod.npts
> 13 tests in matplotlib.toolkits.basemap.pyproj.Proj.__new__
> 18 tests in matplotlib.toolkits.basemap.pyproj.transform
> **********************************************************************
> 1 items had failures:
> 6 of 15 in matplotlib.toolkits.basemap.pyproj.Geod.__new__
> 52 tests in 15 items.
> 46 passed and 6 failed.
> ***Test Failed*** 6 failures.
>
>
> 2007年9月12日, Jeff Whitaker < js...@fa... 
> <mailto:js...@fa...>>:
>
> David Huard wrote:
> > Hi, the pyproj package seems to cause a problem in the polarmap
> > example of the basemap toolkit.
> >
> > Thanks,
> >
> > david
> >
> > huardda@allard:~/svnrepos/toolkits/basemap/examples$ python
> polarmaps.py
> > min/max etopo20 data:
> > -9026.625 6228.8125
> > plotting North Polar Lambert Azimuthal Equal Area basemap ...
> > plotting North Polar Stereographic basemap ...
> > plotting North Polar Azimuthal Equidistant basemap ...
> > Traceback (most recent call last):
> > File "polarmaps.py", line 51, in ?
> > resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/basemap.py",
>
> > line 1111, in __init__
> >
> >
> az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[0]),math.radians(latsnew[0]),radians=True)
> > File
> >
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
>
> > line 478, in inv
> > _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> > File "_geod.pyx", line 123, in _geod.Geod._inv
> > ValueError: undefined inverse geodesic (may be an antipodal point)
> >
> >
> >
> > This is from a fresh SVN version of both matplotlib and basemap.
> > Linux, Ubuntu edgy, Xeon-64.
>
> David: Odd - I can't reproduce that on my mac. Can you try this
>
> >>> from matplotlib.toolkits.basemap import pyproj
> >>> pyproj.test()
>
> and let me know if any of the tests fail?
>
> -Jeff
>
> --
> Jeffrey S. Whitaker Phone : (303)497-6313
> Meteorologist FAX : (303)497-6449
> NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
> <mailto:Jef...@no...>
> 325 Broadway Office : Skaggs Research Cntr 1D-124
> Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
>
>
David: It might be a 64-bit issue. I'll try to find a 64-bit Xeon box 
to debug it, and let you know what I find.
-Jeff
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

On 9/12/07, Michael Droettboom <md...@st...> wrote:
> Thanks for taking the time. Very helpful (and please excuse the mess in
> the code -- I was just trying to get something end-to-end working before
> refining/optimizing/documenting etc...)
I think this is definitely the right approach -- get something that
works in the existing framework and understand where the various
issues are and then try and peel away the stuff that is not ideal. I
looked at the tick labels -- if you just comment out the
transformation offset
trans = trans + Affine2D().translate(0, -1 * self._padPixels)
the tick labels show up too (minus the pad of course), so my guess is
some reference is being lost in the addition....

John Hunter wrote:
> On 9/12/07, Michael Droettboom <md...@st...> wrote:
> I need to spend more time reading through your code before I comment
> further, but I just wanted to make a quick comment vis-a-vis the
> locators and formatters. I commited these changes to your branch, and
> autoscaling is now working there :-) I'll keep poking and learning
> more about what you are doing before commenting on some of your bigger
> questions.
> 
> I made a couple of comments in affine.py as well, prefixed by 'JDH'
Thanks for taking the time. Very helpful (and please excuse the mess in 
the code -- I was just trying to get something end-to-end working before 
refining/optimizing/documenting etc...)
Cheers,
Mike
-- 
Michael Droettboom
Operations and Engineering Division
Space Telescope Science Institute
Operated by AURA for NASA

On 9/12/07, Michael Droettboom <md...@st...> wrote:
> This approach got closer, until I hit the wall that dependencies work at
> an even lower level -- single lazy values get borrowed from one bounding
> box and referenced in another (e.g. Axes.autoscale_view()) Certainly,
> this could be implemented in my new affine-based framework, but then
> we're almost back to square one and have basically re-implemented
> transforms.py/.cpp into something that is probably slower -- though
> perhaps more flexible in that more kinds of transforms could be added
> using only Python. Of course, autoscale_view() (and other instances of
> this) could be rewritten to work differently, but it's hard to know
> where that might end.
The locators do have a reference to the datalim and viewlim intervals,
which is is what they use to compute their autoscale limits and tick
locations, but they return scalars, and autoscale_view simply sets the
new limits with these scalars. So the fact that there is a reference
here is easy to work around. I made a minor change in your code
(ticker.py and axis.py) to illustrate. Instead of relying on the
Interval to pass information from the Axis -> Locator/Formatter, I
simply set the axis instance instead. Then, eg, the Locator can do
 vmin, vmax = self.axis.get_view_interval()
 dmin, dmax = self.axis.get_data_interval()
so there are no confusing intertwined references to deal with, and the
axis can be responsible for knowings its data and view limits, which
seem reasonable. I made these changes just to the MaxNLocator and
ScalarFormatter classes for proof of concept, but it should be trivial
to port to the others. I think in general communicating by scalar
values passed explicitly or through callbacks will make for clearer
code than the deeply nested references we have been using the in the
existingcode.
There are places where one bounding box value is shared with another
(most clearly in sharex and sharey support, eg
 left = self._sharex.viewLim.xmin()
The ability to "share" and axis, eg so changes in pan and zoom on one
are reflected in another, is extremely useful, but a better approach
may be to use callbacks (or something like them) rather than shared,
composited transforms which are updated in place.
I need to spend more time reading through your code before I comment
further, but I just wanted to make a quick comment vis-a-vis the
locators and formatters. I commited these changes to your branch, and
autoscaling is now working there :-) I'll keep poking and learning
more about what you are doing before commenting on some of your bigger
questions.
I made a couple of comments in affine.py as well, prefixed by 'JDH'
JDH

David Huard wrote:
> Hi, the pyproj package seems to cause a problem in the polarmap 
> example of the basemap toolkit.
>
> Thanks,
>
> david
>
> huardda@allard:~/svnrepos/toolkits/basemap/examples$ python polarmaps.py
> min/max etopo20 data:
> -9026.625 6228.8125
> plotting North Polar Lambert Azimuthal Equal Area basemap ...
> plotting North Polar Stereographic basemap ...
> plotting North Polar Azimuthal Equidistant basemap ...
> Traceback (most recent call last):
> File "polarmaps.py", line 51, in ?
> resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/basemap.py", 
> line 1111, in __init__
> 
> az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[0]),math.radians(latsnew[0]),radians=True)
> File 
> "/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py", 
> line 478, in inv
> _Geod._inv(self, inx, iny, inz, ind, radians=radians)
> File "_geod.pyx", line 123, in _geod.Geod._inv
> ValueError: undefined inverse geodesic (may be an antipodal point)
>
>
>
> This is from a fresh SVN version of both matplotlib and basemap. 
> Linux, Ubuntu edgy, Xeon-64.
David: Odd - I can't reproduce that on my mac. Can you try this
 >>> from matplotlib.toolkits.basemap import pyproj
 >>> pyproj.test()
and let me know if any of the tests fail?
-Jeff
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Hi, the pyproj package seems to cause a problem in the polarmap example of
the basemap toolkit.
Thanks,
david
huardda@allard:~/svnrepos/toolkits/basemap/examples$ python polarmaps.py
min/max etopo20 data:
-9026.625 6228.8125
plotting North Polar Lambert Azimuthal Equal Area basemap ...
plotting North Polar Stereographic basemap ...
plotting North Polar Azimuthal Equidistant basemap ...
Traceback (most recent call last):
 File "polarmaps.py", line 51, in ?
 resolution='c',area_thresh=10000.,lat_0=lat_0,lon_0=lon_0_ortho)
 File
"/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/basemap.py",
line 1111, in __init__
 az1,alpha21,dist=gc.inv(lon_0,lat_0,math.radians(lonsnew[0]),
math.radians(latsnew[0]),radians=True)
 File
"/usr/local/lib/python2.4/site-packages/matplotlib/toolkits/basemap/pyproj.py",
line 478, in inv
 _Geod._inv(self, inx, iny, inz, ind, radians=radians)
 File "_geod.pyx", line 123, in _geod.Geod._inv
ValueError: undefined inverse geodesic (may be an antipodal point)
This is from a fresh SVN version of both matplotlib and basemap. Linux,
Ubuntu edgy, Xeon-64.

I should also add -- it's only working with the Agg backend.
John Hunter wrote:
> On 9/12/07, Michael Droettboom <md...@st...> wrote:
> 
>> If you check out r3835 from my branch, simple_plot.py is working, with
>> the exception of things that rely on this really low-level
>> interdependence, e.g. the data limits.)
> 
> I am at 3836 in the transforms branch, but I do not see "pbox".
> Perhaps you forgot to svn add it?
> 
> JDH
-- 
Michael Droettboom
Operations and Engineering Division
Space Telescope Science Institute
Operated by AURA for NASA

Yes. Sorry. It's in r3837 on the branch.
Cheers,
Mike
John Hunter wrote:
> On 9/12/07, Michael Droettboom <md...@st...> wrote:
> 
>> If you check out r3835 from my branch, simple_plot.py is working, with
>> the exception of things that rely on this really low-level
>> interdependence, e.g. the data limits.)
> 
> I am at 3836 in the transforms branch, but I do not see "pbox".
> Perhaps you forgot to svn add it?
> 
> JDH
-- 
Michael Droettboom
Operations and Engineering Division
Space Telescope Science Institute
Operated by AURA for NASA

On 9/12/07, Michael Droettboom <md...@st...> wrote:
> If you check out r3835 from my branch, simple_plot.py is working, with
> the exception of things that rely on this really low-level
> interdependence, e.g. the data limits.)
I am at 3836 in the transforms branch, but I do not see "pbox".
Perhaps you forgot to svn add it?
JDH

[Background: I'm working on refactoring the transforms framework with 
the end goal of making it easier to add new kinds of non-linear 
transforms and projections to matplotlib. I've been talking a bit with 
John Hunter about this -- this question is mainly for John and Ken 
McIvor, though there are probably some other interested parties on this 
list as well.]
I've studied John's mpl1.py and Ken's mpl1_displaypdf.py to try to get a 
sense of where things could go. I appreciate the ideas both of these 
present as clean slates -- however, I think what I'm running into is 
"how to get there from here" in manageable steps.
My first baby step in this large task has been to try to remove 
transforms.py/cpp and replace it with something based on standard 3x3 
affine matrices, using Python/numpy only. The way transforms.py/.cpp 
works now, everything is built around live updates of a tree of 
interdependent bounding boxes and transforms, where a change to a single 
scalar in any object automatically propagates through the tree.
My first thought was to make something out of immutable transforms where 
a transform would have to be calculated from its dependencies 
immediately before drawing, and therefore get rid of these "magical" 
side-effects by not allowing transforms to change in place. Reading 
between the lines, this seems to be what mpl1_displaypdf.py suggests. I 
quickly came to the conclusion that that is perhaps a step too far -- 
matplotlib is very much built around these side-effects and I would hate 
to replace hundreds of lines of well-tested code. On the other hand, 
there is probably a pattern to those changes, and it may be worth the 
effort if others agree it's useful.
My second kick at the can was to build a live-updating tree of 
transforms. This is similar to what I saw in mpl1.py using "changed" 
callbacks so that a change in one transform would affect all transforms 
that depend on it.
[I worry about a pure callback approach because of the likelihood of 
computing many partial values. For example, if 'a' depends on 'b' and 
'c', and I change 'b' then 'c', 'a' will get recomputed twice. Instead, 
I used an "invalidation" technique, where a change in b simply 
invalidates a, and a doesn't get recomputed until it is later requested. 
 This is something we used a lot when I programmed for gaming hardware. 
 The resulting semantics are very similar to using callbacks, however.]
This approach got closer, until I hit the wall that dependencies work at 
an even lower level -- single lazy values get borrowed from one bounding 
box and referenced in another (e.g. Axes.autoscale_view()) Certainly, 
this could be implemented in my new affine-based framework, but then 
we're almost back to square one and have basically re-implemented 
transforms.py/.cpp into something that is probably slower -- though 
perhaps more flexible in that more kinds of transforms could be added 
using only Python. Of course, autoscale_view() (and other instances of 
this) could be rewritten to work differently, but it's hard to know 
where that might end.
(You can see my semi-working sketch of this here:
http://matplotlib.svn.sourceforge.net/viewvc/matplotlib/branches/transforms/lib/matplotlib/affine.py?revision=3835&view=markup
If you check out r3835 from my branch, simple_plot.py is working, with 
the exception of things that rely on this really low-level 
interdependence, e.g. the data limits.)
So, I feel like I'm going in a bit of a circle here, and I might need a 
reality check. I thought I'd better check in and see where you guys 
(who've thought about this a lot longer than I have) see this going. A 
statement of objectives of this part of the task would be helpful. 
(e.g. what's the biggest problem with how transforms work now, and what 
model would be a better fit). John, I know you've mentioned some to me 
before, e.g. the LazyValue concept is quirky and relies on C and the PDF 
stateful transforms model is close, but not quite what we need, etc. I 
feel I have a better sense of the overall code structure now, but you 
guys may have a better "gut" sense of what will fit best.
My next planned step, to move more (affine) transformations to the 
backends to allow the same path data to be transformed in the backend 
without retransmitting/converting the path data each time, doesn't 
actually seem to be dependent on getting the above done. (The existing 
transforms.cpp code already has a way to get a representation as an 
affine matrix). John, I have a note from our phone conversation that 
indicates you thought these two things would be dependent, but I don't 
remember the reason you gave -- could you maybe refresh my memory -- I 
was much less aware of the code structure then.
Sorry if this is sort of an open-ended question... Any pointers or 
impressions, no matter how small, are appreciated.
Cheers,
Mike
-- 
Michael Droettboom
Operations and Engineering Division
Space Telescope Science Institute
Operated by AURA for NASA