-
Notifications
You must be signed in to change notification settings - Fork 4
/
zonalDT.R
84 lines (83 loc) · 3.01 KB
/
zonalDT.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
##' A fast \code{data.table}-based alternative to
##' \code{\link[raster:zonal]{raster::zonal()}}.
##'
##' @title Speedy Zonal Statistics
##' @param x A \code{Raster*} to the totality of whose values
##' \code{fun} should be applied within each zone.
##' @param z A categorical \code{RasterLayer} with codes representing
##' zones.
##' @param fun A name or character string giving the function to be
##' applied to summarize the values by zone. It needs to return a
##' single (or at least a length-one vector). If \code{x} might
##' contain any \code{NA} values, it should be equipped to handle
##' them. For large rasters, this function needs to be one, like
##' \code{sum()} whose value is the same even if carried out in a
##' two-stage application (i.e. first to data subsets and then to
##' the results of those subset applications).
##' @param na.rm Logical. If \code{TRUE}, \code{NA} values in \code{x}
##' are ignored.
##' @return A \code{data.table} with a summary value for each zone.
##' @export
##' @author Joshua O'Brien
##' @examples
##' r <- raster(ncols = 10, nrows = 10)
##' r[] <- runif(ncell(r)) * 1:ncell(r)
##' z <- r
##' z[] <- rep(1:5, each = 20)
##' ## for big files, use a character value rather than a function
##' zonalDT(r, z, "sum")
##'
##' ## for smaller files you can also provide a function
##' zonalDT(r, z, mean)
##' zonalDT(r, z, min)
##'
##' ## multiple layers
##' zonalDT(stack(r, r*10), z, "sum")
zonalDT <- function(x, z, fun = sum, na.rm = TRUE) {
compareRaster(x, z)
fun <- match.fun(fun)
nx <- nlayers(x)
if (canProcessInMemory(x)) {
ll <- vector(mode = "list", length = nx)
names(ll) <- names(x)
DT <- data.table(x = getValues(x),
z = getValues(z))
for (j in 1:nx) {
X <- DT[, fun(.SD[[j]], na.rm = na.rm), by = z]
names(X) <- c("z", names(x)[j])
setkey(X, "z")
ll[[j]] <- X
}
X <- Reduce(merge, ll)
if (na.rm) {
X <- X[complete.cases(X), ]
}
return(X)
} else {
tr <- blockSize(x)
res <- vector(mode = "list", length = tr$n)
for (i in 1:tr$n) {
ll <- vector(mode = "list", length = nx)
names(ll) <- names(x)
DT <- data.table(x = getValues(x, row = tr$row[i],
nrows = tr$nrows[i]),
z = getValues(z, row = tr$row[i],
nrows = tr$nrows[i]))
for (j in 1:nx) {
X <- DT[, fun(.SD[[j]], na.rm = na.rm), by = z]
names(X) <- c("z", names(x)[j])
setkey(X, "z")
ll[[j]] <- X
}
res[[i]] <- Reduce(merge, ll)
}
X <- rbindlist(res)
X <- X[, lapply(.SD, fun, na.rm = na.rm), by = "z",
.SDcols = names(X)[-1]]
setkey(X, "z")
if (na.rm) {
X <- X[complete.cases(X), ]
}
X
}
}