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

Thu Apr 2 18:02:23 CEST 2009

Thanks to Roger (and with many useful comments from Augustin) GRASS just got another big user  -
possible even many more :)

Roger has prepared an update of his package spgrass6 that can now be used to control the GRASS
processing from R (functions: "initGRASS", "execGRASS"). This way one can skip some manual settings
of GRASS (environmental parameters) and just continue with the analysis (this is now really easy).
These are the installation steps (MS Windows machines):

1. Download and install GRASS 6 (e.g. from http://grass.itc.it/grass63/binary/mswindows/native/)
2. Install the updated version of the spgrass6 package (from
http://spatial.nhh.no/R/Devel/spgrass6_0.6-1.zip)
3. Start R and then run the following code:

# obtain the data:
> download.file("http://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)

> library(spgrass6)  
# Location of your GRASS installation:
> loc <- initGRASS("C:/GRASS", home=tempdir())
> loc
# Import the ArcInfo ASCII file to GRASS:
> execGRASS("r.in.gdal", flags="o", parameters=list(input="DEM25m.asc", output="DEM"))
> execGRASS("g.region", parameters=list(rast="DEM"))
> gmeta6()
# extract the drainage network:
> execGRASS("r.watershed", flags="overwrite", parameters=list(elevation="DEM", stream="stream",
threshold=as.integer(50)))
# thin the raster map so it can be converted to vectors:
> execGRASS("r.thin", parameters=list(input="stream", output="streamt"))
# convert to vectors:
> execGRASS("r.to.vect", parameters=list(input="streamt", output="streamt", feature="line"))
> streamt <- readVECT6("streamt")
> plot(streamt)

Here are some screen shots:
http://geomorphometry.org/R.asp 

Tom Hengl
http://spatial-analyst.net 

> -----Original Message-----
> From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
> Sent: Wednesday, March 25, 2009 1:27 PM
> To: Agustin Lobo
> Cc: Tomislav Hengl; carlos.grohmann at gmail.com; r-sig-geo at stat.math.ethz.ch
> Subject: Re: [R-sig-Geo] Running analysis on DEMs using GRASS (from R?)
> 
> On Wed, 25 Mar 2009, Agustin Lobo wrote:
> 
> > Tomislav,
> >
> > You set up your grass location and mapset, import your data
> > etc. outside R. Then, from the grass console, you start R and run
> > an scrip like the one I sent you or tits individual commands. Once you've
> > done
> >
> > require(spgrass6) #Package connecting R<->grass
> > G <- gmeta6() #Read grass environment
> >
> > you can run grass commands using system(), pass results from
> > R to grass and from grass to R (It's done in the script
> > I sent you), import raster files (which I normally
> > don't do: raster files are too large for R in general, I prefer to
> > use grass for processing and then pass data as tables to R for
> > the stats; but have no experience with the new raster package made
> > by Roger
> 
> The raster package on R-forge is by Robert Hijmans and a team of
> developers - very interesting, and often the most active project on the
> whole (very active) site!
> 
> As you write, spgrass6 itself uses system() within its functions. A year
> ago I did some work on an automatic interface generator for GRASS commands
> in R using the interface descriptions provided by GRASS modules in XML.
> Would there be interest in using these to try to automate the writing of
> GRASS commands from R? Would anyone like to help? Should spgrass6 be moved
> from sourceforge CVS to R-forge SVN to make it easier for others to join
> in?
> 
> Roger
> 
> 
> >
> > Actually, I often use R as an scripting language for grass. I used
> > to do scripts in the Bourne shell itself (actually, it was csh), but using R
> > (perhaps python also)
> > is much easier (and powerful).
> >
> > I do not see that many differences between
> > the way you are using saga and the way you would use grass, except that
> > for grass you must have your location, mapset and region set. This
> > is actually very powerful, as (re)defining regions with g.region lets you
> > restrict
> > your analysis to an area of interest, which you can change by
> > repeated calls to g.region that could depend upon a dynamic process.
> > The same can be done by setting a MASK.
> >
> > The script that I sent you is part of a process simulating the spread of an
> > invasive
> > plant. Note that all control is within R. You can do the same with R alone,
> > but
> > if your raster layers are large, it could become very slow (note that there
> > is for() loop
> > in the script, and for() loops have a very low performance in R; as most of
> > the commands run within the loop are grass commands run through system(),
> > there is minimum memory built up).
> >
> > I've used the grass that comes within QGIS in windows sometimes, but I feel
> > really not comfortable
> > doing this type of work outside linux though.
> >
> > You seem experienced with saga. Actually, would like to talk to you on a
> > possible use
> > of sage from within qgis thanks to the fact that both saga and qgis make use
> > of python.
> >
> > Agus
> >
> > Tomislav Hengl wrote:
> >> 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
> >>>
> >>>
> >>>
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >>
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> 
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand at nhh.no