Difference between revisions of "Unleash the power of GRASS GIS at US-IALE 2017"
Wenzeslaus (talk | contribs)
(→Scripting with Python: content from Using GRASS GIS through Python and tangible interfaces (workshop at FOSS4G NA 2016))
Wenzeslaus (talk | contribs)
(→Scripting with R: session info inspired by R statistics/rgrass7 and see also for R and Python)
|Line 345:||Line 345:|
= Scripting with R =
= Scripting with R =
Revision as of 19:55, 28 March 2017
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.
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:
|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.
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.
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
g.region -por in menu Settings - Region - Display region to see current region settings
Computational region concept: A raster with large extent (blue) is displayed as well as another raster with smaller extent (green). The computational region (red) is now set to match the smaller raster, so all the computations are limited to the smaller raster extent even if the input is the larger raster. (Not shown on the image: Also the resolution, not only the extent, matches the resolution of the smaller raster.)
The numeric values of computational region can be checked using:
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:
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
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.
We can explore our study area in 3D view.
- Add elevation_30m and uncheck or remove any other layers.
- 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).
- Adjust the view (perspective, height, vertical exaggeration)
- In Data tab, set Fine mode resolution to 1 and set landuse_2011 as the color of the surface.
- When finished, switch back to 2D view.
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())"
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
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.
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('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.
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 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']) 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.
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.