ok, <br>From your answers, I've tried to play with the SpatialPointsDataFrame object imported from GRASS to fit the gstat krige function requirements. Basic idea seemed to work with the "data" attribute of the imported object, adding the missing coordinates columns : <br>
<br>>sitesR <- readVECT6("grass_sites")<br>>class(sitesR)<br>
> [1] "SpatialPointsDataFrame"<br>> d = attr(sitesR,"data")<br>>d<br><br> site value cat<br>1 Ail 81.80 1<br>2 All 38.50 2<br>3 Arg 103.50 3<br>
4 Avb 52.50 4<br>...<br><br>> sitesR$x <- coordinates(sitesR)[,1]<br>>sitesR$y <- coordinates(sitesR)[,2]<br>> d = attr(sitesR,"data")<br>> class(d)<br>[1] "data.frame"<br>
> d<br> site value cat x y<br>1 Ail 81.80 1 825094.9 6796083<br>2 All 38.50 2 758780.4 6851508<br>3 Arg 103.50 3 818973.3 6796125<br>4 Avb 52.50 4 775136.0 6877271<br>
...<br><br>> coordinates(d) = ~x+y<br>> attr(d, "proj4string") <-CRS(G$proj4)<br><br>>G <- gmeta6()<br>>grd <- gmeta2grd()<br>>slot(grd, "cellsize") <- c(1000, 1000)<br>>data <- data.frame(list(k=rep(1,(floor(G$cols)*floor(G$rows)))))<br>
> mask_SG <- SpatialGridDataFrame(grd, data=data, CRS(G$proj4))<br><br><br>>predictors = "x+y"<br>
>kriging_result = autoKrige(as.formula(paste(column,"~",predictors)), d, mask_SG)<br><br>Result was : <br><br>Error in gstat.formula.predict(d$formula, newdata, na.action = na.action, : <br> NROW(locs) != NROW(X): this should not occur<br>
Warning messages:<br>1: 'newdata' had 7900 rows but variable(s) found have 73 rows <br>2: 'newdata' had 7900 rows but variable(s) found have 73 rows <br><br>1)do you think this approach could be successful?<br>
2)If so, can you help me with this error? <br><br>Mathieu<br><br><div class="gmail_quote">2009/6/26 Edzer Pebesma <span dir="ltr"><<a href="mailto:edzer.pebesma@uni-muenster.de" target="_blank">edzer.pebesma@uni-muenster.de</a>></span><br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div><br>
<br>
Roger Bivand wrote:<br>
> On Thu, 25 Jun 2009, Edzer Pebesma wrote:<br>
><br>
>> Mathieu, that could very well be the case. Note that x+y do not<br>
>> "generally" refer to x- and y-coordinates, but are the actual names of<br>
>> the coordinates (or something else). From your example I cannot find out<br>
>> whether the coordinates are named x and y. In addition, they should have<br>
>> identical names in sitesR and mask_SG.<br>
><br>
> Yes, it is beginning to look like the names in mask_SG, isn't it? Do<br>
> you have a view on whether degree= should be mandatory, or does the<br>
> infrastructure recognise coordinate RHS variables and rescale them<br>
> automatically (I doubt the latter)? Why do users use trend variables<br>
> in this way - the range of values in the (X'X) matrix becomes enormous<br>
> with the usual metre metric for coordinates?<br>
</div>Another question we could ask (ourselves) is why in R we allow users to<br>
give their coordinates the names they like. I remember that for that<br>
sole reason (need to rename) I added the function<br>
<br>
coordnames(x) = c("u", "v")<br>
<br>
to sp. In case they're called "lon" and "lat", they also don't change<br>
name after projection. But maybe I talk too much with people worrying<br>
about semantics, these days. ;-)<br>
--<br>
<font color="#888888">Edzer<br>
</font><div><div></div><div>><br>
> Roger<br>
><br>
>><br>
>> Best regards,<br>
>> --<br>
>> Edzer<br>
>><br>
>> mathieu grelier wrote:<br>
>><br>
>> Dear all,<br>
>> I am trying to achieve universal kriging with GRASS and R through<br>
>> spgrass6<br>
>> 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>
>>> Dear all,<br>
>>> I am trying to achieve universal kriging with GRASS and R through<br>
>>> spgrass6<br>
>>> 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>
>>><br>
>>>> data(meuse)<br>
>>>> coordinates(meuse) = ~x+y<br>
>>>> data(meuse.grid)<br>
>>>> gridded(meuse.grid) = ~x+y<br>
>>>> column = "zinc"<br>
>>>><br>
>>> ##Ordinary kriging :<br>
>>><br>
>>>> predictors = "1"<br>
>>>> autoKrige(as.formula(paste(column,"~", predictors)), meuse,<br>
>>>> meuse.grid)<br>
>>>><br>
>>> ##Universal kriging :<br>
>>><br>
>>>> predictors = "x+y"<br>
>>>> autoKrige(as.formula(paste(column,"~", predictors)), meuse,<br>
>>>> meuse.grid)<br>
>>>><br>
>>><br>
>>> But using spgrass6, it becomes (some steps have been skipped):<br>
>>><br>
>>>> sitesR <- readVECT6("grass_sites")<br>
>>>><br>
>>> ...<br>
>>><br>
>>>> mask_SG <- SpatialGridDataFrame(grd, data=data, CRS(G$proj4))<br>
>>>><br>
>>> ...<br>
>>><br>
>>> ##Ordinary kriging => WORKS:<br>
>>><br>
>>>> predictors = "1"<br>
>>>> kriging_result = autoKrige(as.formula(paste(column,"~",predictors)),<br>
>>>><br>
>>> sitesR, mask_SG)<br>
>>><br>
>>> ##Universal kriging => ERROR:<br>
>>><br>
>>>> predictors = "x+y"<br>
>>>> kriging_result = autoKrige(as.formula(paste(column,"~",predictors)),<br>
>>>><br>
>>> sitesR, mask_SG)<br>
>>><br>
>>>> Error in eval(expr, envir, enclos) : "x" object not found<br>
>>>><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<br>
>>> same<br>
>>> structure ?<br>
>>><br>
>>>> class(meuse)<br>
>>>><br>
>>> [1] "data.frame"<br>
>>><br>
>>>> attributes(meuse)<br>
>>>><br>
>>> $names<br>
>>> [1] "x" "y" "cadmium" "copper" "lead" "zinc" "elev"<br>
>>><br>
>>>> class(sitesR)<br>
>>>><br>
>>> [1] "SpatialPointsDataFrame"<br>
>>><br>
>>>> attributes(sitesR)<br>
>>>><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>
>>><br>
>>><br>
>>> _______________________________________________<br>
>>> grass-stats mailing list<br>
>>> <a href="mailto:grass-stats@lists.osgeo.org" target="_blank">grass-stats@lists.osgeo.org</a><br>
>>> <a href="http://lists.osgeo.org/mailman/listinfo/grass-stats" target="_blank">http://lists.osgeo.org/mailman/listinfo/grass-stats</a><br>
>>><br>
>><br>
>><br>
><br>
<br>
--<br>
</div></div><div><div></div><div>Edzer Pebesma<br>
Institute for Geoinformatics (ifgi), University of Münster<br>
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251<br>
8333081, Fax: +49 251 8339763 <a href="http://ifgi.uni-muenster.de/" target="_blank">http://ifgi.uni-muenster.de/</a><br>
<a href="http://www.springer.com/978-0-387-78170-9" target="_blank">http://www.springer.com/978-0-387-78170-9</a> <a href="mailto:e.pebesma@wwu.de" target="_blank">e.pebesma@wwu.de</a><br>
<br>
</div></div></blockquote></div><br>