GIS/MEA 584:
Mapping and Analysis Using UAS

Introduction to Python scripting in GRASS GIS 7

Python code can be executed in GRASS GIS GUI in the Python tab in the Layer Manager. You can enter any Python code, we can, for example, do a simple calculation: 
import math
5 * math.sin(math.pi/2)
This is an quick way to run small code snippets and test Python code you are currently working on.

You can find your own way of running Python code in GRASS GIS, but for now use the Simple Python Editor which comes with GRASS GIS 7.2 or newer or the interactive Python console in the Python tab where you should only paste one line at a time. For scripts written in an external editor such as Spyder, Atom, or Notepad++, you can use File > Launch script which will run your script in the Console tab and you will be able to run it again by using arrow up and Enter key as in case of other commands.

We are using North Carolina GRASS Location in this assignment. If you already have North Carolina GRASS GIS sample dataset, you don't have to download anything.

Python Scripting Library

The GRASS GIS 7 Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include:
  • run_command: most often used with modules which output raster/vector data and where text output is not expected
  • read_command: used when we are interested in text output
  • parse_command: used with modules producing text output as key=value pair
  • write_command: for modules expecting text input from either standard input or file
In addition, this library provides several wrapper functions for selected general-purpose modules.

Calling GRASS GIS modules

We start by importing GRASS GIS Python Scripting Library: 
import grass.script as gs
First, we set the computational extent and resolution to the raster layer elevation. 
gs.run_command('g.region', raster='elevation')
The run_command() function is used most often. Here, we apply the focal operation average (r.neighbors) to smooth the elevation raster layer. 
gs.run_command('r.neighbors', input='elevation', output='elev_smoothed', method='average', flags='c')
Note that the syntax is similar to the syntax used in the manuals, command line, or Bash; one of the bigger differences is that the flags are specified in a parameter instead of using dashes.

Calling GRASS GIS modules with textual input or output

Textual output from modules can be captured using the read_command() function. First we get the current computational region settings: 
print(gs.read_command('g.region', flags='p'))
Then we get univariate raster statistics and we add a g flag to what we would use in the command line. 
print(gs.read_command('r.univar', map='elev_smoothed', flags='g'))
Certain modules can produce output in key-value format which is enabled by the g flag. The parse_command() function automatically parses this output and returns a dictionary. In the following example, we call g.proj to get the projection parameters of the actual location. 
gs.parse_command('g.proj', flags='g')
For comparison, below is the same example, but using the read_command() function. 
print(gs.read_command('g.proj', flags='g'))
Certain modules require the text input be in a file or provided as standard input. Using the write_command() function we can conveniently pass the string as the standard input to the module. Here, we are creating a new vector with one point with the v.in.ascii module. Note that the stdin parameter is not used as a module parameter, but its content is passed as the standard input to the subprocess. 
gs.write_command('v.in.ascii', input='-', stdin='%s|%s' % (635818.8, 221342.4), output='view_point')

Convenient wrapper functions

Some modules have wrapper functions to simplify the most common tasks. We can obtain the information about the vector layer which we just created with the v.info wrapper. 
gs.vector_info('view_point')
It is also possible to retrieve the raster layer history (r.support) and layer information (r.info) or to query (r.what) raster layer pixel values. The pixel (cell) values are obtained using the raster_what() function: 
gs.raster_what('elevation', [[635818.8, 221342.4], [635710, 221330]])
As another example, the r.mapcalc wrapper for raster algebra allows using a long expressions. 
gs.mapcalc("elev_strip = if(elevation > 100 && elevation < 125, elevation, null())")
print(gs.read_command('r.univar', map='elev_strip', flags='g'))
The g.region wrapper is a convenient way to retrieve the current region settings (i.e., computational region). It returns a dictionary with values converted to appropriate types (floats and ints). 
region = gs.region()
print region
# cell area in map units (in projected Locations)
region['nsres'] * region['ewres']
We can list data stored in a GRASS GIS location with g.list wrappers. With this function, the map layers are grouped by mapsets (in this example, raster layers): 
gs.list_grouped(['raster'])
Here is an example of a different g.list wrapper which structures the output as list of pairs (name, mapset). We also obtain current mapset with g.gisenv wrapper. 
current_mapset = gs.gisenv()['MAPSET']
gs.list_pairs('raster', mapset=current_mapset)

PyGRASS

Standard scripting with GRASS GIS modules in Python may seem cumbersome in certain cases namely when you have to do to things like iterating over vector features or raster rows/columns. PyGRASS library introduces several objects that allow to interact directly with the data, individual features and cells, using the underlying C API of GRASS GIS libraries.

Working with vector data

Import the necessary classes: 
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Point
Create an instance of a new vector map: 
points = VectorTopo('points')
Open the map in write mode: 
points.open(mode='w')
Create some vector geometry features, like two points: 
point1 = Point(635818.8, 221342.4)
point2 = Point(633627.7, 227050.2)
Add the above two points to the new vector map with: 
points.write(point1)
points.write(point2)
Finally close the vector map: 
points.close()
Now do the same thing using the context manager syntax provided by with, and creating and setting also the attribute table: 
# Define the columns of the new vector map
cols = [(u'cat',       'INTEGER PRIMARY KEY'),
        (u'name',      'TEXT')]

with VectorTopo('points', mode='w', tab_cols=cols, overwrite=True) as points:
    # save the point, category (identifier) and the attribute
    points.write(point1, cat=1, attrs=('pub', ))
    points.write(point2, cat=2, attrs=('restaurant', ))
    # save the changes to the database
    points.table.conn.commit()
Note: We don't need to close the vector map because it is already closed by the context manager (used through the with statement).

Reading an existing vector map: 

with VectorTopo('points', mode='r') as points:
    for point in points:
        print(point, point.attrs['name'])
The above content is derived from workshop How to write a Python GRASS GIS 7 addon presented at FOSS4G Europe 2015. The relevant parts which contain more details are: