https://github.com/cran/fields
Tip revision: 04279c16ce718b82025aae97f6fcd30ba3a8b1c5 authored by Doug Nychka on 05 September 2011, 20:18:33 UTC
version 6.6
version 6.6
Tip revision: 04279c1
stationary.taper.cov.R
# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"stationary.taper.cov" <- function(x1, x2, Covariance = "Exponential",
Taper = "Wendland", Dist.args = NULL, Taper.args = NULL,
theta = 1, V = NULL, C = NA, marginal = FALSE, spam.format = TRUE,
verbose = FALSE, ...) {
# get covariance function arguments from call
Cov.args <- list(...)
# coerce x1 and x2 to matrices
if (is.data.frame(x1))
x1 <- as.matrix(x1)
if (!is.matrix(x1))
x1 <- matrix(c(x1), ncol = 1)
if (missing(x2))
x2 <- x1
if (is.data.frame(x2))
x2 <- as.matrix(x1)
if (!is.matrix(x2))
x2 <- matrix(c(x2), ncol = 1)
d <- ncol(x1)
n1 <- nrow(x1)
n2 <- nrow(x2)
# Default taper arguments that are particular to the Wendland.
# Make sure dimension argument is added.
if (Taper == "Wendland") {
if (is.null(Taper.args)) {
Taper.args <- list(theta = 1, k = 2, dimension = ncol(x1))
}
if (is.null(Taper.args$dimension)) {
Taper.args$dimension <- ncol(x1)
}
}
#
# Add in general defaults for taper arguments if not Wendland
# theta = 1.0 is the default range for the taper.
if (is.null(Taper.args)) {
Taper.args <- list(theta = 1)
}
#
# separate out a single scalar transformation and a
# more complicated scaling and rotation.
# this is done partly to have the use of great circle distance make sense
# by applying the scaling _after_ finding the distance.
#
# flag for great circle distance
great.circle <- ifelse(is.null(Dist.args$method), FALSE,
Dist.args$method == "greatcircle")
# check form of theta
if (length(theta) > 1) {
stop("theta as a matrix has been depreciated, use the V argument")
}
#
# following now treats V as a full matrix for scaling and rotation.
#
if (!is.null(V)) {
# try to catch error of mixing great circle distance with a
# linear scaling of coordinates.
if (theta != 1) {
stop("can't specify both theta and V!")
}
if (great.circle) {
stop("Can not mix great circle distance\nwith general scaling (V argument or vecotr of theta's)")
}
x1 <- x1 %*% t(solve(V))
x2 <- x2 %*% t(solve(V))
}
#
# locations are now scaled and rotated correctly
# copy taper range
if (great.circle) {
# set the delta cutoff to be in scale of angular latitude.
# figure out if scale is in miles or kilometers
miles <- ifelse(is.null(Dist.args$miles), TRUE, Dist.args$miles)
delta <- (180/pi) * Taper.args$theta/ifelse(miles, 3963.34,
6378.388)
}
else {
delta <- Taper.args$theta
}
if (length(delta) > 1) {
stop("taper range must be a scalar")
}
#NOTE tapering is applied to the _scaled_ locations.
# now apply covariance function to pairwise distance matrix, or multiply
# by C vector or just find marginal variance
if (!marginal) {
# find nearest neighbor distances based on taper threshhold.
# This is hardwired to 'nearest.dist' function from spam.
# note that delta is taken from the taper range not theta or V
sM <- do.call("nearest.dist", c(list(x1, x2, delta = delta,
upper = NULL), Dist.args))
# sM@entries are the pairwise distances up to distance taper.range.
# apply covariance and taper to these.
# note rescaling by theta and taper ranges.
sM@entries <- do.call(Covariance, c(list(d = sM@entries/theta),
Cov.args)) * do.call(Taper, c(list(d = sM@entries),
Taper.args))
# if verbose print out each component separately
if (verbose) {
print(sM@entries/theta)
print(do.call(Covariance, c(list(d = sM@entries/theta),
Cov.args)))
print(do.call(Taper, c(list(d = sM@entries), Taper.args)))
}
if (is.na(C[1])) {
# decide whether to return sM in spam sparse form or as a full matrix
if (spam.format) {
return(sM)
}
else {
return(as.matrix(sM))
}
}
else {
# other option is to sparse multiply cross covariance by C
return(sM %*% C)
}
}
else {
# find marginal variance and return a vector.
sigma2 <- do.call(Covariance, c(list(d = 0), Cov.args)) *
do.call(Taper, c(list(d = 0), Taper.args))
return(rep(sigma2, nrow(x1)))
}
# should not get here!
}