I am trying to write a code in R that use gstat library in order to create an interpolation. I have already read the gstat manual and based on some examples on internet I had managed to write this code (this is only a part):
g <- gstat(id="tec", formula=TEC ~ 1, data=data) ##I create an object
v <- variogram(g) # plot the empirical variogram
plot(v)
mod<-vgm(sill=var(data$TEC),model="Sph",range=200,nugget=200) #create the variogram model
v.fit <- fit.variogram(v, model=mod,fit.method=1) #fit the empirical variogram
Theor_variogram=plot(variogram(g),v.fit,main="WLS Model") #plot the theoretical variogram
plot(Theor_variogram)
## Kriging interpolation
p <- predict.gstat(g, model=v.fit, newdata=predGrid)
My problem is that, when I run the last command (predict) instead of getting a result with ordinary kriging interpolation, I get one with inverse distance weighted (IDW). I read in the gstat manual that: "When no variograms are specified, inverse distance weighted interpolation is the default action. When variograms are specified the default prediction method is ordinary kriging."
But, as you can see in my code, I specify both the empirical and theoretical variogram. Do you know why I keep getting IDW instead of ordinary kriging? Can it be related with the type of coordinates that I have? If for example I have coordinates close to each other, or if the region of interest is too big? Any help would be really useful.
Thanks in advance Dimitris
You need to include the creation of the gstat object, not in de prediction phase:
g <- gstat(id="tec", formula=TEC ~ 1, data=data, model = v.fit)
However, I would recommend using the standard interface for gstat
using krige
. This combines the building of the gstat
object and the prediction into one functions. Very rarely do you need to build gstat
objects yourself. For example:
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y
m <- vgm(.59, "Sph", 874, .04)
# OK:
x <- krige(log(zinc)~1, meuse, meuse.grid, model = m)
You could also use the automap
package (which I am the author of) and let the variogram model be automatically fitted to the data. For example using the meuse
dataset:
library(automap)
kr = autoKrige(log(zinc)~1, meuse, meuse.grid)
This will automatically build a sample variogram, and fit a variogram model to that sample variogram.