Working with GRASS GIS: A guided tour

display difference map

Basic overview of this session

In this session, we will go through the most commonly used spatial analysis techniques in GRASS GIS. The session is a hands-on with a common exercise which covers basic raster processing capabilities of GRASS GIS.

Course content

Data

Data download

For this exercise, we will use some GeoTIFF files and the GRASS GIS location "ecad_ll"

Workflow overview

How do we work with spatial data in GRASS GIS?

Assuming we already have our maps in the GRASSDBASE (i.e., grassdata/location/mapset), the first step is to import the datasets to GRASS GIS. We will import the elevation data required for the exercise, and the rest of the data is already imported to the GRASS GIS mapset ecad50.

Here is a list of major GRASS GIS modules we will be using in this exercise and the links to their manuals.

  1. Start page of the GRASS GIS manual
  2. Raster processing modules
  3. Vector processing modules
  4. Imagery processing modules
  5. GRASS GIS addons
  6. General commands to start GRASS GIS, set region, list data: grass72, g.region, g.list
  7. Import raster, vector data, import and project on the fly: r.in.gdal, v.in.ogr, r.import, v.import
  8. Raster map info, map calculator, univariate statistics, aggregation: r.info, r.mapcalc, r.univar, r.series,
  9. Vector map info, data base select, univariate statistics, vector to raster: v.info, v.db.select, v.univar, v.to.rast,

Spatial data analysis with GRASS GIS

Download the shared data. Copy the given GRASS GIS location to this path: $HOME/grassdata. There are two mapsets in the provided location "ecad_ll":

  1. PERMANENT
  2. ecad50 - Where all the ECAD datasets are stored

We will now create a new mapset "user1" in the same location "ecad_ll" to do all our analysis. To make use of the ecad dataset, we will add the mapset "ecad50" in the search path. We can run the commands from the command line or use the main GUI and copy the corresponding commands for future replication of the workflow. Note the "Copy" button in the GUI of each module. Maybe you might want to get help about the options and flags of the different commands, e.g. r.colors --help

        
grass72 $HOME/grassdata/ecad_ll/user1/ -c --gui
g.list rast 
g.mapsets -p
g.mapsets mapset=ecad50 operation=add
g.mapsets -p
g.list rast
        
        
# Launch GRASS GIS, -c creates new mapset user1, --gui lauches the gui
grass72 $HOME/grassdata/ecad_ll/user1/ -c --gui
# List all the raster maps in all the mapsets in the search path
g.list type=rast 
# List all the mapsets in the search path
g.mapsets -p
# Add the mapset ecad50 to the search path
g.mapsets mapset=ecad50 operation=add
# List all the mapsets in the search path
g.mapsets -p
# List all the raster maps in all the mapsets in the search path
g.list type=rast
        
    

We will now import the elevation map, provided as a GeoTIFF to the current mapset. Then set the color palette to elevation

        
cd $HOME/ecad5_geotiffs_LL
r.in.gdal input=elev_0.25deg_reg_v6.0.tif output=elev_ecad
r.colors map=elev_ecad color=elevation
        
        
# Change dirctory to the one with TIFF files
cd $HOME/ecad5_geotiffs_LL
# Import the elevation raster map in GeoTIFF format into GRASS GIS
r.in.gdal input=elev_0.25deg_reg_v6.0.tif output=elev_ecad
# Set the color palette of the elevation raster map to 'elevation'
r.colors map=elev_ecad color=elevation
        
    

Let us now display the elevation map and add decorations.

 
# Open a monitor
d.mon wx0

# Display a raster map
d.rast map=elev_ecad

# Display a vector map
d.vect map=country_boundaries type=boundary

# Add raster legend
d.legend -t -s -b raster=elev_ecad title=ELEVATION title_fontsize=20 font=sans fontsize=18

# Add North arrow
d.northarrow style=1b text_color=black

# Add text
d.text -b text="Resampled elevation map" color=black bgcolor=229:229:229 align=cc font=sans size=8
 

display DEM
Elevation in meters

Data aggregation in GRASS GIS

Let us check how to aggregate multiple maps at different time periods based on different statistical methods

