#····································································
# kriging.R (npsp package)
#····································································
# np.kriging(object, ...)
# np.kriging.default(object, svm, lp.resid, ngrid, ...)
# np.kriging.np.geo(object, ngrid, ...)
# kriging.simple(x, y, newx, svm)
# kriging.simple.solve(x, newx, svm)
#
# (c) R. Fernandez-Casal
#
# NOTE: Press Ctrl + Shift + O to show document outline in RStudio
#····································································
# PENDENTE:
# - spam = TRUE
# - newx = FALSE -> CV? chol2inv.spam
#····································································
#····································································
# np.kriging(object, ...) ----
#····································································
#' Nonparametric (residual) kriging
#'
#' Compute simple kriging or residual kriging predictions
#' (and also the corresponding simple kriging standard errors
# *and cross-validation measures, not yet...*
#' ). Currently, only global (residual) simple kriging is implemented.
#' @aliases kriging
#' @param object object used to select a method:
#' local polynomial estimate of the trend (class \code{\link{locpol.bin}})
#' or nonparametric geostatistical model (class extending \code{\link{np.geo}}).
#' @param ... further arguments passed to or from other methods.
#' @export
#····································································
np.kriging <- function(object, ...) {
UseMethod("np.kriging")
} # S3 generic function np.kriging
#' @rdname np.kriging
#' @method np.kriging default
#' @param svm semivariogram model (of class extending \code{\link{svarmod}}).
# or covariance matrix of the data.
#' @param ngrid number of grid nodes in each dimension.
#' @param lp.resid residuals (defaults to \code{residuals(object)}).
#' @param intermediate logical, determines whether the intermediate computations
#' are included in the output (component \code{kriging}; see Value).
#' These calculations can be reused, e.g. for bootstrap.
#' @return
#' \code{np.kriging()}, and \code{kriging.simple()} when \code{newx} defines
#' gridded data (extends \code{grid.par} or \code{data.grid} classes),
#' returns an S3 object of class \code{krig.grid} (kriging results + grid par.).
#' A \code{\link{data.grid}} object with the additional (some optional) components:
#' \item{kpred}{vector or array (dimension \code{$grid$n}) with the kriging predictions. }
#' \item{ksd}{vector or array with the kriging standard deviations. }
#' \item{kriging}{(if requested) a list with 4 components:
#' \itemize{
#' \item{\code{lambda} matrix of kriging weights (columns correspond with predictions
#' and rows with data)).}
#' \item{\code{cov.est} (estimated) covariance matrix of the data.}
#' \item{\code{chol} Cholesky factorization of \code{cov.est}.}
#' \item{\code{cov.pred} matrix of (estimated) covariances between data (rows)
#' and predictions (columns).}
#' }}
#' When \code{newx} is a matrix of coordinates (where each row is a prediction location),
#' \code{kriging.simple()} returns a list with the previous components (\code{kpred}, \code{ksd}
#' and, if requested, \code{kriging}).
#' @examples
#' geomod <- np.fitgeo(aquifer[,1:2], aquifer$head)
#' krig.grid <- np.kriging(geomod, ngrid = c(96, 96)) # 9216 locations
#' old.par <- par(mfrow = c(1,2))
#' simage(krig.grid, 'kpred', main = 'Kriging predictions',
#' xlab = "Longitude", ylab = "Latitude", reset = FALSE )
#' simage(krig.grid, 'ksd', main = 'Kriging sd', xlab = "Longitude",
#' ylab = "Latitude" , col = hot.colors(256), reset = FALSE)
#' par(old.par)
#' @seealso \code{\link{np.fitgeo}}, \code{\link{locpol}}, \code{\link{np.svar}}.
#' @export
#····································································
np.kriging.default <- function(object, svm, lp.resid = NULL, ngrid = object$grid$n,
intermediate = FALSE, ...) {
stopifnot(inherits(object, "locpol.bin" ))
if (!is.null(object$mask) && is.null(object$window)) # Interpolar mask?
warning("A resized 'object' can not be masked (`ngrid != object$grid$n`).")
if (lp.hd <- any(ngrid != object$grid$n)) {
bin <- with(object$data, binning(x, y, nbin = ngrid, window = object$window))
object <- locpol(bin, h = object$locpol$h)
lp.resid <- residuals(object)
}
if (is.null(lp.resid))
lp.resid <- if(inherits(svm, "fitsvar")) svm$esv$data$y else residuals(object)
krig.grid <- kriging.simple(x = object$data$x, y = lp.resid, newx = object,
svm = svm, intermediate = intermediate)
krig.grid$kpred <- object$est + krig.grid$kpred
krig.grid$trend <- object$est
if (!is.null(object$mask)) is.na(krig.grid$ksd) <- !object$mask
# krig.grid$lp <- object
# krig.grid$residuals <- lp.resid
return(krig.grid)
}
#' @rdname np.kriging
#' @method np.kriging np.geo
#' @export
#····································································
np.kriging.np.geo <- function(object, ngrid = object$grid$n,
intermediate = FALSE, ...) {
np.kriging.default(object, object$svm, object$residuals, ngrid = ngrid,
intermediate = intermediate, ...)
}
#' @rdname np.kriging
#' @inheritParams locpol.default
#' @param x vector/matrix with data locations
#' (each component/row is an observation location).
#' @param newx vector/matrix with the (irregular) locations to predict
#' (each component/row is a prediction location).
#' or an object extending \code{\link{grid.par}}-\code{\link{class}}
#' (\code{\link{data.grid}}).
#' @export
#····································································
# res <- krig.grid$kriging
kriging.simple <- function(x, y, newx, svm, intermediate = FALSE) {
x <- as.matrix(x)
if ( !identical(length(y), nrow(x)) )
stop("arguments 'y' and 'x' do not have the same length.")
gridded <- inherits(newx, "grid.par") || inherits(newx, "data.grid")
if (gridded) {
grid <- if(inherits(newx, "grid.par")) newx else newx$grid
newx <- coords(newx)
}
res <- .kriging.simple.solve(x = x, newx = newx, svm = svm)
kpred <- drop(y %*% res$lambda)
ksd2 <- with(res, cov.est[1,1] - colSums( lambda * cov.pred ))
# ksd2 <- with(res, cov.est[1,1] - spam::colSums.spam( lambda * cov.pred ))
# ksd2 <- with(res, spam::diag(cov.est) - spam::colSums.spam( lambda * cov.pred ))
if (gridded) {
result <- data.grid(kpred = kpred, ksd = sqrt(ksd2), grid = grid)
oldClass(result) <- c("krig.grid", oldClass(result))
} else {
result <- list(kpred = kpred, ksd = sqrt(ksd2))
}
if (intermediate) result$kriging <- res
return(result)
}
#' @rdname npsp-internals
#' @inheritParams kriging.simple
#' @keywords internal
#····································································
.kriging.simple.solve <- function(x, newx, svm) {
# .DNRM2_R <- npsp:::.DNRM2_R
# Evitar recalcular matriz de covarianza y factorizacion
# spam = TRUE
# newx = FALSE -> CV? chol2inv.spam
if (!inherits(svm, "svarmod"))
stop("'svm' must be a semivariogram model ('svarmod' class).")
x <- as.matrix(x)
newx <- as.matrix(newx)
if(ncol(x) != ncol(newx))
stop("arguments 'x' and 'newx' have incompatible dimensions")
# Covariance matrix of the data
# cov.est <- varcov(svm, coords = x)
cov.est <- spam::as.spam(varcov(svm, coords = x))
# Covariances between prediction grid points and data
dist2 <- as.numeric(.DNRM2_R(newx, x))
cov.pred <- matrix(covar(svm, dist2, discretize = TRUE), nrow = nrow(x))
# cov.pred <- spam::as.spam(matrix(covar(svm, dist2, discretize = TRUE), nrow = nrow(x)))
# Kriging system
# library(spam)
if (!requireNamespace("spam")) stop("'spam' package must be installed...")
# chol <- spam::chol.spam(spam::as.spam(cov.est))
chol <- spam::chol.spam(cov.est)
lambda <- spam::solve.spam(chol, spam::as.spam(cov.pred))
# lambda <- spam::solve.spam(chol, cov.pred)
# res <- solve.spam(chol, as.spam(cov.pred))
return(list(lambda = lambda, cov.est = cov.est, chol = chol, cov.pred = cov.pred))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.