matplotlib-users Mailing List for matplotlib

On Tue, 13 Jul 2004, John Hunter wrote:

John:

>
> OK, made some headway here.
>
> When profiling for performace, I often turn off the GUI so my numbers
> aren't influenced by the user interface
>
> > /usr/local/lib/python/profile.py plotmap.py -dAgg > prof.out
>
> Here are the profile results sorted by cumulative time for your script
> - only the most expensive functions are shown
>
>   7.590  proj.py:43(__call__)
>   7.080  os.py:611(popen2)
>   7.050  popen2.py:31(__init__)
>   7.050  popen2.py:143(popen2)
>   2.610  matlab.py:1305(savefig)
>
> On my system, over half the time is spent running your popen process.
> Can't help you with this one :-).  To get better numbers and focus on
> matplotlib, I cached the xy points that the popen process is
> generating to a pickled list.

OK - looks like I have to work on the proj module. I probably should make 
it a C extension that calls the proj library directly instead of 
communicating with the proj application via pipes.

>
> You're using pcolor to generate the colormap.  pcolor uses
> collections, which is already pretty fast (agg implements collection
> drawing in extension code, no python loops needed to make a pcolor).
> However, you can get better performace still from imshow with the data
> limits set.  This will approximately double the performance of the
> image part of the map.  I had to make some changes to matplotlib (see
> below) because image origin wasn't playing nicely with image extent -
> this also fixes Andrew's bug.
>
> For the line part of the map, I extended the
> matplotlib.collections.LineCollection class to handle a sequence of
> lines, where each line is defined by a list of tuples (x0,y0), (x1,
> y1), ...  Thus all of your lines are handled by a single object,
> rather than having 1800+ separate line objects created in plot.
> Again, no python loops required.
>
> In the current form, the code takes about 1.15s to run on my machine
> and is about 30x faster than the original code you posted which
> includes the data loading part.  Nonetheless, the matplotlib part is
> much faster too, as you'll see when you interact with the data.
>
> I'm including a link to a matplotlib snapshot below which includes the
> required changes.  As a bonus, you can try out the new navigation
> toolbar (in progress, only works w/ gtk and gtkagg).  It includes a
> view limits stack, hand pan, and zoom to rectangle.  Much nicer for
> map navigation.  And with the changes, you can actually interact with
> your data with reasonable performance.  I need to add some more
> features to the toolbar, but give it a test drive and let me know if
> you have suggestions.  Set 'toolbar: toolbar2' in matplotlibrc

Excellent!  I'll do some more work on it and repost a (hopefully) faster 
version.

Thanks for the quick response.

-Jeff


-- 
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/CDC  R/CDC1        Email  : Jef...@no...
325 Broadway                Web    : www.cdc.noaa.gov/~jsw
Boulder, CO, USA 80303-3328 Office : Skaggs Research Cntr 1D-124

OK, made some headway here.

When profiling for performace, I often turn off the GUI so my numbers
aren't influenced by the user interface

 > /usr/local/lib/python/profile.py plotmap.py -dAgg > prof.out
 
Here are the profile results sorted by cumulative time for your script
- only the most expensive functions are shown

   7.590  proj.py:43(__call__)
   7.080  os.py:611(popen2)
   7.050  popen2.py:31(__init__)
   7.050  popen2.py:143(popen2)
   2.610  matlab.py:1305(savefig)

On my system, over half the time is spent running your popen process.
Can't help you with this one :-).  To get better numbers and focus on
matplotlib, I cached the xy points that the popen process is
generating to a pickled list.

You're using pcolor to generate the colormap.  pcolor uses
collections, which is already pretty fast (agg implements collection
drawing in extension code, no python loops needed to make a pcolor).
However, you can get better performace still from imshow with the data
limits set.  This will approximately double the performance of the
image part of the map.  I had to make some changes to matplotlib (see
below) because image origin wasn't playing nicely with image extent -
this also fixes Andrew's bug.

For the line part of the map, I extended the
matplotlib.collections.LineCollection class to handle a sequence of
lines, where each line is defined by a list of tuples (x0,y0), (x1,
y1), ...  Thus all of your lines are handled by a single object,
rather than having 1800+ separate line objects created in plot.
Again, no python loops required.

