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

[Roger Bivand]

On Wed, 25 Jan 2006, ivan marchesini wrote:

> Dear Roger...
> we have tried to use your suggestions...
> we have tried using ours data...
> all seems to work fine, but at the end (writeRast6sp) we have a strange
> error regarding the dimension of the cells:
> 
> error in writeAsciGrid (.........)
> Asciigrid does not support grids with non square cells
> 
> but this is not clear because we have used Gmeta6() to collect the
> location data of grass and before we have set "g.region res=20".
> 
> now I have seen that "g.region -p" gives me a resolution of
> ns=20.0003 (or similar)
> eo=19.99998 (or similar)
> 
> do you think is possible that this is the problem???

Yes, Arc ASCII grids have to have exactly square cells. I guess you can 
also force it by using:

cellsize=c(20, 20)

in the GridTopology() call.

Roger

> now we are trying to use a perfect integer region (ns=20 and eo=20)
> I will let you know the results....
> 
> best wishes
> 
> 
> 
> Il giorno lun, 23/01/2006 alle 21.09 +0100, Roger Bivand ha scritto:
> > 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