R/kriging.R

In geostats: An Introduction to Statistics for Geoscientists

.snr <- function(lsnr){
    snr <- rep(0,3)
    snr[1] <- exp(lsnr[1]) # sill
    snr[2] <- snr[1]*exp(lsnr[2])/(1+exp(lsnr[2])) # nugget
    snr[3] <- exp(lsnr[3]) # range
    snr
}
.lsnr <- function(snr){
    lsnr <- rep(0,3)
    lsnr[1] <- log(snr[1])
    lsnr[2] <- log(snr[2]) - log(snr[1]-snr[2])
    lsnr[3] <- log(snr[3])
    lsnr
}

semivarmod <- function(h,lsnr,model='spherical'){
    snr <- .snr(lsnr)
    s <- snr[1]
    n <- snr[2]
    r <- snr[3]
    out <- rep(0,length(h))
    if (model=='spherical'){
        i <- (h<r) & (h>0)
        out[i] <- n + (s-n)*(1.5*h[i]/r-0.5*(h[i]/r)^3)
        out[h>=r] <- s
    } else if (model=='exponential'){
        i <- (h>0)
        out[i] <- s + (n-s)*exp(-h[i]/r)
    } else if (model=='gaussian'){
        i <- (h>0)
        out[i] <- s + (n-s) * exp(-(h[i]/r)^2)
    } else {
        stop('Invalid semivariogram model')
    }
    out
}

#' @title semivariogram
#' @description Plots the semivariance of spatial data against
#'     inter-sample distance, and fits a spherical equation to it.
#' @param x numerical vector
#' @param y numerical vector of the same length as \code{x}
#' @param z numerical vector of the same length as \code{x}
#' @param bw (optional) the bin width of the semivariance search
#'     algorithm
#' @param nb (optional) the maximum number of bins to evaluate
#' @param plot logical. If \code{FALSE}, suppresses the graphical
#'     output
#' @param fit logical. If \code{TRUE}, returns the sill, nugget and
#'     range.
#' @param model the parametric model to fit to the empirical
#'     semivariogram (only used if \code{fit=TRUE}).
#' @param ... optional arguments to be passed on to the generic
#'     \code{plot} function
#' @return returns a list with the estimated semivariances at
#'     different distances for the data, and (if \code{fit=TRUE}), a
#'     vector with the sill, nugget and range.
#' @examples
#' data(meuse,package='geostats')
#' semivariogram(x=meuse$x,y=meuse$y,z=log(meuse$cadmium))
#' @export
semivariogram <- function(x,y,z,bw=NULL,nb=13,plot=TRUE,fit=TRUE,
                          model=c('spherical','exponential','gaussian'),
                          ...){
    d <- as.matrix(stats::dist(cbind(x,y)))
    if (is.null(bw)){
        N <- length(z)
        bw <- stats::quantile(d,probs=min(1,2/N))
    }
    sv <- rep(NA,nb)
    for (i in 1:nb){
        ij <- which((d>=((i-1)*bw)) & d<(i*bw),arr.ind=TRUE)
        if (any(ij)) sv[i] <- mean((z[ij[,1]]-z[ij[,2]])^2)/2
    }
    h <- bw*((1:nb)-1/2) # lag
    if (plot){
        plot(h,sv,type='p',xlab='lag',ylab=expression(gamma(lag)),
             ylim=c(0,max(sv,na.rm=TRUE)),...)
    }
    out <- list()
    out$model = model[1]
    out$h <- h
    out$sv <- sv
    if (fit){
        misfit <- function(lsnr,h,sv,model){
            svp <- semivarmod(h,lsnr=lsnr,model=model)
            sum((sv-svp)^2,na.rm=TRUE)
        }
        lsnr <- stats::optim(par=c(log(max(sv,na.rm=TRUE)),-2,log(max(h))),
                             fn=misfit,h=h,sv=sv,model=model[1])$par
        if (plot){
            x <- seq(from=.Machine$double.xmin,to=max(h),length.out=20)
            y <- semivarmod(h=x,lsnr=lsnr,model=model[1])
            graphics::lines(x,y)
        }
        out$snr <- .snr(lsnr)
        names(out$snr) <- c('sill','nugget','range')
    }
    invisible(out)
}

#' @title kriging
#' @description Ordinary kriging interpolation of spatial
#'     data. Implements a simple version of ordinary kriging that uses
#'     all the data in a training set to predict the z-value of some
#'     test data, using a semivariogram model generated by the
#'     \code{\link{semivariogram}} function.
#' @param x numerical vector of training data
#' @param y numerical vector of the same length as \code{x}
#' @param z numerical vector of the same length as \code{x}
#' @param xi scalar or vector with the \code{x}-coordinates of the
#'     points at which the \code{z}-values are to be evaluated.
#' @param yi scalar or vector with the \code{y}-coordinates of the
#'     points at which the \code{z}-values are to be evaluated.
#' @param svm output of the \code{\link{semivariogram}} function, a
#'     3-element vector with the sill, nugget and range of the
#'     semivariogram fit.
#' @param grid logical. If \code{TRUE}, evaluates the kriging
#'     interpolator along a regular grid of values defined by
#'     \code{xi} and \code{yi}.
#' @param err logical. If \code{TRUE}, returns the variance of the
#'     kriging estimate.
#' @return either a vector (if \code{grid=FALSE}) or a matrix (if
#'     \code{grid=TRUE}) of kriging interpolations. In the latter
#'     case, values that are more than 10\% out of the data range are
#'     given \code{NA} values.
#' @examples
#' data(meuse,package='geostats')
#' x <- meuse$x
#' y <- meuse$y
#' z <- log(meuse$cadmium)
#' svm <- semivariogram(x=x,y=y,z=z)
#' kriging(x=x,y=y,z=z,xi=179850,yi=331650,svm=svm,grid=TRUE)
#' @export
kriging <- function(x,y,z,xi,yi,svm,grid=FALSE,err=FALSE){
    # training data
    lsnr <- .lsnr(svm$snr)
    d <- as.matrix(stats::dist(cbind(x,y)))
    m <- semivarmod(h=d,lsnr=lsnr,model=svm$model)
    N <- length(x)
    M <- matrix(m,N,N)
    W <- rbind(cbind(M,1),1)
    W[N+1,N+1] <- 0
    out <- NA*xi
    Y <- c(z,0)
    if (grid){
        xyi <- expand.grid(xi,yi)
        Xi <- xyi[,1]
        Yi <- xyi[,2]
    } else {
        Xi <- xi
        Yi <- yi
    }
    vzi <- NA*Xi
    if (grid){
        buffer <- 0.05
        dx <- buffer*(max(x)-min(x))
        dy <- buffer*(max(y)-min(y))
        xy <- cbind(c(x-dx,x-dx,x+dx,x+dx),
                    c(y-dy,y+dy,y-dy,y+dy))
        i <- chull(xy)
        good <- which(inside(cbind(Xi,Yi),xy[i,,drop=FALSE]))
    } else {
        good <- 1:length(vzi)
    }
    for (i in good){
        h <- sqrt((Xi[i]-x)^2+(Yi[i]-y)^2)
        B <- c(semivarmod(h=h,lsnr=lsnr,model=svm$model),1)
        L <- MASS::ginv(W) %*% B
        if (err)
            vzi[i] <- t(B) %*% L
        else
            vzi[i] <- Y %*% L
    }
    if (grid){
        Ni <- length(xi)
        out <- matrix(vzi,Ni,Ni)
    } else {
        out <- vzi
    }
    out
}