GRASS Python Scripting Library

The code in lib/python/ provides grass.script in order to support GRASS scripts written in Python. The scripts directory of GRASS 7 contains a series of examples actually provided to the end users (while the script in GRASS 6 are shell scripts).

Python Scripting Library code details:

• for GRASS 6: core.py, db.py, raster.py, vector.py, setup.py, array.py task.py
• for GRASS 7: core.py, db.py, raster.py, vector.py, setup.py, array.py task.py

Examples

Display example

Example of Python script, which is processed by g.parser:

#!/usr/bin/env python
#
############################################################################
#
# MODULE:      d.shadedmap
# AUTHOR(S):   Unknown; updated to GRASS 5.7 by Michael Barton
#              Converted to Python by Glynn Clements
# PURPOSE:     Uses d.his to drape a color raster over a shaded relief map
# COPYRIGHT:   (C) 2004,2008,2009 by the GRASS Development Team
#
#              This program is free software under the GNU General Public
#              License (>=v2). Read the file COPYING that comes with GRASS
#              for details.
#
#############################################################################

#%Module
#% description: Drapes a color raster over a shaded relief map using d.his
#%End
#%option
#% key: reliefmap
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of shaded relief or aspect map
#% required : yes
#%end
#%option
#% key: drapemap
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of raster to drape over relief map
#% required : yes
#%end
#%option
#% key: brighten
#% type: integer
#% description: Percent to brighten
#% options: -99-99
#% answer: 0
#%end

import sys
from grass.script import core as grass

def main():
    drape_map = options['drapemap']
    relief_map = options['reliefmap']
    brighten = options['brighten']
    ret = grass.run_command("d.his", h_map = drape_map,  i_map = relief_map, brighten = brighten)
    sys.exit(ret)

if __name__ == "__main__":
    options, flags = grass.parser()
    main()

Parsing the options and flags

grass.parser() is an interface to g.parser, and allows to parse the options and flags passed to your script on the command line. It is to be called at the top-level:

if __name__ == "__main__":
    options, flags = grass.parser()
    main()

Global variables "options" and "flags" are Python dictionaries containing the options/flags values, keyed by lower-case option/flag names. The values in "options" are strings, those in "flags" are Python booleans. All those variables have to be previously declared in the header of your script.

>>> options, flags = grass.parser()
>>> options
{'input': 'my_map', 'output': 'map_out', 'option1': '21.472', 'option2': ''}
>>> flags
{'c': True, 'm': False}

Example for embedding r.mapcalc (map algebra)

grass.mapcalc() accepts a template string followed by keyword arguments for the substitutions, e.g. (code snippets):

grass.mapcalc("${out} = ${rast1} + ${rast2}",
              out = options['output'],
              rast1 = options['raster1'],
              rast2 = options['raster2'])

Best practice: first copy all of the options[] into separate variables at the beginning of main(), i.e.:

def main():
    output = options['output']
    raster1 = options['raster1']
    raster2 = options['raster2']
 
    ...
 
    grass.mapcalc("${out} = ${rast1} + ${rast2}",
                  out = output,
                  rast1 = raster1,
                  rast2 = raster2)

Example for parsing raster category labels

How to obtain the text labels

    # dump cats to file to avoid "too many argument" problem:
    p = grass.pipe_command('r.category', map = rastertmp, fs = ';', quiet = True)
    cats = []
    for line in p.stdout:
        cats.append(line.rstrip('\r\n').split(';')[0])
    p.wait()

    number = len(cats)
    if number < 1:
        grass.fatal(_("No categories found in raster map"))

Example for parsing category numbers

Q: How to obtain the number of cells of a certain category?

A: It is recommended to use pipe_command() and parse the output, e.g.:

       p = grass.pipe_command('r.stats',flags='c',input='map')
       result = {}
       for line in p.stdout:
           val,count = line.strip().split()
           result[int(val)] = int(count)
       p.wait()

Example for getting the region's number of rows and columns

Q: How to obtain the number of rows and columns of the current region?

A: It is recommended to use the "grass.region()" function which will create a dictionary with values for extents and resolution, e.g.:

