Unit 11 - PyGRASS scripting

Let’s start with a Python script created by Graphical modeler in Unit 10. Before saving the script remove the lines below to avoid generating GUI dialogs when launching the script (will be explained in Unit 12 - Script User Interface).

#%module
#% description: NDVI model version 2
#%end

Then run a script from Layer Manager grass-script-load Launch user-defined script main toolbar.

Note

Before starting a script GRASS can ask you to add script directory path into GRASS_ADDON_PATH. It can be useful if you will run script(s) from this directory more often. Then you don’t need to define full path to the scripts, its name will be enough.

../_images/addon-path.png

Fig. 77 Register script directory into GRASS Addon Path.

After selecting a script to run, the command output will be printed in Layer Manager Console tab.

../_images/script-output.png

Fig. 78 Script output printed in Layer Manager Console.

Python script generated by Graphical Modeler is based on GRASS Scripting Library. As a first step the script will be rewritten into PyGRASS syntax.

PyGRASS syntax

Open exported Python script by your favorite editor or if your do not have any just use GRASS-integrated Python editor grass-python Open a simple Python code editor.

../_images/editor.png

Fig. 79 Simple GRASS Python code editor in action.

Load script by grass-open Open and replace every occurence of core.run_command function by PyGRASS equivalent. PyGRASS allows calling GRASS modules similarly as GRASS Scripting Library does (see Unit 10 - Python intro). The module caller is represented by Module class. In contrast to GRASS Scripting Library which defines several routines to run module (core.run_command, core.read_command, or core.feed_command) in PyGRASS there is only one caller technique.

Replace all occurrence of core.run_command function by Module caller, see example below.

from grass.script import run_command

run_command("v.overlay",
             overwrite = True,
             ainput = "jena_boundary@PERMANENT",
             alayer = "1",
             atype = "auto",
             binput = "MaskFeature@PERMANENT",
             blayer = "1",
             btype = "area",
             operator = "not",
             output = "region_mask",
             olayer = "1,0,0",
             snap = 1e-8)

by

from grass.pygrass.modules import Module

Module("v.overlay",
       overwrite = True,
       ainput = "jena_boundary@PERMANENT",
       alayer = "1",
       atype = "auto",
       binput = "MaskFeature@PERMANENT",
       blayer = "1",
       btype = "area",
       operator = "not",
       output = "region_mask",
       olayer = "1,0,0",
       snap = 1e-8)

Warning

There are some caveats. Mupliple options given as a string in GRASS Scripting Library must be given as a list of strings in PyGRASS, see v.clean example below.

run_command("v.clean",
            ...
            type = "point,line,boundary,centroid,area,face,kernel",
            ...
Module("v.clean",
       ...
       type = ["point","line","boundary","centroid","area","face","kernel"],
       ...

In the next step the script will be improved by printing NDVI value statistics (be aware of indentation), see Unit 10.

    ret = Module('r.univar', flags='g', map='ndvi', stdout_=PIPE)
    stats = parse_key_val(ret.outputs.stdout)
    print ('-' * 80)
    print ('NDVI value statistics')
    print ('-' * 80)
    print ('NDVI min value: {0:.4f}'.format(float(stats['min'])))
    print ('NDVI max value: {0:.4f}'.format(float(stats['max'])))
    print ('NDVI mean value: {0:.4f}'.format(float(stats['mean'])))

Note

Import relevant functions

from subprocess import PIPE

from grass.script import parser, parse_key_val

Launch script by grass-execute Run and check out an output in Layer Manager Console tab.

../_images/run-script.svg

Fig. 80 Run script from Python editor.

Statistics

Also NDVI classes statistics could be reported. Area size can be computed by v.report.

    print ('-' * 80)
    print ('NDVI class statistics')
    print ('-' * 80)
    ret = Module('v.report', map='ndvi_vector', option='area', stdout_=PIPE)
    for line in ret.outputs.stdout.splitlines()[1:]: # skip first line (cat|label|area)
        # parse line (eg. 1||2712850)
        data = line.split('|')
        cat = data[0]
        area = float(data[-1])
        print ('NDVI class {0}: {1:.1f} ha'.format(cat, area/1e4)) 

Output of v.report module need to be parsed. Unfortunately the command does not offer a shell script output similarly to r.univar. Python fuctions like splitlines() and split() can be used.

At the end NDVI zonal statistics for each class can be computed:

  • zonal statistics can be computed by v.rast.stats and stored in attribute table
  • attributes can be printed by v.db.select

Todo

Can be simplified since this statistics is already reported by v.report

    # v.to.rast: use -c flag for updating statistics if exists
    Module('v.rast.stats', flags='c', map='ndvi_vector', raster='ndvi',
           column_prefix='ndvi', method=['minimum','maximum','average'])
    # v.db.select: don't print column names (-c)
    ret = Module('v.db.select', flags='c', map='ndvi_vector', separator='comma', stdout_=PIPE)
    for line in ret.outputs.stdout.splitlines():
        # parse line (eg. 1,,-0.433962264150943,0.740350877192983,0.051388909449992)
        cat,label,min,max,mean = line.split(',')
        print ('NDVI class {0}: {1:.4f} (min) {2:.4f} (max) {3:.4f} (mean)'.format(
        cat, float(min), float(max), float(mean)))

Example of script output below.

--------------------------------------------------------------------------------
NDVI value statistics
--------------------------------------------------------------------------------
NDVI min value: -0.6094
NDVI max value: 0.9997
NDVI mean value: 0.6485
--------------------------------------------------------------------------------
NDVI class statistics
--------------------------------------------------------------------------------
NDVI class 1: 271.3 ha
NDVI class 2: 2438.7 ha
NDVI class 3: 7561.0 ha
--------------------------------------------------------------------------------
NDVI class 1: -0.4340 (min) 0.7404 (max) 0.0514 (mean)
NDVI class 2: -0.1738 (min) 0.8547 (max) 0.3262 (mean)
NDVI class 3: -0.6094 (min) 0.9997 (max) 0.7740 (mean)

Tip

In order to simplify testing and increase code readability our code could be split into two functions: compute() and stats().

def main(options, flags):
    compute(options)
    stats(options)

    return 0

Sample script to download: ndvi-v2.py