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

[Agustin Lobo]

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
>
>