[GRASS-dev] matplotlib example script

Glynn Clements glynn at gclements.plus.com
Thu Jul 24 15:24:56 EDT 2008


Michael Barton wrote:

> I've attached an example script and output for using the MatPlotLib  
> Python library in GRASS.

One of the Apple Mail developers needs a dose of the clue-stick. The
script (and image) were attached to the HTML alternative of the
multipart/alternative section, so they don't exist for anyone using
the plain-text alternative.

> This script replicates the functionality of  
> d.histogram, but makes a much nicer looking plot. It also has more  
> options, although I've included only a very few of the formatting  
> possibilities here. The script sends the output from r.stats into 9  
> lines of code to produce the histogram with formatting options. Most  
> of it is done by 4 lines. This can be run from the command line and  
> does not require a GUI to be installed. It writes to a PNG file, but  
> could easily be set to create a TclTk or wxPython interactive window  
> as a user-selectable option.

Before we start using MatPlotLib extensively, we need to figure out
how to handle output format selection so that all script can work with
any backend, and the scripts can be used alongside d.* commands
without the caller needing to distingiush between them.

Or, at least, we need to figure out an interface that will allow us to
do that later without having to modify large numbers of scripts. IOW,
somthing like a grass.plot module which hides the details behind
begin() and end() methods, so the scripts themselves don't reference
specific format types.

A couple of brief comments on the code:

> if os.environ['PYTHONPATH'] != None:
>     os.environ['PYTHONPATH'] = '/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages:'+os.environ['PYTHONPATH']

This doesn't belong here.

>     mapstats = grass.read_command('r.stats', input = options['input'], fs = ',',
>                                     nsteps = '255', flags = 'Acn')
>     histlist = mapstats.splitlines()
>     for pair in histlist:

Ideally, read_command() shouldn't used for anything which might create
a lot of output.

For processing "bulk" output, iterate over stdout, with e.g.:

	mapstats_p = grass.start_command('r.stats', input = options['input'], fs = ',',
	         	                  nsteps = '255', flags = 'Acn', stdout = subprocess.PIPE)
	for pair in mapstats_p.stdout:
	        try:
        	    pairlist = pair.strip().split(',')
	...

	mapstats_p.wait()

[Note the added strip(); the "for line in file" idiom doesn't strip the
line terminator from the line.]

As that's likely to be quite common, I have added grass.pipe_command(),
which is like start_command() but adds the "stdout = subprocess.PIPE"
automatically, so you can use:

	mapstats_p = grass.pipe_command('r.stats', input = options['input'], fs = ',',
	         	                  nsteps = '255', flags = 'Acn')
	for pair in mapstats_p.stdout:
	...

It returns the Popen object rather than the pipe itself, as you need
the Popen object for the wait() call (even if you don't want the exit
code, the child process will remain as a zombie until you call wait()).

>             for i in range(cells):
>                 plotlist.append(val)

Ouch. This is going to cause problems (i.e. fail) for large maps (even
for moderately-sized maps, it will be slow). As it stands, this script
*isn't* a substitute for d.histogram.

As axes.hist() just calls numpy.histogram() then uses axes.bar() (or
axes.barh()) to plot the result, you may as well just use axes.bar()
directly on the value,count pairs emitted by r.stats.

If you specifically want to bin the data, it would be better to either
add an option to r.stats, or to bin the data in the script as it is
being read, rather than read an entire raster map into a Python list.

>         if len(fcolor.split(':')) == 3:
>             #rgb color
>             r = float(fcolor.split(':')[0])
>             g = float(fcolor.split(':')[1])
>             b = float(fcolor.split(':')[2])
>             #translate to mpl rgb format
>             r = r / 255
>             g = g / 255
>             b = b / 255
>             fcolor = (r,g,b)

This should probably be made into a library function; it would also
need the named colours (lib/gis/named_colr.c).

>     if options['linecolor'] != None:   
>         lcolor = options['linecolor']
>     else: fcolor = 'k'

The else clause should presumably be lcolor. Does this not also need to
handle R:G:B?

-- 
Glynn Clements <glynn at gclements.plus.com>


More information about the grass-dev mailing list