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

[ivan marchesini]

Dear Roger, dear grass/R users...
we would like to execute a kriging process of points data only inside a
well defined area.
we have realised a raster map that we would use as MASK inside R.
how can we impose to the krige command to perform only into the points
contained in the MASk???
we have imported the raster map in R using readFLOAT6sp and the
resulting map (named "structure") presents some points having value=11
and a lot of points having value=NA.
then we have created a

mask_struttura<-SpatialGridDataFrame(grd,data=structure,,proj4string=CRS(G$proj4))

later we have used the gstat krige command to perform the prediction,
defining the parameter newdata=mask_structure
but we obbtain this error:
NROWS(locs) != NROWS(X)

using simply all the region to perform the prediction 
mask_SG <- SpatialGridDataFrame(grd, data=list(k=rep(1, G$cols*G
$rows)),proj4string=CRS(G$proj4))
as you suggested, all is right...  
where are we wronging???  :-(
thank you 

Ivan

PS: we have started a wiki for the grass/R kriging problem in the grass
wiki.. as requested from someone.. 

 








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