[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