[STATSGRASS] raster cross-stats

[Roger Bivand]

Markus!

On Sun, 30 Dec 2001, Markus Neteler wrote:

> In fact I want to run
>  R CMD BATCH 
> 
> I have a small text file (Rbatchfile):
> library(GRASS)
> G <- gmeta()
> mapname <- "$1"
> map <- rast.get(G, mapname)
> map.frame <- data.frame(east(G), north(G), map$mapname)
> summary(map.frame)
> 
> Then I run the batch job:
> 
> R BATCH Rbatchfile aspect
> 
> cat aspect   # for some reason stored into this filename
> [...]
> Running in /usr/local/share/grassdata/spearfish/user1 
> > G <- gmeta()
> > mapname <- "$1"
> > map <- rast.get(G, mapname)
> Error in rast.get(G, mapname) : raster map: $1 not found
> Execution halted
> 
> When hard-coding the map name in the text file, everything runs
> well. But I don't get it flexible.
> 

The short answer is Sys.getenv().

Here I started R, did:

> data(cars)
> write.table(cars, "cars.txt")
> q()

to generate a test file. The batch file then refers to environment
variables which I assigned the file names to be used:

$ cat Rbatchfile.in
cars <- read.table(Sys.getenv("R_TABLE_IN"))
sink(Sys.getenv("THIS_R_OUT"))
cor(cars$dist, cars$speed)
summary(lm(dist ~ speed, cars))
sink()

$ echo $R_TABLE_IN
cars.txt
$ echo $THIS_R_OUT
R.out.1

So the batch file found them, and wrote the results to the output file.

$ R CMD BATCH Rbatchfile.in
$ cat R.out.1
[1] 0.8068949

Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.069  -9.525  -2.272   9.215  43.201 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.5791     6.7584  -2.601   0.0123 *  
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-Squared: 0.6511,     Adjusted R-squared: 0.6438 
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12 

Even works in a loop in a shell script, replacing the strings needed in R.

For your case that would be:

library(GRASS)
G <- gmeta()
mapname <- Sys.getenv("THIS_RASTER")
map <- rast.get(G, mapname)
map.frame <- data.frame(east(G), north(G), map$mapname)
sink(Sys.getenv("THIS_OUTFILE"))
summary(map.frame)
sink()

for shell variables $THIS_RASTER and $THIS_OUTFILE. I don't think it will
fly for numbered arguments to the R CMD BATCH script.

Very best wishes for the New Year to you and others on this list!

Roger

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no
and: Department of Geography and Regional Development, University of
Gdansk, al. Mar. J. Pilsudskiego 46, PL-81 378 Gdynia, Poland.