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.ST.R
VgmAverage = function(ret) {
# take care of missing rows...
gamma = apply(do.call(cbind, lapply(ret, function(x) x$gamma)), 1, mean,
na.rm = TRUE)
dist = apply(do.call(cbind, lapply(ret, function(x) x$dist)), 1, mean,
na.rm = TRUE)
np = apply(do.call(cbind, lapply(ret, function(x) x$np)), 1, sum,
na.rm = TRUE)
v = data.frame(np = np, dist = dist, gamma = gamma)
class(v) = class(ret[[1]])
attr(v, "boundaries") = attr(ret[[1]], "boundaries")
v
}
StVgmLag = function(formula, data, dt, pseudo, ...) {
d = dim(data)
ret = vector("list", d[2] - dt)
if (dt == 0) {
for (i in 1:d[2])
ret[[i]] = variogram(formula, data[,i],...)
} else {
for (i in 1:(d[2] - dt)) {
obj = gstat(NULL, paste("D", i, sep=""), formula, data[,i])
obj = gstat(obj, paste("D", i+dt, sep=""), formula, data[,i+dt])
ret[[i]] = variogram(obj, cross = "ONLY", pseudo = pseudo, ...)
}
}
VgmAverage(ret)
}
variogramST = function(formula, locations, data, ..., tlags = 0:15,
progress = TRUE, pseudo = TRUE) {
if (missing(data))
data = locations
stopifnot(is(data, "STFDF") || is(data, "STSDF"))
it = index(data@time)
if (is.regular(zoo(1:length(it), it), strict = TRUE)) {
twidth = diff(it)[1]
tlags = tlags[tlags <= min(max(tlags), length(unique(it)) - 1)]
} else {
warning("strictly irregular time steps were assumed to be regular")
twidth = mean(diff(it))
}
ret = vector("list", length(tlags))
obj = NULL
t = twidth * tlags
if (progress)
pb = txtProgressBar(style = 3, max = length(tlags))
for (dt in seq(along = tlags)) {
ret[[dt]] = StVgmLag(formula, data, tlags[dt], pseudo = pseudo, ...)
ret[[dt]]$id = paste("lag", dt, sep="")
if (progress)
setTxtProgressBar(pb, dt)
}
if (progress)
close(pb)
# add time lag:
v = do.call(rbind, ret)
v$timelag = rep(t, sapply(ret, nrow))
if (is(t, "yearmon"))
class(v$timelag) = "yearmon"
b = attr(ret[[2]], "boundaries")
ix = findInterval(v$dist, b)
spacelags = b[1:(length(b)-1)] + diff(b)/2
v$spacelag = spacelags[ix]
if (isTRUE(!is.projected(data)))
attr(v$spacelag, "units") = "km"
class(v) = c("StVariogram", "data.frame")
v
}
plot.StVariogram = function(x, ..., col = bpy.colors(), xlab, ylab, map = TRUE,
convertMonths = FALSE) {
lst = list(...)
if (!is.null(lst$col.regions))
col = lst$col.regions
if (is(x$timelag, "yearmon")) {
if (convertMonths) {
x$timelag = as.numeric(x$timelag) * 12
attr(x$timelag, "units") = "months"
} else
attr(x$timelag, "units") = "years"
}
if (missing(xlab)) {
xlab = "distance"
u = attr(x$spacelag, "units")
if (!is.null(u))
xlab = paste(xlab, " (", u, ")", sep="")
}
if (missing(ylab)) {
ylab = "time lag"
u = attr(x$timelag, "units")
if (!is.null(u))
ylab = paste(ylab, " (", u, ")", sep="")
}
if (!is.null(x$model)) {
v0 = rbind(x, x)
v0$what = c(rep("sample", nrow(x)), rep("model", nrow(x)))
v0$gamma = c(x$gamma, x$model)
x = v0
}
if (map) {
if (!is.null(x$model))
f = gamma ~ spacelag + timelag | what
else
f = gamma ~ spacelag + timelag
levelplot(f, x, xlab = xlab, ylab = ylab, col.regions = col, ...)
} else {
if (!is.null(x$model))
f = gamma ~ dist | what
else
f = gamma ~ dist
xyplot(f, x, groups = x$id, type='b',
auto.key = list(space = "right"), xlab = xlab, ...)
}
}
variogram.ST00 = function(formula, locations, data, ...,
nt = 15, twidth, tcutoff, dX) {
# Sat Jul 16 19:51:37 CEST 2011
if (missing(data))
data = locations
stopifnot(is(data, "STIDF"))
# da = as(data, "STIDF")
it = index(data@time)
if (is.regular(it, strict = TRUE)) {
if (missing(twidth))
twidth = diff(it)[1]
nt = min(nt, length(unique(it))/2)
} else {
if (!(missing(twidth) && missing(tcutoff)))
stop("time is not strictly regular: specify twidth and/or tcutoff")
if (!missing(twidth) && !missing(tcutoff)) {
nt = tcutoff %/% twidth
if (tcutoff %% twidth > 0)
nt = nt+1
}
}
if (missing(tcutoff) && !missing(twidth))
tcutoff = nt * twidth
if (missing(twidth) && !missing(tcutoff))
twidth = tcutoff / nt
# now we have all three; continue with twidth and tcutoff.
#ret = vector("list", 2*nt+1)
ret = vector("list", nt)
obj = NULL
da = data
df = as(da, "data.frame")
dsp = addAttrToGeom(da@sp, df) # now Spatial, with time as attribute
dsp$RES = residuals(lm(formula, df))
t = twidth * (0:(nt-1))
for (dt in 1:nt) {
id = paste("B", dt-1, sep="")
df0 = as(dsp, "SpatialPointsDataFrame")
df0$TIME = df0$time - t[dt] # Back i-1 steps
obj = gstat(obj, id, RES~TIME, df0)
}
if (missing(dX))
dX = twidth/2
v = variogram(obj, dX = dX, cross = "ST", progress = TRUE, ...)
# add time lag:
tlag = substring(as.character(v$id), first = 5)
tlag[tlag == ""] = "0"
v$timelag = t[as.numeric(tlag)+1]
b = attr(v, "boundaries")
ix = findInterval(v$dist, b)
spacelags = b[1:(length(b)-1)] + diff(b)/2
v$spacelag = spacelags[ix]
if (isTRUE(!is.projected(data)))
attr(v$spacelag, "units") = "km"
class(v) = c("StVariogram", "data.frame")
v
}