[STATSGRASS] writing grass format from R

[Roger Bivand]

On Wed, 11 Jun 2003, Edzer J. Pebesma wrote:

> Roger Bivand wrote:
> 
> >I would think that the output of gstat in R could be moved back to the
> >GRASS database using the R/GRASS interface, and then accessed from
> >Mapserver. I don't have access to a computer with the necessary software
> >now, but will post a suggestion on Thursday morning - the R/GRASS
> >interface expects a vector ordered in the way that GRASS rasters are
> >stored.
> >
> >  
> >
> Gstat outputs a vector in the order of the prediction locations input 
> (newdata).  So, if
> this comes from a GRASS import, everything should work fine; however I'd 
> have to
> check what happens if part of this vector consists of NA's (missing values).
> 

It turns out to be fairly simple:

> library(GRASS)
> library(gstat)
> data(utm.maas)
> example(gmeta)
> names(utm.maas)
 [1] "east"    "north"   "x"       "y"       "elev"    "d.river" "Cd"     
 [8] "Cu"      "Pb"      "Zn"      "LOI"     "Fldf"    "Soil"   
> bmd142 <- vgm(158232, "Exp", 361.6, 24803.4)
> print(bmd142)
  model    psill range
1   Nug  24803.4   0.0
2   Exp 158232.0 361.6
> res <- krige(Zn~1, ~east+north, utm.maas, 
newdata=data.frame(east=east(G), north=north(G)), model=bmd142)
[using ordinary kriging]
> str(res)
> plot(G, res$var1.pred*maasmask)
> rast.put(G, "gstatOK", res$var1.pred*maasmask)
> retrieve <- rast.get(G, "gstatOK")
> str(retrieve)
List of 1
 $ gstatOK: num [1:60088] NA NA NA NA NA NA NA NA NA NA ...
> plot(G, retrieve$gstatOK)

the key being using newdata=data.frame() with the same coordinate location 
names as in the coordinate data frame. vgm() is here initialised with the 
values from Burrough & McDonnell p. 142. You can even try this outside 
GRASS, because the examples (if run outside GRASS) create a temporary 
LOCATION and MAPSET to try things out on. This will most likely replace 
the sgeostat dependency in the interface package, especially when I get 
round to writing vignettes.

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