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

[Roger Bivand]

On Mon, 23 Jan 2006, Roger Bivand wrote:

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

This is a simpler and (I think) more correct version of the gstat code:

library(gstat)
cvgm <- variogram(zinc ~ 1, locations=meuse, width=100, cutoff=1000)
efitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100,
  nugget=1))
OK_pred <- krige(zinc ~ 1, locations=meuse, newdata=mask_SG,
  model=efitted) # make the kriging prediction
names(OK_pred)
writeRast6sp(OK_pred, "OK_pred", zcol="var1.pred", NODATA=-9999)

Because gstat now uses these classes directly, the arguments can be 
simplified.

> 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