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

From GRASS-Wiki
Jump to navigation Jump to search
Line 283: Line 283:

   r.forestfrag input=forest output=fragmentation window=15
   r.forestfrag input=forest output=fragmentation window=15
Report the distribution of the fragmentation categories:
  r.report map=fragmentation units=k,p

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

Revision as of 03:04, 5 April 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

For the following examples we will use raster landuse96_28m from our North Carolina dataset, where patches represent different land cover. We will use module r.neighbors and addons r.diversity and {{addonCmd|r.forestfrag}. Install the addon first:

 g.extension r.diversity
 g.extension r.forestfrag

Then close GRASS GUI and open it again.


First, we compute richness (number of unique classes) using r.neighbors. By using moving window analysis, we create a raster representing spatially variable richness. The window size is in cells.

 g.region raster=landuse96_28m
 r.neighbors input=landuse96_28m output=richness method=diversity size=15

Landscape indices

Addon r.diversity computes landscape indices using moving window. It is based on r.li modules for landscape structure analysis. In this example, we compute Simpson, Shannon, and Renyi diversity indices with a range of window sizes:

 r.diversity input=landuse96_28m prefix=index alpha=0.8 size=9,15,21 method=simpson,shannon,renyi

This generates 9 rasters with names like index_simpson_size_9. We can add them to Map Display using Add multiple raster or vector layers in Layer manager toolbar (top).

For viewing, we should set the same color table for rasters maps of the same index.

 r.colors map=index_shannon_size_21,index_shannon_size_15,index_shannon_size_9 color=viridis
 r.colors map=index_renyi_size_21_alpha_0.8,index_renyi_size_15_alpha_0.8,index_renyi_size_15_alpha_0.8 color=viridis
 # we use grey1.0 color ramp because simpson is from 0-1
 r.colors map=index_simpson_size_21,index_simpson_size_15,index_simpson_size_9 color=grey1.0

Forest fragmentation

Forest fragmentation computed with addon r.forestfrag

Module r.forestfrag computes the forest fragmentation following the methodology proposed by Riitters et al. (2000).

First mark all cells which are forest as 1 and everything else as zero:

 # first set region
 g.region raster=landclass96
 # list classes:
 r.category map=landclass96
 # select classes 
 r.mapcalc "forest = if(landclass96 == 5, 1, 0)"

Use the new forest presence raster map to compute the forest fragmentation index with window size 15:

 r.forestfrag input=forest output=fragmentation window=15

Report the distribution of the fragmentation categories:

 r.report map=fragmentation units=k,p

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:


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:


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.


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'], 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)


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.

Creating a new GRASS GIS Location and importing data

For the following example with R, we need to first create a new GRASS Location and import new data. We will use PRISM temperature data and vector boundaries of US states. We will create new Location based on EPSG code 4269 (NAD83 datum).

  1. Start GRASS GIS
  2. Select New in the left part of the welcome screen to start Location Wizard
  3. In the wizard, type name of the new Location, for example PRISM, press Next
  4. Choose Select EPSG code of spatial reference system, press Next
  5. Type 4269 in the EPSG code field, press Next
  6. Select 1 from datum transformation dialog
  7. Review PROJ.4 definition and press "Finish"

Mapset PERMANENT is automatically created, so you can start GRASS session.

Download and unzip the following datasets into a folder:

Set current working directory in GRASS GIS to simplify finding the extracted files. In menu Settings - GRASS working environment - Change working directory set the directory where you have the data. Alternatively, you can also type cd in the GUI command cosole and select the directory.

Now we can import them into GRASS GIS:

 r.import input=PRISM_tmean_30yr_normal_4kmM2_annual_bil.bil output=temp_mean
 v.import input=cb_2015_us_state_20m.shp output=boundaries
 r.import input=US_elevation.tif output=elevation

We don't specify full path to the files here thanks to setting the working directory above. If you use dialog, you can browse to the files using the Browse button.

Now display the data and set the appropriate color ramp for the temperature:

 r.colors map=temp_mean@PERMANENT color=celsius

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.

We will run R within GRASS GIS session (the first way). Launch R inside GRASS GIS and install rgrass7 and rgdal package


We can execute GRASS modules using execGRASS function:

execGRASS("g.region", raster="temp_mean", flags="p")

We will analyze the relationship between temperature and elevation and latitude.

First we will generate raster of latitude values:

execGRASS("r.mapcalc", expression="latitude = y()")

Note: in projected coordinate systems you can use r.latlong.

Then, we will generate random points and sample the datasets.

execGRASS("v.random", output="samples", npoints=1000)
# this will restrict sampling to the boundaries of USA
# we are overwriting vector samples, so we need to use overwrite flag
execGRASS("v.random", output="samples", npoints=1000, restrict="boundaries", flags=c("overwrite"))
# create attribute table
execGRASS("v.db.addtable", map="samples", columns=c("elevation double precision", "latitude double precision", "temp double precision"))
# sample individual rasters
execGRASS("v.what.rast", map="samples", raster="temp_mean", column="temp")
execGRASS("v.what.rast", map="samples", raster="latitude", column="latitude")
execGRASS("v.what.rast", map="samples", raster="elevation", column="elevation")

Note: there is an addon module r.sample.category which simplifies this process by sampling multiple rasters.

Now open GRASS GIS attribute table manager to inspect the sampled values or use v.db.select to list the values. Explore the dataset in R:

samples <- readVECT("samples")

Compute multivariate linear model:

linmodel <- lm(temp ~ elevation + latitude, samples)

Predict temperature using this model:

maps <- readRAST(c("elevation", "latitude"))
maps$temp_model <- predict(linmodel, newdata=maps)
spplot(maps, "temp_model")
# write modeled temperature to GRASS raster and set color ramp
writeRAST(maps, "temp_model", zcol="temp_model")
execGRASS("r.colors", map="temp_model", color="celsius")

Compare simple linear model to real data:

execGRASS("r.mapcalc", expression="diff = temp_mean - temp_model")
execGRASS("r.colors", map="diff", color="differences")

In GRASS GUI, add layers temp_mean and temp_model, select them and go to File - Map Swipe to compare visually the modeled and real temperatures.

See also