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