#!/usr/bin/env python
#-*- coding:utf-8 -*-
#
############################################################################
#
# MODULE:       g.region.resolution
# AUTHOR(S):    based on a post at GRASS-USER mailing list [1]               
# PURPOSE:	Parses "g.region -g", prints out number of rows, cols
# COPYLEFT:     ;-)
# COMMENT:      ...a lot of comments to be easy-to-read for/by beginners
#
#############################################################################
#
#%Module
#% description: Print number of rows, cols of current geographic region
#% keywords: region
#%end

# importing required modules
import sys # the sys module [2]
from grass.script import core as grass # the core module [3]

# information about imported modules can be obtained using the dir() function
# e.g.: dir(sys)

# define the "main" function: get number of rows, cols of region
def main():
    
    # #######################################################################
    # the following commented code works but is kept only for learning purposes
     
    ## assigning the output of the command "g.region -g" in a string called "return_rows_x_cols"
    # return_rows_x_cols = grass.read_command('g.region', flags = 'g')
    
    ## parsing arguments of interest (rows, cols) in a dictionary named "rows_x_cols"
    # rows_x_cols = grass.parse_key_val(return_rows_x_cols)
    
    ## selectively print rows, cols from the dictionary "rows_x_cols"
    # print 'rows=%d \ncols=%d' % (int(rows_x_cols['rows']), int(rows_x_cols['cols']))
    
    # #######################################################################
    
    # faster/ easier way: use of the "grass.region()" function
    gregion = grass.region()
    rows = gregion['rows']
    cols = gregion['cols']
    
    # print rows, cols properly formated 
    print 'rows=%d \ncols=%d' % (rows, cols)

# this "if" condition instructs execution of code contained in this script, *only* if the script is being executed directly 
if __name__ == "__main__": # this allows the script to be used as a module in other scripts or as a standalone script
    options, flags = grass.parser() #
    sys.exit(main()) #

# Links
# [1] http://n2.nabble.com/Getting-rows-cols-of-a-region-in-a-script-tp2787474p2787509.html
# [2] http://www.python.org/doc/2.5.2/lib/module-sys.html
# [3] http://download.osgeo.org/grass/grass6_progman/pythonlib.html#pythonCore

Managing mapsets

To check if a certain mapset exists in the active location, use:

       grass.script.mapsets(False)

... returns a list of mapsets in the current location.

r.mapcalc example

Example of Python script, which is processed by g.parser:

The shell script line:

  r.mapcalc "MASK = if(($cloudResampName < 0.01000),1,null())"

would be written like this:

       import grass.script as grass

       ...

       grass.mapcalc("MASK=if(($cloudResampName < 0.01000),1,null())",
                     cloudResampName = cloudResampName)

The first argument to the mapcalc function is a template (see the Python library documentation for string.Template). Any keyword arguments (other than quiet, verbose or overwrite) specify substitutions.

Using output from GRASS modules in the script

Sometimes you need to use the output of a module for the next step. There are dedicated functions to obtain the result of, for example, a statistical analysis.

Example: get the range of a raster map and use it in r.mapcalc. Here you can use grass.script.raster_info(), e.g.:

       import grass.script as grass

       max = grass.raster_info(inmap)['max']
       grass.mapcalc("$outmap = $inmap / $max",
                     inmap = inmap, outmap = outmap, max = max)

Calling a GRASS module in Python

Imagine, you wanted to execute this command in Python:

  r.profile -g input=mymap output=newfile profile=12244.256,-295112.597,12128.012,-295293.77

