[STATSGRASS] unfilled matrices with R and gstat in grass 6.2.0beta3

[Roger Bivand]

On Thu, 15 Feb 2007 roger at spinn.net wrote:

> 
> Roger, 
> 
> Thanks for your help.  Most of my applications use data with an external 
> trend, so I normally use universal kriging.  My belief from reading the help 
> for "krige" and experimenting with it a little is that in order to use 
> universal kriging the independent variables that appear in the formula must 
> also appear in the "locations" and "newdata" arguments to "krige". 
> 
> We typically use "east" and "north" as the coordinate names.  Much of the 
> work in "read" is there to make sure the SpatialPointsDataFrame object used 
> as the "locations" argument for "krige" has data with those names.  I use a 
> SpatialPixelsDataFrame object as the "newdata" argument to "krige" because I 
> can name the independent variables in a SpatialPixelsDataFrame.  I have not 
> found a way to name the independent variables in a SpatialGridDataFrame. 
> 
> If there is some way to use universal kriging without going through all 
> those manipulations it would simplify my work a lot. 
> 

My understanding is that using the degree= argument to gstat(), then 
predict() from the gstat object, will include the point coordinates. You 
can check against surf.ls in the spatial package:

Grb <- gstat(formula=rb_elevati ~ 1, data=rb, degree=2)
predict(Grb, rb)

library(spatial)
x <- coordinates(rb)[, 1]
y <- coordinates(rb)[, 2]
print(predict(surf.ls(2, x, y, z=rb$rb_elevati), x, y))

If I recall correctly, adding a model= argument to the gstat() object 
definition will get you UK with your chosen trend.

Hope this helps,

Roger

> 
> Roger Miller 
> 
> 
> Roger Bivand writes: 
> 
> > On Wed, 14 Feb 2007, Roger Miller wrote: 
> > 
> >> I realized just recently that my reply on this thread went to Edzer
> >> rather than back to the list. 
> >> 
> >> Thanks for your replies.  The scripts are attached.  Both the wrong and
> >> the right (but odd) results resulted from running the same script with
> >> different data sets and different Grass region settings. 
> >> 
> >> The scripts are run in sequence "startup" followed by "read" and finally
> >> "map".
> > 
> > OK. Replace "startup" by calling gmeta2grd(), see the example on the 
> > gmeta6 help page. That gets the SpatialGridDataFrame object, if you need a 
> > SpatialPixelsDataFrame object, coerce using  
> > 
> > as(obj, "SpatialPixelsDataFrame") 
> > 
> > I have no idea what you are doing in "read", all you need to do is call 
> > readVECT6(), the rest is I think superfluous (with the possible exception 
> > of overwriting the WKT-style CRS object with the GRASS-style). The CRS 
> > issue can be handled as: 
> > 
> > proj4string(rb) <- CRS(G$proj4) 
> > 
> > I would not have written "map" anyway. If you show your users how to use 
> > the command line directly, they'll be able to fish themselves from then 
> > on, and variogram fitting does need interaction.  
> > 
> > Please try to get this working _using the command line only_ until it 
> > works, then simplify it until it cannot be reasonably simplified further. 
> > Then you will have what you need to build a TUI/GUI (a tcltk one for Rcmdr 
> > came past earlier this week ...) if it is worth the trouble (IMHO it 
> > usually isn't). 
> > 
> > Roger 
> > 
> >> 
> >> I also attached the script "tile" that I used to construct the 16 maps
> >> that I could patch together to get the correct result. 
> >> 
> >> The workspace is rather large, but it can be downloaded from 
> >> 
> >> http://members.spinn.net/~roger 
> >> 
> >> If that doesn't work then let me know and I'll try something else. 
> >> 
> >> 
> >> Roger Miller 
> >> 
> >> 
> > 
> > -- 
> > 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 
> > 
> >  
> > 
> >  
> > 
> 

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