[GRASSLIST:9996] kriging and grass6 - SOLVED

[ivan marchesini]

Ok Roger...
we have tried with a regular mesh grid (sqare cells) and all is gone
well...
Thank you very much!!!!
I hope all these e-mails will help someone else with grass,R and
geostatistics

Ivan



Il giorno mer, 25/01/2006 alle 10.47 +0100, Roger Bivand ha scritto:
> 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
> > > > > > 
> > > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> 
-- 
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