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

[Roger Miller]

Roger,

I learned several things.  The scripts are working fine now.  Thanks.

Roger Miller

On Thu, 2007-02-15 at 21:24 +0100, Roger Bivand wrote:
> 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 
> > > 
> > >  
> > > 
> > >  
> > > 
> > 
>