ok <br>First I think everyone here has its own domain of competence, so please don't answer as if I was at school, I'm not sure all these classes are so obvious for all GRASS users. Please clearly display the R requirements to post in the list if they exist<br>
I just try to implement UK in a GRASS process for a larger mapping system architecture. I think in some way it could be feasible : it should not be for me to understand how.<br><br>Also always consider improving your doc and/or your design : it shouldn't be normal for a non-specialist to spend so much time to understand how to manipulate the objects to achieve UK with GRASS.<br>
<br>Obviously I don't have time to study the R internals and I've just tried to quickly reproduce the working R code with meuse dataset, which use a dataframe as input data.<br>By the way, the rest of the code is yours, you personally gave me at the time. <br>
I understood that this is a colunm name problem and now I'm no longer sure what I want to do is possible without modification in gstat.<br>But I thought maybe there was a workaround?<br><br>Now I see it makes no sense but the initial R example seems to work with a data.frame as input data and it induced me in error.<br>
Best regards.<br><br>Mathieu<br><br><br><div class="gmail_quote">2009/6/29 Roger Bivand <span dir="ltr"><<a href="mailto:Roger.Bivand@nhh.no">Roger.Bivand@nhh.no</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 class="im">On Mon, 29 Jun 2009, mathieu grelier wrote:<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
ok,<br>
>From your answers, I've tried to play with the SpatialPointsDataFrame object<br>
imported from GRASS to fit the gstat krige function requirements. Basic idea<br>
seemed to work with the "data" attribute of the imported object, adding the<br>
missing coordinates columns :<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
sitesR <- readVECT6("grass_sites")<br>
class(sitesR)<br>
[1] "SpatialPointsDataFrame"<br>
d = attr(sitesR,"data")<br>
d<br>
</blockquote>
<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>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
sitesR$x <- coordinates(sitesR)[,1]<br>
sitesR$y <- coordinates(sitesR)[,2]<br>
d = attr(sitesR,"data")<br>
</blockquote></blockquote>
<br></div>
**NO**, there are no attr() in an S4 class. Just do the obvious thing:<br>
<br>
d <- as(sitesR, "data.frame")<br>
<br>
**BUT** you shouldn't be doing this at all, the data= argument in gstat can take a SpatialPointsDataFrame.<div class="im"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
class(d)<br>
</blockquote>
[1] "data.frame"<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
d<br>
</blockquote>
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>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
coordinates(d) = ~x+y<br>
</blockquote></blockquote>
<br></div>
(Check, it has lost its x and y columns again, hasn't it?)<br>
str(d)<div class="im"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
attr(d, "proj4string") <-CRS(G$proj4)<br>
</blockquote></blockquote>
<br></div>
**NO**!!! These are S4 classes, no attr()!!!<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
G <- gmeta6()<br>
grd <- gmeta2grd()<br>
</blockquote></blockquote>
<br>
**PLEASE** don't do this! You cannot set one slot without resetting the others if you want to cover the GRASS region. If you want something different, use g.region in GRASS and just gmeta2grd().<div class="im"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
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>
</blockquote>
<br>
</blockquote>
<br></div>
What is this supposed to be? What we've told you already is that mask_SG needs to be a Spatial*DataFrame with "x" and "y" columns, why are you not doing it?<div class="im"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
predictors = "x+y"<br>
kriging_result = autoKrige(as.formula(paste(column,"~",predictors)), d,<br>
</blockquote>
mask_SG)<br>
<br>
Result was :<br>
<br>
Error in gstat.formula.predict(d$formula, newdata, na.action = na.action, :<br>
<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>
</blockquote>
<br></div>
No traceback()? Why? Isn't it obvious that it has found the x and y values in d (73 long)?<div class="im"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
1)do you think this approach could be successful?<br>
2)If so, can you help me with this error?<br>
</blockquote>
<br></div>
Try doing everything in very small steps manually. Use gstat not automap, to avoid the extra layer of uncertainty. Simply try to set up a gstat() object for trend surface and predict from it:<br>
<br>
g_tr1a <- gstat("tr1a", <whatever column was> ~ 1, d, degree=1)<br>
pr_tr1a <- predict(g_tr1a, mask_SG)<br>
<br>
Once that (untried) works, go further, not until. I still suspect that this is to do with the colnames of the coordinates, (x, y) in one case and (coords1, coords2) or some such in the other.<br>
<br>
Please do study the gstat help pages and the examples, they will help, but you are jumping to far too many conclusions.<br><font color="#888888">
<br>
Roger</font><div><div></div><div class="h5"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
Mathieu<br>
<br>
2009/6/26 Edzer Pebesma <<a href="mailto:edzer.pebesma@uni-muenster.de" target="_blank">edzer.pebesma@uni-muenster.de</a>><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
<br>
Roger Bivand wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
On Thu, 25 Jun 2009, Edzer Pebesma wrote:<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
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>
</blockquote>
<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>
</blockquote>
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>
Edzer<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
Roger<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<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>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
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>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
data(meuse)<br>
coordinates(meuse) = ~x+y<br>
data(meuse.grid)<br>
gridded(meuse.grid) = ~x+y<br>
column = "zinc"<br>
<br>
</blockquote>
##Ordinary kriging :<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
predictors = "1"<br>
autoKrige(as.formula(paste(column,"~", predictors)), meuse,<br>
meuse.grid)<br>
<br>
</blockquote>
##Universal kriging :<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
predictors = "x+y"<br>
autoKrige(as.formula(paste(column,"~", predictors)), meuse,<br>
meuse.grid)<br>
<br>
</blockquote>
<br>
But using spgrass6, it becomes (some steps have been skipped):<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
sitesR <- readVECT6("grass_sites")<br>
<br>
</blockquote>
...<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
mask_SG <- SpatialGridDataFrame(grd, data=data, CRS(G$proj4))<br>
<br>
</blockquote>
...<br>
<br>
##Ordinary kriging => WORKS:<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
predictors = "1"<br>
kriging_result = autoKrige(as.formula(paste(column,"~",predictors)),<br>
<br>
</blockquote>
sitesR, mask_SG)<br>
<br>
##Universal kriging => ERROR:<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
predictors = "x+y"<br>
kriging_result = autoKrige(as.formula(paste(column,"~",predictors)),<br>
<br>
</blockquote>
sitesR, mask_SG)<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Error in eval(expr, envir, enclos) : "x" object not found<br>
<br>
</blockquote>
<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>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
class(meuse)<br>
<br>
</blockquote>
[1] "data.frame"<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
attributes(meuse)<br>
<br>
</blockquote>
$names<br>
[1] "x" "y" "cadmium" "copper" "lead" "zinc" "elev"<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
class(sitesR)<br>
<br>
</blockquote>
[1] "SpatialPointsDataFrame"<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
attributes(sitesR)<br>
<br>
</blockquote>
$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>
</blockquote></blockquote></blockquote>
------------------------------------------------------------------------<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<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>
</blockquote>
<br>
<br>
</blockquote>
<br>
</blockquote>
<br>
--<br>
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>
<br>
</blockquote>
<br>
</blockquote>
<br></div></div><div><div></div><div class="h5">
-- <br>
Roger Bivand<br>
Economic Geography Section, Department of Economics, Norwegian School of<br>
Economics and Business Administration, Helleveien 30, N-5045 Bergen,<br>
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43<br>
e-mail: <a href="mailto:Roger.Bivand@nhh.no" target="_blank">Roger.Bivand@nhh.no</a><br>
</div></div></blockquote></div><br>