PyGRASS is a library originally developed during the Google Summer of Code 2012. The PyGRASS library adds two main functionalities:
For further discussion about the implementation ideas and performance are presented in the article:
Standard scripting with GRASS modules in Python may sometime seem discouraging especially when you have to do conceptually simple things like: iterate over vector features or raster rows/columns. With the standard scripting approach you would have to use several workaround such as:
Most of GRASS GIS is written in C programming language. Using the C API, all this is much more simple since you can directly work on GRASS GIS data and do just what you need to do. However, you perhaps want to stick to Python. Here, the PyGRASS library introduced several objects that allow to interact directly with the data using the underlying C API of GRASS GIS. Three main areas are covered so far:
Import the necessary classes:
In [ ]:
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Point
Create an instance of the vector map, with:
In [ ]:
view_points = VectorTopo('view_points')
Open the map in write mode, with:
In [ ]:
view_points.open(mode='w')
Create some vector geometry features, like two points:
In [ ]:
point1 = Point(635818.8, 221342.4)
point2 = Point(633627.7, 227050.2)
Add the above two points to the new vector map with:
In [ ]:
view_points.write(point1)
view_points.write(point2)
Finally close the vector map, with:
In [ ]:
view_points.close()
Now do the same thing using the context manager syntax provided by with, and creating and setting also the attribute table, with:
In [ ]:
# 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.
In [ ]:
with VectorTopo('view_points', mode='r') as points:
for point in points:
print(point, point.attrs['name'])
A class with similar behaviour has been created for raster maps, import this class with:
In [ ]:
from grass.pygrass.raster import RasterRow
Instanziate a raster class with:
In [ ]:
elev = RasterRow('elevation', mapset='PERMANENT')
Similar to the vector approach shown above, we can open a raster map with:
In [ ]:
elev.open(mode='r')
Now we have row by row access to the whole raster map. Below we show an example to access the first row:
In [ ]:
row = elev[0]
In [ ]:
row
Since row instances inherit from the numpy arrays, you can treat them as numpy arrays:
In [ ]:
row * 2
In [ ]:
row.shape
To convert the whole raster map to a numpy array, and then back to a GRASS GIS raster map, we can define two functions like:
In [ ]:
import numpy as np
from grass.pygrass.raster.buffer import Buffer
from grass.pygrass.gis.region import Region
def raster2numpy(rastname, mapset=''):
"""Return a numpy array from a raster map"""
with RasterRow(rastname, mapset=mapset, mode='r') as rast:
return np.array(rast)
def numpy2raster(array, mtype, rastname, overwrite=False):
"""Save a numpy array to a raster map"""
reg = Region()
if (reg.rows, reg.cols) != array.shape:
msg = "Region and array are different: %r != %r"
raise TypeError(msg % ((reg.rows, reg.cols), array.shape))
with RasterRow(rastname, mode='w', mtype=mtype, overwrite=overwrite) as new:
newrow = Buffer((array.shape[1],), mtype=mtype)
for row in array:
newrow[:] = row[:]
new.put_row(newrow)
In [ ]:
elev = raster2numpy('elevation', mapset='PERMANENT')
In [ ]:
elev
In [ ]:
numpy2raster(elev * 2, mtype='FCELL', rastname='doubled_elev', overwrite=True)
In [ ]:
# import a function to render the results
from render import view
from grass.pygrass.modules.shortcuts import raster as r
# set the colors of the input and output maps to a predefined elevation color table
r.colors(map='elevation@PERMANENT,doubled_elev', color='elevation')
In [ ]:
view(['elevation@PERMANENT'])
In [ ]:
view(['doubled_elev'])
Import the Module class from the pygrass library, with:
In [ ]:
from grass.pygrass.modules import Module
Execute a GRASS GIS module with (here: viewshed calculation):
In [ ]:
viewshed = Module('r.viewshed', input='elevation', output='viewshed2', coordinates=(635818.8, 221342.4), flags='c', overwrite=True)
We can play with the python instance asking for inputs and outputs parameters and flags
In [ ]:
viewshed.inputs.input
In [ ]:
viewshed.inputs.coordinates
In [ ]:
viewshed.outputs.output
In [ ]:
viewshed.flags.c
We can change the parameters values, parameterwise with:
In [ ]:
viewshed.inputs.coordinates = [(635825, 221350)]
Or re-execute the module passing only the parameter(s) that need(s) to be different:
In [ ]:
viewshed(output='viewshed_newpoint')
We can retrieve the command string that will be execute in both bash and Python style:
In [ ]:
viewshed.get_bash()
In [ ]:
viewshed.get_python()
The run_ parameter defines if the module must be execute or not, the default value is True.
The finish_ parameter defines if the module has to wait the end of the process, default is True
The env_ parameter defines the environmental variable that will be set and used by the process.
The parameters stdin_, stdout_ and stderr_ are used to pass or capture the input or the output of the GRASS GIS modules.
In [ ]:
import subprocess as sub
# Another way to use the Module class objects is through the shortcuts
# import the shortcuts for the general modules
from grass.pygrass.modules.shortcuts import general as g
Now we can call and capture the stdout of a GRASS GIS modules, with:
In [ ]:
greg = g.region(flags='p', stdout_=sub.PIPE)
print(greg.outputs.stdout)