Geographic Resources Analysis Support System (GRASS) is the pioneer member of the is a very powerful geographic information system (GIS) software for geospatial data analysis, image processing, spatiotemporal modeling, and visualization (https://grass.osgeo.org/). Similarly, R is a free statistical software, and has now over 15,000 extension packages developed independently. Combining their capabilities will greatly improve a user’s processing power.
In this module, the students will learn to
use GRASS GIS spatial data processing within R and RStudio.
translate GRASS GIS commands into their R equivalent.
learn that you can use R without detailed knowledge in R
programming, as long as you know how to use a specific extension
package. In this study, we will use rgrass extension
package.
Please follow the steps below to complete GRASS GIS task within R and RStudio IDE:
The initGRASS command from rgrass library
either opens (if existing) or creates a GRASS database. You need to
supply the GISDbase folder C:/session10, location
loc and mapset PERMANENT.
library("rgrass")
initGRASS("C:/Program Files/GRASS GIS 8.2",
gisDbase="C:/session10", location="loc",
mapset="PERMANENT", override=TRUE)If the GRASS database is newly created, you need to set the Coordinate Reference System (CRS) first. This tutorial uses UTM 51N.
execGRASS("g.proj", flags = "c", epsg = 32651)We will import raster and vector
layers on the newly created GRASS database, and then and change the
environmental settings using g.region
If the GRASS database is newly created, you need to import raster and vector datasets, and setup its working environment.
The next command will work only if the digital elevation model (DEM) dem_region8.tif is stored in C:/session10/. Otherwise, download the image in your local folder and place its appropriate folder location.
#GRASS: r.in.gdal input=D:/Jun_turnover_files/DEM/dem_region8.tif output=dem8 --o
execGRASS("r.in.gdal",
input="C:/session10/dem_region8.tif",
output="dem8",
flag="overwrite")
# check if raster dem8 is properly imported
execGRASS("g.list",
type="raster", pattern="dem8") # GRASS: g.list type=rasteras the AOI region settings of the GRASS database
# GRASS equivalent: g.region -p
execGRASS("g.region", flag="p")# GRASS equivalent: g.region raster=dem8 -p
execGRASS("g.region", raster="dem8", flag="p")The next command will work only if the vector layer samarisland.geojson is saved in C:/session10/ folder. It contains a boundary delineation of the 3 provinces of Samar.
# GRASS: v.in.ogr input=D:/gisclass/samarisland.geojson output=samarisland
execGRASS("v.in.ogr",
input="C:/session10/samarisland.geojson",
output="samarisland",
flag="overwrite")
# check if properly imported
# GRASS: g.list type=vector
execGRASS("g.list", type="vector", pattern="samarisland") # g.list type=vector# check the original column names
# GRASS: db.columns table=samarisland
execGRASS("db.columns", table="samarisland")
# delete unnecessary columns
# GRASS: v.db.dropcolumn map=samarisland columns=ID_0,ISO,ID_1,NAME_0,VARNAME_1,NL_NAME_1,HASC_1,CC_1,TYPE_1,ENGTYPE_1,VALIDFR_1,VALIDTO_1,REMARKS_1
execGRASS("v.db.dropcolumn",
map="samarisland",
columns="ID_0,ISO,ID_1,NAME_0,VARNAME_1,NL_NAME_1,HASC_1,CC_1,TYPE_1,ENGTYPE_1,VALIDFR_1,VALIDTO_1,REMARKS_1")
# check the remaining column names
execGRASS("db.columns", table="samarisland")Add a new column and assign a value of 1. It will be used in creating a raster layer with a value of 1 for all the 3 samar provinces. We will use this raster to mask or limit the the analysis within the terrestrial boundaries of the 3 provinces of Samar. Neighboring provinces and the sea will be excluded in the analysis as long as a mask is set.
# "add an column name imask"
# v.db.addcolumn map=samarisland columns="imask integer"
execGRASS("v.db.addcolumn",
map="samarisland",
columns="imask integer")
# assign a value of 1 to imask column
# GRASS: v.db.update map=samarisland column=imask value=1
execGRASS("v.db.update",
map="samarisland",
column="imask",
value="1")To enhance processing speed and use lesser file storage, we will
reduce our area of interest using g.region settings based from a smaller
extent of a raster. To do this, we have to create a new raster covering
only the 3 Samar provinces. Then, use that smaller raster compared to
the original dem8 raster, as the new basis for the g.region
settings.
Use v.to.rast (https://grass.osgeo.org/grass82/manuals/v.to.rast.html) to convert from vector to raster. The code below will create a raster with a value of 1 for the 3 provinces of Samar, the rest will have a NULL value.
# vector to raster convertion
#v.to.rast input=samarisland type=area output=samisland use=attr attribute_column=imask value=1 --o
execGRASS("v.to.rast",
input="samarisland",
type="area",
output="samisland",
use="attr",
attribute_column="imask",
flag="overwrite")
# check if raster layers samisland is created
execGRASS("g.list", type="raster", pattern="samisland")We will reduce the extent of our region area of interest to the
extent samisland raster, which is smaller than the extent
of dem8 used earlier. We do this by changing the g.region
area of interest to the non-zero extent of the samarisland
raster. The code below will reduce the extent of the area of interest in
the succeeding operations The GRASS command g.region -p
will display their difference before and after running the
g.region zoom=samisland code.
# GRASS: g.region -p
execGRASS("g.region", flag="p")# GRASS:g.region zoom=samislsamislandand
execGRASS("g.region", zoom="samisland")# GRASS: g.region -p
execGRASS("g.region", flag="p")We can further limit the processing area to the non-zero part of the
input mask, which can either be a raster or
vector layer. We will limit processing only within the
island of Samar using the r.mask module. Data processing
will be within the non-zero pixels of the input raster
samisland mask, while the sea and the neighboring provinces
are ignored.
#GRASS: r.mask raster=samisland
execGRASS("r.mask", raster="samisland",
flag="overwrite")Run execGRASS("r.mask", flag="r") to remove the mask
#Delineate a watershed basin boundary
#GRASS: r.watershed elevation=dem8 threshold=50000 stream=d.stream --o
execGRASS("r.watershed",
elevation="dem8",
threshold=50000,
stream="rstream",
flag="overwrite")If there is no error, the next command will show the
rstream raster
execGRASS("g.list", type="raster", pattern="rstream")The r.thin command will make stream width into 1 pixel
wide only in prepation for vector convertion
#r.thin --o input=d.stream output=stream_thin
execGRASS("r.thin",input="rstream",
output="rstream_thin",
flag="overwrite")#r.to.vect --o type=line input=stream_thin output=vstream
execGRASS("r.to.vect",
type="line",
input="rstream_thin",
output="vstream",
flag="overwrite")Any area with 8 degrees or less slope will be considered as flat.
#r.geomorphon elevation=dem8 forms=rgeomorph flat=5
execGRASS("r.geomorphon",
elevation="dem8",
forms="rgeomorph",
flat=8,
flag="overwrite")
#display landform categories
#r.category rgeomorph
execGRASS("r.category", map="rgeomorph")The r.watershed module above delineates a watershed basin boundary
and stream network from a digital elevation model (DEM). Use dem8 as the
input DEM, threshold or a minimum 300,000 pixels (resolution is
approcimately 30x30 meter) to be considered a watershed. Any
sub-watershed that has less than 300,000 pixels are ignored, leading to
NULL pixels along the coastal areas.
#r.watershed elevation=dem8 basin=rbasin threshold=300000 --o
execGRASS("r.watershed",
elevation="dem8",
threshold=300000,
basin="rbasin",
flag="overwrite")#r.recode input=rbasin output=rbasin4 rules=C:/session10/rc.txt
execGRASS("r.recode",
input="rbasin",
output="rbasin2",
rules="C:/session10/rc.txt",
flag="overwrite")The rc.txt contains:
2:2:1:1
4:4:2:2
6:6:3:3
8:8:4:4
12:16:5:5
18:26:6:6
28:32:7:7
34:34:8:8
36:36:9:9
38:38:10:10
40:40:11:11
42:42:12:12execGRASS("r.report",
map="rbasin2",
units="h",
output="C:/session10/statreport.txt")#r.to.vect type=area input=rbasin2 output=vbasin --o -sv
execGRASS("r.to.vect",
type="area",
input="rbasin2",
output="vbasin",
flag=c("overwrite", "s", "v"))
execGRASS("g.list",type="vector", pattern="vbasin")# add a hectare column
# execGRASS("v.db.dropcolumn",
# map="vbasin",
# columns="hectare")
execGRASS("v.db.addcolumn",
map="vbasin",
columns="hectare double precision")
execGRASS("db.columns", table="vbasin")
# compute area in hectare
#v.to.db --o map=vwatershed2 layer=1 type=boundary option=area columns=hectare unit=hectares
execGRASS("v.to.db",
map="vbasin",
type="boundary",
option="area",
columns="hectare",
units="hectares",
flag="overwrite")
# display column names of the attribute table
execGRASS("db.columns", table="vbasin")#execGRASS("db.describe", table="vbasin")
execGRASS("v.db.select",map="vbasin",
columns=c("cat","hectare"),
format="csv")#r.cross input=rbasin,rgeomorph output=crossbasmorp --o
execGRASS("r.cross",
input=c("rbasin2","rgeomorph"),
output="crossbasmorp",
flag="overwrite")execGRASS("g.list", type="raster",pattern="crossbasmorp")#r.report map=crossbasmorp units=h output=C:/session10/crosstab.txt
execGRASS("r.report", map="crossbasmorp", units="h",
output="C:/session10/crosstab.txt",
flag="overwrite")This tutorial may be your first use of R, though not yet that
extensive since we used only the two (2) rgrass package
functions initGRASS and execGRASS.
initGRASS simply opened an invisible GRASS GIS interface
within R. On the other hand, execGRASS run GRASS GIS
commands with almost exactly the same syntax as running in the GRASS
console or command prompt.
In summary, we were able to:
Though we are able to demonstrate that we can use R and RStudio as the integrated development environment (IDE) for GRASS, in the next sessions, we shall discover the power of R both in spatial and nonspatial analysis of geospatial data.