In the current form, the code takes about 1.15s to run on my machine
and is about 30x faster than the original code you posted which
includes the data loading part.  Nonetheless, the matplotlib part is
much faster too, as you'll see when you interact with the data.

I'm including a link to a matplotlib snapshot below which includes the
required changes.  As a bonus, you can try out the new navigation
toolbar (in progress, only works w/ gtk and gtkagg).  It includes a
view limits stack, hand pan, and zoom to rectangle.  Much nicer for
map navigation.  And with the changes, you can actually interact with
your data with reasonable performance.  I need to add some more
features to the toolbar, but give it a test drive and let me know if
you have suggestions.  Set 'toolbar: toolbar2' in matplotlibrc

Thanks much for the fink package - I'll continue to point OS X users
to your site!

JDH

### matplotlib shapshot:
    http://nitace.bsd.uchicago.edu:8080/files/share/matplotlib-0.60.3a.tar.gz

### Here is a link to the xy point data:
    http://nitace.bsd.uchicago.edu:8080/files/share/xypts.pickled 

### Here is the code snippet I used to generate it

xypts = []
for line in wcl:
    splitline = line.split()
    if splitline[0] == '#':
        segnum = segnum + 1
	if segnum > 1:
            # convert lon,lat to map coordinates x,y
            xys = zip(*proj(Numeric.array(lons),Numeric.array(lats)))
            xypts.append( xys)
pickle.dump(xypts, file('xypts.pickled', 'w') )


### Here is the modified plotmap script which uses imshow and line
### collections

import time
from matplotlib.matlab import *
from matplotlib.collections import LineCollection
from proj import Proj
import Numeric, cPickle

# set the default params for imshow
rc('image', origin='lower', cmap='jet')

ax = subplot(111)
nx = 349; ny = 277
dx = 32463.41; dy = 32463.41
xmax = (nx-1)*dx; ymax = (ny-1)*dy # size of domain to plot

datin = open('topodata.pickle','rb')
C = cPickle.load(datin)
# use imshow rather than pcolor for speed
im = ax.imshow(C, interpolation='nearest',
               extent=(0, xmax, 0, ymax))

xypts = cPickle.load(file('xypts.pickled', 'r') )
# all args are sequences, length 1 in case of linewidths and
# antialiased
collection = LineCollection(segments     = xypts,
                            colors       = ( (0,0,0,1), ), # black
                            linewidths   = (1.5,),
                            antialiaseds = (0,),  # turn off aa for speed
                            )

ax.add_collection(collection)

# you have to manually update the datalim; this is a bit ugly so I'll
# work on the interface
xs = [ x for thisline in xypts for x,y in thisline ]
ys = [ y for thisline in xypts for x,y in thisline ]
minx, maxx = min(xs), max(xs)
miny, maxy = min(ys), max(ys)
ax.update_datalim( ((minx, miny), (maxx, maxy)))


axis([0, xmax, 0, ymax])
ax.set_xticks([]) # no ticks
ax.set_yticks([])
title('40 km Topography - Lambert Conformal Conic Projection')
#savefig('test4')
show()

>>>>> "Andrew" == Andrew Straw <str...@as...> writes:

  Andrew> OK, I discovered a rather easy demonstration of what looks
  Andrew> to be a OK, (possibly the same) bug (adding 'ylim' to the
  Andrew> image_origin demo):

    from matplotlib.matlab import *
    x = arange(100.0); x.shape = 10,10

    subplot(211)
    title('blue should be up')
    imshow(x, origin='upper', interpolation='nearest')
    set(gca(),'ylim',(0,23))

    subplot(212)
    title('blue should be down')
    imshow(x, origin='lower', interpolation='nearest')
    set(gca(),'ylim',(0,23))

    #savefig('image_origin2')
    show()

I agree it looks like a bug.  Just to make sure we are on the same
page: what do you thing the correct behavior should be here?  Do you
think subplot (211) is correct?  

The question is: how should image origin interact with the axes
ylimits?  In this example, should both images fill the space from 0-10
on the y axes, regardless of the ylimit setting?  In the origin upper
case, should the blue strip be at 10 and the red strip at zero?  In
the origin lower case, should the blue be at zero and the red at 10?

Inquiring minds want to know.

Thanks,
JDH

