A simple example of Ordinary Kriging For Envinronmental Interpolation in R.

Feature: automatically finding the semivariogram using automap package

Data was taken from the lecture: Nicolas Christon. Ordinary kriging using geoR and gstat. Statistics C173/C273. Available online at :http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/c173c273_lec11_w11.pdf

Packages requirement are sp and automap. They can be downloaded by using install.packages command or Tools/Install Packages…

library(sp)
## Warning: package 'sp' was built under R version 3.2.5
library(automap)
## Warning: package 'automap' was built under R version 3.2.5

Known sampled data of X,Y coordinate with corresponding of Z value

X = c(61,63,64,68,71,73,75)
Y = c(139,140,129,128,140,141,128)
Z = c(477,696,227,646,606,791,783)

Unknown sampled data of X, Y coordinate

X1 = 65; Y1 = 137

Question is that what is the Z1 value corresponding to coordinate (X1,Y1). We create data frame for known sampled data

knowndt = data.frame(X,Y,Z)

Here’s create coordinates for known sampled data

coordinates(knowndt)  <- ~ X + Y

Here’s create data frame for unknown sampled data

unknowndt = data.frame(X1,Y1)

Here’s create coordinate for known sampled data

coordinates(unknowndt)  <- ~ X1 + Y1

Case 1: set the nugget, range,and sill as the same with those in the lecture

krigemod<- autoKrige(Z ~ 1, knowndt, unknowndt,
                     model = "Exp",  #use an exponential model. 
                     fix.values = c(0,3.33,10)) #C0,h,C1
## [using ordinary kriging]
predZ = krigemod$krige_output@data$var1.pred
predZ
## [1] 592.7587

The result of Z1(predZ) is the same as in the lecture note In spite of exponential model, Automap package supports other models such as Sph(shperical), Gau(gaussian), Mat(A model of the Matern familiy), and Ste Matern (M. Stein’s parameterization)

Case 2: dont set parameter, let autoKrige function automatically finds the suitable parameters

krigemod<- autoKrige(Z ~ 1, knowndt, unknowndt,
                     model = "Exp",  #use an exponential model. 
                     fix.values = c(NA,NA,NA)) #C0,h,C1
## [using ordinary kriging]
predZ = krigemod$krige_output@data$var1.pred
predZ
## [1] 603.7142

the result of Z1(predZ) in this case is 604, an difference of 1.8% compared with the result in the lecture note. It should be an acceptable result.