Nothing
.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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.