Search code examples
rkriginggstat

Variogram fit using gstat package in R


The following code is for predicting v value at three locations with the kriging method in R using gstat package.

library(gstat);
library(sp);
walk470 <- read.table("D:/kriging/walk470.txt",header=T)
attach(walk470)
coordinates(walk470) = ~x+y
walk.var1 <- variogram(v ~ x+y,data=walk470,width=10)
plot(walk.var1,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 5")
model1.out <- fit.variogram(walk.var1,vgm(70000,"Sph",40,20000))
plot(walk.var1, model=model1.out,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 10")
predpts <- matrix(c(60,190,225,50,110,185),ncol=2,byrow=T)
predpts.g <- data.frame(x=predpts[,1],y=predpts[,2])
coordinates(predpts.g) <- ~x+y
g <- gstat(NULL,"new.v",v~1,data=walk470,model=model1.out)
three.pred <- predict(g,predpts.g)
print(three.pred)

I want to know why for fitting variogram model, we need to provide the sill, nugget and range value beforehand using vgm() method. From kriging theory, I thought we have to calculate these values by minimizing the WLS objective function.

~Regards, Chandan


Solution

  • Since one of your parameters is the range of a spherical variogram model, the optimization problem solved by WLS is nonlinear. Nonlinear optimization typically requires initial values (see e.g. ?optim), and that is what you pass to the initial vgm call. The actual value passed is also somewhat important: if it is far outside a range of reasonable values (nearly zero or VERY large), the fit will not succeed.

    If you would only fit sill and nugget (passing fit.ranges=TRUE to fit.variogram), the problem would be linear and could in principle be done without initial value. In that case, arbitrary values far outside the data range will still work.