Skip to content

Variogram Error

2 messages · Saman Monfared, Edzer Pebesma

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





-- 
Saman Monfared
Msc, Department of Statistics, Shiraz University,
Shiraz 71454, Iran
Email: Samanmonfared1 at gmail.com

Tel: +98 917 5305167
#
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: