[R-sig-Geo] Running analysis on DEMs using GRASS (from R?)
Agustin Lobo
alobolistas at gmail.com
Wed Mar 25 12:27:57 CET 2009
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).
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
>
>
More information about the R-sig-Geo
mailing list