Let us create the following aggregated maps:

        
g.region rast=precip_0_25deg_1951_1980_monthly_sums.1
g.region -p
g.list rast pattern="precip_0_25deg_1951_1980_monthly_sums.*"
r.series in=`g.list rast pattern="precip_0_25deg_1951_1980_monthly_sums.*" sep=","` out=precip_0_25deg.1951.1980.sum method=sum
r.colors precip_0_25deg.1951.1980.sum col=precipitation
g.list rast pattern="precip_0_25deg.1981_2010.*"
r.series in=`g.list rast pattern="precip_0_25deg.1981_2010.*" sep=","` out=precip_0_25deg.1981.2010.sum method=sum
r.colors precip_0_25deg.1981.2010.sum col=precipitation
r.series in=`g.list rast pattern="tmean_0_25deg.1981_2010*" sep=","` out=tmean_0_25deg.1981_2010.avg method=average
r.colors tmean_0_25deg.1981_2010.avg col=celsius
        
        
# Set the computational region
g.region rast=precip_0_25deg_1951_1980_monthly_sums.1
# Display the region
g.region -p
# list all monthly maps for the period 1951-1980 needed for the annual precipitation aggregation
g.list rast pattern="precip_0_25deg_1951_1980_monthly_sums.*"
# Now use the r.series command to create annual precip for the period 1951 to 1980
r.series in=`g.list rast pattern="precip_0_25deg_1951_1980_monthly_sums.*" sep=","` out=precip_0_25deg.1951.1980.sum method=sum
# Set the color
r.colors precip_0_25deg.1951.1980.sum col=precipitation
# list all monthly maps for the period 1981-2010 needed for the annual precipitation aggregation
g.list rast pattern="precip_0_25deg.1981_2010.*"
# Now use the r.series command to create annual precip for the period 1981 to 2010
r.series in=`g.list rast pattern="precip_0_25deg.1981_2010.*" sep=","` out=precip_0_25deg.1981.2010.sum method=sum
# Set the color
r.colors precip_0_25deg.1981.2010.sum col=precipitation
# Aggregate the temperature maps average annual temperature
r.series in=`g.list rast pattern="tmean_0_25deg.1981_2010*" sep=","` out=tmean_0_25deg.1981_2010.avg method=average
# Set the color
r.colors tmean_0_25deg.1981_2010.avg col=celsius
        
 

display annual mean temp
Annual Mean temperature (1981 - 2010)

Task:Now try aggregating the mean temperature maps for the period 1951 to 1980 to create annual mean temperature for the same period. Display the annual precipitation and temperature maps with legend and compare the legends using WxGUI or d.mon. For more details on r.series check here.

Univariate statistics in GRASS GIS

Now lets compute univariate statistics of the developed annual maps over the entire area and also over a particular country.

        
r.univar tmean_0_25deg.1981_2010.avg -e -g
r.univar precip_0_25deg.1981.2010.sum -e -g
v.db.select country_boundaries|head -5
v.db.select country_boundaries|grep Netherlands
r.mask vect=country_boundaries cats=118
r.univar tmean_0_25deg.1981_2010.avg -e -g
r.univar precip_0_25deg.1981.2010.sum -e -g
r.mask -r
        
        
# Compute extended univariate statistics
r.univar tmean_0_25deg.1981_2010.avg -e -g
# Compute extended univariate statistics
r.univar precip_0_25deg.1981.2010.sum -e -g
# Get a feel of the vector data and the attributes
v.db.select country_boundaries|head -5
# Find the attributes related to The Netherlands
v.db.select country_boundaries|grep Netherlands
# Mask out all the area other than The Netherlands
r.mask vect=country_boundaries cats=118
# Repeat the univariate statistics for The Netherlands
r.univar tmean_0_25deg.1981_2010.avg -e -g
# Repeat the univariate statistics for The Netherlands
r.univar precip_0_25deg.1981.2010.sum -e -g
# Remove the mask
r.mask -r
        
 

display univar stats
Univariate statistics of annual precipitation (1981 to 2010)

Raster map algebra using map calculator

Let us compute a difference map from the two 30 years aggregated precipitation maps.

        
g.region rast=precip_0_25deg.1981_2010.06.sum -p
r.mapcalc "diff_1951_1981_2010 = precip_0_25deg.1951.1980.sum - precip_0_25deg.1981.2010.sum"
r.colors -n diff_1951_1981_2010 color=differences
        
        
# we want to be sure to have the computational region set			
g.region rast=precip_0_25deg.1981_2010.06.sum -p
# Map algebra using map calculator
r.mapcalc "diff_1951_1981_2010 = precip_0_25deg.1951.1980.sum - precip_0_25deg.1981.2010.sum"
# apply inverse “differences” color table
r.colors raster=diff_1951_1981_2010 color=differences
        
 

Task:Visualize the difference map with the legend and compute univariate statistics over The Netherlands.

display difference map
Difference between precipitation (1951-1980) and (1981-2010)

Zonal statistics in GRASS GIS

We want to extract a table with basic statistics over all the countries in Europe.

Main steps are:

        
v.info -c country_boundaries
v.db.select country_boundaries column="NAME"
g.region rast=precip_0_25deg.1981.2010.sum -p
v.to.rast country_boundaries out=country_boundaries use=cat labelcolumn="NAME"
r.category country_boundaries
r.univar -t precip_0_25deg.1981.2010.sum zones=country_boundaries sep="," output=countries_precip_1981_2010_sum.csv
r.univar -t diff_1951_1981_2010 zones=country_boundaries sep="," output=countries_precip_1951_1980_vs_1981_2010_sum.csv
        
        
# check the available column (or use wxGUI attribute manager)
v.info -c country_boundaries
# look at country names
v.db.select country_boundaries column="NAME"
# we want to be sure to have the computational region set	
g.region rast=precip_0_25deg.1981.2010.sum -p
# for zonal statistics, we convert the vector polyons to raster model
v.to.rast country_boundaries out=country_boundaries use=cat labelcolumn="NAME"
# check raster polygons
r.category country_boundaries
# get the statistics for each country as a table
r.univar -t precip_0_25deg.1981.2010.sum zones=country_boundaries sep="," output=countries_precip_1981_2010_sum.csv
r.univar -t diff_1951_1981_2010 zones=country_boundaries sep="," output=countries_precip_1951_1980_vs_1981_2010_sum.csv
        
 

