[GRASSLIST:9997] Re: R, kriging and grass6

[Roger Bivand]

On Mon, 23 Jan 2006, Kirk R. Wythers wrote:

> Roger,
> 
> Thanks for putting this together. I do have a very basic question  
> about getting the vector points from a GRASS vector to R. Can you  
> roll this description back a couple of steps and start with a GRASS  
> vector? Is this example based on the maas data?

Yes, this uses the data in the sp package for R:

library(sp)
library(maptools)
data(meuse)
coordinates(meuse) = c("x", "y")
writePointsShape(meuse, "meuse")
data(meuse.grid)
coordinates(meuse.grid) = c("x", "y")
meuse.grid_sp <- SpatialPixelsDataFrame(meuse.grid, data=as(meuse.grid, 
  "data.frame"))
writeAsciiGrid(meuse.grid_sp, "soil", fn="soil.aai")
meuse.grid_sp$ffreq1 <- as.integer(meuse.grid_sp$ffreq)
writeAsciiGrid(meuse.grid_sp, "ffreq1", fn="ffreq1.aai")

then in GRASS in an arbitrary location:

r.in.gdal input=soil.aai output=soil location=meuse

and moving to location meuse (more or less):

g.mapset PERMANENT
g.proj -c proj4='+proj=stere +lat_0=52.15616055555555 
  +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 
  +ellps=bessel +units=m'
# watch line breaks
g.mapset rsb
# some fiddling to get the projection data right
# then some more to get ewres and nsres to 40
r.in.gdal -o input=ffreq1.aai output=ffreq1
v.in.ogr -o dsn=. layer=meuse output=meuse

then R and the kriging code.

Hope this helps,

Roger

> 
> 
>   Kirk
> On Jan 23, 2006, at 2:09 PM, Roger Bivand wrote:
> 
> > This is the dense version, feedback welcome!
> >
> > library(spgrass6) # load package running within GRASS 6
> > meuse <- getSites6sp("meuse") # get vector points as
> > SpatialPointsDataFrame
> > class(meuse)
> > G <- gmeta6() # get region
> > grd <- GridTopology(cellcentre.offset=c(G$west+(G$ewres/2),
> >   G$south+(G$nsres/2)), cellsize=c(G$ewres, G$nsres),
> >   cells.dim=c(G$cols, G$rows))
> > mask_SG <- SpatialGridDataFrame(grd, data=list(k=rep(1, G$cols*G 
> > $rows)),
> >   proj4string=CRS(G$proj4)) # create a SpatialGridDataFrame
> > class(mask_SG)
> > cvgm <- variogram(zinc~1, loc=~x+y, data=meuse, width=100,  
> > cutoff=1000)
> > efitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100,
> >   nugget=1))
> > OK_pred <- krige(id="OK_pred", formula = zinc ~ 1, locations = ~ x  
> > + y,
> >   data = as(meuse, "data.frame"), newdata=mask_SG, model=efitted)
> >   # make the kriging prediction
> > names(OK_pred)
> > writeRast6sp(OK_pred, "OK_pred", zcol = "OK_pred.pred", NODATA=-9999)
> >
> 
> 

-- 
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