All arguments except the first (which is a flag) are keyword arguments, i.e. arg = val. For the flag, use flags = 'g' (note that "-g" would be the negative of a Python variable named "g"!). So:

       grass.run_command(
               'r.profile',
               input = input_map,
               output = output_file,
               profile = [12244.256,-295112.597,12128.012,-295293.77]

or:

               profile = [(12244.256,-295112.597),(12128.012,-295293.77)]

i.e. you need to provide the keyword, and the argument must be a valid Python expression. Function run_command() etc accept lists and tuples.

Differences between run_command() and read_command():

• run_command() executes the command and waits for it to terminate; it doesn't redirect any of the standard streams.
• read_command() executes the command with stdout redirected to a pipe, and reads everything written to it. Once the command terminates, it returns the data written to stdout as a string.

How to retrieve error messages from read_command():

None of the existing *_command functions redirect stderr. You can do so with e.g.:

def read2_command(*args, **kwargs):
   kwargs['stdout'] = grass.PIPE
   kwargs['stderr'] = grass.PIPE
   ps = grass.start_command(*args, **kwargs)
   return ps.communicate()

This behaves like read_command() except that it returns a tuple of (stdout,stderr) rather than just stdout.

Percentage output for progress of computation

A) Within a Python script, the grass.script.core.percent() module method wraps the g.message -p command.

B) If you call a GRASS command within the Python code, you have to parse the output by setting GRASS_MESSAGE_FORMAT=gui in the environment when running the command and read from the command's stderr; e.g.

       import grass.script as grass
       env = os.environ.copy()
       env['GRASS_MESSAGE_FORMAT'] = 'gui'
       p = grass.start_command(..., stderr = grass.PIPE, env = env)
       # read from p.stderr
       p.wait()

If you need to capture both stdout and stderr, you need to use threads, select, or non-blocking I/O to consume data from both streams as it is generated in order to avoid deadlock.

ALTERNATIVE:

Redirect both stdout and stderr to the same pipe (and hope that the normal output doesn't include anything which will be mistaken for progress/error/etc messages):

       p = grass.start_command(..., stdout = grass.PIPE, stderr = grass.STDOUT, env = env)

NULL data management

How to analyse if there are only NULL cells in a map:

If a map contains only null cells, its minimum and maximum will be "NULL":

       $ r.mapcalc 'foo = null()'
       $ r.info -r foo
       min=NULL
       max=NULL

Using the Python API, The 'min' and 'max' values in the result of the raster_info() function will be None.

Counting cells

Counting cells is far more expensive than simply determining whether there are any non-null cells. Counting cells requires reading the entire map, while the r.info approach only needs to read the metadata files.

If you do need to count cells, r.stats is likely to be more efficient than {{cmd|r.univar}.

A count loop:

       while grass.raster_info(inmap)['max'] is not None:
           ...

Path to GISDBASE

In order to a avoid hardcoded paths to GRASS mapset files like the SQLite DB file, you can get the GISDBASE variable from the environment:

       import grass.script as grass
       import os.path

       env = grass.gisenv()

       gisdbase = env['GISDBASE']
       location = env['LOCATION_NAME']
       mapset = env['MAPSET']

       path = os.path.join(gisdbase, location, mapset, 'sqlite.db')

Use Python reserved keyword

Question: r.resamp.bspline uses 'lambda' as a command line parameter name, but when you try to use it with grass.run_command() you get an error as lambda is a python reserved keyword. How to work around that?

Answer: Prepend an underscore to the name, i.e.:

       grass.run_command('r.resamp.bspline', _lambda = ...)

Controlling the PNG display driver

Code fragment to control the pngdriver in Python:

import os
import sys
from grass.script import core as grass
def main():
       os.environ['GRASS_PNGFILE'] = filename
       os.environ['GRASS_WIDTH'] = str(width)
       os.environ['GRASS_HEIGHT'] = str(height)
       grass.run_command('d.his', i='elevation_shade', h='elevation')

Sophisticated cleanup procedure

Scripts which create several temporary files need a more sophisticated cleanup procedure which deletes all the tmp maps which have been created. This procedure should also work if the script stops (e.g due to an error).

Solution: Define a list of map names which starts out empty and has names appended to it as the names are generated. Code fragment:

       tmp_rast = []

       def cleanup():
           for rast in tmp_rast:
               grass.run_command("g.remove",
                                 rast = rast,
                                 quiet = True)

       def main():
           ...
           while ...:
               next_rast = ...
               tmp_rast.append(next_rast)
               ...

       if __name__ == "__main__":
           options, flags = grass.parser()
           atexit.register(cleanup)
           sys.exit(main())