Dear all,<br>I am trying to achieve universal kriging with GRASS and R through spgrass6 package.<br>Maybe you could help me to find the right R syntax.<br>As in the gstat doc, the polynomial I try to use for now is "x+y"<br>
<br>
In R, with automap package, this works :<br>>data(meuse)<br>
>coordinates(meuse) = ~x+y<br>
>data(meuse.grid)<br><div>
>gridded(meuse.grid) = ~x+y<br></div>
>column = "zinc"<br>##Ordinary kriging : <br>>predictors = "1"<br>
>autoKrige(as.formula(paste(column,"~", predictors)), meuse, meuse.grid)<br>##Universal kriging :<br>
>predictors = "x+y"<br>
>autoKrige(as.formula(paste(column,"~", predictors)), meuse, meuse.grid)<br><br>But using spgrass6, it becomes (some steps have been skipped):<br>>sitesR <- readVECT6("grass_sites")<br>...<br>
>mask_SG <- SpatialGridDataFrame(grd, data=data, CRS(G$proj4))<br>...<br><br>##Ordinary kriging => WORKS: <br>>predictors = "1"<br>>kriging_result = autoKrige(as.formula(paste(column,"~",predictors)), sitesR, mask_SG)<br>
<br>##Universal kriging => ERROR: <br>
>predictors = "x+y"<br>
>kriging_result = autoKrige(as.formula(paste(column,"~",predictors)), sitesR, mask_SG)<br>>Error in eval(expr, envir, enclos) : "x" object not found<br><br>Why can't we use the x an y columns here?<br>
Maybe it is because the meuse and sites R data.frames don't have the same structure ?<br>> class(meuse)<br>[1] "data.frame"<br>> attributes(meuse)<br>$names<br> [1] "x" "y" "cadmium" "copper" "lead" "zinc" "elev" <br>
>class(sitesR)<br>[1] "SpatialPointsDataFrame"<br>> attributes(sitesR)<br>$bbox<br>...<br>$proj4string<br>...<br>$coords<br>...<br>$data<br> site value cat x y<br>1Ai 81.80 1 825094.9 6796083<br>
2 Al 38.50 2 758780.4 6851508<br>3 Ar 103.50 3 818973.3 6796125<br>4 Av 52.50 4 775136.0 6877271<br>...<br><br><br><br>