On Tue, 13 Jul 2004, John Hunter wrote:

>
> I am having trouble running this - proj.py is making the system call
> 'stdin,stdout=os.popen2(cmd,mode='b')' where cmd is
>
> proj +R=63712000 +lon_0=-107 +proj=lcc +lat_1=50 +lat_2=50 -b
>
> but I don't have proj on my system.  Am I missing something?  Is this
> a tool I can easily install?

John:  Yes, it's easy to install.  You can grab it from 
http://proj.maptools.org.  A simple ./configure; make; make install worked 
for me.

>
> Also, while I have your ear, what version of matplotlib are you
> working with and what backend do you typically use?

I have 0.60.2 and usually use the default gtkagg backend.  I work mostly 
on OS X (I'm actually the fink package maintainer).

>
> Nice map, by the way!
>

Thanks - very nice software!

-Jeff

-- 
Jeffrey S. Whitaker         Phone : (303)497-6313
NOAA/OAR/CDC  R/CDC1        FAX   : (303)497-6449
325 Broadway                Web   : http://www.cdc.noaa.gov/~jsw
Boulder, CO, USA 80305-3328 Office: Skaggs Research Cntr 1D-124

>>>>> "Jeff" == Jeff Whitaker <js...@fa...> writes:

    Jeff> Hi:

    Jeff> I've hacked together some scripts that let me plot data on
    Jeff> top of world maps with different map projections in
    Jeff> matplotlib.  It's quite slow, but it works!  There's an
    Jeff> example at
    Jeff> http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/plotmap.py
    Jeff> (the plot should look like
    Jeff> http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/plotmap.png).
    Jeff> To run the example you'll need the data files
    Jeff> http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/wcl.txt
    Jeff> (world coastlines) and
    Jeff> http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/topodata.pickle
    Jeff> (data to plot), as well as a module I wrote to compute
    Jeff> cartographic transformations using proj4
    Jeff> http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/python/proj.py.


    Jeff> If anyone knows of a faster way to do this, I'd love to hear
    Jeff> about it.

I am having trouble running this - proj.py is making the system call
'stdin,stdout=os.popen2(cmd,mode='b')' where cmd is

 proj +R=63712000 +lon_0=-107 +proj=lcc +lat_1=50 +lat_2=50 -b

but I don't have proj on my system.  Am I missing something?  Is this
a tool I can easily install?

Also, while I have your ear, what version of matplotlib are you
working with and what backend do you typically use?

Nice map, by the way!

JDH

Hi:

I've hacked together some scripts that let me plot data on top of world 
maps with different map projections in matplotlib.  It's quite slow, but 
it works!  There's an example at 
http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/plotmap.py (the plot 
should look like 
http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/plotmap.png).  To run 
the example you'll need the data files 
http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/wcl.txt (world 
coastlines) and 
http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/topodata.pickle (data to 
plot), as well as a module I wrote to compute cartographic transformations 
using proj4 
http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/python/proj.py.


If anyone knows of a faster way to do this, I'd love to hear about it.

-Jeff




-- 
Jeffrey S. Whitaker         Phone : (303)497-6313
NOAA/OAR/CDC  R/CDC1        FAX   : (303)497-6449
325 Broadway                Web   : http://www.cdc.noaa.gov/~jsw
Boulder, CO, USA 80305-3328 Office: Skaggs Research Cntr 1D-124

John Hunter wrote:

>
>    Andrew> John+Perry+whoever: the only matplotlib bug I found this
>    Andrew> time is worked around by the line
>
>    >> rcParams['image.origin'] = 'upper' # 'lower': nav toolbar
>    >> problem -ADS
>
>    Andrew> We should fix that at some point.
>
>What was the bug?  I didn't notice any problems with lower.
>  
>
OK, I discovered a rather easy demonstration of what looks to be a 
(possibly the same) bug (adding 'ylim' to the image_origin demo):

from matplotlib.matlab import *

x = arange(100.0); x.shape = 10,10


subplot(211)
title('blue should be up')
imshow(x, origin='upper', interpolation='nearest')
set(gca(),'ylim',(0,23))

subplot(212)
title('blue should be down')
imshow(x, origin='lower', interpolation='nearest')
set(gca(),'ylim',(0,23))

#savefig('image_origin2')
show()