swh:1:snp:33a53053e50f7abe7d281cc0c803be827debf4a3
Tip revision: 48d0b1457c4d7f60524adbe61de8b574c2eb3cc6 authored by Edzer Pebesma on 22 June 2012, 17:21:36 UTC
version 1.0-12
version 1.0-12
Tip revision: 48d0b14
variogram.default.R
# $Id: variogram.default.q,v 1.29 2009-11-02 21:33:17 edzer Exp $
"variogram.default" <-
function(object, locations, X, cutoff, width = cutoff/15.0, alpha = 0,
beta = 0, tol.hor = 90/length(alpha), tol.ver = 90/length(beta),
cressie = FALSE, dX = numeric(0), boundaries = numeric(0),
cloud = FALSE, trend.beta = NULL, debug.level = 1, cross = TRUE,
grid, map = FALSE, g = NULL, ..., projected = TRUE, lambda = 1.0,
verbose = FALSE, covariogram = FALSE, PR = FALSE, pseudo = -1)
{
id1 = id2 = 0
ret = NULL
if (missing(cutoff)) {
if (is.logical(map) && map == TRUE)
stop("for variogram maps, supply at least a cutoff")
cutoff = numeric(0)
}
else if (width <= 0)
stop("argument width should be positive")
if (missing(width))
width = numeric(0)
if (cloud == TRUE)
width = 0
if (is.logical(map) && map == TRUE) {
cells.dim = length(seq(-cutoff, cutoff, by = width))
grid.topology = GridTopology(rep(-cutoff, 2), rep(width, 2), rep(cells.dim, 2))
map = SpatialGrid(grid = grid.topology)
}
if (any(is.na(boundaries)))
stop("no NA values allowed in boundaries")
.Call("gstat_init", as.integer(debug.level))
if (!is.null(g) && !is.null(g$set))
gstat.load.set(g$set)
id.names = NULL
if (is.list(object) && is.list(locations)) {
nvars = length(object)
if (!is.null(names(object)))
id.names = names(object)
for (i in 1:nvars) {
if (missing(X))
Xloc = rep(1, length(object[[i]]))
else Xloc = X[[i]]
t.beta = numeric(0)
if (!is.null(trend.beta) && length(trend.beta) > 0)
t.beta = trend.beta[[i]]
else t.beta = numeric(0)
if (missing(grid) || !is.list(grid)) {
if (!is.null(g) && gridded(g$data[i]$data))
grd = unlist(gridparameters(g$data[i]$data))
else
grd = numeric(0)
} else
grd = grid[[i]]
.Call("gstat_new_data", as.double(object[[i]]),
as.double(locations[[i]]), as.double(Xloc),
as.integer(1), as.double(t.beta), as.integer(-1),
as.integer(0), as.double(-1), as.integer(1),
double(0), grd, as.integer(0), as.integer(projected),
as.integer(0), as.double(lambda), as.integer(0))
if (!is.null(g) && !is.null(g$model[[id.names[i]]]))
load.variogram.model(g$model[[id.names[i]]], c(i - 1, i - 1))
}
} else
stop("argument object and locations should be lists")
if (inherits(map, "SpatialGrid"))
map = as.double(unlist(gridparameters(as(map, "SpatialGrid"))))
pos = 0
ids = NULL
bnd = NULL
is.direct = NULL
if (PR) {
if (any(object[[1]] <= 0))
warning("pairwise relative variogram assumes non-negative data")
covariogram = 2
}
if (cross == "ONLY") {
stopifnot(nvars >= 2)
id.range = 2:nvars
} else if (cross == TRUE)
id.range = nvars:1
else
id.range = 1:nvars
for (id1 in id.range) {
if (cross == "ST")
id2range = 1
else if (cross == "ONLY")
id2range = 1:(id1-1)
else
id2range = ifelse(cross, 1, id1):id1
for (id2 in id2range) {
if (verbose) cat(".")
if (is.null(id.names))
id = ifelse(id1 == id2, paste(id1), cross.name(id2, id1))
else id = ifelse(id1 == id2, paste(id.names[id1]),
cross.name(id.names[id2], id.names[id1]))
for (a in alpha) {
for (b in beta) {
direction = as.numeric(c(a, b, tol.hor, tol.ver))
# boundaries nees to be passed into the c code only ONCE:
if (is.null(bnd))
bnd = boundaries
else # NOT the first time we're here, so...
bnd = numeric(0)
ret.call = .Call("gstat_variogram",
as.integer(c(id1 - 1, id2 - 1)),
as.numeric(cutoff), as.numeric(width),
as.numeric(direction), as.integer(cressie),
as.numeric(dX), as.numeric(bnd), map,
as.integer(covariogram),
as.integer(pseudo))
if (is.logical(map) && map == FALSE) {
np = ret.call[[1]]
sel = np > 0
n.dir = length(sel[sel])
if (n.dir > 0) {
dist = ret.call[[2]]
gamma = ret.call[[3]]
ret_cw = ret.call[[4]] # Sat Jul 16 16:19:59 CEST 2011
dir.a = rep(a, n.dir)
dir.b = rep(b, n.dir)
ids = c(ids, rep(id, n.dir))
is.direct = c(is.direct, id1 == id2)
df = data.frame(np = np[sel], dist = dist[sel],
gamma = gamma[sel], dir.hor = dir.a, dir.ver = dir.b)
if (pos > 0)
ret[(pos + 1):(pos + n.dir), ] = df
else ret = df
pos = pos + n.dir
}
} else {
if (is.null(ret)) {
ret = data.frame(ret.call[[1]], ret.call[[2]],
ret.call[[4]], ret.call[[3]])
names = c("dx", "dy", id, paste("np", id, sep="."))
} else {
ret = data.frame(ret, ret.call[[4]], ret.call[[3]])
names = c(names, id, paste("np", id, sep="."))
}
names(ret) = names
}
}
}
}
}
if (verbose) cat("\n")
.Call("gstat_exit", NULL)
if (is.logical(map) && map == FALSE) {
if (!is.null(ids)) {
ret$id = factor(ids, levels = unique(ids))
attr(ret, "direct") = data.frame(id = unique(ids), is.direct = is.direct)
if (cloud) {
class(ret) = c("variogramCloud", "data.frame")
attr(ret, ".BigInt") = 2^(4 * .Machine$sizeof.long)
} else {
class(ret) = c("gstatVariogram", "data.frame")
if (length(boundaries) == 0) {
cutoff = ret_cw[1]
width = ret_cw[2]
boundaries = seq(0, cutoff, width)
if (!identical(max(boundaries), cutoff))
boundaries = c(boundaries, cutoff)
}
attr(ret, "boundaries") = boundaries
attr(ret, "pseudo") = ret_cw[3]
row.names(ret) = NULL
}
}
} else {
coordinates(ret) = c("dx", "dy")
gridded(ret) = TRUE
ret = list(map = ret)
class(ret) = c("variogramMap", "list")
}
if (!is.null(ret)) {
if (PR)
attr(ret, "what") = "pairwise relative semivariance"
else if (covariogram)
attr(ret, "what") = "covariance"
else if (cressie)
attr(ret, "what") = "Cressie's semivariance"
else
attr(ret, "what") = "semivariance"
}
ret
}