R/ksline.R
In rundel/geoR: Analysis of Geostatistical Data
"ksline" <-
function (geodata, coords=geodata$coords, data=geodata$data,
locations, borders = NULL,
cov.model = "matern",
cov.pars = stop("covariance parameters (sigmasq and phi) needed"),
kappa = 0.5, nugget = 0, micro.scale = 0,
lambda = 1, m0 = "ok", nwin = "full",
n.samples.backtransform = 500,
trend = 1, d = 2, ktedata = NULL, ktelocations = NULL,
aniso.pars = NULL, signal = FALSE, dist.epsilon = 1e-10,
messages)
{
if(missing(messages))
messages.screen <- as.logical(ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages")))
else messages.screen <- messages
##
## selecting locations inside the borders
##
locations <- .check.locations(locations)
## Checking for 1D prediction
if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
krige1d <- TRUE
else krige1d <- FALSE
if(!is.null(borders)){
locations <- locations.inside(locations, borders, bound=TRUE)
if(messages.screen)
cat("ksline: results will be returned only for prediction locations inside the borders\n")
}
call.fc <- match.call()
cov.model <- match.arg(cov.model, choices = .geoR.cov.models)
if(abs(lambda-1) > 0.001) {
if(messages.screen)
cat("ksline: Box-Cox data transformation performed.\n")
if(abs(lambda) < 0.001)
data <- log(data)
else data <- ((data^lambda) - 1)/lambda
}
coords <- as.matrix(coords)
locations <- as.matrix(locations)
dimnames(coords) <- list(NULL, NULL)
dimnames(locations) <- list(NULL, NULL)
if(!is.null(ktedata) & !is.null(ktelocations) & m0 != "kte"){
cat("ksline: external variable (covariate) provided. Kriging ste to KTE\n")
m0 <- "kte"
}
##
## anisotropy correction
##
if(!is.null(aniso.pars)) {
if(length(aniso.pars) != 2 | mode(aniso.pars) != "numeric")
stop("anisotropy parameters must be provided as a numeric vector with two elements: the rotation angle (in radians) and the anisotropy ratio (a number greater than 1)")
if(messages.screen)
cat("ksline: anisotropy correction performed\n")
coords.c <- coords.aniso(coords = coords, aniso.pars = aniso.pars)
locations.c <- coords.aniso(coords = locations, aniso.pars = aniso.pars)
}
else {
coords.c <- coords
locations.c <- locations
}
## 2. Preparing KTE matrices #####
##
if(m0 == "kte") {
ktedata <- as.matrix(ktedata)
ktelocations <- as.matrix(ktelocations)
dimnames(ktedata) <- list(NULL, NULL)
dimnames(ktelocations) <- list(NULL, NULL)
}
n <- length(data)
ni <- length(locations[, 1])
tausq <- nugget
sigmasq <- cov.pars[1]
phi <- cov.pars[2]
if(nwin == "full") {
est <- rep(0, ni)
dif <- rep(0, ni)
kvar <- rep(0, ni)
sumw <- rep(0, ni)
wofmean <- rep(0, ni)
iv <- varcov.spatial(coords = coords.c, cov.model = cov.model,
kappa = kappa, nugget = nugget, cov.pars = cov.pars,
inv = TRUE, det = FALSE, func.inv = "cholesky")$inverse
av <- mean(data)
sd <- sqrt(var(data))
one <- rep(1, n)
tone <- t(one)
toneiv <- crossprod(one, iv)
den <- solve(toneiv %*% one)
ml <- den %*% toneiv %*% data
kmsd <- sqrt(den)
means <- c(average = av, stdev = sd, kmean = ml, kmsd = kmsd)
if(m0 != "kt") {
mktlocations <- "Constant trend"
beta <- ml
}
else {
mktlocations <- rep(0, ni)
if(m0 == "kt" & trend == 1) {
if(d == 1) {
xmat <- cbind(rep(1, n), coords[, 2])
xmati <- cbind(rep(1, ni), locations[, 2])
}
else {
xmat <- cbind(rep(1, n), coords[, 1], coords[, 2])
xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2])
}
iviv <- solve(crossprod(xmat,iv) %*% xmat)
txiv <- crossprod(xmat,iv)
beta <- iviv %*% txiv %*% data
mkt <- xmat %*% beta
}
if(m0 == "kt" & trend == 2) {
if(d == 1) {
xmat <- cbind(rep(1, n), coords[, 2], (coords[, 2])^2)
xmati <- cbind(rep(1, ni), locations[, 2], (locations[,
2])^2)
}
else {
xmat <- cbind(rep(1, n), coords[, 1], coords[, 2],
(coords[, 1])^2, (coords[, 2])^2, coords[, 1] * coords[,
2])
xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2], (locations[, 1])^2, (locations[, 2])^2, locations[,
1] * locations[, 2])
}
iviv <- solve(crossprod(xmat,iv) %*% xmat)
txiv <- crossprod(xmat,iv)
beta <- iviv %*% txiv %*% data
mkt <- xmat %*% beta
}
}
if(m0 != "kte")
mktelocations <- "No external trend"
else {
if(m0 == "kte") {
mktelocations <- rep(0, ni)
xmat <- cbind(rep(1, n), ktedata)
xmati <- cbind(rep(1, ni), ktelocations)
iviv <- solve(crossprod(xmat,iv) %*% xmat)
txiv <- crossprod(xmat,iv)
beta <- iviv %*% txiv %*% data
mkte <- xmat %*% beta
}
}
for (i in 1:ni) {
if(messages.screen) {
if(ni < 11)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
else {
if(ni < 101 & (i%%10 == 1))
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(ni > 100 & i%%100 == 1)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(i == ni)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
}
}
coords0 <- cbind((coords.c[, 1] - locations.c[i, 1]), (coords.c[, 2] - locations.c[i,
2]))
dm0 <- sqrt(coords0[, 1]^2 + coords0[, 2]^2)
v0 <- cov.spatial(obj = dm0, cov.model = cov.model,
kappa = kappa, cov.pars = cov.pars)
v0[dm0 < dist.epsilon] <- micro.scale + sigmasq
tv0 <- t(v0)
v0iv <- crossprod(v0, iv)
v0ivv0 <- v0iv %*% v0
skw <- crossprod(v0,iv)
wofmean[i] <- 1 - sum(skw)
##
## 4.2.1 Simple kriging with known mean
##
if(mode(m0) == "numeric") {
dif[i] <- skw %*% (data - m0)
est[i] <- m0 + dif[i]
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0
else kvar[i] <- tausq + sigmasq - v0ivv0
sumw[i] <- sum(skw)
}
##
## 4.2.2 Simple kriging with data average mean
##
if(m0 == "av") {
dif[i] <- skw %*% (data - av)
est[i] <- av + dif[i]
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0
else kvar[i] <- tausq + sigmasq - v0ivv0
sumw[i] <- sum(((tone/n) + skw - ((skw %*% one %*%
tone)/n)))
}
##
## 4.2.3 Ordinary kriging (or SK with G.L.S. mean)
##
if(m0 == "ok") {
dif[i] <- skw %*% (data - ml)
est[i] <- ml + dif[i]
redu <- as.vector(1 - toneiv %*% v0)
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0 + (redu %*%
den %*% redu)
else kvar[i] <- tausq + sigmasq - v0ivv0 + (
redu %*% den %*% redu)
sumw[i] <- sum((den %*% one + tv0 - v0iv %*%
one %*% den %*% tone) %*% iv)
}
##
## 4.2.4 Universal Kriging (or Kriging with trend model)
##
if(m0 == "kt") {
dif[i] <- skw %*% (data - mkt)
est[i] <- xmati[i, ] %*% beta + dif[i]
redu <- as.vector(xmati[i, ]) - as.vector(
txiv %*% v0)
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0 + (redu %*%
iviv %*% redu)
else kvar[i] <- tausq + sigmasq - v0ivv0 + (
redu %*% iviv %*% redu)
sumw[i] <- sum(skw + xmati[i, ] %*% iviv %*%
txiv - skw %*% xmat %*% iviv %*% txiv)
mktlocations[i] <- xmati[i, ] %*% beta
}
##
## 4.2.5 Kriging with external trend
##
if(m0 == "kte") {
dif[i] <- skw %*% (data - mkte)
est[i] <- xmati[i, ] %*% beta + dif[i]
redu <- as.vector(xmati[i, ]) - as.vector(
txiv %*% v0)
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0 - (redu %*%
iviv %*% redu)
else kvar[i] <- tausq + sigmasq - v0ivv0 + (
redu %*% iviv %*% redu)
sumw[i] <- sum(skw + xmati[i, ] %*% iviv %*%
txiv - skw %*% xmat %*% iviv %*% txiv)
mktelocations[i] <- xmati[i, ] %*% beta
}
NULL
}
message <- "Kriging performed using global neighbourhood"
if(messages.screen)
cat(paste(message,"\n"))
results <- list(predict = est, krige.var = kvar, dif = dif, summary = means,
ktrend = mktlocations, ktetrend = mktelocations, beta = beta,
wofmean = wofmean)
}
else {
nwin <- min(n, nwin)
avwin <- rep(0, ni)
sdwin <- rep(0, ni)
mlwin <- rep(0, ni)
kmsdwin <- rep(0, ni)
estwin <- rep(0, ni)
difwin <- rep(0, ni)
kvarwin <- rep(0, ni)
sumwwin <- rep(0, ni)
wofmean <- rep(0, ni)
if(m0 != "kt")
mkt <- "Constant position trend"
else mkt <- rep(0, ni)
if(m0 != "kte")
mkte <- "No external trend"
else mkte <- rep(0, ni)
if(m0 != "kt" & m0 != "kte")
betawin <- "No polynomial or external trend"
if(m0 == "kt") {
if(trend == 1) {
if(d == 1)
xmati <- cbind(rep(1, ni), locations[, 2])
else xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2])
}
if(trend == 2) {
if(d == 1)
xmati <- cbind(rep(1, ni), locations[, 2], locations[,
2]^2)
else xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2], (locations[, 1])^2, (locations[, 2])^2, locations[, 1] *
locations[, 2])
}
betawin <- matrix(0, nrow = (ncol(xmati) * ni), ncol = ncol(xmati))
}
if(m0 == "kte") {
xmati <- cbind(rep(1, ni), ktelocations)
if(is.vector(ktedata) == TRUE)
betawin <- matrix(0, nrow = (2 * ni), ncol = 2)
else betawin <- matrix(0, nrow = ((ncol(ktedata) +
1) * ni), ncol = (ncol(ktedata) + 1))
}
for (i in 1:ni) {
temp.win <- .ksline.aux.1(coords = coords, coords.c = coords.c,
data = data, n = n,
locations = locations[i, ],
locations.c = locations.c[i, ],
cov.pars = cov.pars, nugget = nugget,
cov.model = cov.model, kappa = kappa, m0 = m0,
nwin = nwin, trend = trend, d = d, ktedata =
ktedata, ktelocations = ktelocations,
micro.scale = micro.scale,
location.number = i, xmati = xmati[i, ],
mkte = NULL, mkt = NULL, betawin = NULL,
signal = signal, dist.epsilon = dist.epsilon)
avwin[i] <- temp.win$avwin
sdwin[i] <- temp.win$sdwin
mlwin[i] <- temp.win$mlwin
kmsdwin[i] <- temp.win$kmsdwin
estwin[i] <- temp.win$estwin
difwin[i] <- temp.win$difwin
kvarwin[i] <- temp.win$kvarwin
sumwwin[i] <- temp.win$sumwwin
wofmean[i] <- temp.win$wofmean
if(m0 == "kt")
mkt[i] <- temp.win$mkt
if(m0 == "kte")
mkte[i] <- temp.win$mkte
if(m0 == "kt" | m0 == "kte")
betawin[i, ] <- temp.win$betawin
if(messages.screen) {
if(ni < 11)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
else {
if(ni < 101 & (i%%10 == 1))
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(ni > 100 & i%%100 == 1)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(i == ni)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
}
}
}
message <- "kriging performed in moving neighbourhood"
if(messages.screen)
cat(paste(message,"\n"))
results <- list(predict = estwin, krige.var = kvarwin, dif = difwin,
avtrend = avwin, sd = sdwin, oktrend = mlwin, oksd = kmsdwin,
ktrend = mkt, ktetrend = mkte, beta = betawin, sumw = sumwwin,
wofmean = wofmean)
}
if(abs(lambda - 1) > 0.001) {
if(messages.screen)
cat("Back-transforming the predicted mean and variance.\n")
if(abs(lambda) < 0.001) {
predict.transf <- results$predict
results$predict <- exp(predict.transf) - 0.5 * results$krige.var
results$krige.var <- (exp(2 * predict.transf - results$krige.var)) * (exp(results$krige.var) - 1)
}
if(lambda > 0.001) {
if(messages.screen)
cat("Back-transformation by simulating from the normal predictive distribution\n")
##ap.warn <- options()$warn
##options(warn = -1)
temp.data <- suppressWarnings(matrix(rnorm(ni * n.samples.backtransform,
mean = results$predict,
sd = sqrt(results$krige.var)),
nrow = ni))
##options(warn = ap.warn)
temp.data[(results$krige.var == 0), ] <- results$predict[(results$krige.var == 0)]
temp.data[temp.data < -1/lambda] <- -1/lambda
temp.data <- ((temp.data * lambda) + 1)^(1/lambda)
### temp.data[is.na(temp.data)] <- Inf
results$predict <- as.vector(rowMeans(temp.data))
results$krige.var <- as.vector(apply(temp.data, 1, var))
}
if(lambda < -0.001) {
cat("Resulting distribution has no mean for lambda < 0 - back transformation not performed\n"
)
}
}
results$locations <- locations
results$message <- message
results$call <- call.fc
attr(results, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
attr(results, "prediction.locations") <- call.fc$locations
if(!is.null(call.fc$borders))
attr(results, "borders") <- call.fc$borders
oldClass(results) <- c("kriging")
return(invisible(results))
}
".ksline.aux.1" <-
function (coords, coords.c, data, n, locations, locations.c, cov.pars,
nugget, cov.model, kappa,
m0, nwin, trend, d, ktedata, ktelocations, mbased,
micro.scale, location.number,
xmati, mkte, mkt, betawin, signal, dist.epsilon)
{
i <- location.number
sigmasq <- cov.pars[1]
phi <- cov.pars[2]
tausq <- nugget
coords0 <- cbind((coords.c[, 1] - locations.c[1]), (coords.c[, 2] -
locations.c[2]))
dm0 <- sqrt(coords0[, 1]^2 + coords0[, 2]^2)
coordswin <- coords[order(dm0)[1:nwin], ]
coordswin.c <- coords.c[order(dm0)[1:nwin], ]
datawin <- data[order(dm0)[1:nwin]]
ivwin <- varcov.spatial(coords = coordswin.c, cov.model = cov.model,
kappa = kappa, nugget = nugget, cov.pars = cov.pars, inv = TRUE,
det = FALSE, func.inv = "cholesky", only.decomposition = FALSE)$inverse
avwin <- mean(datawin)
sdwin <- sqrt(var(datawin))
onewin <- rep(1, nwin)
toneivwin <- crossprod(onewin, ivwin)
denwin <- solve(toneivwin %*% onewin)
mlwin <- denwin %*% toneivwin %*% datawin
kmsdwin <- sqrt(denwin)
coords0win <- cbind((coordswin[, 1] - locations[1]), (coordswin[, 2] -
locations[2]))
coords0win.c <- cbind((coordswin.c[, 1] - locations.c[1]), (coordswin.c[
, 2] - locations.c[2]))
dm0win <- sqrt(coords0win.c[, 1]^2 + coords0win.c[, 2]^2)
v0win <- cov.spatial(obj = dm0win, cov.model = cov.model, kappa = kappa,
cov.pars = cov.pars)
v0win[dm0win < dist.epsilon] <- micro.scale + sigmasq
skwwin <- crossprod(v0win, ivwin)
wofmean <- 1 - sum(skwwin)
if(m0 == "kt" & trend == 1) {
if(d == 1)
xmatwin <- cbind(rep(1, nwin), coordswin[, 2])
else xmatwin <- cbind(rep(1, nwin), coordswin[, 1], coordswin[
, 2])
txivwin <- crossprod(xmatwin, ivwin)
ivivwin <- solve(txivwin %*% xmatwin)
betawin <- ivivwin %*% txivwin %*% datawin
mktwin <- xmatwin %*% betawin
mkt <- xmati %*% betawin
}
if(m0 == "kt" & trend == 2) {
if(d == 1)
xmatwin <- cbind(rep(1, nwin), coordswin[, 2], (
coordswin[, 2])^2)
else xmatwin <- cbind(rep(1, nwin), coordswin[, 1], (coordswin[
, 1])^2, coordswin[, 2], (coordswin[, 2])^
2, coordswin[, 1] * coordswin[, 2])
xmatwin.cent <- xmatwin
xmatwin.cent[, 2] <- xmatwin.cent[, 2] - mean(xmatwin[, 2])
xmatwin.cent[, 3] <- xmatwin.cent[, 3] - mean(xmatwin[, 3])
ivivwin <- solve(crossprod(xmatwin.cent, ivwin) %*%
xmatwin.cent)
txivwin <- crossprod(xmatwin.cent, ivwin)
betawin <- ivivwin %*% txivwin %*% datawin
betawin <- mean(datawin) - crossprod(betawin, c(0, mean(xmatwin[
, 2]), mean(xmatwin[, 3])))
mktwin <- xmatwin %*% betawin
mkt <- xmati %*% betawin
}
if(m0 == "kte") {
if(is.vector(ktedata))
ktedatawin <- ktedata[order(dm0)[1:nwin]]
else ktedatawin <- ktedata[order(dm0)[1:nwin], ]
xmatwin <- cbind(rep(1, nwin), ktedatawin)
ivivwin <- solve(crossprod(xmatwin, ivwin) %*% xmatwin)
txivwin <- crossprod(xmatwin, ivwin)
betawin <- ivivwin %*% txivwin %*% datawin
mktewin <- xmatwin %*% betawin
mkte <- xmati %*% betawin
}
##
## Simple kriging with know mean
##
if(mode(m0) == "numeric") {
difwin <- skwwin %*% (data - m0)
# estwin <- m0win + difwin
estwin <- m0 + difwin
if(signal)
kvarwin <- sigmasq - crossprod(v0win, ivwin) %*% v0win
else kvarwin <- tausq + sigmasq - crossprod(v0win, ivwin) %*%
v0win
sumwwin <- sum(skwwin)
}
##
## 4.2.2 Simple kriging with data average mean
##
if(m0 == "av") {
difwin <- skwwin %*% (datawin - avwin)
estwin <- avwin + difwin
if(signal)
kvarwin <- sigmasq - crossprod(v0win, ivwin) %*% v0win
else kvarwin <- tausq + sigmasq - crossprod(v0win, ivwin) %*%
v0win
sumwwin <- sum(((t(onewin)/nwin) + skwwin - ((skwwin %*% onewin %*%
t(onewin))/n)))
}
##
## Ordinary kriging (or SK with G.L.S. mean)
##
if(m0 == "ok") {
difwin <- skwwin %*% (datawin - mlwin)
estwin <- mlwin + difwin
redu <- as.vector(1 - toneivwin %*% v0win)
if(signal)
kvarwin <- sigmasq - v0win %*% ivwin %*% v0win + (
redu %*% denwin %*% redu)
else kvarwin <- tausq + sigmasq - v0win %*% ivwin %*% v0win +
(redu %*% denwin %*% redu)
sumwwin <- sum((denwin %*% onewin + t(v0win) - crossprod(v0win,
ivwin) %*% onewin %*% denwin %*% t(onewin)) %*% ivwin)
}
##
## Universal Kriging (or Kriging with trend model)
##
if(m0 == "kt") {
difwin <- skwwin %*% (datawin - mktwin)
estwin <- mkt + difwin
xmati <- as.vector(xmati)
redu <- as.vector(xmati) - as.vector(txivwin %*% v0win)
if(signal)
kvarwin <- sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
else kvarwin <- tausq + sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
sumwwin <- sum(skwwin + xmati %*% ivivwin %*% txivwin - skwwin %*%
xmatwin %*% ivivwin %*% txivwin)
}
##
## Kriging with external trend
##
if(m0 == "kte") {
difwin <- skwwin %*% (datawin - mktewin)
estwin <- mkte + difwin
xmati <- as.vector(xmati)
redu <- as.vector(xmati) - as.vector(txivwin %*% v0win)
if(signal)
kvarwin <- sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
else kvarwin <- tausq + sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
sumwwin <- sum(skwwin + xmati %*% ivivwin %*% txivwin - skwwin %*%
xmatwin %*% ivivwin %*% txivwin)
}
##
##
##
results <- list(avwin = avwin, sdwin = sdwin, mlwin = mlwin, kmsdwin =
kmsdwin, wofmean = wofmean, betawin = betawin, mkt = mkt, mkte
= mkte, difwin = difwin, estwin = estwin, kvarwin = kvarwin,
sumwwin = sumwwin)
return(results)
}
rundel/geoR documentation built on May 18, 2019, 11:28 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
"ksline" <-
function (geodata, coords=geodata$coords, data=geodata$data,
locations, borders = NULL,
cov.model = "matern",
cov.pars = stop("covariance parameters (sigmasq and phi) needed"),
kappa = 0.5, nugget = 0, micro.scale = 0,
lambda = 1, m0 = "ok", nwin = "full",
n.samples.backtransform = 500,
trend = 1, d = 2, ktedata = NULL, ktelocations = NULL,
aniso.pars = NULL, signal = FALSE, dist.epsilon = 1e-10,
messages)
{
if(missing(messages))
messages.screen <- as.logical(ifelse(is.null(getOption("geoR.messages")), TRUE, getOption("geoR.messages")))
else messages.screen <- messages
##
## selecting locations inside the borders
##
locations <- .check.locations(locations)
## Checking for 1D prediction
if(length(unique(locations[,1])) == 1 | length(unique(locations[,2])) == 1)
krige1d <- TRUE
else krige1d <- FALSE
if(!is.null(borders)){
locations <- locations.inside(locations, borders, bound=TRUE)
if(messages.screen)
cat("ksline: results will be returned only for prediction locations inside the borders\n")
}
call.fc <- match.call()
cov.model <- match.arg(cov.model, choices = .geoR.cov.models)
if(abs(lambda-1) > 0.001) {
if(messages.screen)
cat("ksline: Box-Cox data transformation performed.\n")
if(abs(lambda) < 0.001)
data <- log(data)
else data <- ((data^lambda) - 1)/lambda
}
coords <- as.matrix(coords)
locations <- as.matrix(locations)
dimnames(coords) <- list(NULL, NULL)
dimnames(locations) <- list(NULL, NULL)
if(!is.null(ktedata) & !is.null(ktelocations) & m0 != "kte"){
cat("ksline: external variable (covariate) provided. Kriging ste to KTE\n")
m0 <- "kte"
}
##
## anisotropy correction
##
if(!is.null(aniso.pars)) {
if(length(aniso.pars) != 2 | mode(aniso.pars) != "numeric")
stop("anisotropy parameters must be provided as a numeric vector with two elements: the rotation angle (in radians) and the anisotropy ratio (a number greater than 1)")
if(messages.screen)
cat("ksline: anisotropy correction performed\n")
coords.c <- coords.aniso(coords = coords, aniso.pars = aniso.pars)
locations.c <- coords.aniso(coords = locations, aniso.pars = aniso.pars)
}
else {
coords.c <- coords
locations.c <- locations
}
## 2. Preparing KTE matrices #####
##
if(m0 == "kte") {
ktedata <- as.matrix(ktedata)
ktelocations <- as.matrix(ktelocations)
dimnames(ktedata) <- list(NULL, NULL)
dimnames(ktelocations) <- list(NULL, NULL)
}
n <- length(data)
ni <- length(locations[, 1])
tausq <- nugget
sigmasq <- cov.pars[1]
phi <- cov.pars[2]
if(nwin == "full") {
est <- rep(0, ni)
dif <- rep(0, ni)
kvar <- rep(0, ni)
sumw <- rep(0, ni)
wofmean <- rep(0, ni)
iv <- varcov.spatial(coords = coords.c, cov.model = cov.model,
kappa = kappa, nugget = nugget, cov.pars = cov.pars,
inv = TRUE, det = FALSE, func.inv = "cholesky")$inverse
av <- mean(data)
sd <- sqrt(var(data))
one <- rep(1, n)
tone <- t(one)
toneiv <- crossprod(one, iv)
den <- solve(toneiv %*% one)
ml <- den %*% toneiv %*% data
kmsd <- sqrt(den)
means <- c(average = av, stdev = sd, kmean = ml, kmsd = kmsd)
if(m0 != "kt") {
mktlocations <- "Constant trend"
beta <- ml
}
else {
mktlocations <- rep(0, ni)
if(m0 == "kt" & trend == 1) {
if(d == 1) {
xmat <- cbind(rep(1, n), coords[, 2])
xmati <- cbind(rep(1, ni), locations[, 2])
}
else {
xmat <- cbind(rep(1, n), coords[, 1], coords[, 2])
xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2])
}
iviv <- solve(crossprod(xmat,iv) %*% xmat)
txiv <- crossprod(xmat,iv)
beta <- iviv %*% txiv %*% data
mkt <- xmat %*% beta
}
if(m0 == "kt" & trend == 2) {
if(d == 1) {
xmat <- cbind(rep(1, n), coords[, 2], (coords[, 2])^2)
xmati <- cbind(rep(1, ni), locations[, 2], (locations[,
2])^2)
}
else {
xmat <- cbind(rep(1, n), coords[, 1], coords[, 2],
(coords[, 1])^2, (coords[, 2])^2, coords[, 1] * coords[,
2])
xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2], (locations[, 1])^2, (locations[, 2])^2, locations[,
1] * locations[, 2])
}
iviv <- solve(crossprod(xmat,iv) %*% xmat)
txiv <- crossprod(xmat,iv)
beta <- iviv %*% txiv %*% data
mkt <- xmat %*% beta
}
}
if(m0 != "kte")
mktelocations <- "No external trend"
else {
if(m0 == "kte") {
mktelocations <- rep(0, ni)
xmat <- cbind(rep(1, n), ktedata)
xmati <- cbind(rep(1, ni), ktelocations)
iviv <- solve(crossprod(xmat,iv) %*% xmat)
txiv <- crossprod(xmat,iv)
beta <- iviv %*% txiv %*% data
mkte <- xmat %*% beta
}
}
for (i in 1:ni) {
if(messages.screen) {
if(ni < 11)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
else {
if(ni < 101 & (i%%10 == 1))
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(ni > 100 & i%%100 == 1)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(i == ni)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
}
}
coords0 <- cbind((coords.c[, 1] - locations.c[i, 1]), (coords.c[, 2] - locations.c[i,
2]))
dm0 <- sqrt(coords0[, 1]^2 + coords0[, 2]^2)
v0 <- cov.spatial(obj = dm0, cov.model = cov.model,
kappa = kappa, cov.pars = cov.pars)
v0[dm0 < dist.epsilon] <- micro.scale + sigmasq
tv0 <- t(v0)
v0iv <- crossprod(v0, iv)
v0ivv0 <- v0iv %*% v0
skw <- crossprod(v0,iv)
wofmean[i] <- 1 - sum(skw)
##
## 4.2.1 Simple kriging with known mean
##
if(mode(m0) == "numeric") {
dif[i] <- skw %*% (data - m0)
est[i] <- m0 + dif[i]
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0
else kvar[i] <- tausq + sigmasq - v0ivv0
sumw[i] <- sum(skw)
}
##
## 4.2.2 Simple kriging with data average mean
##
if(m0 == "av") {
dif[i] <- skw %*% (data - av)
est[i] <- av + dif[i]
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0
else kvar[i] <- tausq + sigmasq - v0ivv0
sumw[i] <- sum(((tone/n) + skw - ((skw %*% one %*%
tone)/n)))
}
##
## 4.2.3 Ordinary kriging (or SK with G.L.S. mean)
##
if(m0 == "ok") {
dif[i] <- skw %*% (data - ml)
est[i] <- ml + dif[i]
redu <- as.vector(1 - toneiv %*% v0)
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0 + (redu %*%
den %*% redu)
else kvar[i] <- tausq + sigmasq - v0ivv0 + (
redu %*% den %*% redu)
sumw[i] <- sum((den %*% one + tv0 - v0iv %*%
one %*% den %*% tone) %*% iv)
}
##
## 4.2.4 Universal Kriging (or Kriging with trend model)
##
if(m0 == "kt") {
dif[i] <- skw %*% (data - mkt)
est[i] <- xmati[i, ] %*% beta + dif[i]
redu <- as.vector(xmati[i, ]) - as.vector(
txiv %*% v0)
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0 + (redu %*%
iviv %*% redu)
else kvar[i] <- tausq + sigmasq - v0ivv0 + (
redu %*% iviv %*% redu)
sumw[i] <- sum(skw + xmati[i, ] %*% iviv %*%
txiv - skw %*% xmat %*% iviv %*% txiv)
mktlocations[i] <- xmati[i, ] %*% beta
}
##
## 4.2.5 Kriging with external trend
##
if(m0 == "kte") {
dif[i] <- skw %*% (data - mkte)
est[i] <- xmati[i, ] %*% beta + dif[i]
redu <- as.vector(xmati[i, ]) - as.vector(
txiv %*% v0)
if(signal == TRUE)
kvar[i] <- sigmasq - v0ivv0 - (redu %*%
iviv %*% redu)
else kvar[i] <- tausq + sigmasq - v0ivv0 + (
redu %*% iviv %*% redu)
sumw[i] <- sum(skw + xmati[i, ] %*% iviv %*%
txiv - skw %*% xmat %*% iviv %*% txiv)
mktelocations[i] <- xmati[i, ] %*% beta
}
NULL
}
message <- "Kriging performed using global neighbourhood"
if(messages.screen)
cat(paste(message,"\n"))
results <- list(predict = est, krige.var = kvar, dif = dif, summary = means,
ktrend = mktlocations, ktetrend = mktelocations, beta = beta,
wofmean = wofmean)
}
else {
nwin <- min(n, nwin)
avwin <- rep(0, ni)
sdwin <- rep(0, ni)
mlwin <- rep(0, ni)
kmsdwin <- rep(0, ni)
estwin <- rep(0, ni)
difwin <- rep(0, ni)
kvarwin <- rep(0, ni)
sumwwin <- rep(0, ni)
wofmean <- rep(0, ni)
if(m0 != "kt")
mkt <- "Constant position trend"
else mkt <- rep(0, ni)
if(m0 != "kte")
mkte <- "No external trend"
else mkte <- rep(0, ni)
if(m0 != "kt" & m0 != "kte")
betawin <- "No polynomial or external trend"
if(m0 == "kt") {
if(trend == 1) {
if(d == 1)
xmati <- cbind(rep(1, ni), locations[, 2])
else xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2])
}
if(trend == 2) {
if(d == 1)
xmati <- cbind(rep(1, ni), locations[, 2], locations[,
2]^2)
else xmati <- cbind(rep(1, ni), locations[, 1], locations[,
2], (locations[, 1])^2, (locations[, 2])^2, locations[, 1] *
locations[, 2])
}
betawin <- matrix(0, nrow = (ncol(xmati) * ni), ncol = ncol(xmati))
}
if(m0 == "kte") {
xmati <- cbind(rep(1, ni), ktelocations)
if(is.vector(ktedata) == TRUE)
betawin <- matrix(0, nrow = (2 * ni), ncol = 2)
else betawin <- matrix(0, nrow = ((ncol(ktedata) +
1) * ni), ncol = (ncol(ktedata) + 1))
}
for (i in 1:ni) {
temp.win <- .ksline.aux.1(coords = coords, coords.c = coords.c,
data = data, n = n,
locations = locations[i, ],
locations.c = locations.c[i, ],
cov.pars = cov.pars, nugget = nugget,
cov.model = cov.model, kappa = kappa, m0 = m0,
nwin = nwin, trend = trend, d = d, ktedata =
ktedata, ktelocations = ktelocations,
micro.scale = micro.scale,
location.number = i, xmati = xmati[i, ],
mkte = NULL, mkt = NULL, betawin = NULL,
signal = signal, dist.epsilon = dist.epsilon)
avwin[i] <- temp.win$avwin
sdwin[i] <- temp.win$sdwin
mlwin[i] <- temp.win$mlwin
kmsdwin[i] <- temp.win$kmsdwin
estwin[i] <- temp.win$estwin
difwin[i] <- temp.win$difwin
kvarwin[i] <- temp.win$kvarwin
sumwwin[i] <- temp.win$sumwwin
wofmean[i] <- temp.win$wofmean
if(m0 == "kt")
mkt[i] <- temp.win$mkt
if(m0 == "kte")
mkte[i] <- temp.win$mkte
if(m0 == "kt" | m0 == "kte")
betawin[i, ] <- temp.win$betawin
if(messages.screen) {
if(ni < 11)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
else {
if(ni < 101 & (i%%10 == 1))
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(ni > 100 & i%%100 == 1)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
if(i == ni)
cat(paste("ksline: kriging location: ", i, "out of",
ni, "\n"))
}
}
}
message <- "kriging performed in moving neighbourhood"
if(messages.screen)
cat(paste(message,"\n"))
results <- list(predict = estwin, krige.var = kvarwin, dif = difwin,
avtrend = avwin, sd = sdwin, oktrend = mlwin, oksd = kmsdwin,
ktrend = mkt, ktetrend = mkte, beta = betawin, sumw = sumwwin,
wofmean = wofmean)
}
if(abs(lambda - 1) > 0.001) {
if(messages.screen)
cat("Back-transforming the predicted mean and variance.\n")
if(abs(lambda) < 0.001) {
predict.transf <- results$predict
results$predict <- exp(predict.transf) - 0.5 * results$krige.var
results$krige.var <- (exp(2 * predict.transf - results$krige.var)) * (exp(results$krige.var) - 1)
}
if(lambda > 0.001) {
if(messages.screen)
cat("Back-transformation by simulating from the normal predictive distribution\n")
##ap.warn <- options()$warn
##options(warn = -1)
temp.data <- suppressWarnings(matrix(rnorm(ni * n.samples.backtransform,
mean = results$predict,
sd = sqrt(results$krige.var)),
nrow = ni))
##options(warn = ap.warn)
temp.data[(results$krige.var == 0), ] <- results$predict[(results$krige.var == 0)]
temp.data[temp.data < -1/lambda] <- -1/lambda
temp.data <- ((temp.data * lambda) + 1)^(1/lambda)
### temp.data[is.na(temp.data)] <- Inf
results$predict <- as.vector(rowMeans(temp.data))
results$krige.var <- as.vector(apply(temp.data, 1, var))
}
if(lambda < -0.001) {
cat("Resulting distribution has no mean for lambda < 0 - back transformation not performed\n"
)
}
}
results$locations <- locations
results$message <- message
results$call <- call.fc
attr(results, 'sp.dim') <- ifelse(krige1d, "1d", "2d")
attr(results, "prediction.locations") <- call.fc$locations
if(!is.null(call.fc$borders))
attr(results, "borders") <- call.fc$borders
oldClass(results) <- c("kriging")
return(invisible(results))
}
".ksline.aux.1" <-
function (coords, coords.c, data, n, locations, locations.c, cov.pars,
nugget, cov.model, kappa,
m0, nwin, trend, d, ktedata, ktelocations, mbased,
micro.scale, location.number,
xmati, mkte, mkt, betawin, signal, dist.epsilon)
{
i <- location.number
sigmasq <- cov.pars[1]
phi <- cov.pars[2]
tausq <- nugget
coords0 <- cbind((coords.c[, 1] - locations.c[1]), (coords.c[, 2] -
locations.c[2]))
dm0 <- sqrt(coords0[, 1]^2 + coords0[, 2]^2)
coordswin <- coords[order(dm0)[1:nwin], ]
coordswin.c <- coords.c[order(dm0)[1:nwin], ]
datawin <- data[order(dm0)[1:nwin]]
ivwin <- varcov.spatial(coords = coordswin.c, cov.model = cov.model,
kappa = kappa, nugget = nugget, cov.pars = cov.pars, inv = TRUE,
det = FALSE, func.inv = "cholesky", only.decomposition = FALSE)$inverse
avwin <- mean(datawin)
sdwin <- sqrt(var(datawin))
onewin <- rep(1, nwin)
toneivwin <- crossprod(onewin, ivwin)
denwin <- solve(toneivwin %*% onewin)
mlwin <- denwin %*% toneivwin %*% datawin
kmsdwin <- sqrt(denwin)
coords0win <- cbind((coordswin[, 1] - locations[1]), (coordswin[, 2] -
locations[2]))
coords0win.c <- cbind((coordswin.c[, 1] - locations.c[1]), (coordswin.c[
, 2] - locations.c[2]))
dm0win <- sqrt(coords0win.c[, 1]^2 + coords0win.c[, 2]^2)
v0win <- cov.spatial(obj = dm0win, cov.model = cov.model, kappa = kappa,
cov.pars = cov.pars)
v0win[dm0win < dist.epsilon] <- micro.scale + sigmasq
skwwin <- crossprod(v0win, ivwin)
wofmean <- 1 - sum(skwwin)
if(m0 == "kt" & trend == 1) {
if(d == 1)
xmatwin <- cbind(rep(1, nwin), coordswin[, 2])
else xmatwin <- cbind(rep(1, nwin), coordswin[, 1], coordswin[
, 2])
txivwin <- crossprod(xmatwin, ivwin)
ivivwin <- solve(txivwin %*% xmatwin)
betawin <- ivivwin %*% txivwin %*% datawin
mktwin <- xmatwin %*% betawin
mkt <- xmati %*% betawin
}
if(m0 == "kt" & trend == 2) {
if(d == 1)
xmatwin <- cbind(rep(1, nwin), coordswin[, 2], (
coordswin[, 2])^2)
else xmatwin <- cbind(rep(1, nwin), coordswin[, 1], (coordswin[
, 1])^2, coordswin[, 2], (coordswin[, 2])^
2, coordswin[, 1] * coordswin[, 2])
xmatwin.cent <- xmatwin
xmatwin.cent[, 2] <- xmatwin.cent[, 2] - mean(xmatwin[, 2])
xmatwin.cent[, 3] <- xmatwin.cent[, 3] - mean(xmatwin[, 3])
ivivwin <- solve(crossprod(xmatwin.cent, ivwin) %*%
xmatwin.cent)
txivwin <- crossprod(xmatwin.cent, ivwin)
betawin <- ivivwin %*% txivwin %*% datawin
betawin <- mean(datawin) - crossprod(betawin, c(0, mean(xmatwin[
, 2]), mean(xmatwin[, 3])))
mktwin <- xmatwin %*% betawin
mkt <- xmati %*% betawin
}
if(m0 == "kte") {
if(is.vector(ktedata))
ktedatawin <- ktedata[order(dm0)[1:nwin]]
else ktedatawin <- ktedata[order(dm0)[1:nwin], ]
xmatwin <- cbind(rep(1, nwin), ktedatawin)
ivivwin <- solve(crossprod(xmatwin, ivwin) %*% xmatwin)
txivwin <- crossprod(xmatwin, ivwin)
betawin <- ivivwin %*% txivwin %*% datawin
mktewin <- xmatwin %*% betawin
mkte <- xmati %*% betawin
}
##
## Simple kriging with know mean
##
if(mode(m0) == "numeric") {
difwin <- skwwin %*% (data - m0)
# estwin <- m0win + difwin
estwin <- m0 + difwin
if(signal)
kvarwin <- sigmasq - crossprod(v0win, ivwin) %*% v0win
else kvarwin <- tausq + sigmasq - crossprod(v0win, ivwin) %*%
v0win
sumwwin <- sum(skwwin)
}
##
## 4.2.2 Simple kriging with data average mean
##
if(m0 == "av") {
difwin <- skwwin %*% (datawin - avwin)
estwin <- avwin + difwin
if(signal)
kvarwin <- sigmasq - crossprod(v0win, ivwin) %*% v0win
else kvarwin <- tausq + sigmasq - crossprod(v0win, ivwin) %*%
v0win
sumwwin <- sum(((t(onewin)/nwin) + skwwin - ((skwwin %*% onewin %*%
t(onewin))/n)))
}
##
## Ordinary kriging (or SK with G.L.S. mean)
##
if(m0 == "ok") {
difwin <- skwwin %*% (datawin - mlwin)
estwin <- mlwin + difwin
redu <- as.vector(1 - toneivwin %*% v0win)
if(signal)
kvarwin <- sigmasq - v0win %*% ivwin %*% v0win + (
redu %*% denwin %*% redu)
else kvarwin <- tausq + sigmasq - v0win %*% ivwin %*% v0win +
(redu %*% denwin %*% redu)
sumwwin <- sum((denwin %*% onewin + t(v0win) - crossprod(v0win,
ivwin) %*% onewin %*% denwin %*% t(onewin)) %*% ivwin)
}
##
## Universal Kriging (or Kriging with trend model)
##
if(m0 == "kt") {
difwin <- skwwin %*% (datawin - mktwin)
estwin <- mkt + difwin
xmati <- as.vector(xmati)
redu <- as.vector(xmati) - as.vector(txivwin %*% v0win)
if(signal)
kvarwin <- sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
else kvarwin <- tausq + sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
sumwwin <- sum(skwwin + xmati %*% ivivwin %*% txivwin - skwwin %*%
xmatwin %*% ivivwin %*% txivwin)
}
##
## Kriging with external trend
##
if(m0 == "kte") {
difwin <- skwwin %*% (datawin - mktewin)
estwin <- mkte + difwin
xmati <- as.vector(xmati)
redu <- as.vector(xmati) - as.vector(txivwin %*% v0win)
if(signal)
kvarwin <- sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
else kvarwin <- tausq + sigmasq - (v0win %*% ivwin %*% v0win) + (redu %*% ivivwin %*% redu)
sumwwin <- sum(skwwin + xmati %*% ivivwin %*% txivwin - skwwin %*%
xmatwin %*% ivivwin %*% txivwin)
}
##
##
##
results <- list(avwin = avwin, sdwin = sdwin, mlwin = mlwin, kmsdwin =
kmsdwin, wofmean = wofmean, betawin = betawin, mkt = mkt, mkte
= mkte, difwin = difwin, estwin = estwin, kvarwin = kvarwin,
sumwwin = sumwwin)
return(results)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.