View source for Using GRASS GIS through Python and tangible interfaces (workshop at FOSS4G NA 2016)

{{ToModify}}
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 [https://geospatial.ncsu.edu/osgeorel/tangible-landscape.html 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.
[https://geospatial.ncsu.edu/osgeorel/tangible-landscape.html 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''': [https://grass.osgeo.org GRASS GIS 7]

'''Data''': Download [https://grass.osgeo.org/sampledata/north_carolina/nc_spm_08_grass7.zip complete North Carolina sample dataset].

= GRASS GIS introduction =
You need to create a GRASS database we will use for the tutorial. Please download the [https://grass.osgeo.org/sampledata/north_carolina/nc_spm_08_grass7.zip dataset] for the workshop, noting where the files are located on your local directory. 
Now, create (unless you already have it) a directory named <tt>grassdata</tt> (GRASS database) in your home folder (or Documents), unzip the downloaded data into this directory. You should now have a Location <tt>nc_spm_08_grass7</tt> in grassdata. Start GRASS GIS in Mapset <tt>user1</tt>.

= GRASS GIS Python API =

== Python Scripting library ==
[[File:GRASS GUI Python shell.png|thumbnail|right|Python shell in GRASS GIS GUI]]
The GRASS GIS 7 Python Scripting Library provides functions to call GRASS modules within scripts as subprocesses. The most often used functions include:
* '''[https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.run_command run_command]''': most often used with modules which output raster/vector data where text output is not expected
* '''[https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.read_command read_command]''': used when we are interested in text output
* '''[https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.parse_command parse_command]''': used with modules producing text output as key=value pair
* '''[https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.write_command 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 <tt>python myscript.py</tt> 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:
<source lang="python">
import grass.script as gscript
</source>

Before running any GRASS raster modules, you need to set the computational region using {{cmd|g.region}}. In this example, we set the computational extent and resolution to the raster layer <tt>elevation</tt>.
<source lang="python">
gscript.run_command('g.region', raster='elevation')
</source>

The run_command() function is the most commonly used one. Here, we apply the focal operation average ({{cmd|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.
<source lang="python">
gscript.run_command('r.neighbors', input='elevation', output='elev_smoothed', method='average', flags='c')
</source>
If we run the Python commands from GUI Python console, we can use <tt>AddLayer</tt> to add the newly created layer:
<source lang="python">
AddLayer('elev_smoothed')
</source>

=== Calling GRASS GIS modules with textual input or output ===
Textual output from modules can be captured using the read_command() function.
<source lang="python">
gscript.read_command('g.region', flags='p')
</source>
<source lang="python">
gscript.read_command('r.univar', map='elev_smoothed', flags='g')
</source>

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 {{cmd|g.proj}} to display the projection parameters of the actual location.
<source lang="python">
gscript.parse_command('g.proj', flags='g')
</source>

For comparison, below is the same example, but using the read_command() function.
<source lang="python">
gscript.read_command('g.proj', flags='g')
</source>

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 {{cmd|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.
<source lang="python">
gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (635818, 221342), output='point')
</source>

If we run the Python commands from GUI Python console, we can use <tt>AddLayer</tt> to add the newly created layer:
<source lang="python">
AddLayer('point')
</source>

=== Convenient wrapper functions ===
Some modules have wrapper functions to simplify frequent tasks.
For example we can obtain the information about a raster layer with [https://grass.osgeo.org/grass70/manuals/libpython/script.html?highlight=mapcalc#script.raster.raster_info raster_info] which is a wrapper of {{cmd|r.info}},
or a vector layer with [https://grass.osgeo.org/grass70/manuals/libpython/script.html?script.vector.vector_info vector_info].
<source lang="python">
gscript.raster_info('elevation')
gscript.vector_info('point')
</source>

Another example is using [https://grass.osgeo.org/grass70/manuals/libpython/script.html?highlight=mapcalc#script.raster.mapcalc r.mapcalc wrapper] for raster algebra:
<source lang="python">
gscript.mapcalc("elev_strip = if(elevation > 100 && elevation < 125, elevation, null())")
gscript.read_command('r.univar', map='elev_strip', flags='g')
</source>

Function [https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.region 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). 
<source lang="python">
region = gscript.region()
print region
# cell area in map units (in projected Locations)
region['nsres'] * region['ewres']
</source>

We can list data stored in a GRASS GIS location with {{cmd|g.list}} wrappers. With [https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.list_grouped list_grouped], the map layers are grouped by mapsets (in this example, raster layers):
<source lang="python">
gscript.list_grouped(type=['raster'])
gscript.list_grouped(type=['raster'], pattern="landuse*")
</source>

Here is an example of a different {{cmd|g.list}} wrapper [https://grass.osgeo.org/grass70/manuals/libpython/script.html#script.core.list_pairs list_pairs] which structures the output as list of pairs (name, mapset). We obtain current mapset with {{cmd|g.gisenv}} wrapper.
<source lang="python">
current_mapset = gscript.gisenv()['MAPSET']
gscript.list_pairs('raster', mapset=current_mapset)
</source>

=== Exercise ===

== PyGRASS ==
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:
''[http://www.mdpi.com/2220-9964/2/1/201 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 [https://grass.osgeo.org/grass70/manuals/libpython/pygrass_vector.html manual page]) ===
Create a new vector map. Import the necessary classes:
<source lang="python">
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Point
</source>

Create an instance of a vector map:
<source lang="python">
my_points = VectorTopo('my_points')
</source>

Open the map in write mode:
<source lang="python">
my_points.open(mode='w')
</source>

Create some vector geometry features, like two points:
<source lang="python">
point1 = Point(635818.8, 221342.4)
point2 = Point(633627.7, 227050.2)
</source>

Add the above two points to the new vector map:
<source lang="python">
my_points.write(point1)
my_points.write(point2)
</source>

Finally close the vector map:
<source lang="python">
my_points.close()
</source>

Display the newly created vector (from GUI Python Shell):
<source lang="python">
AddLayer('my_points')
</source>

Now do the same thing using the context manager syntax and set also the attribute table:
<source lang="python">
# 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()
</source>
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:
<source lang="python">
with VectorTopo('my_points', mode='r') as points:
    for point in points:
        print(point, point.attrs['name'])
</source>

= Tangible Landscape =
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 [http://www.springer.com/us/book/9783319257730 book], [https://www.youtube.com/channel/UCc37pVh-WE46Xkqeq-KZQsA Youtube channel], and [https://plus.google.com/112857255995703553144/photos Google+ photos].

The software is free and open source and you can [https://github.com/ncsu-osgeorel/grass-tangeoms download it] and use it in your own setup of Tangible Landscape.


== Example: computing slope ==
<source lang="python">
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 = 'elevation'
    env = None
    run_slope(scanned_elev=elevation, env=env)
</source>

[[Category: Tutorial]]
[[Category: Python]]