1# Author: Robert J. Hijmans
2# Date : March 2009
3# Version 0.9
4# Licence GPL v3
5
6
7
8setMethod('zonal', signature(x='RasterLayer', z='RasterLayer'),
9	function(x, z, fun='mean', digits=0, na.rm=TRUE, ...) {
10
11		# backward compatibility
12		if (!is.null(list(...)$stat)) {
13			stop('argument "stat" was replaced by "fun"')
14		}
15
16		compareRaster(c(x, z))
17		stopifnot(hasValues(z))
18		stopifnot(hasValues(x))
19
20		layernames <- names(x)
21
22		if (canProcessInMemory(x, 3)) {
23			inmem <- TRUE
24		} else {
25			inmem <- FALSE
26		}
27
28
29		if (inmem) {
30			pb <- pbCreate(2, label='zonal', ...)
31			if (isTRUE(try(fun == 'count', silent=TRUE))) {
32				func <- function(x, na.rm) {
33					if (na.rm) {
34						length(stats::na.omit(x))
35					} else {
36						length(x)
37					}
38				}
39			} else {
40				func <- match.fun(fun)
41			}
42			x <- getValues(x)
43			z <- round(getValues(z), digits=digits)
44			pb <- pbStep(pb, 1)
45			alltab <- tapply(x, z, FUN=func, na.rm=na.rm)
46
47			if (is.array(alltab)) { # multiple numbers
48				id <- as.numeric(dimnames(alltab)[[1]])
49				alltab <- matrix(unlist(alltab, use.names = FALSE), nrow=dim(alltab), byrow=TRUE)
50				alltab <- cbind(id, alltab)
51			} else {
52				alltab <- cbind(as.numeric(names(alltab)), alltab)
53			}
54			pb <- pbStep(pb, 2)
55			colnames(alltab)[1] <- 'zone'
56			d <- dim(alltab)[2]
57			if (d==2) {
58				if (is.character(fun)) {
59					colnames(alltab)[2] <- fun[1]
60				} else {
61					colnames(alltab)[2] <- 'value'
62				}
63			} else {
64				colnames(alltab)[2:d] <- paste0('value_', 1:(d-1))
65			}
66
67		} else {
68
69			if (class(fun)[1] != 'character') {
70				stop("RasterLayers cannot be processed in memory.\n You can use fun='sum', 'mean', 'sd', 'min', 'max', or 'count' but not a function")
71			}
72			if (! fun %in% c('sum', 'mean', 'sd', 'min', 'max', 'count')) {
73				stop("fun can be 'sum', 'mean', 'sd', 'min', 'max', or 'count'")
74			}
75			sdtab <- FALSE
76			counts <- FALSE
77			if (fun == 'count') {
78				func1 <- function(x, na.rm) {
79					if (na.rm) {
80						length(stats::na.omit(x))
81					} else {
82						length(x)
83					}
84				}
85				func2 <- sum
86			} else {
87				func1 <- func2 <- match.fun(fun)
88			}
89			if ( fun == 'mean' | fun == 'sd') {
90				func1 <- func2 <- sum
91				counts <- TRUE
92				if (fun == 'sd') {
93					sdtab <- TRUE
94				}
95			}
96
97			alltab <- array(dim=0)
98			sqtab <- cnttab <- alltab
99
100			tr <- blockSize(x, n=2)
101			pb <- pbCreate(tr$n, label='zonal', ...)
102
103			#nc <- nlayers(x)
104			#nc1 <- nc + 1
105			#nc2 <- 2:nc1
106			#nc2 <- 2
107			x <- readStart(x, ...)
108			z <- readStart(z, ...)
109
110			for (i in 1:tr$n) {
111				d <- cbind(getValues(x, row=tr$row[i], nrows=tr$nrows[i]))
112				Z <- round(getValues(z, row=tr$row[i], nrows=tr$nrows[i]), digits=digits)
113				#cat(i, '\n')
114				#utils::flush.console()
115
116				a <- tapply(d, Z, FUN=func1, na.rm=na.rm)
117				a <- cbind(as.numeric(names(a)), a)
118				alltab <- rbind(alltab, a)
119				if (counts) {
120					if (na.rm) {
121						a <- tapply(d, Z, FUN=function(x)length(stats::na.omit(x)))
122						a <- cbind(as.numeric(names(a)), a)
123						cnttab <- rbind(cnttab, a)
124						if (sdtab) {
125							a <- tapply( d^2, Z, FUN=function(x)sum(stats::na.omit(x)))
126							a <- cbind(as.numeric(names(a)), a)
127							sqtab <- rbind(sqtab, a)
128						}
129					} else {
130						a <- tapply(d, Z, FUN=length)
131						a <- cbind(as.numeric(names(a)), a)
132						cnttab <- rbind(cnttab, a)
133						if (sdtab) {
134							a <- tapply(d^2, Z, FUN=sum)
135							a <- cbind(as.numeric(names(a)), a)
136							sqtab <- rbind(sqtab, a)
137						}
138					}
139				}
140				if (length(alltab) > 10000) {
141					alltab <- tapply(alltab[,2], alltab[,1], FUN=func2, na.rm=na.rm)
142					alltab <- cbind(as.numeric(names(alltab)), alltab)
143					if (counts) {
144						cnttab <- tapply(cnttab[,2], cnttab[,1], FUN=sum, na.rm=na.rm)
145						cnttab <- cbind(as.numeric(names(cnttab)), cnttab)
146						if (sdtab) {
147							sqtab <- tapply(sqtab[,2], sqtab[,1], FUN=sum, na.rm=na.rm)
148							sqtab <- cbind(as.numeric(names(sqtab)), sqtab)
149						}
150					}
151				}
152				pbStep(pb, i)
153			}
154			x <- readStop(x)
155			z <- readStop(z)
156
157			alltab <- tapply(alltab[,2], alltab[,1], FUN=func2, na.rm=na.rm)
158			alltab <- cbind(as.numeric(names(alltab)), alltab)
159			if (counts) {
160				cnttab <- tapply(cnttab[,2], cnttab[,1], FUN=sum)
161				cnttab <- cbind(as.numeric(names(cnttab)), cnttab)
162				alltab[,2] <- alltab[,2] / cnttab[,2]
163				if (sdtab) {
164					sqtab <- tapply(sqtab[,2], sqtab[,1], FUN=sum, na.rm=na.rm)
165					sqtab <- cbind(as.numeric(names(sqtab)), sqtab)
166					alltab[,2] <- sqrt(( (sqtab[,2] / cnttab[,2]) - (alltab[,2])^2 ) * (cnttab[,2]/(cnttab[,2]-1)))
167				}
168
169			}
170			colnames(alltab)[1] <- 'zone'
171			if (is.character(fun)) {
172				colnames(alltab)[2] <- fun
173			} else {
174				colnames(alltab)[2] <- 'value'
175			}
176		}
177		#alltab <- as.matrix(alltab)
178		pbClose(pb)
179		return(alltab)
180	}
181)
182
183#zonal(r, z, 'sd')
184
185
186
187
188setMethod('zonal', signature(x='RasterStackBrick', z='RasterLayer'),
189	function(x, z, fun='mean', digits=0, na.rm=TRUE, ...) {
190
191		# backward compatibility
192		if (!is.null(list(...)$stat)) {
193			stop('argument "stat" was replaced by "fun"')
194		}
195
196		compareRaster(c(x, z))
197		stopifnot(hasValues(z))
198		stopifnot(hasValues(x))
199
200		layernames <- names(x)
201
202		if (canProcessInMemory(x, 3)) {
203			inmem <- TRUE
204		} else {
205			inmem <- FALSE
206		}
207
208		if (inmem) {
209			pb <- pbCreate(2, label='zonal', ...)
210			if (isTRUE(try(fun == 'count', silent=TRUE))) {
211				func <- function(x, na.rm) {
212					if (na.rm) {
213						length(stats::na.omit(x))
214					} else {
215						length(x)
216					}
217				}
218			} else {
219				func <- match.fun(fun)
220			}
221
222			x <- getValues(x)
223			x <- cbind(x, round(getValues(z), digits=digits))
224			pb <- pbStep(pb, 1)
225			alltab <- aggregate(x[,1:(ncol(x)-1)], by=list(x[,ncol(x)]), FUN=func, na.rm=na.rm)
226			fun <- 'value'
227			pb <- pbStep(pb, 2)
228
229		} else {
230
231			if (class(fun)[1] != 'character') {
232				stop("RasterLayers cannot be processed in memory.\n You can use fun='sum', 'mean', 'sd', 'min', 'max', or 'count' but not a function")
233			}
234			if (! fun %in% c('sum', 'mean', 'sd', 'min', 'max', 'count')) {
235				stop("fun can be 'sum', 'mean', 'sd', 'min', 'max', or 'count'")
236			}
237			sdtab <- FALSE
238			counts <- FALSE
239
240			if (fun == 'count') {
241				func1 <- function(x, na.rm) {
242					if (na.rm) {
243						length(stats::na.omit(x))
244					} else {
245						length(x)
246					}
247				}
248				func2 <- sum
249			} else {
250				func1 <- func2 <- match.fun(fun)
251			}
252			if ( fun == 'mean' | fun == 'sd') {
253				func1 <- func2 <- sum
254				counts <- TRUE
255				if (fun == 'sd') {
256					sdtab <- TRUE
257				}
258			}
259
260			alltab <- array(dim=0)
261			sqtab <- cnttab <- alltab
262
263			tr <- blockSize(x, n=2)
264			pb <- pbCreate(tr$n, label='zonal', ...)
265
266			nc <- nlayers(x)
267			nc1 <- nc + 1
268			nc2 <- 2:nc1
269
270			# for a RasterStack it would be more efficient to loop over the layers
271			x <- readStart(x, ...)
272			z <- readStart(z, ...)
273
274			for (i in 1:tr$n) {
275				d <- cbind(getValues(x, row=tr$row[i], nrows=tr$nrows[i]),
276					 round(getValues(z, row=tr$row[i], nrows=tr$nrows[i]), digits=digits))
277				#cat(i, '\n')
278				#utils::flush.console()
279				alltab <- rbind(alltab, aggregate(d[,1:nc], by=list(d[,nc1]), FUN=func1, na.rm=na.rm))
280				if (counts) {
281					if (na.rm) {
282						cnttab <- rbind(cnttab, aggregate(d[,1:nc], by=list(d[,nc1]), FUN=function(x)length(stats::na.omit(x))))
283						if (sdtab) {
284							sqtab <- rbind(sqtab, aggregate( (d[,1:nc])^2, by=list(d[,nc1]), FUN=function(x)sum(stats::na.omit(x))))
285						}
286					} else {
287						cnttab <- rbind(cnttab, aggregate(d[,1:nc], by=list(d[,nc1]), FUN=length))
288						if (sdtab) {
289							sqtab <- rbind(sqtab, aggregate( (d[,1:nc])^2, by=list(d[,nc]), FUN=sum))
290						}
291					}
292				}
293				if (length(alltab) > 10000) {
294					alltab <- aggregate(alltab[,nc2], by=list(alltab[,1]), FUN=func2, na.rm=na.rm)
295					if (counts) {
296						cnttab <- aggregate(cnttab[,nc2], by=list(cnttab[,1]), FUN=sum, na.rm=na.rm)
297						if (sdtab) {
298							sqtab <- aggregate(sqtab[,nc2], by=list(sqtab[,1]), FUN=sum, na.rm=na.rm)
299						}
300					}
301				}
302				pbStep(pb, i)
303			}
304			x <- readStop(x)
305			z <- readStop(z)
306
307			alltab <- aggregate(alltab[,nc2], by=list(alltab[,1]), FUN=func2, na.rm=na.rm)
308			if (counts) {
309				cnttab <- aggregate(cnttab[,nc2], by=list(cnttab[,1]), FUN=sum)
310				alltab[,nc2] <- alltab[,nc2] / cnttab[,nc2]
311				if (sdtab) {
312					sqtab <- aggregate(sqtab[,nc2], by=list(sqtab[,1]), FUN=sum, na.rm=na.rm)
313					alltab[,nc2] <- sqrt(( (sqtab[,nc2] / cnttab[,nc2]) - (alltab[nc2])^2 ) * (cnttab[,nc2]/(cnttab[,nc2]-1)))
314				}
315
316			}
317
318		}
319
320		alltab <- as.matrix(alltab)
321		colnames(alltab)[1] <- 'zone'
322		if (ncol(alltab) > 2) {
323			colnames(alltab)[2:ncol(alltab)] <- layernames
324		} else {
325			colnames(alltab)[2] <- fun[1]
326		}
327		pbClose(pb)
328
329		return(alltab)
330	}
331)
332
333#zonal(r, z, 'sd')
334
335
336