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

[roger at spinn.net]

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. 


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