[GRASS-dev] matplotlib example script

[Michael Barton]

Thanks for the input Glynn. A few responses below.


On Jul 24, 2008, at 12:24 PM, Glynn Clements wrote:

>
> 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.

Hmmm. I might be able to help that, but it's hard to remember to  
switch the way the program works between people here who expect html- 
like and plain text.

>
>
>> 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.

I agree. One possibility is to use the ps backend if you are indeed  
planning to switch the output to postscript. In this example, I just  
wanted to show that TclTk, wxPython, Qt or other GUI is not required.

>
>
> 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.

Forgot to take this out. This is an issue that my system is having and  
not only does it not belong in this script, it doesn't help inside  
this script.


>
>
>>    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.]
>

The strip() is probably unnecessary; I have a habit of worrying about  
extraneous spaces causing formating problems.

> 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()).

Thanks.

>
>
>>            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.

Well, I was worried about that too. However, even with the 10m DEM in  
spearfish, the main time is spent in running r.stats. The rest of the  
code takes no appreciable time to run. I should try it on a big Terra/ 
ASTER or Landsat ETM file and see how long it takes.

The fastest thing probably is to just read the GRASS histogram file.  
But this doesn't permit changing bin sizes or other plotting options.

>
>
> 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.

I agree about r.stats. Wouldn't binning in the script take as long as  
the current method?


>
>
>>        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).

Yes. And it might be one in MatPlotLib. It has some color function,  
but I haven't explored them yet. Do the grass color names have html  
equivalents? If so, we don't need /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?

You are right.

Thanks again
Michael