I have a point dataset which I'd like to Krige, ideally using an open-source software package. If possible, I'd also like to choose the semi-variogram model during the process to improve the estimation.
Open-Source Kriging – Methods for Kriging in Open-Source GIS
interpolationkrigingopen-source-gisspatial statistics
Related Solutions
I think there's a difference in the variogram between the model fitted here with gstat and the one in the linked question.
In the linked question, the model is expressed in terms of the covariance as C:
C(h) = σ² = 4, if |h| = 0
C(h) = (σ² - a)exp(-3|h|/r) = (4 - 1)exp(-3|h|/6), if |h| > 0
In R,
C = function(h){3*exp(-3*h/6)}
And the variogram from that covariance is σ²-C(h):
VC = function(h){4-C(h)}
Compare that with the variogram from gstat:
VCgstat = vgm(psill=3, nugget=1, range=6, "Exp")
plot(variogramLine(VCgstat,maxdist=50), type="l")
lines(0:50, VC(0:50), col="red")
I suspect this difference is a difference in the definition of "range" for the model. I don't think there is a clear agreement on the meaning of "range" for exponential variograms.
A little experimenting with the range parameter of the gstat model shows it is exactly equal to the VC function when range=2:
VCgstat2 = vgm(psill=3, nugget=1, range=2, "Exp")
And if I krige using that I get close agreement with the linked post:
> krige(z~1, p, SpatialPoints(cbind(1,1)), vgm(psill=3, nugget=1, range=2, "Exp"))
[using ordinary kriging]
coordinates var1.pred var1.var
1 (1, 1) 6.887426 3.136307
I'm not sure where this change has come from - possibly there's a multiply by sigma-squared minus 1 that shouldn't be there - but you should always carefully read (and then check!) the meaning of all parameters.
Best Answer
Depending on which Kriging type you want to apply, there are different packages to choose from:
Ordinary Kriging
The most common version is implemented for example in:
Simple Kriging
Simple Kriging uses the average of the entire data set while Ordinary Kriging uses a local average. Therefore, Simple Kriging can be less accurate, but it generally produces "smoother" results. It's implemented in:
Universal Kriging
Universal Kriging allows for consideration of drift in data. Implementations are included in:
Other Kriging Types
GRASS v.krige also supports Block Kriging.
HPGL implements a big number of less known Kriging methods (check the manual for more information on those):
SAGA offers different versions of both Ordinary and Universal Kriging.
Gstat krige additionally supports Block and Point Kriging.