Difference between revisions of "Unleash the power of GRASS GIS at US-IALE 2017"

From GRASS-Wiki
Jump to navigation Jump to search
(→‎Scripting with R: session info inspired by R statistics/rgrass7 and see also for R and Python)
(lower headings level by one, follows Wikipedia standard)
Line 4: Line 4:
This is material for [[GRASS GIS at US-IALE 2017 Annual Meeting|US-IALE 2017]] workshop ''Unleash the power of GRASS GIS'' held in Baltimore April 11, 2017. This workshop introduces GRASS GIS and processing capabilities relevant to landscape ecology.  
This is material for [[GRASS GIS at US-IALE 2017 Annual Meeting|US-IALE 2017]] workshop ''Unleash the power of GRASS GIS'' held in Baltimore April 11, 2017. This workshop introduces GRASS GIS and processing capabilities relevant to landscape ecology.  


= GRASS GIS introduction =
== GRASS GIS introduction ==
Here we provide an overview of the GRASS GIS project: [https://grass.osgeo.org grass.osgeo.org] that might be helpful to review if you are a first time user. For this exercise it's not necessary to have a full understanding of how to use GRASS GIS. However, you will need to know how to place your data in the correct GRASS GIS database directory, as well as some basic GRASS functionality. Here we introduce main concepts necessary for running the tutorial:
Here we provide an overview of the GRASS GIS project: [https://grass.osgeo.org grass.osgeo.org] that might be helpful to review if you are a first time user. For this exercise it's not necessary to have a full understanding of how to use GRASS GIS. However, you will need to know how to place your data in the correct GRASS GIS database directory, as well as some basic GRASS functionality. Here we introduce main concepts necessary for running the tutorial:
== Setting up GRASS for the tutorial ==
== Setting up GRASS for the tutorial ==
GRASS uses unique database terminology and structure ([https://grass.osgeo.org/grass71/manuals/grass_database.html GRASS database]) that are important to understand for the set up of this tutorial, as you will need to place the required data (e.g. Location) in a specific GRASS database. In the following we review important terminology and give step by step directions on how to download and place you data in the correct location.
GRASS uses unique database terminology and structure ([https://grass.osgeo.org/grass71/manuals/grass_database.html GRASS database]) that are important to understand for the set up of this tutorial, as you will need to place the required data (e.g. Location) in a specific GRASS database. In the following we review important terminology and give step by step directions on how to download and place you data in the correct location.


=== Structure of the GRASS GIS Spatial Database ===
==== Structure of the GRASS GIS Spatial Database ====


* A '''GRASS GIS Spatial Database''' (''GRASS database'') consists of directory with specific Locations (projects) where data (data layers/maps) are stored.
* A '''GRASS GIS Spatial Database''' (''GRASS database'') consists of directory with specific Locations (projects) where data (data layers/maps) are stored.
Line 17: Line 17:
[[File:Grass database.png]]
[[File:Grass database.png]]


=== Creating a GRASS database for the tutorial ===
==== Creating a GRASS database for the tutorial ====
You need to create a GRASS database with the Mapset that we will use for the tutorial before we can run the FUTURES model. Please download the [http://fatra.cnr.ncsu.edu/futures/futures_Asheville_ncspm.zip GRASS Location] 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>futures_ncspm</tt> in <tt>grassdata</tt>.
You need to create a GRASS database with the Mapset that we will use for the tutorial before we can run the FUTURES model. Please download the [http://fatra.cnr.ncsu.edu/futures/futures_Asheville_ncspm.zip GRASS Location] 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>futures_ncspm</tt> in <tt>grassdata</tt>.


<center><gallery perrow=1 widths=500 heights=250>Image:GRASS FUTURES startup.png|GRASS GIS 7.0.3 startup dialog with downloaded Location and Mapsets for FUTURES workshop</gallery></center>
<center><gallery perrow=1 widths=500 heights=250>Image:GRASS FUTURES startup.png|GRASS GIS 7.0.3 startup dialog with downloaded Location and Mapsets for FUTURES workshop</gallery></center>


=== Displaying and exploring data ===
==== Displaying and exploring data ====


Now that we have the data in the correct GRASS database, we can launch the Graphical User Interface (GUI) in Mapset ''practice1''.
Now that we have the data in the correct GRASS database, we can launch the Graphical User Interface (GUI) in Mapset ''practice1''.
Line 35: Line 35:
</gallery></center>
</gallery></center>


== GRASS GIS modules ==
=== GRASS GIS modules ===


One of the advantages of GRASS is the diversity and number of modules that let you analyze all manner of spatial and temporal. GRASS GIS has over [https://grass.osgeo.org/grass70/manuals/full_index.html 500 different modules] in the core distribution and over [https://grass.osgeo.org/grass70/manuals/addons/ 230 addon modules] that can be used to prepare and analyze data layers.
One of the advantages of GRASS is the diversity and number of modules that let you analyze all manner of spatial and temporal. GRASS GIS has over [https://grass.osgeo.org/grass70/manuals/full_index.html 500 different modules] in the core distribution and over [https://grass.osgeo.org/grass70/manuals/addons/ 230 addon modules] that can be used to prepare and analyze data layers.
Line 67: Line 67:
The name of the module helps to understand its function, for example ''v.in.lidar'' starts with ''v'' so it deals with vector maps, the name follows with ''in'' which indicates that the module is for importing the data into GRASS GIS Spatial Database and finally ''lidar'' indicates that it deals with lidar point clouds.
The name of the module helps to understand its function, for example ''v.in.lidar'' starts with ''v'' so it deals with vector maps, the name follows with ''in'' which indicates that the module is for importing the data into GRASS GIS Spatial Database and finally ''lidar'' indicates that it deals with lidar point clouds.


=== Finding and running a module ===
==== Finding and running a module ====


To find a module for your analysis, type the term into the search box into the ''Modules'' tab in the ''Layer Manager'', then keep pressing Enter until you find your module.  
To find a module for your analysis, type the term into the search box into the ''Modules'' tab in the ''Layer Manager'', then keep pressing Enter until you find your module.  
Line 78: Line 78:
</gallery></center>
</gallery></center>


=== Running a module as a command ===
==== Running a module as a command ====


If you already know the name of the module, you can just use it in the command line. The GUI offers a Command console tab with command line specifically build for running GRASS GIS modules. If you type module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get GUI dialog for the module.
If you already know the name of the module, you can just use it in the command line. The GUI offers a Command console tab with command line specifically build for running GRASS GIS modules. If you type module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get GUI dialog for the module.
Line 88: Line 88:
You can use the command line to run also whole commands for example when you get a command, i.e. module and list of parameters, in the instructions.
You can use the command line to run also whole commands for example when you get a command, i.e. module and list of parameters, in the instructions.


=== Command line vs. GUI interface===
==== Command line vs. GUI interface====
GRASS modules can be executed either through a GUI or command line interface. The GUI offers a user-friendly approach to executing modules where the user can navigate to data layers that they would like to analyze and modify processing options with simple check boxes. The GUI also offers an easily accessible manual on how to execute a model. The command line interface allows users to execute a module using command prompts specific to that module. This is handy when you are running similar analyses with minor modification or are familiar with the module commands for quick efficient processing. In this workshop we provide module prompts that can be copy and pasted into the command line for our workflow, but you can use both GUI and command line depending on personal preference. Look how {{cmd|wxGUI.modules|version=71|desc=GUI and command line interface}} represent the same tool.<br />''Task:'' compute aspect (orientation) from provided digital elevation model using module {{cmd|r.slope.aspect}} using both module dialog and command line.
GRASS modules can be executed either through a GUI or command line interface. The GUI offers a user-friendly approach to executing modules where the user can navigate to data layers that they would like to analyze and modify processing options with simple check boxes. The GUI also offers an easily accessible manual on how to execute a model. The command line interface allows users to execute a module using command prompts specific to that module. This is handy when you are running similar analyses with minor modification or are familiar with the module commands for quick efficient processing. In this workshop we provide module prompts that can be copy and pasted into the command line for our workflow, but you can use both GUI and command line depending on personal preference. Look how {{cmd|wxGUI.modules|version=71|desc=GUI and command line interface}} represent the same tool.<br />''Task:'' compute aspect (orientation) from provided digital elevation model using module {{cmd|r.slope.aspect}} using both module dialog and command line.
* How to find modules? Modules are organized by their functionality in wxGUI menu, or we can search for them in Search modules tab. If we already know which module to use, we can just type it in the wxGUI command console.
* How to find modules? Modules are organized by their functionality in wxGUI menu, or we can search for them in Search modules tab. If we already know which module to use, we can just type it in the wxGUI command console.


=== Module parameters ===
==== Module parameters ====


<center><gallery perrow=3 widths=380 heights=190>
<center><gallery perrow=3 widths=380 heights=190>
Line 104: Line 104:
Conversely, you can fill the GUI dialog parameter by parameter when you have the command.
Conversely, you can fill the GUI dialog parameter by parameter when you have the command.


== Computational region ==
=== Computational region ===


Before we use a module to compute a new raster map, we must set properly computational region. All raster computations will be performed in the specified extent and with the given resolution.
Before we use a module to compute a new raster map, we must set properly computational region. All raster computations will be performed in the specified extent and with the given resolution.
Line 151: Line 151:
  g.region vector=lakes res=10 -a -p
  g.region vector=lakes res=10 -a -p


== Running modules ==
=== Running modules ===


Find the module for computing slope and aspect in menu or the module tree under ''Raster → Terrain analysis → Slope'' and aspect or simply run {{cmd|r.slope.aspect}}.
Find the module for computing slope and aspect in menu or the module tree under ''Raster → Terrain analysis → Slope'' and aspect or simply run {{cmd|r.slope.aspect}}.
Line 162: Line 162:
</gallery></center>
</gallery></center>


== 3D view ==
=== 3D view ===
We can explore our study area in 3D view.
We can explore our study area in 3D view.
# Add <tt>elevation_30m</tt> and uncheck or remove any other layers.
# Add <tt>elevation_30m</tt> and uncheck or remove any other layers.
Line 171: Line 171:
[[File:WxGUI 3Dview.png|700px|thumb|center| Land cover 2011 draped over elevation]]
[[File:WxGUI 3Dview.png|700px|thumb|center| Land cover 2011 draped over elevation]]


= Raster and vector analysis =
== Raster and vector analysis ==


== Distance from forest edge ==
=== Distance from forest edge ===


Use raster map algebra to extract just the given forest class (here 5) from the land classification raster:
Use raster map algebra to extract just the given forest class (here 5) from the land classification raster:
Line 194: Line 194:
</gallery></center>
</gallery></center>


== Point statistics ==
=== Point statistics ===


=== Importing a Shapefile ===
==== Importing a Shapefile ====


Download sample Shapefile [http://ncsu-geoforall-lab.github.io/grass-intro-workshop/files/points_of_interest.zip points_of_interest.zip] and unzip it.
Download sample Shapefile [http://ncsu-geoforall-lab.github.io/grass-intro-workshop/files/points_of_interest.zip points_of_interest.zip] and unzip it.
Line 204: Line 204:
  v.in.ogr input=/path/to/points_of_interest.shp output=points_of_interest
  v.in.ogr input=/path/to/points_of_interest.shp output=points_of_interest


=== Generating a hexagonal grid ===
==== Generating a hexagonal grid ====


To compute point density in a hexagonal grid for the vector map points_of_interest use the vector map itself to set extent of the computational region. The resolution is based on the desired size of hexagons.
To compute point density in a hexagonal grid for the vector map points_of_interest use the vector map itself to set extent of the computational region. The resolution is based on the desired size of hexagons.
Line 214: Line 214:
  v.mkgrid map=hexagons -h
  v.mkgrid map=hexagons -h


=== Computing statistics of points in polygons ===
==== Computing statistics of points in polygons ====


The following counts the number of points per hexagon using the v.vect.stats module.
The following counts the number of points per hexagon using the v.vect.stats module.
Line 228: Line 228:
</gallery></center>
</gallery></center>


= Landscape structure analysis =
== Landscape structure analysis ==


= Lidar data processing =
== Lidar data processing ==


= Spatio-temporal data handling and visualization =
== Spatio-temporal data handling and visualization ==


= Scripting with Python =
== Scripting with Python ==


The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python editor integrated in GRASS GIS accessible from the toolbar or the Python tab in the Layer Manager. Another option is to write the Python code in your favorite plain text editor like Notepad++ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu File -> Launch script.
The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python editor integrated in GRASS GIS accessible from the toolbar or the Python tab in the Layer Manager. Another option is to write the Python code in your favorite plain text editor like Notepad++ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu File -> Launch script.
Line 253: Line 253:
Besides, this library provides several wrapper functions for often called modules.
Besides, this library provides several wrapper functions for often called modules.


=== Calling GRASS GIS 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.
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.


Line 277: Line 277:
</source>
</source>


=== Calling GRASS GIS modules with textual input or output ===
==== Calling GRASS GIS modules with textual input or output ====
Textual output from modules can be captured using the read_command() function.
Textual output from modules can be captured using the read_command() function.
<source lang="python">
<source lang="python">
Line 306: Line 306:
</source>
</source>


=== Convenient wrapper functions ===
==== Convenient wrapper functions ====
Some modules have wrapper functions to simplify frequent tasks.
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}},
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}},
Line 341: Line 341:
</source>
</source>


=== Exercise ===
==== Exercise ====
Export all raster layers from your mapset with a name prefix "elev_*" as GeoTiff (see {{cmd|r.out.gdal}}). Don't forget to set the current region ({{cmd|g.region}}) for each map in order to match the individual exported raster layer extents and resolutions since they may differ from each other.
Export all raster layers from your mapset with a name prefix "elev_*" as GeoTiff (see {{cmd|r.out.gdal}}). Don't forget to set the current region ({{cmd|g.region}}) for each map in order to match the individual exported raster layer extents and resolutions since they may differ from each other.


= Scripting with R =
== Scripting with R ==


Using R and GRASS GIS together can be done in two ways:
Using R and GRASS GIS together can be done in two ways:
Line 355: Line 355:
** Use the <code>initGRASS()</code> function to start GRASS GIS inside R.
** Use the <code>initGRASS()</code> function to start GRASS GIS inside R.


= See also =
== See also ==


* [[GRASS and Python]]
* [[GRASS and Python]]

Revision as of 19:21, 30 March 2017

Us-iale logo.jpg
Grassgis logo colorlogo text whitebg.png

This is material for US-IALE 2017 workshop Unleash the power of GRASS GIS held in Baltimore April 11, 2017. This workshop introduces GRASS GIS and processing capabilities relevant to landscape ecology.

GRASS GIS introduction

Here we provide an overview of the GRASS GIS project: grass.osgeo.org that might be helpful to review if you are a first time user. For this exercise it's not necessary to have a full understanding of how to use GRASS GIS. However, you will need to know how to place your data in the correct GRASS GIS database directory, as well as some basic GRASS functionality. Here we introduce main concepts necessary for running the tutorial:

Setting up GRASS for the tutorial

GRASS uses unique database terminology and structure (GRASS database) that are important to understand for the set up of this tutorial, as you will need to place the required data (e.g. Location) in a specific GRASS database. In the following we review important terminology and give step by step directions on how to download and place you data in the correct location.

Structure of the GRASS GIS Spatial Database

  • A GRASS GIS Spatial Database (GRASS database) consists of directory with specific Locations (projects) where data (data layers/maps) are stored.
  • Location is a directory with data related to one geographic location or a project. All data within one Location has the same coordinate reference system.
  • Mapset is a collection of maps within Location, containing data related to a specific task, user or a smaller project.

Grass database.png

Creating a GRASS database for the tutorial

You need to create a GRASS database with the Mapset that we will use for the tutorial before we can run the FUTURES model. Please download the GRASS Location 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 futures_ncspm in grassdata.

Displaying and exploring data

Now that we have the data in the correct GRASS database, we can launch the Graphical User Interface (GUI) in Mapset practice1.

The GUI interface allows you to display raster, vector data as well as navigate through zooming in and out. More advanced exploration and visualization is also possible using, e.g., queries and adding legend. The screenshots below depicts how you can add different map layers (left) and display the metadata of your data layers.

GRASS GIS modules

One of the advantages of GRASS is the diversity and number of modules that let you analyze all manner of spatial and temporal. GRASS GIS has over 500 different modules in the core distribution and over 230 addon modules that can be used to prepare and analyze data layers.

GRASS functionality is available through modules (tools, functions). Modules respect the following naming conventions:

Prefix Function Example
r.* raster processing r.mapcalc: map algebra
v.* vector processing v.clean: topological cleaning
i.* imagery processing i.segment: object recognition
db.* database management db.select: select values from table
r3.* 3D raster processing r3.stats: 3D raster statistics
t.* temporal data processing t.rast.aggregate: temporal aggregation
g.* general data management g.rename: renames map
d.* display d.rast: display raster map

These are the main groups of modules. There is few more for specific purposes. Note also that some modules have multiple dots in their names. This often suggests further grouping. For example, modules staring with v.net. deal with vector network analysis.

The name of the module helps to understand its function, for example v.in.lidar starts with v so it deals with vector maps, the name follows with in which indicates that the module is for importing the data into GRASS GIS Spatial Database and finally lidar indicates that it deals with lidar point clouds.

Finding and running a module

To find a module for your analysis, type the term into the search box into the Modules tab in the Layer Manager, then keep pressing Enter until you find your module.

Alternatively, you can just browse through the module tree in the Modules tab. You can also browse through the main menu. For example, to find information about a raster map, use: Raster → Reports and statistics → Basic raster metadata.

Running a module as a command

If you already know the name of the module, you can just use it in the command line. The GUI offers a Command console tab with command line specifically build for running GRASS GIS modules. If you type module name there, you will get suggestions for automatic completion of the name. After pressing Enter, you will get GUI dialog for the module.

You can use the command line to run also whole commands for example when you get a command, i.e. module and list of parameters, in the instructions.

Command line vs. GUI interface

GRASS modules can be executed either through a GUI or command line interface. The GUI offers a user-friendly approach to executing modules where the user can navigate to data layers that they would like to analyze and modify processing options with simple check boxes. The GUI also offers an easily accessible manual on how to execute a model. The command line interface allows users to execute a module using command prompts specific to that module. This is handy when you are running similar analyses with minor modification or are familiar with the module commands for quick efficient processing. In this workshop we provide module prompts that can be copy and pasted into the command line for our workflow, but you can use both GUI and command line depending on personal preference. Look how GUI and command line interface represent the same tool.
Task: compute aspect (orientation) from provided digital elevation model using module r.slope.aspect using both module dialog and command line.

  • How to find modules? Modules are organized by their functionality in wxGUI menu, or we can search for them in Search modules tab. If we already know which module to use, we can just type it in the wxGUI command console.

Module parameters

The same analysis can be done using the following command:

r.neighbors -c input=elevation output=elev_smooth size=5

Conversely, you can fill the GUI dialog parameter by parameter when you have the command.

Computational region

Before we use a module to compute a new raster map, we must set properly computational region. All raster computations will be performed in the specified extent and with the given resolution.

Computational region is an important raster concept in GRASS GIS. In GRASS a computational region can be set, subsetting larger extent data for quicker testing of analysis or analysis of specific regions based on administrative units. We provide a few points to keep in mind when using the computational region function:

  • defined by region extent and raster resolution
  • applies to all raster operations
  • persists between GRASS sessions, can be different for different mapsets
  • advantages: keeps your results consistent, avoid clipping, for computationally demanding tasks set region to smaller extent, check your result is good and then set the computational region to the entire study area and rerun analysis
  • run g.region -p or in menu Settings - Region - Display region to see current region settings

The numeric values of computational region can be checked using:

g.region -p

After executing the command you will get something like this:

north:      220750
south:      220000
west:       638300
east:       639000
nsres:      1
ewres:      1
rows:       750
cols:       700
cells:      525000

Computational region can be set also using a vector map. In that case, only extent is set (as vector map does not have any resolution - at least not in the way raster map does). In GUI, this can be done in the same way as for the raster map. In the command line, it looks like this:

g.region vector=lakes

Resolution can be set separately using the res parameter of the g.region module. The units are the units of the current location, in our case meters. This can be done in the Resolution tab of the g.region dialog or in the command line in the following way (using also the -a flag to print the new values):

g.region res=3 -p

The new resolution may be slightly modified in this case to fit into the extent which we are not changing. However, often we want the resolution to be the exact value we provide and we are fine with a slight modification of the extent. That's what -a flag is for.

The following example command will use the extent from the vector named lakes, use resolution 10, modify the extent to align it to this 10 meter resolution, and print the values of this new computational region settings:

g.region vector=lakes res=10 -a -p

Running modules

Find the module for computing slope and aspect in menu or the module tree under Raster → Terrain analysis → Slope and aspect or simply run r.slope.aspect.

3D view

We can explore our study area in 3D view.

  1. Add elevation_30m and uncheck or remove any other layers.
  2. Zoom to an area around Asheville and in Map Display select Various zoom options - Set computational region extent from display. Switch to 3D view (in the right corner on Map Display).
  3. Adjust the view (perspective, height, vertical exaggeration)
  4. In Data tab, set Fine mode resolution to 1 and set landuse_2011 as the color of the surface.
  5. When finished, switch back to 2D view.
Land cover 2011 draped over elevation

Raster and vector analysis

Distance from forest edge

Use raster map algebra to extract just the given forest class (here 5) from the land classification raster:

r.mapcalc "forest = if(landclass96 == 5, 1, null())"

The if() function we used has three parameters with the following syntax:

if(condition, value used when it is true, value used when it is false)

Then we used operator == which evaluates as true when both sides are equal. Finally we used null() function which represents NULL (no data) value.

Now we can get distance to the edge of the forest using r.grow.distance module which computes distances to areas with values in areas without values (with NULLs) or the other way around. By default it would give us distance to the edge of the forest from outside of the forest, but we are now using the -n flag to obtain distance to the edge from within the forest itself:

r.grow.distance -n input=forest distance=distance

Point statistics

Importing a Shapefile

Download sample Shapefile points_of_interest.zip and unzip it.

Import the file using v.in.ogr module. Note that you need to specify the full path to the file.

v.in.ogr input=/path/to/points_of_interest.shp output=points_of_interest

Generating a hexagonal grid

To compute point density in a hexagonal grid for the vector map points_of_interest use the vector map itself to set extent of the computational region. The resolution is based on the desired size of hexagons.

g.region vector=points_of_interest res=2000 -pa

Although computation region is usually not used in vector processing, the hexagonal grid is created as a vector map based on the previously selected extent and size of the grid.

v.mkgrid map=hexagons -h

Computing statistics of points in polygons

The following counts the number of points per hexagon using the v.vect.stats module.

v.vect.stats points=points_of_interest areas=hexagons count_column=count

The last command sets the vector map color table to viridis based on the count column. Use color table ryb if you have GRASS GIS 7.0.

v.colors map=hexagons use=attr column=count color=viridis

Landscape structure analysis

Lidar data processing

Spatio-temporal data handling and visualization

Scripting with Python

The simplest way to execute the Python code which uses GRASS GIS packages is to use Simple Python editor integrated in GRASS GIS accessible from the toolbar or the Python tab in the Layer Manager. Another option is to write the Python code in your favorite plain text editor like Notepad++ (note that Python editors are plain text editors). Then run the script in GRASS GIS using the main menu File -> Launch script.



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 g.region. In this example, 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')

If we run the Python commands from GUI Python console, we can use AddLayer to add the newly created layer:

AddLayer('elev_smoothed')

Calling GRASS GIS modules with textual input or output

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

gscript.read_command('g.region', flags='p')
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.

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, 221342), output='point')

If we run the Python commands from GUI Python console, we can use AddLayer to add the newly created layer:

AddLayer('point')

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 r.info, or a vector layer with vector_info.

gscript.raster_info('elevation')
gscript.vector_info('point')

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 g.list wrappers. With list_grouped, the map layers are grouped by mapsets (in this example, raster layers):

gscript.list_grouped(type=['raster'])
gscript.list_grouped(type=['raster'], pattern="landuse*")

Here is an example of a different g.list wrapper list_pairs 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)

Exercise

Export all raster layers from your mapset with a name prefix "elev_*" as GeoTiff (see r.out.gdal). Don't forget to set the current region (g.region) for each map in order to match the individual exported raster layer extents and resolutions since they may differ from each other.

Scripting with R

Using R and GRASS GIS together can be done in two ways:

  • Using R within GRASS GIS session, i.e. you start R (or RStudio) from the GRASS GIS command line.
    • You work with data in GRASS GIS Spatial Database using GRASS GIS
    • Do not use the initGRASS() function (GRASS GIS is running already).
  • Using GRASS GIS within a R session, i.e. you connect to a GRASS GIS Spatial Database from within R (or RStudio).
    • You put data into GRASS GIS Spatial Database just to perform the GRASS GIS computations.
    • Use the initGRASS() function to start GRASS GIS inside R.

See also