[R-sig-Geo] Running analysis on DEMs using GRASS (from R?)

Tomislav Hengl T.Hengl at uva.nl
Wed Mar 25 11:24:10 CET 2009

Dear Augustin,

Thank you for your message and examples shown (flavored with a bits of Portuguese). I am aware of
the spgrass6 package. I also know that there are plenty of examples for the R GRASS connection e.g.
http://casoilresource.lawr.ucdavis.edu/drupal/node/438, but I would still like to see how an R
script would look like for the case study shown below
("http:/www.geomorphometry.org/data/DEM25m.zip"; extraction of channel network). 

I am still confused whether GRASS can be run from R (or is it only the other way around?)? how to
import the data into GRASS (e.g. ArcInfo ASCII)? how to set-up the parameters of the GRASS
environment (working directory, projection system etc)? how to control the processing from R (e.g.
using "system(...)"), and then read the results of analysis back to R? 

Please excuse my ignorance (and the fact that I am a MS Windows user). I am just a beginner with
GRASS, so it could be that things are much simpler than they seem (it could also be that they are
also much more complicated).

Many thanks,

Tom Hengl

-----Original Message-----
From: Agustin Lobo [mailto:alobolistas at gmail.com] 
Sent: Wednesday, March 25, 2009 10:46 AM
To: Tomislav Hengl
Cc: r-sig-geo at stat.math.ethz.ch; carlos.grohmann at gmail.com
Subject: Re: [R-sig-Geo] Running analysis on DEMs using GRASS (from R?)


This would be part of an R script using grass:

"brachspread" <- function(itmax=1,a=-1)
#Do the next once per session
#Execute grass command: import tiff file into grass
#manchas (importado de ****)
#noocup (importado de no_ocupable1.tif, 3==no ocupable)
#(senderos importado de senderos_raster10m.tif)
#trailbuf importado de senderos
#r.buffer -z --o in=senderos out=trailbuf distance=10 units=meters
#camino:1; contiguo:2; resto, NULL

require(spgrass6) #Package connecting R<->grass
G <- gmeta6() #Read grass environment
system("g.region rast=manchas")

system("r.mapcalc 'habqual = if(trailbuf==1,0,trailbuf)'")
system("r.mapcalc 'habqual = if(habqual==2,0.8,habqual)'")
system("r.mapcalc 'habqual = if(isnull(habqual),0.5,habqual)'")

#Execute grass command: set region to raster
system("g.region rast=manchas")
system("r.mapcalc 'manchasit = manchas'")

#inicio del bucle de itersciones tiempo
it <-0
system("r.mapcalc 'prob = 0.000'")
for(it in 1:itmax) {
system("r.clump --o in=manchasit out=manchasit.clmp")
#Execute grass command: create stats file
system("r.stats -cl in=manchasit.clmp out=delme.txt")
#Read stats file into R:
mistat <- read.table("delme.txt",header=F,comment.char="*")
system("rm delme.txt") #delete garbage
#Write reclass file for grass
write.table(mistat, file="delme.txt", row.names=F, col.names=F ,sep="=")
system("r.reclass --o in=manchasit.clmp out=manchas.sze rules=delme.txt")
system("rm delme.txt")


Tomislav Hengl wrote:
> Dear list,
> I am trying to test some DEM analysis functions in GRASS GIS (under MS Windows machines). This is
> first contact with GRASS scripting (I have successfully installed GRASS 6.4 for windows from
> http://grass.itc.it/grass63/binary/mswindows/native/).
> For example, with SAGA GIS, I can completely control the analysis from R and use SAGA as external
> application to run the processing, e.g. to extract the drainage networks from a DEM:
>> library(rgdal)
>> library(RSAGA)
> # obtain the data:
> download.file("http://www.geomorphometry.org/data/DEM25m.zip",
> destfile=paste(getwd(),"DEM25m.zip",sep="/"))
> fname <- zip.file.extract(file="DEM25m.asc", zipname="DEM25m.zip")
> file.copy(fname, "./DEM25m.asc", overwrite=TRUE)
> # load to SAGA format:
>> rsaga.esri.to.sgrd(in.grids="dem25m.asc", out.sgrds="dem25m.sgrd", in.path=getwd())
> # Fill the spurious sinks:
>> rsaga.geoprocessor(lib="ta_preprocessor", module=2, param=list(DEM="dem25m.sgrd",
> RESULT="dem25m_f.sgrd", MINSLOPE=0.05))
> # extract the channel network:
>> rsaga.geoprocessor(lib="ta_channels", module=0, param=list(ELEVATION="dem25m_f.sgrd",
> CHNLNTWRK="channel_ntwrk.sgrd", CHNLROUTE="channel_route.sgrd", SHAPES="channels.shp",
> INIT_GRID="dem25m_f.sgrd", DIV_CELLS=3, MINLEN=40))
> # read back to R:
>> channels <- readOGR("channels.shp", "channels")
>> spplot(channels["LENGTH"], col.regions=bpy.colors())
> How could I run the same processing using GRASS? Can you control GRASS from R at all? (I would be
> bit disappointed with GRASS if I am not able to run the analysis from R).
> Many thanks for your time and suggestions in advance!
> Tom Hengl
> http://home.medewerker.uva.nl/t.hengl/
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

More information about the R-sig-Geo mailing list