[R-sig-Geo] Variogram Error

[Edzer Pebesma]

Saman,

Ben Graeler managed to reproduce the same effect by:

library(gstat)
data(meuse)
coordinates(meuse) <- ~x+y

obj = gstat(NULL, "D1", zinc~1, meuse[1,], set = list(zero_dist = 3))
obj = gstat(obj, "D2", zinc~1, meuse[2,])
variogram(obj, cross = "ONLY", pseudo = T, boundaries=c(0,100))

Essentially, for ST variograms, for each time slice the mean was removed
before computing the (pseudo) cross variogram. This will usually have
little effect, unless you have very few, or as in this (and your) case,
only one observation per time slice.

We changed the code now such that the mean is not removed. Updates are
committed r-forge; the next version on CRAN will include this update.




On 04/14/2013 06:10 PM, Saman Monfared wrote:
> Hi,
> I am trying  to fit a st variogram .The values of variogram$gamma
> in all st lags are 0. I don't no why this occurs!
> 
> Data is attached.
> 
> my program is:
> rm(list=ls())
> library(gstat)
> library(spacetime)
> library(maptools)
> library(RColorBrewer)
> library(maps)
> library(lattice)
> data<-read.table("cancer.rate.txt",header=TRUE)
> cancer.loc<-read.table("bladder.loc.txt",header=TRUE)
> pts<-cbind(cancer.loc$x,cancer.loc$y)
> pts = SpatialPoints(pts)
> data$time = ISOdate(data$year, data$mounth, data$day)
> ind<-as.matrix(cbind(data$indexs,data$indext))
> xx = STSDF(pts,data$time,data,ind)
> vv<-variogram(Rate~1,xx,width=200000,cutoff=3000000,tlags=1:8)
> plot(vv,map=T)
> wireframe(gamma~spacelag+timelag,vv,col="red",drape=T,outer =T,
> scales=list(arrows=FALSE))
> separableModel<-vgmST("separable",space=vgm(.05,"Gau",10000, 0.1),
> time =vgm(.02,"Gau", 100, 0),sill=.2,nugget=0)
> v.f<-fit.StVariogram(vv,separableModel,method="L-BFGS-B")
> v.f
> wireframe(model~spacelag+timelag,variogramSurface(v.f,vv),drape=T,
> aspect = c(1,1),panel.aspect =1,scales=list(arrows=F))
> plot(vv,v.f)
> grd<-read.table("farsgrid.txt",header=T)
> grd = SpatialPixels(SpatialPoints(cbind(grd$x,grd$y)))
> plot(grd)
> n =3
> library(xts)
> tgrd = seq(max(index(xx)+31622400*2),max(index(xx)+3*31622400),length=2)
> pred.grd = STF(grd, tgrd)
> plot(pred.grd )
> cancer.ST = krigeST(Rate~1,xx,pred.grd ,separableModel)
> stplot(cancer.ST)
> 
> 
> 
> 
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de