Unit 11 - PyGRASS scripting

Let’s start with a Python script created by Graphical modeler in Unit 10 saved into file.

Important

Before saving the script we will remove the lines below to avoid generating GUI dialogs when launching the script.

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

Then we can try to run a script from Layer Manager ^{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.

Fig. 78 Register script directory into GRASS Addon Path.

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

Fig. 79 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 scripting

Open Python script by your favorite editor or if your do not have any use GRASS integrated Python editor _{Open a simple Python code editor}.

Fig. 80 Simple GRASS Python code editor in action.

Load script by ^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 = "oslo@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 = "oslo@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.

from subprocess import PIPE
from grass.script import parse_key_val

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'])))

Launch script by ^Run and check out an output in Layer Manager Console tab.

Fig. 81 Run script from Python editor.

Statistics

Also NDVI classes statistics could be reported. Area size can be easily 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 shell script output similarly to r.univar. We will implement our own parsing technique based on Python functions like splitlines() and split().

Than also 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

# 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.9914
NDVI max value: 0.9982
NDVI mean value: 0.5650
--------------------------------------------------------------------------------
NDVI class statistics
--------------------------------------------------------------------------------
NDVI class 1: 3588.8 ha
NDVI class 2: 3864.8 ha
NDVI class 3: 15828.2 ha
NDVI class 1: -0.9914 (min) 0.9710 (max) -0.2770 (mean)
NDVI class 2: -0.9286 (min) 0.9487 (max) 0.2796 (mean)
NDVI class 3: -0.9545 (min) 0.9982 (max) 0.8258 (mean)

Tip

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

def main():
    compute()
    stats()

    return 0

Sample script to download: ndvi-v2.py