https://github.com/cran/gstat
Tip revision: cc0e158f1e0ce7dce51ae4c6b8e63147c5df3986 authored by Edzer J. Pebesma on 19 October 2004, 03:25:34 UTC
version 0.9-16
version 0.9-16
Tip revision: cc0e158
predict.gstat.R
"predict.gstat" <-
function (object, newdata, block = numeric(0), nsim = 0, indicators = FALSE,
BLUE = FALSE, debug.level = 1, mask, na.action = na.pass, ...)
{
if (missing(object) || length(object$data) < 1)
stop("no data available")
if (!inherits(object, "gstat"))
stop("first argument should be of class gstat")
if (extends(class(newdata), "SpatialDataFrame") && require(sp))
use.sdf = TRUE
else
use.sdf = FALSE
.Call("gstat_init", as.integer(debug.level), PACKAGE = "gstat")
if (!missing(mask)) {
cat("argument mask is deprecated:")
stop("use a missing value pattern in newdata instead")
}
nvars = length(object$data)
new.X = NULL
for (i in 1:length(object$data)) {
name = names(object$data)[i]
d = object$data[[i]]
if (d$nmax == Inf)
nmax = as.integer(-1)
else nmax = as.integer(d$nmax)
nmin = as.integer(max(0, d$nmin))
if (d$maxdist == Inf)
maxdist = as.numeric(-1)
else maxdist = d$maxdist
if (d$dummy) {
tr = terms(d$locations)
if (is.null(d$beta) || length(d$beta) == 0)
stop("dummy data should have beta defined")
if (d$degree != 0)
stop("dummy data cannot have non-zero degree arg; use formula")
loc.dim = length(attr(tr, "term.labels"))
.Call("gstat_new_dummy_data", as.integer(loc.dim),
as.integer(d$has.intercept), as.double(d$beta),
nmax, nmin, maxdist, as.integer(d$vfn)
, PACKAGE = "gstat"
)
} else {
if (is.null(d$weights))
w = numeric(0)
else
w = d$weights
raw = gstat.formula(d$formula, d$locations, d$data)
.Call("gstat_new_data", as.double(raw$y), as.double(raw$locations),
as.double(raw$X), as.integer(raw$has.intercept),
as.double(d$beta), nmax, nmin, maxdist, as.integer(d$vfn),
as.numeric(w), double(0), as.integer(d$degree)
, PACKAGE = "gstat"
)
}
if (!is.null(object$model[[name]]))
load.variogram.model(object$model[[name]], c(i - 1, i - 1))
raw = gstat.formula.predict(d$formula, d$locations, newdata,
na.action = na.action)
if (is.null(new.X))
new.X = raw$X
else new.X = cbind(new.X, raw$X)
if (i > 1) {
for (j in 1:(i - 1)) {
cross = cross.name(names(object$data)[j], name)
if (!is.null(object$model[[cross]]))
load.variogram.model(object$model[[cross]],
c(i - 1, j - 1))
}
}
}
if (!is.null(object$set))
gstat.load.set(object$set)
if (!is.null(object$merge))
gstat.load.merge(object)
if (!is.null(dim(block))) { # i.e., block is data.frame or matrix
block = data.matrix(block) # converts to numeric
block.cols = ncol(block)
} else {
block = as.numeric(block) # make sure it's not integer
block.cols = numeric(0)
}
# handle NA's in the parts of newdata used:
valid.pattern = NULL
if (any(is.na(raw$locations)) || any(is.na(new.X))) {
valid.pattern = !(apply(cbind(raw$locations, new.X), 1,
function(x) any(is.na(x))))
raw$locations.all = raw$locations
raw$locations = as.matrix(raw$locations[valid.pattern, ])
new.X = as.matrix(new.X[valid.pattern, ])
}
if (nsim) {
if (indicators == TRUE)
nsim = -abs(nsim)
# random path: randomly permute row indices
perm = sample(seq(along = new.X[, 1]))
ret = .Call("gstat_predict", as.integer(nrow(as.matrix(new.X))),
as.vector(raw$locations[perm, ]), as.vector(new.X[perm,]),
as.integer(block.cols), as.vector(block),
as.integer(nsim), as.integer(BLUE)
, PACKAGE = "gstat"
)[[1]]
ret = data.frame(cbind(raw$locations,
matrix(ret[order(perm),], nrow(as.matrix(new.X)),
max(2, abs(nsim) * nvars))))
}
else {
ret = .Call("gstat_predict", as.integer(nrow(as.matrix(new.X))),
as.vector(raw$locations), as.vector(new.X), as.integer(block.cols),
as.vector(block), as.integer(nsim), as.integer(BLUE)
, PACKAGE = "gstat"
)[[1]]
ret = data.frame(cbind(raw$locations, ret))
}
.Call("gstat_exit", NULL, PACKAGE = "gstat")
if (!is.null(valid.pattern) && any(valid.pattern)) {
ret.all = data.frame(matrix(NA, length(valid.pattern), ncol(ret)))
ret.all[, 1:ncol(raw$locations.all)] = raw$locations.all
ret.all[valid.pattern, ] = ret
ret = ret.all
}
if (abs(nsim) > 0) {
names.vars = names(object$data)
if (length(names.vars) > 1)
names.vars = paste(rep(names.vars, each = abs(nsim)),
paste("sim", 1:abs(nsim), sep = ""), sep = ".")
else names.vars = paste("sim", 1:abs(nsim), sep = "")
if (abs(nsim) == 1)
ret = ret[, 1:(ncol(ret) - 1)]
} else
names.vars = create.gstat.names(names(object$data))
names(ret) = c(dimnames(raw$locations)[[2]], names.vars)
if (use.sdf) {
coordinates(ret) = dimnames(raw$locations)[[2]]
gridded(ret) = gridded(newdata)
}
return(ret)
}
# call with: create.gstat.names(names(object$data))
# creates the names of the output columns in case of (multivariable) prediction
create.gstat.names <- function(ids, names.sep = ".") {
nvars = length(ids)
names.vars = character(nvars * 2 + nvars * (nvars - 1)/2)
pos = 1
for (i in 1:length(ids)) {
name = ids[i]
names.vars[1 + (i - 1) * 2] = paste(name, "pred", sep = names.sep)
names.vars[2 + (i - 1) * 2] = paste(name, "var", sep = names.sep)
if (i > 1) {
for (j in 1:(i - 1)) {
cross = paste(ids[j], name, sep = names.sep)
names.vars[nvars * 2 + pos] = paste("cov", cross,
sep = names.sep)
pos = pos + 1
}
}
}
return(names.vars)
}