Cluster analysis in GRASS GIS

Let us now do climatic characterisation of Europe using unsupervised classification

Main steps are:

        
g.list rast pattern="*_0_25deg.1981_2010.*.*"
i.group group=ecad_1981_2010 subgroup=ecad_1981_2010 input=`g.list rast pattern="*_0_25deg.1981_2010.*.*" sep=","`
i.cluster group=ecad_1981_2010 subgroup=ecad_1981_2010 sig=ecad_1981_2010 classes=5 separation=0.5
i.maxlik group=ecad_1981_2010 subgroup=ecad_1981_2010 sig=ecad_1981_2010 output=ecad_1981_2010.class rej=ecad_1981_2010.class.rej
        
        
# pick the variables of interest
g.list rast pattern="*_0_25deg.1981_2010.*.*"
# add maps to an imagery group for easier management
i.group group=ecad_1981_2010 subgroup=ecad_1981_2010 input=`g.list rast pattern="*_0_25deg.1981_2010.*.*" sep=","`
# statistics for unsupervised classification
i.cluster group=ecad_1981_2010 subgroup=ecad_1981_2010 sig=ecad_1981_2010 classes=5 separation=0.5
# Maximum Likelihood unsupervised classification
i.maxlik group=ecad_1981_2010 subgroup=ecad_1981_2010 sig=ecad_1981_2010 output=ecad_1981_2010.class rej=ecad_1981_2010.class.rej
        
 

display climatic zones
Climatic Zones of Europe

Check the manual of i.cluster and i.maxlik.

Check if the result make sense comparing to the Kottek, M. et al., 2006

Linear models in GRASS GIS

Let us now look at the relation ship between mean temperature, elevation and latitude.

        
r.mapcalc "latitude = y()"
r.regression.line mapx=elev_ecad mapy=tmean_0_25deg.1981_2010.avg
r.mapcalc "tmean_model_line = 6.143929 + 0.002434 * elev_ecad"
r.mapcalc "tmean_model_diff_line = tmean_model_line - tmean_0_25deg.1981_2010.avg"
r.colors tmean_model_diff_line color=differences
r.regression.multi mapx=elev_ecad,latitude mapy=tmean_0_25deg.1981_2010.avg
r.mapcalc "tmean_model_multi = 40.722804 + -0.004234 * elev_ecad + -0.629130 * latitude"
r.mapcalc "tmean_model_diff_multi = tmean_model_multi - tmean_0_25deg.1981_2010.avg"
r.colors tmean_model_diff_multi color=differences
        
        
# Create a latitude map
r.mapcalc "latitude = y()"
# linear regression model
r.regression.line mapx=elev_ecad mapy=tmean_0_25deg.1981_2010.avg
# Apply the model using the model variables
r.mapcalc "tmean_model_line = 6.143929 + 0.002434 * elev_ecad"
# Calculate a difference map
r.mapcalc "tmean_model_diff_line = tmean_model_line - tmean_0_25deg.1981_2010.avg"
# Apply differences color palette
r.colors tmean_model_diff_line color=differences
# Multiple linear regression model
r.regression.multi mapx=elev_ecad,latitude mapy=tmean_0_25deg.1981_2010.avg
# Apply the model using the model variables
r.mapcalc "tmean_model_multi = 40.722804 + -0.004234 * elev_ecad + -0.629130 * latitude"
# Calculate a difference map
r.mapcalc "tmean_model_diff_multi = tmean_model_multi - tmean_0_25deg.1981_2010.avg"
# Apply differences color palette
r.colors tmean_model_diff_multi color=differences
        
 

display regression output
Multiple linear regression model

Visualize the model outputs and difference maps. Please check the manual of r.regression.line and r.regression.multi.

To use R functionalities inside GRASS GIS making use of the datasets in the mapset refer here

Other (very) useful links

References


This tutorial is adapted from Geostat exercises prepared by Markus Neteler and related content can be found in his blog .

Last changed: 02 Nov 2017

GRASS GIS manual main index | Topics index | Keywords Index | Full index | Raster index | Vector index | Temporal index |

Creative Commons License
Hands-on to GIS and RS with GRASS GIS by Veronica Andreo, Sajid Pareeth and Paulo van Breugel is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License - Thanks to Vaclav Petras for the style.