[GRASS-stats] universal kriging with spgrass6

mathieu grelier greliermathieu at gmail.com
Thu Jun 25 09:17:34 EDT 2009


Dear all,
I am trying to achieve universal kriging with GRASS and R through spgrass6
package.
Maybe you could help me to find the right R syntax.
As in the gstat doc, the polynomial I try to use for now is "x+y"

In R, with automap package, this works :
>data(meuse)
>coordinates(meuse) = ~x+y
>data(meuse.grid)
>gridded(meuse.grid) = ~x+y
>column = "zinc"
##Ordinary kriging :
>predictors = "1"
>autoKrige(as.formula(paste(column,"~", predictors)), meuse, meuse.grid)
##Universal kriging :
>predictors = "x+y"
>autoKrige(as.formula(paste(column,"~", predictors)), meuse, meuse.grid)

But using spgrass6, it becomes (some steps have been skipped):
>sitesR <- readVECT6("grass_sites")
...
>mask_SG <- SpatialGridDataFrame(grd, data=data, CRS(G$proj4))
...

##Ordinary kriging => WORKS:
>predictors = "1"
>kriging_result = autoKrige(as.formula(paste(column,"~",predictors)),
sitesR, mask_SG)

##Universal kriging => ERROR:
>predictors = "x+y"
>kriging_result = autoKrige(as.formula(paste(column,"~",predictors)),
sitesR, mask_SG)
>Error in eval(expr, envir, enclos) : "x" object not found

Why can't we use the x an y columns here?
Maybe it is because the meuse and sites R data.frames don't have the same
structure ?
> class(meuse)
[1] "data.frame"
> attributes(meuse)
$names
 [1] "x"       "y"       "cadmium" "copper"  "lead"    "zinc"    "elev"
>class(sitesR)
[1] "SpatialPointsDataFrame"
> attributes(sitesR)
$bbox
...
$proj4string
...
$coords
...
$data
           site value cat        x       y
1Ai      81.80   1 825094.9 6796083
2     Al      38.50   2 758780.4 6851508
3     Ar     103.50   3 818973.3 6796125
4           Av      52.50   4 775136.0 6877271
...
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-stats/attachments/20090625/19d14c3f/attachment.html


More information about the grass-stats mailing list