Python scripting for GRASS

Example: Flooding simulation

Instead of creating our first Python script for GRASS from scratch we will re-use model created in the previous section. The Graphic Modeler allows on-the-fly conversion to simple Python script (tab Python editor).

Fig. 57 Model as Python script.

We save this script to a new file by clicking on Save as button. Python scripts can be launched afterwards from menu File ‣ Launch script or from the toolbar .

We open this file in GRASS Python editor or in your favourite editor. GRASS Python editor can be launched from Layer Manager (tab Python, button Simple editor), or from the main toolbar .

Fig. 58 Launching GRASS Python editor from Layer Manager.

Our goal is to modify the script to be capable performing flooding simulation multiple times based on start and end water level and step value.

One of the advantages is that the resultant script behaves as normal GRASS module. If user doesn’t enter all required parameters GRASS generates GUI dialog similarly to other GRASS commands.

Fig. 59 Generated GUI dialog of our Python lake script.

Resultant script to download: lake.py

#!/usr/bin/env python
#
##############################################################################
#
# MODULE:       model
#
# AUTHOR(S):    landa
#
# PURPOSE:      Script generated by wxGUI Graphical Modeler.
#
# DATE:         Mon Jul  4 18:02:51 2016
#
##############################################################################

#%module
#% description: Performs flooding simulation.
#%end

#%option 
#% key: start_level
#% description: Start water level in meters
#% required: yes
#% answer: 990
#% type: integer
#% end
#%option 
#% key: end_level
#% description: End water level in meters
#% required: yes
#% answer: 995
#% type: integer
#% end
#%option 
#% key: step
#% description: Level step
#% required: yes
#% answer: 1
#% type: integer
#% end

import sys
import os
import atexit

from grass.script import parser, run_command, message

def cleanup():
    pass

def main():
    start_level = int(options['start_level'])
    end_level = int(options['end_level'])
    step = int(options['step'])

    for level in range(start_level, end_level+1, step):
        message("Computing water level: {}".format(level))
        compute(level)

def compute(level):
    run_command("r.lake",
                overwrite = True,
                quiet = True,
                elevation = "dem37",
                water_level = level,
                lake = 'lake_{}'.format(level),
                coordinates = '531963.147664,5626869.62523')

    run_command("r.report",
                flags = 'hn',
                overwrite = True,
                map = 'lake_{}'.format(level),
                units = "h",
                null_value = "*",
                page_length = 0,
                page_width = 79,
                nsteps = 1)

    run_command("r.mapcalc",
                overwrite = True,
                expression = "lake_mask = if ( lake_{}, 1, null() )".format(level))

    run_command("r.to.vect",
                flags = 'vt',
                overwrite = True,
                quiet = True,
                input = "lake_mask",
                output = 'lake_{}'.format(level),
                type = "area",
                column = "value")

    run_command("v.overlay",
                flags = 't',
                overwrite = True,
                quiet = True,
                ainput = "Cesta",
                alayer = "1",
                atype = "auto",
                binput = 'lake_{}'.format(level),
                blayer = "1",
                btype = "area",
                operator = "and",
                output = 'road_lake_{}'.format(level),
                olayer = "1,0,0",
                snap = 1e-8)

    run_command("v.to.db",
                flags = 'pc',
                quiet = True,
                map = 'road_lake_{}'.format(level),
                layer = "1",
                type = "point,line,boundary,centroid",
                option = "length",
                columns = "dummy",
                query_layer = "1",
                separator = "pipe")


    return 0

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

Real stuff: DEM from Lidar data

Script to download: isprs.py

Tip

The script could be easily parallelized by ParallelModuleQueue using pyGRASS.

#!/usr/bin/env python

#%module
#% description: Imports lidar and orthophoto tile, create mask and interpolate DEM.
#%end

#%option
#% key: tile
#% description: Tile name (eg. TANV37)
#% multiple: yes
#% required: yes
#%end
#%option
#% key: resolution
#% description: Output DEM resolution (in meters)
#% type: double
#% answer: 1.0
#%end

import os
import sys
import atexit
import grass.script as grass

def cleanup():
    grass.run_command('g.remove', type='raster', pattern='{}.*'.format(options['tile']))
    grass.run_command('r.mask', flags='r')
    
def import_lidar(tile):
    print('Importing lidar data...')
    grass.run_command('v.in.lidar', flags='ot', input='pr_{}_5g.laz'.format(tile),
                      output='pr_{}_5g'.format(tile))

def import_orthophoto(tile):
    print('Importing orthophoto...')
    grass.run_command('r.import', flags='o', input='{}.tif'.format(tile), output=tile)
    grass.run_command('g.region', raster='{}.1'.format(tile))
    grass.run_command('r.composite', red='{}.1'.format(tile), green='{}.2'.format(tile),
                      blue='{}.3'.format(tile), output='{}.rgb'.format(tile))
    grass.run_command('r.neighbors', input='{}.rgb'.format(tile), output='{}.mode'.format(tile), method='mode')
    grass.run_command('r.mapcalc',
                      expression="{tile} = if ( isnull( {tile}.1 + {tile}.2 + {tile}.3 ), {tile}.mode, {tile}.rgb )".format(tile=tile))
    grass.run_command('r.colors', map=tile, raster='{}.rgb'.format(tile))

def create_dem(tile, resolution):
    print('Creating DEM...')
    grass.run_command('r.mask', raster=tile)
    grass.run_command('g.region', res=resolution, flags='a')
    grass.run_command('v.surf.rst', input='pr_{}_5g'.format(tile), elevation='dem_{}'.format(tile),
                      slope='slope_{}'.format(tile), npmin=80, tension=20, smooth=1)

def main():
    os.environ['GRASS_OVERWRITE'] = '1'
    os.environ['GRASS_VERBOSE'] = '0'
    
    for tile in options['tile'].split(','):
        print('Processing tile {}...'.format(tile))
        import_lidar(tile)
        import_orthophoto(tile)
        create_dem(tile, options['resolution'])
        
    return 0

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