Kriging In Splus

Boris & Sandra Kozintsev, Benjamin Kedem, Department of Mathematics, University of Maryland, 1999.

This web page contains algorithms for Ordinary and Trans-Gaussian Kriging (see Cressie (1993), Statistics for Spatial Data, Wiley: New York, p.123, 137), and their Splus implementation.

See also: Generating Gaussian Random Fields.

Ordinary Kriging

Trans-Gaussian Kriging

An Splus Implementation

The ordinary kriging is a particular case of the trans-Gaussian kriging when the transformation function is a transaltion, phi(x)=x+c. Therefore we implement the second algorithm only.

The covariance function r is defined as function(x, y, theta1, theta2, tau). For example, for the Spherical covariance function

the corresponding Splus code is
   sphercov <- function (x, y, theta1, theta2, tau) {  
      d <- sqrt (x^2 + y^2)  
      ifelse (d < theta1,
             (1- 1.5*(d/theta1)+0.5*(d/theta1)^3) / tau,
Note that the covariance function must be able to accept vectors for x and y.

The transformation function, its inverse, and its two derivatives are defined as function(x, lambda). For example, for the Box-Cox model,

the corresponding Splus code is
  phiInv <- function (x, lambda)
     if (lambda==0) log(x) else (x^lambda-1)/lambda
  phi <- function(x, lambda)
     if (lambda==0) exp(x) else (x*lambda+1)^(1/lambda)
  phiPrime <- function (x, lambda)
     if (lambda==0) exp(x) else (x*lambda+1)^(1/lambda-1)
  phiDouble <- function (x, lambda)
     if (lambda==0) exp(x) else
        lambda * (1/lambda - 1) * (lambda * x + 1)^(1/lambda-2)
Note that selecting lambda=1 in this model means using the ordinary kriging.

To perform kriging, the observed data must stored in the form of an n x 3 matrix, where each row represents an (x, y, z) point.

Then, the following function performs the trans-Gaussian kriging. The x- and y-coordinates of the locations at which the prediction is requested are passed to the function in the parameters sx and sy respectively.

   kriging <- function(sx, sy, data, covfun, theta1, theta2, tau, lambda=1) {
     if (length(dim(data)) != 2 | dim(data)[2] != 3)
                   stop("\nData must be an Nx3 matrix!")
     lensx <- length (sx)
     if (length(sy) != lensx) stop("\nsx and sy must be same length!")
     n <- dim(data)[1]
     x <- c(data[,rep(1:2, n)])
     y <- c(t(data[,c(rep(1,n), rep(2,n))]))
     dif <- matrix(t(matrix(x - y, ncol=2*n)), ncol=2, byrow=T)
     C <- matrix (covfun(dif[,1], dif[,2], theta1, theta2, tau), nrow=n)
     Cinv <- solve(C)
     sumCinv <- sum(Cinv)
     Z <- c(phiInv(data[,3], lambda))
     mu <- sum (Cinv %*% Z)/sumCinv
     prediction <- numeric(lensx)
     leftend <- numeric(lensx)
     rightend <- numeric(lensx)
     for (i in 1:lensx) {
       vectorc <- covfun(data[,1]-sx[i], data[,2]-sy[i], theta1, theta2, tau)
       m <- (1 - sum (Cinv %*% vectorc)) / sumCinv
       vectorlambda <- c((vectorc + m) %*% Cinv)
       sigma2 <- covfun(0, 0, theta1, theta2, tau) -
                 sum(vectorlambda*vectorc) + m
       prediction[i] <- phi(sum(vectorlambda*Z), lambda) +
                         phiDouble(mu, lambda) * (sigma2/2 - m)
       leftend[i] <- prediction[i] - 1.96*phiPrime(mu, lambda)*sqrt (sigma2)
       rightend[i] <- prediction[i] + 1.96*phiPrime(mu, lambda)*sqrt (sigma2)
     list (prediction=prediction, leftend=leftend, rightend=rightend)
The function returns a list of point predictors and left and right ends of the 95% prediction intervals.

Suppose that you have a text file with x, y, z data. For example, run our on-line field generation script (another window will open), and "Save As" the file under the grey-level picture on the left into the directory from which you started Splus, under a name field.txt.

To perform the ordinary kriging prediction at location (4, 6) assuming a Spherical(4, 1, 1) covariance model, type
  data <- matrix(scan(file="field.txt", skip=1), byrow=T, ncol=3)  
  kriging(4, 6, data, sphercov, 4, 1, 1)

The trans-Gaussian kriging under the Box-Cox model requires positive data. We can shift the observations up by 10, and perform the trans-Gaussian kriging with lambda=0.5:
  data[,3] <- data[,3] + 10  
  kriging(4, 6, data, sphercov, 4, 1, 1, 0.5)  

Prediction on the Grid

The kriging function above can predict several points at a time by accepting vectors for sx and sy. Therefore, to predict on the grid, one only has to create the vectors with x- and y- coordinates of the points of the grid.

Start by enumerating the projections of the grid on the X and Y axis, for example
  x <- seq (from=0, to=5, by=0.5)  
  y <- seq (from=0, to=3, by=0.5)  

Now create a grid as the cartesian product of x and y by typing
  grid <- cbind(c(matrix (rep (y, length(x)), ncol=length(y), byrow=T)),
                rep (x, length(y)))

To predict on the grid and save the point predictors to a variable called surface, type
  surface <- kriging (grid[,1], grid[,2], data,
                                sphercov, 4, 1, 1)$prediction

To view the prediction surface, type
  persp (x, y, matrix(surface, ncol=length(y)))

To compare the prediction to the original (observed) surface, type
  persp (matrix(data[,3], ncol=4))

Saving Prediction Results to a Text File

To save the results of the kriging in the example above, type
  write (t(cbind(grid, surface)), file = "pred3.txt", ncolumns=3)
This creates a file pred3.txt in the directory from which you started Splus with the data in the (x, y, z) form.

If you want to output Z values only, type
  write (t(surface), file = "predZ.txt")