[GRASS-stats] universal kriging with spgrass6

[Roger Bivand]

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?

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
>>
>
>

-- 
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