-
Notifications
You must be signed in to change notification settings - Fork 29
/
raster.change.R
184 lines (179 loc) · 7.1 KB
/
raster.change.R
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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
#' @title Raster change between two nominal rasters
#' @description
#' Compares two categorical rasters with a variety of statistical options
#'
#' @param x A terra SpatRaster
#' @param y A terra SpatRaster for comparison to x
#' @param s Integer or matrix for defining Kernel,
#' must be odd but not necessarily square
#' @param stat Statistic to use in comparison, please see details for
#' options.
#' @param ... Additional arguments passed to terra::focalPairs
#'
#' @details
#' This function provides a various statistics for comparing two classified maps.
#' Valid options are:
#' * kappa - Cohen's Kappa
#' * t.test - Two-tailed paired t-test
#' * cor - Persons Correlation
#' * entropy - Delta entropy
#' * cross-entropy - Cross-entropy loss function
#' * divergence - Kullback-Leibler divergence (relative entropy)
#'
#' Kappa and t-test values < 0 are reported as 0. For a weighted kappa, a matrix must
#' be provided that correspond to the pairwise weights for all values in both rasters.
#' Delta entropy is derived by calculating Shannon's on each focal window then
#' differencing them (e(x) - e(y)). The s argument can be a single scalar, defining
#' a symmetrical kernel, two scalers defining the dimensions of the kernel eg., c(3,5)
#' or a matrix defining the kernel say, resulting from terra::focalMat
#'
#' @return
#' A terra SpatRaster layer containing one of the following layers:
#' * kappa - Kappa or Weighted Kappa statistic (if stat = "kappa")
#' * correlation - Paired t.test statistic (if stat = "cor")
#' * entropy - Local entropy (if stat = "entropy")
#' * divergence - Kullback-Leibler divergence (if stat = "divergence")
#' * cross.entropy - Local Cross-entropy (if stat = "cross.entropy")
#' * t.test - Paired t.test statistic (if stat = "t.test")
#' * p.value - p-value of the paired t.test statistic (if stat = "t.test")
#'
#' @author Jeffrey S. Evans <jeffrey_evans@@tnc.org>
#'
#' @references
#' Cohen, J. (1960). A coefficient of agreement for nominal scales. Educational
#' and Psychological Measurement, 20:37-46
#'
#' McHugh M.L. (2012) Interrater reliability: the kappa statistic.
#' Biochemia medica, 22(3):276–282.
#'
#' Kullback, S., R.A. Leibler (1951). On information and sufficiency. Annals of
#' Mathematical Statistics. 22(1):79–86
#'
#' @examples
#' \donttest{
#' library(sf)
#' library(terra)
#'
#' e <- ext(179407.8, 181087.9, 331134.4, 332332.1)
#' r1 <- rast(e, resolution=20)
#' r1[] <- sample(1:5, ncell(r1), replace=TRUE)
#' r2 <- rast(e, resolution=20)
#' r2[] <- sample(1:5, ncell(r2), replace=TRUE)
#'
#' d = 5 # kernel
#' ( r.kappa <- raster.change(r1, r2, s = d) )
#' ( r.ttest <- raster.change(r1, r2, s = d, stat="t.test") )
#' ( r.ent <- raster.change(r1, r2, s = d, stat="entropy") )
#' ( r.cor <- raster.change(r1, r2, s = d, stat="cor") )
#' ( r.ce <- raster.change(r1, r2, s = d, stat = "cross-entropy") )
#' ( r.kl <- raster.change(r1, r2, s = d, stat = "divergence") )
#'
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(3,2))
#' plot(r.kappa, main="Kappa")
#' plot(r.ttest[[1]], main="Paired t-test")
#' plot(r.ent, main="Delta Entropy")
#' plot(r.cor, main="Rank Correlation")
#' plot(r.kl, main="Kullback-Leibler")
#' plot(r.ce, main="cross-entropy")
#' par(opar)
#' }
#'
#' @md
#' @export raster.change
raster.change <- function(x, y, s = 3, stat = c("kappa", "t.test",
"cor", "entropy", "cross-entropy",
"divergence"), ...) {
stat = stat[1]
if (!inherits(x, "SpatRaster"))
stop(deparse(substitute(x)), " Must be a terra SpatRaster object")
if (!inherits(y, "SpatRaster"))
stop(deparse(substitute(y)), " Must be a terra SpatRaster object")
if(!any((dim(x)[1:2] == dim(y)[1:2])))
stop("Rasters dimensions do not match")
if(!terra::ext(x) == terra::ext(y))
stop("Rasters extents do not match")
if(stat == "wkappa")
stop("Sorry; weighted kappa is not yet implemented")
if(!any( (stat %in% c("kappa", "wkappa", "t.test", "cor",
"entropy", "cross-entropy", "divergence"))) )
stop("Not a valid option for evaluation statistic")
if(!inherits(s, "matrix")) {
if(any(round(s) %% 2 == 0))
stop("Dimensions for kernel must be an odd number")
if(length(s) < 2) s = c(s,s)
s <- matrix(1, s[1], s[2])
} else {
if(any(dim(s) %% 2 == 0))
stop("Dimensions for kernel must be an odd number")
}
cohens <- function(x, y) {
x <- factor(x)
y <- factor(y)
lvl <- unique(c(levels(x), levels(y)))
x <- factor(x, levels = lvl)
y <- factor(y, levels = lvl)
x <- table(x, y)
d <- diag(x)
n <- sum(x)
nc <- ncol(x)
colFreqs <- colSums(x) / n
rowFreqs <- rowSums(x) / n
return( (sum(d)/n - crossprod(colFreqs, rowFreqs)[1])/
(1 - crossprod(colFreqs, rowFreqs)[1]) )
}
divergence <- function(x, y, type = c("Kullback-Leibler", "cross-entropy")) {
type = type[1]
if(!is.vector(x) | !is.vector(y))
stop("x and y must be numeric of character vectors")
if(any(type %in% c("Kullback-Leibler", "cross-entropy")==FALSE))
stop("Not a valid option for statistic type")
q <- table(x) / sum(table(x)) # observed or approximated
p <- table(y) / sum(table(y)) # estimated or probability
classes <- intersect(names(q), names(p))
p <- p[which(names(p) %in% classes)]
q <- q[which(names(q) %in% classes)]
if(type == "cross-entropy") {
return( -sum( q, log(p) ) )
} else if(type == "Kullback-Leibler") {
return( sum( p * log(p / q) ) )
}
}
ent <- function(x,y) {
return((-sum(prop.table(x) * log(prop.table(x)))) -
(-sum(prop.table(y) * log(prop.table(y)))))
}
ttest <- function(x,y) {
tt <- tryCatch(stats::t.test(x,y),
error = function(e) {})
if(is.null(tt)) {
tt <- c(NA,NA)
} else {
tt <- c(tt$statistic, tt$p.value)
}
return(tt)
}
if(stat == "kappa") {
r <- terra::focalPairs(c(x,y), w=s, function(x, y) { cohens(x, y) }, ...)
}
if(stat == "divergence") {
r <- terra::focalPairs(c(x,y), w=s, function(x, y) {
divergence(x, y, type = "Kullback-Leibler") }, ...)
}
if(stat == "cross-entropy") {
r <- terra::focalPairs(c(x,y), w=s, function(x, y) { divergence(x, y,
type = "cross-entropy") }, ...)
}
if(stat == "cor") {
r <- terra::focalPairs(c(x,y), w=s, function(x, y){ stats::cor(x, y,
method = "spearman") }, ...)
}
if(stat == "entropy") {
r <- terra::focalPairs(c(x,y), w=s, function(x, y) { ent(x, y) }, ...)
}
if(stat == "t.test") {
r <- terra::focalPairs(c(x,y), w=s, function(x, y){ ttest(x,y) }, ...)
names(r) <- c("t", "p.value")
}
return(r)
}