Using GRASS GIS through Python and tangible interfaces (workshop at FOSS4G NA 2016)
This is material for FOSS4G NA 2016 workshop Using GRASS GIS through Python and tangible interfaces held in Raleigh May 2, 2016 - 09:00 to 13:00.
Learn about scripting, graphical and tangible (!) interfaces for GRASS GIS, the powerful desktop GIS and geoprocessing backend. We will start with the Python interface and finish with Tangible Landscape, a new tangible interface for GRASS GIS. Python is the primary scripting language for GRASS GIS. We will demonstrate how to use Python to automate your geoprocessing workflows with GRASS GIS modules and develop custom algorithms using a Pythonic interface to access low level GRASS GIS library functions. We will also review several tips and tricks for parallelization. Tangible Landscape is an example of how the GRASS GIS Python API can be used to build new, cutting edge tools and advanced applications. Tangible Landscape is a collaborative 3D sketching tool which couples a 3D scanner, a projector and a physical 3D model with GRASS GIS. The workshop will be a truly hands-on experience – you will play with Tangible Landscape, using your hands to shape a sand model and drive geospatial processes.
Software: GRASS GIS 7
Data: Download complete North Carolina sample dataset.
GRASS GIS introduction
You need to create a GRASS database we will use for the tutorial. Please download the dataset for the workshop, noting where the files are located on your local directory. Now, create (unless you already have it) a directory named grassdata (GRASS database) in your home folder (or Documents), unzip the downloaded data into this directory. You should now have a Location nc_spm_08_grass7 in grassdata. Start GRASS GIS in Mapset user1.
GRASS GIS Python API
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
Besides, this library provides several wrapper functions for often called modules.
Calling GRASS GIS modules
We will use GRASS GUI Python Shell to run the commands. For longer scripts, you can create a text file, save it into your current working directory and run it with python myscript.py from the GUI command console or terminal.
Tip: When copying Python code snippets to GUI Python shell, right click at the position and select Paste Plus in the context menu. Otherwise multiline code snippets won't work.
We start by importing GRASS GIS Python Scripting Library:
import grass.script as gscript
Before running any GRASS raster modules, you need to set the computational region using. In this example, we set the computational extent and resolution to the raster layer elevation.
The run_command() function is the most commonly used one. Here, we apply the focal operation average () 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')
If we run the Python commands from GUI Python console, we can use AddLayer to add the newly created layer:
Calling GRASS GIS modules with textual input or output
Textual output from modules can be captured using the read_command() function.
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 callto display the projection parameters of the actual location.
For comparison, below is the same example, but using the read_command() function.
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. 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, 221342), output='point')
If we run the Python commands from GUI Python console, we can use AddLayer to add the newly created layer:
Convenient wrapper functions
Some modules have wrapper functions to simplify frequent tasks. For example we can obtain the information about a raster layer with raster_info which is a wrapper of , or a vector layer with vector_info.
Another example is using r.mapcalc wrapper for raster algebra:
gscript.mapcalc("elev_strip = if(elevation > 100 && elevation < 125, elevation, null())") gscript.read_command('r.univar', map='elev_strip', flags='g')
Function region 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 list_grouped, the map layers are grouped by mapsets (in this example, raster layers):wrappers. With
gscript.list_grouped(type=['raster']) gscript.list_grouped(type=['raster'], pattern="landuse*")
Here is an example of a different list_pairs which structures the output as list of pairs (name, mapset). We obtain current mapset with wrapper.wrapper
current_mapset = gscript.gisenv()['MAPSET'] gscript.list_pairs('raster', mapset=current_mapset)
PyGRASS is a library originally developed during the Google Summer of Code 2012. PyGRASS library adds two main functionalities:
- Python interface through the ctypes binding of the C API of GRASS, to read and write natively GRASS GIS 7 data structures,
- GRASS GIS module interface using objects to check the parameters and execute the respective modules.
For further discussion about the implementation ideas and performance are presented in the article: Zambelli, P.; Gebbert, S.; Ciolli, M. Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS). ISPRS Int. J. Geo-Inf. 2013, 2, 201-219.
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. Using the C API (most of GRASS GIS is written in C), 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.
Working with vector data (see manual page)
Create a new vector map. Import the necessary classes:
from grass.pygrass.vector import VectorTopo from grass.pygrass.vector.geometry import Point
Create an instance of a vector map:
my_points = VectorTopo('my_points')
Open the map in write mode:
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:
Finally close the vector map:
Display the newly created vector (from GUI Python Shell):
Now do the same thing using the context manager syntax and set also the attribute table:
# Define the columns of the new vector map cols = [(u'cat', 'INTEGER PRIMARY KEY'), (u'name', 'TEXT')] with VectorTopo('my_points', mode='w', tab_cols=cols, overwrite=True) as my_points: # save the point and the attribute my_points.write(point1, ('pub', )) my_points.write(point2, ('restaurant', )) # save the changes to the database my_points.table.conn.commit()
Note: we don't need to close the vector map because it is already closed by the context manager.
Go to Layer Manager, right click on 'my_points' and Show attribute data.
Read an existing vector map:
with VectorTopo('my_points', mode='r') as points: for point in points: print(point, point.attrs['name'])
Tangible Landscape is a collaborative modeling environment for analysis of terrain changes. We couple a scanner, projector and a physical 3D model with GRASS GIS. We can analyze the impact of terrain changes by capturing the changes on the model, bringing them into the GIS, performing desired analysis or simulation and projecting the results back on the model in real-time. Tangible Landscape, as an easy-to-use 3D sketching tool, enables rapid design and scenarios testing for people with different backgrounds and computer knowledge, as well as support for decision-making process.
Tangible Landscape has been used for a variety of applications supported by the extensive set of geospatial analysis and modeling tools available in GRASS GIS. We have explored how dune breaches affect the extent of coastal flooding, the impact of different building configurations on cast shadows and solar energy potential, and the effectiveness of various landscape designs for controlling runoff and erosion. Have a look at our book, Youtube channel, and Google+ photos.
The software is free and open source and you can download it and use it in your own setup of Tangible Landscape.
Basic example how we can compute topographic slope in degrees:
import grass.script as gscript def run_slope(scanned_elev, env, **kwargs): gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', env=env) if __name__ == '__main__': elevation = 'elev_lid792_1m' env = None run_slope(scanned_elev=elevation, env=env)
We can compute slope in percent and reclassify the values into discrete intervals:
def run_slope(scanned_elev, env, **kwargs): gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', format='percent', env=env) # reclassify using rules passed as a string to standard input, 0:2:1 means reclassify interval 0 to 2 percent of slope to category 1 rules = ['0:2:1', '2:5:2', '5:8:3', '8:15:4', '15:30:5', '30:*:6'] gscript.write_command('r.recode', input='slope', output='slope_class', rules='-', stdin='\n'.join[rules], env=env) # set new color table: green - yellow - red gscript.run_command('r.colors', map='slope_class', color='gyr', env=env)
Task: Compute topographic aspect (slope orientation) usingand reclassify it into 8 main directions.
Explore concavity/convexity of terrain with curvatures computed by:
def run_curvatures(scanned_elev, env, **kwargs): gscript.run_command('r.param.scale', input=scanned_elev, output='profile_curv', method='profc', size=10, env=env) gscript.run_command('r.param.scale', input=scanned_elev, output='tangential_curv', method='crosc', size=10, env=env)
Change window size and see how it influences the result.
There are 2 different methods for landform (ridge, valley, ...) identification in GRASS GIS. First one is implemented inand uses curvatures. The second is implemented in an addon which uses multiscale line-of-sight approach. This addon can be installed from GUI command line:
or follow this tutorial.
def run_curvatures(scanned_elev, env, **kwargs): gscript.run_command('r.param.scale', input=scanned_elev, output='landforms1', method='feature', size=10, env=env) gscript.run_command('r.geomorphon', dem=scanned_elev, forms='landforms2', search=16, skip=6, env=env)
Building model based on difference raster
Modulefills a lake to a target water level from a given start point or seed raster. The resulting raster map contains cells with values representing lake depth and NULL for all other cells beyond the lake.
Example showing basic usage:
import grass.script as gscript def run_lake(scanned_elev, env, **kwargs): coordinates = [638830, 220150] gscript.run_command('r.lake', elevation=scanned_elev, lake='output_lake', coordinates=coordinates, water_level=115, env=env) if __name__ == '__main__': elevation = 'elev_lid792_1m' env = None run_lake(scanned_elev=elevation, env=env)
Create a lake where water level is relative to the altitude of the seed cell. Use function raster_what to obtain the elevation value of the DEM and create a lake with water level being for example 5 m higher.
def run_lake(scanned_elev, env, **kwargs): coordinates = [638830, 220150] res = gscript.raster_what(map=scanned_elev, coord=[coordinates]) elev_value = res[scanned_elev]['value'] gscript.run_command('r.lake', elevation=scanned_elev, lake='output_lake', coordinates=coordinates, water_level=elev_value + 5, env=env)
Flow accumulation and watersheds
We can derive watersheds and flow accumulation using modulein one command:
def run_hydro(scanned_elev, env, **kwargs): gscript.run_command('r.watershed', elevation=scanned_elev, accumulation='flow_accum', basin='watersheds', threshold=1000, env=env)
Compute average slope value for each watershed. To do that use zonal statistics using:
def run_watershed_slope(scanned_elev, env, **kwargs): gscript.run_command('r.watershed', elevation=scanned_elev, accumulation='flow_accum', basin='watersheds', threshold=1000, env=env) gscript.run_command('r.slope', elevation=scanned_elev, slope='slope', env=env) gscript.run_command('r.stats.zonal', base='watersheds', cover='slope', method='average', output='watersheds_slope', env=env) gscript.run_command('r.colors', map='watersheds_slope', color='bgyr', env=env)