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

[ivan marchesini]

Thank you very much Roger....
I will try this afternoon to use your hints,...
I'll keep you informed about the results!!!!
thank you very much!!!!!!
ivan



Il giorno lun, 23/01/2006 alle 21.36 +0100, Roger Bivand ha scritto:
> 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
> 
> 
-- 
Ivan Marchesini
Department of Civil and Environmental Engineering
University of Perugia
Via G. Duranti 93/a 
06125
Perugia (Italy)
e-mail: marchesini at unipg.it
        ivan.marchesini at gmail.com
tel: +39(0)755853760
fax: +39(0)755853756