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

[Tomislav Hengl]

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
http://home.medewerker.uva.nl/t.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?)

Tomislav,

This would be part of an R script using grass:

"brachspread" <- function(itmax=1,a=-1)
{
#Do the next once per session
#attach("/media/mifat32/Rutils/utiles.rda")
#Execute grass command: import tiff file into grass
#Data:
#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")

etc

Tomislav Hengl wrote:
> Dear list,
>
> I am trying to test some DEM analysis functions in GRASS GIS (under MS Windows machines). This is
my
> 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
a
> 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
>
>