[GRASS-stats] universal kriging with spgrass6
Roger Bivand
Roger.Bivand at nhh.no
Mon Jun 29 06:26:30 EDT 2009
On Mon, 29 Jun 2009, mathieu grelier wrote:
> 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")
**NO**, there are no attr() in an S4 class. Just do the obvious thing:
d <- as(sitesR, "data.frame")
**BUT** you shouldn't be doing this at all, the data= argument in gstat
can take a SpatialPointsDataFrame.
>> 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
(Check, it has lost its x and y columns again, hasn't it?)
str(d)
>> attr(d, "proj4string") <-CRS(G$proj4)
**NO**!!! These are S4 classes, no attr()!!!
>
>> G <- gmeta6()
>> grd <- gmeta2grd()
**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().
>> 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))
>
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?
>
>> 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
>
No traceback()? Why? Isn't it obvious that it has found the x and y values
in d (73 long)?
> 1)do you think this approach could be successful?
> 2)If so, can you help me with this error?
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:
g_tr1a <- gstat("tr1a", <whatever column was> ~ 1, d, degree=1)
pr_tr1a <- predict(g_tr1a, mask_SG)
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.
Please do study the gstat help pages and the examples, they will help, but
you are jumping to far too many conclusions.
Roger
>
> 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
>>
>>
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the grass-stats
mailing list