1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
MNF = function(Sigma.Noise, Sigma, x) {
e = eigen(t(solve(t(Sigma), t(Sigma.Noise)))) # Sigma.Noise %*% Sigma^-1
e$sdev = sqrt(e$values)
e$rotation = e$vectors
e$x = x %*% e$rotation
colnames(e$x) = paste0("MNF", 1:ncol(e$x))
e$scale = FALSE
e$center = FALSE
e$vectors = NULL
class(e) = c("mnf", "prcomp")
e
}
mnf = function(x, ...) UseMethod("mnf")
mnf.matrix = function(x, ..., Sigma.Noise, use = "complete.obs") {
if (nrow(x) <= ncol(x)) # do diff
warning("matrix rank deficient: covariance matrices will be singular")
Sigma = cov(x, use = use)
if (missing(Sigma.Noise))
Sigma.Noise = 0.5 * cov(apply(x, 2, diff), use = use)
MNF(Sigma.Noise, Sigma, x)
}
mnf.mts = function(x, ..., use = "complete.obs") {
ret = NextMethod(as.matrix(x))
attr(ret$x, "tsp") = attr(x, "tsp")
class(ret$x) = class(x)
ret
}
mnf.zoo = function(x, ..., use = "complete.obs") {
ret = mnf(as.matrix(x), ..., use = use)
ret$x = zoo(ret$x, index(x))
ret
}
mnf.SpatialPixelsDataFrame = function(x, ..., use = "complete.obs") {
mnf.SpatialGridDataFrame(x, ..., use = use) # gets $x right
}
mnf.SpatialGridDataFrame = function(x, ..., Sigma.Noise, use = "complete.obs") {
if (missing(Sigma.Noise)) {
get.dx = function(x) as.vector(x[-1,] - x[-nrow(x),]) # Delta_h
get.dy = function(x) as.vector(x[,-1] - x[,-ncol(x)]) # Delta_v
a = as(x, "array")
cov.dx = 0.5 * cov(apply(a, 3, get.dx), ..., use = use)
cov.dy = 0.5 * cov(apply(a, 3, get.dy), ..., use = use)
Sigma.Noise = (cov.dx + cov.dy)/2.0
}
ret = MNF(Sigma.Noise, cov(x@data, ..., use = use), as.matrix(x@data))
# put data into original structure:
x@data = as.data.frame(ret$x)
ret$x = x
ret
}
mnf.RasterStack = function(x, ..., use = "complete.obs")
mnf.SpatialGridDataFrame(as(x, "SpatialGridDataFrame"), ..., use = use)
mnf.RasterBrick = function(x, ..., use = "complete.obs")
mnf.SpatialGridDataFrame(as(x, "SpatialGridDataFrame"), ..., use = use)
mnf.STSDF = function(x, ..., use = "complete.obs", mode = "temporal") {
NextMethod(as(x, "STFDF"))
}
mnf.STFDF = function(x, ..., use = "complete.obs", mode = "temporal") {
if (dim(x)[3] != 1)
stop("select a single attribute")
if (mode == "temporal") {
ret = mnf(as(x, "zoo"), use = use, ...)
x@data[[1]] = as.vector(t(ret$x))
row.names(x@sp) = colnames(ret$x)
ret$x = x
} else if (mode == "spatial") {
sp = as(x, "Spatial")
stopifnot(gridded(sp))
ret = mnf(sp, use = use, ...)
} else
stop("unknown mode")
ret
}
predict.mnf = function(object, newdata, ...) {
if (missing(newdata))
object$x
else
scale(newdata, object$center, object$scale) %*% object$vectors
}
|