[STATSGRASS] speed of GRASS/R interface

[Agustin Lobo]

Carlos,

Just as an example, the attached R script which uses
the following grass shell script stats_aux.scr:

r.statistics base=$1 cover=$2 method=average output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.m1
r.statistics base=$1 cover=$2 method=variance output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.v1
r.statistics base=$1 cover=$2 method=sum output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.s1
#
r.statistics base=$1 cover=$3 method=average output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.m2
r.statistics base=$1 cover=$3 method=variance output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.v2
#
r.statistics base=$1 cover=$4 method=average output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.m3
r.statistics base=$1 cover=$4 method=variance output=delme --overwrite
cp /home/alobo/PERU/seg1/2001u/cats/delme delme.v3
#
g.remove rast=delme


The script calculates a table of stratified
statistics (one reference raster as strata and 3 rasters for the
data). I.e., You can calculate the average, sd etc of 3 bands for each
class on a given landcover map. There are probably other
(better) ways of doing this, probably using r.stats instead of 
r.statistics, but this could be a guide. The idea is
use r.stats and get that input into R instead of the
actual rasters. For example
r.stats -c in=miraster out=delme.txt

and then get delme.txt into R with read.table().For stratified stats:

r.stats -c in=reference,band out=delme.txt



Agus



Carlos "Guâno" Grohmann escribió:
> Thanks for the quick reply, Augustin.
> 
> I am doing a multi-scale study of surface roughness. So I have 180 
> rasters, created using 3 different methods, 5 different resolutions, and 
> with 12 moving-window sizes. In R, I am importing the data to produce 
> histograms, density plots and so on. More as a summary statistics tool 
> than really analysis. In R I can plot density curves for various maps in 
> one graphic, and see what is the behaviour of the morphometric variable 
> when resolution or moving-window size changes.
> 
> 
> Carlos
> 
> 
> 
> 
> On 10/2/07, *Agustin Lobo* <Agustin.Lobo at ija.csic.es 
> <mailto:Agustin.Lobo at ija.csic.es>> wrote:
> 
>     It makes sense, R puts everything in memory.
>     In general, unless you have small raster layers,
>     the best is to use the gis for the processing, get
>     intermediate data out of the raster and process these data in R.
>     Then back to gis for displaying if required.
>     In other words, no way of using R instead of r.mapcalc
>     unless the raster is very small.
> 
>     If you can let us know what you want to do with the raster
>     we can try to be more specific on how a combined use of grass/r could
>     solve the problem (hopefully).
> 
>     Agus
> 
>     Carlos "Guâno" Grohmann escribió:
>      > I've been thinking about the time it takes to import a raster
>     from GRASS
>      > to R. I have some rasters with about 2000 rows and 2500 columns.
>      > importing starts fine, creates the header, then goes to 100%, and it
>      > hangs there for about 10 minutes, before get ready for the next
>     command.
>      > Any way of decrease this waiting time?
>      >
>      > cheers
>      >
>      > Carlos
>      >
>      > --
>      > +-----------------------------------------------------------+
>      >               Carlos Henrique Grohmann - Guano
>      >   Visiting Researcher at Kingston University London - UK
>      >   Geologist M.Sc  - Doctorate Student at IGc-USP - Brazil
>      > Linux User #89721  - carlos dot grohmann at gmail dot com
>      > +-----------------------------------------------------------+
>      > _________________
>      > "Good morning, doctors. I have taken the liberty of removing
>     Windows 95
>      > from my hard drive."
>      > --The winning entry in a "What were HAL's first words" contest
>     judged by
>      > 2001: A SPACE ODYSSEY creator Arthur C. Clarke
>      >
>      > Can't stop the signal.
>      >
>      >
>      >
>     ------------------------------------------------------------------------
>      >
>      > _______________________________________________
>      > statsgrass mailing list
>      > statsgrass at grass.itc.it <mailto:statsgrass at grass.itc.it>
>      > http://grass.itc.it/mailman/listinfo/statsgrass
>      >
> 
>     --
>     Dr. Agustin Lobo
>     Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
>     LLuis Sole Sabaris s/n
>     08028 Barcelona
>     Spain
>     Tel. 34 934095410
>     Fax. 34 934110012
>     email: Agustin.Lobo at ija.csic.es <mailto:Agustin.Lobo at ija.csic.es>
>     http://www.ija.csic.es/gt/obster
> 
> 
> 
> 
> -- 
> +-----------------------------------------------------------+
>               Carlos Henrique Grohmann - Guano  
>   Visiting Researcher at Kingston University London - UK
>   Geologist M.Sc  - Doctorate Student at IGc-USP - Brazil
> Linux User #89721  - carlos dot grohmann at gmail dot com
> +-----------------------------------------------------------+
> _________________
> "Good morning, doctors. I have taken the liberty of removing Windows 95 
> from my hard drive."
> --The winning entry in a "What were HAL's first words" contest judged by 
> 2001: A SPACE ODYSSEY creator Arthur C. Clarke
> 
> Can't stop the signal.
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> statsgrass mailing list
> statsgrass at grass.itc.it
> http://grass.itc.it/mailman/listinfo/statsgrass
> 

-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster
-------------- next part --------------
"sgr.stats" <- function(nom.ref="t1pc1e3c", nompc1="t1pc1", nompc2="t1pc2",nompc3="t1pc3")
{
print("Executing sgr.stats...")
require(spgrass6)
G <- gmeta6()
comm <- paste("/home/alobo/PERU/stats_aux.scr", nom.ref, nompc1, nompc2, nompc3)
print(comm)
print("Grass...")
system(comm)
print("back to R...")
print(paste("    command executed in grass:",comm))

#Calculates the stat file for imorm
m1 <- read.table("delme.m1", skip=4,sep=":")
m1[m1=="*"] <- NA
m2 <- read.table("delme.m2", skip=4,sep=":")
m2[m2=="*"] <- NA
m3 <- read.table("delme.m3", skip=4,sep=":")
m3[m3=="*"] <- NA
v1 <- read.table("delme.v1", skip=4,sep=":")
v1[v1=="*"] <- NA
v2 <- read.table("delme.v2", skip=4,sep=":")
v2[v2=="*"] <- NA
v3 <- read.table("delme.v3", skip=4,sep=":")
v3[v3=="*"] <- NA
v1[,2] <- sqrt(v1[,2])
v2[,2] <- sqrt(v2[,2])
v3[,2] <- sqrt(v3[,2])
s1 <- read.table("delme.s1", skip=4,sep=":")
s1[s1=="*"] <- NA
n <- s1[,2]/m1[,2]
lin <- colu <- n*0
cv1 <- v1[,2]*100/m1[,2]
cv2 <- v2[,2]*100/m2[,2]
cv3 <- v3[,2]*100/m3[,2]
res <- cbind(m1[,1],lin, colu,n,
                              m1[,2],m2[,2],m3[,2],
                              v1[,2],v2[,2],v3[,2],
                              cv1,cv2,cv3)
dimnames(res)[[2]] <- c("etiq","lin", "col", "num",
                    "media 1", "media 2", "media 3",
                    "sd 1", "sd 2", "sd 3",
                    "cv 1", "cv 2", "cv 3")
system (paste("rm delme.m*", "delme.s*", "delme.v*"))
res
}