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

[Roger Miller]

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

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

-------------- next part --------------
#map
items<-names(rb)
cat("\nAvailable variables:\n")
for(n in 1:length(items)) {
  cat(paste(items[n],"\n"))
}
cat("\nEnter the formula: ")
f<-readLines(n=1)
if(length(f)>0) form<-formula(f)
rm(items,f)
#
var<-variogram(form,locations=rb)
rng0<-(sum(var$dist*var$np)/sum(var$np))
ngt0<-min(var$gamma)
sll0<-mean(var$gamma)
#
mod<-character()
while ( length(mod)==0 ) {
  cat("\nEnter the model mnemonic or \"list\" to display the options: ")
  mod<-readLines(n=1)
#
  if ( length(mod)==0 ) break
  if ( mod=="list" ) {
    items<-as.character(vgm()$long)
    for (n in 1:length(items)) {
      item<-items[n]
      cat(paste(item,"\n"))
    }    
    rm(n,item,items)
    mod<-character()
  }
}
#
if ( mod=="" ) {
  cat("\nWorking...")
  map<-krige(form,rb,g.grid)
} else {
  if ( mod=="Nug" ) {
    mdl<-vgm(psill=sll0,model=mod,range=0)
  } else {
    cat("Enter \"1\" to include a nugget effect:  ")
    if(readLines(n=1)=="1") {
      mdl<-vgm(psill=sll0,model=mod,nugget=ngt0,range=rng0)
    } else {
      mdl<-vgm(psill=sll0,model=mod,range=rng0)
    }
  }    
  cat("\nWorking...")
  varfit<-fit.variogram(var,model=mdl)
  map<-krige(form,rb,g.grid,model=varfit)
  rm(mdl)
}
rm (n,sll0,rng0,ngt0)
#rm (rb,var,varfit)
-------------- next part --------------
#read
cat("Enter the name of the Grass vector file with data: ")
fname<-readLines(n=1)
rb<-readVECT6(fname)
rbdf<-data.frame(rb)
if(charmatch("east",names(rbdf),nomatch=0)==0) {
  n<-charmatch("coords.x1",names(rbdf),nomatch=0)
  if(n>0) names(rbdf)[n]<-"east"
}
if(charmatch("north",names(rbdf),nomatch=0)==0) {
  n<-charmatch("coords.x2",names(rbdf),nomatch=0)
  if(n>0) names(rbdf)[n]<-"north"
}
rb<-SpatialPointsDataFrame(coordinates(rb),rbdf,proj4str=CRS(G$proj4))
rm(fname,rbdf)
-------------- next part --------------
#startup
library(spgrass6)
library(gstat)
G<-gmeta6()
x<-numeric()
x[1:(G$cols*G$rows)]<-(seq(G$w+G$ewres/2,G$e-G$ewres/2,G$ewres))
y<-numeric()
y[1:(G$cols*G$rows)]<-(seq(G$s+G$nsres/2,G$n-G$nsres/2,G$nsres))
#g.grid<-SpatialPoints(data.frame(east=x,north=y),
#                      proj4string=CRS(G$proj4)
#                      )
g.grid<-SpatialPixelsDataFrame(data.frame(east=x,north=y),
                               data.frame(rep(1,G$cols*G$rows)),
                               proj4string=CRS(G$proj4)
                               )
rm(x,y)
-------------- next part --------------
#tile
library(spgrass6)
library(gstat)
#
system("g.region rast=courson_consensus_rb_surfer")
G<-gmeta6()
source("read")
form<-formula("dbl_3~east+north")
mdl<-vgm(model="Exp",range=1000,psill=10000)
var<-variogram(form,rb)
varfit<-fit.variogram(var,mdl)
xv<-numeric()
yv<-numeric()
x0<-310295
y0<-3975671
cols<-251
rows<-171
inc<-30
#
y<-y0
for( i in 1:4 ) { 
  x<-x0
  y1<-(y+rows*inc)
  for (j in 1:4) {
    x1<-(x+cols*inc)
    system(sprintf("g.region w=%d s=%d e=%d n=%d",x,y,x1,y1))
    G<-gmeta6()
    xv[1:(G$cols*G$rows)]<-(seq(G$w+G$ewres/2,G$e-G$ewres/2,G$ewres))
    yv[1:(G$cols*G$rows)]<-(seq(G$s+G$nsres/2,G$n-G$nsres/2,G$nsres))
    g.grid<-SpatialPixelsDataFrame(data.frame(east=xv,north=yv),
                                   data.frame(rep(1,G$cols*G$rows)),
                                   proj4string=CRS(G$proj4)
                                   )
    map<-krige(form,rb,g.grid,model=varfit) 
    writeRAST6(map,sprintf("testmap%d%d",i,j),zcol="var1.pred")
    x<-x1
  }
  y<-y1
}
rm(x0,y0,cols,rows,inc,form,mdl,var,varfit,yv,xv,y,x,y1,x1,map)
system("g.region rast=courson_consensus_rb_surfer")