GIS595/MEA792: UAV/lidar Data Analytics

This is an unmaintained course material, please see current material at:

Introduction to Python scripting in GRASS GIS 7

Python code can be executed in GRASS GIS GUI in the Python shell 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.

For this class, we will use more advanced method called IPython Notebook which gives us ability to write code, run code and also write nicely formated additional explanatory text and create a report or notes in this way.

To avoid installing IPython Notebook on all machines, we will use OSGeo Live machine from NCSU VLC in the class. Go to VLC reservation page and reserve machine called OSGeo Live. Wait for it to boot and you will get instructions how to connect to it using remote desktop method.

Once you are in OSGeo Live, start GRASS GIS from the menu at the bottom left, submenu Geospatial, Desktop GIS. In GRASS GIS, go to the system command line (not the one in GUI) and start ipython notebook (note the space).

We are using North Carolina GRASS Location in this assignment. if you already have North Carolina GRASS GIS sample dataset or you are using OSGeo Live, 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 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 often called modules.

Calling GRASS GIS modules

We start by importing GRASS GIS Python Scripting Library:

import grass.script as gscript

First, we set the computational extent and resolution to the raster layer elevation.

gscript.run_command('g.region', raster='elevation')

The run_command() function is the most commonly used one. Here, we apply the focal operation average (r.neighbors) to smooth the elevation raster layer. Note that the syntax is similar to bash syntax, just the flags are specified in a parameter.

gscript.run_command('r.neighbors', input='elevation', output='elev_smoothed', method='average', flags='c')

Calling GRASS GIS modules with textual input or output

Textual output from modules can be captured using the read_command() function.

print(gscript.read_command('g.region', flags='p'))

print(gscript.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 this example, we call g.proj to display the projection parameters of the actual location.

gscript.parse_command('g.proj', flags='g')

For comparison, below is the same example, but using the read_command() function.

print(gscript.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 to the module. Here, we are creating a new vector with one point with v.in.ascii. Note that stdin parameter is not used as a module parameter, but its content is passed as standard input to the subprocess.

gscript.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 frequent tasks. We can obtain the information about the vector layer which we just created with the v.info wrapper.

gscript.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.

gscript.raster_what('elevation', [[635818.8, 221342.4], [635710, 221330]])

As another example, the r.mapcalc wrapper for raster algebra allows using a long expressions.

gscript.mapcalc("elev_strip = if(elevation > 100 && elevation < 125, elevation, null())")
print(gscript.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 = gscript.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):

gscript.list_grouped(['raster'])

Here is an example of a different g.list wrapper which structures the output as list of pairs (name, mapset). We obtain current mapset with g.gisenv wrapper.

current_mapset = gscript.gisenv()['MAPSET']
gscript.list_pairs('raster', mapset=current_mapset)

PyGRASS

Standard scripting with GRASS modules in Python may sometime seem cumbersome especially when you have to do conceptually simple things like iterating over vector features or raster rows/columns. PyGRASS library introduced several objects that allow to interact directly with the data using the underlying C API of GRASS GIS.

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:

view_points = VectorTopo('view_points')

Open the map in write mode:

view_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:

view_points.write(point1)
view_points.write(point2)

Finally close the vector map:

view_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('view_points', mode='w', tab_cols=cols, overwrite=True) as view_points:
    # save the point and the attribute
    view_points.write(point1, ('pub', ))
    view_points.write(point2, ('restaurant', ))
    # save the changes to the database
    view_points.table.conn.commit()

Note: we don't need to close the vector map because it is already closed by the context manager.

Reading an existing vector map:

with VectorTopo('view_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 are: