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

[ivan marchesini]

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