[GRASS-stats] universal kriging with spgrass6
mathieu grelier
greliermathieu at gmail.com
Mon Jun 29 06:03:46 EDT 2009
ok,
>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 :
>sitesR <- readVECT6("grass_sites")
>class(sitesR)
> [1] "SpatialPointsDataFrame"
> d = attr(sitesR,"data")
>d
site value cat
1 Ail 81.80 1
2 All 38.50 2
3 Arg 103.50 3
4 Avb 52.50 4
...
> sitesR$x <- coordinates(sitesR)[,1]
>sitesR$y <- coordinates(sitesR)[,2]
> d = attr(sitesR,"data")
> class(d)
[1] "data.frame"
> d
site value cat x y
1 Ail 81.80 1 825094.9 6796083
2 All 38.50 2 758780.4 6851508
3 Arg 103.50 3 818973.3 6796125
4 Avb 52.50 4 775136.0 6877271
...
> coordinates(d) = ~x+y
> attr(d, "proj4string") <-CRS(G$proj4)
>G <- gmeta6()
>grd <- gmeta2grd()
>slot(grd, "cellsize") <- c(1000, 1000)
>data <- data.frame(list(k=rep(1,(floor(G$cols)*floor(G$rows)))))
> mask_SG <- SpatialGridDataFrame(grd, data=data, CRS(G$proj4))
>predictors = "x+y"
>kriging_result = autoKrige(as.formula(paste(column,"~",predictors)), d,
mask_SG)
Result was :
Error in gstat.formula.predict(d$formula, newdata, na.action = na.action, :
NROW(locs) != NROW(X): this should not occur
Warning messages:
1: 'newdata' had 7900 rows but variable(s) found have 73 rows
2: 'newdata' had 7900 rows but variable(s) found have 73 rows
1)do you think this approach could be successful?
2)If so, can you help me with this error?
Mathieu
2009/6/26 Edzer Pebesma <edzer.pebesma at uni-muenster.de>
>
>
> Roger Bivand wrote:
> > On Thu, 25 Jun 2009, Edzer Pebesma wrote:
> >
> >> Mathieu, that could very well be the case. Note that x+y do not
> >> "generally" refer to x- and y-coordinates, but are the actual names of
> >> the coordinates (or something else). From your example I cannot find out
> >> whether the coordinates are named x and y. In addition, they should have
> >> identical names in sitesR and mask_SG.
> >
> > Yes, it is beginning to look like the names in mask_SG, isn't it? Do
> > you have a view on whether degree= should be mandatory, or does the
> > infrastructure recognise coordinate RHS variables and rescale them
> > automatically (I doubt the latter)? Why do users use trend variables
> > in this way - the range of values in the (X'X) matrix becomes enormous
> > with the usual metre metric for coordinates?
> Another question we could ask (ourselves) is why in R we allow users to
> give their coordinates the names they like. I remember that for that
> sole reason (need to rename) I added the function
>
> coordnames(x) = c("u", "v")
>
> to sp. In case they're called "lon" and "lat", they also don't change
> name after projection. But maybe I talk too much with people worrying
> about semantics, these days. ;-)
> --
> Edzer
> >
> > Roger
> >
> >>
> >> Best regards,
> >> --
> >> Edzer
> >>
> >> mathieu grelier wrote:
> >>
> >> 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"
> >>
> >>> 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
> >>> ...
> >>>
> >>>
> >>>
> ------------------------------------------------------------------------
> >>>
> >>>
> >>> _______________________________________________
> >>> grass-stats mailing list
> >>> grass-stats at lists.osgeo.org
> >>> http://lists.osgeo.org/mailman/listinfo/grass-stats
> >>>
> >>
> >>
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
> http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-stats/attachments/20090629/296cc3dc/attachment.html
More information about the grass-stats
mailing list