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

[Roger Bivand]

On Mon, 23 Jan 2006, ivan marchesini wrote:

> Dear Roger...
> thank you for your previous answer... 
> can you fast describe me how to create a SpatialGrid object starting
> from the data returned from Gmeta6()???
> only few hints...  :-)
> I only need this to perform the kriging process, because I have already
> generated the variogram and I only need the grid over which to perform
> the prediction.....
> thank you very much...
> ivan

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


> 
> 
> Il giorno lun, 23/01/2006 alle 13.14 +0100, Roger Bivand ha scritto:
> > On Mon, 23 Jan 2006, ivan marchesini wrote:
> > 
> > > Hi to all...
> > > we are triyng to use kriging using R and grass6,
> > > We have seen a lot of documentation about R and GRASS5 but not so many
> > > about grass6...
> > > the problem is that we are trying to do an interpolation using the
> > > functions gmeta6 and kriging...
> > > if we are not wrong, the R library GRASS is developed for grass5 and
> > > inside there is the function krige.G that use the location data acquired
> > > using the gmeta() function...
> > > the problem is that we have found the gmeta6() function but not the
> > > krige.G6 function and obvously the krige.G doesn't work using the
> > > location data of grass6 (obtained using gmeta6)  
> > > 
> > > how can we simply solve this problem....
> > 
> > krige.G was a nasty hack, and should best be forgotten. The direct route 
> > is to use the sp classes in R and implemented for GRASS in spgrass6 - see 
> > GRASS News for the GRASS side and R News for the R side, combined with a 
> > proper R geostats package. Of those available, gstat is tightly bound to 
> > sp classes, so as far as I know, the only bit that needs doing by hand is 
> > to create a GridTopology object, or a SpatialGrid object, from the data 
> > returned by gmeta6() to pass to the kriging prediction function.
> > 
> > If you like, we can iterate to a working example from passing the vector 
> > points to R from GRASS, doing the modelling and kriging predictions, and 
> > passing the prediction rasters back to GRASS, but I'd like input from 
> > users to make the description sound.
> > 
> > Best wishes,
> > 
> > Roger
> > 
> > > 
> > > thank you
> > > 
> > > Ivan
> > > 
> > > 
> > > 
> > > 
> > 
> 

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