# 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
"vgram" <- function(loc, y, id = NULL, d = NULL, lon.lat = FALSE,
dmax = NULL, N = NULL, breaks = NULL) {
# coerce to matrix
y <- cbind(y)
# if nearest neighbor indices are missing create all possible pairs.
if (is.null(id)) {
n <- nrow(loc)
ind <- rep(1:n, n) > rep(1:n, rep(n, n))
id <- cbind(rep(1:n, n), rep(1:n, rep(n, n)))[ind, ]
}
# if distances are missing calculate these
if (is.null(d)) {
loc <- as.matrix(loc)
if (lon.lat) {
d <- rdist.earth(loc)[id]
}
else {
d <- rdist(loc, loc)[id]
}
}
#
# calculating variogram will average over columns if y is a matrix
#
vg <- 0.5 * rowMeans(cbind((y[id[, 1], ] - y[id[, 2], ])^2),
na.rm = TRUE)
#
#information for returned object
#
call <- match.call()
if (is.null(dmax)) {
dmax <- max(d)
}
od <- order(d)
d <- d[od]
vg <- vg[od]
ind <- d <= dmax & !is.na(vg)
## add a binned variogram if breaks are supplied
out <- list(d = d[ind], vgram = vg[ind], call = call)
if (!is.null(breaks) | !is.null(N)) {
out <- c(out, stats.bin(d[ind], vg[ind], N = N, breaks = breaks))
}
out
}