1# Author: Robert J. Hijmans
2# Date :  March 2016
3# Version 1.0
4# Licence GPL v3
5
6# 'whiches' functions based on code by Data Munger:
7# https://stackoverflow.com/questions/36117678/r-raster-how-to-record-ties-using-which-max/36120244#36120244
8.whiches <- function(i, fun=min, na.rm=TRUE) {
9	w <- getOption('warn')
10	on.exit(options('warn'= w))
11	options('warn'=-1)
12	m <- which(i == fun(i, na.rm=na.rm))
13    sum(m * 10^(rev(seq_along(m)) - 1))
14}
15
16
17
18setMethod("whiches.min", "RasterStackBrick",
19	function(x) {
20
21		whichesMin <- function(i) {
22			m <- which(i == min(i, na.rm=TRUE))
23			sum(m * 10^(rev(seq_along(m)) - 1))
24		}
25
26		r <- raster(x)
27		nl <- nlayers(x)
28		if (nl > 9) {
29			stop('you can use only use this function for an object with less than 10 layers')
30		}
31
32		if (canProcessInMemory(x)) {
33			x <- values(x)
34			d <- dim(x)
35			i <- .rowSums(is.na(x), d[1], d[2]) < nl
36			y <- rep(NA, nrow(x))
37			if (sum(i) > 0) {
38				y[i] <- apply(x[i,], 1, whichesMin)
39			}
40			return( setValues(r, y) )
41		} else {
42			tr <- blockSize(x)
43			x <- readStart(x)
44			out <- raster(x)
45			out <- writeStart(out, '')
46			for (i in 1:tr$n) {
47				v <- getValues(x, row=tr$row[i], nrows=tr$nrows[i])
48				d <- dim(v)
49				j <- .rowSums(is.na(v), d[1], d[2]) < nl
50				y <- rep(NA, nrow(v))
51				if (sum(j) > 0) {
52					y[j] <- apply(v[j,], 1, whichesMin)
53				}
54				out <- writeValues(out, y, tr$row[i])
55			}
56			out <- writeStop(out)
57			x <- readStop(x)
58			return(out)
59		}
60	}
61)
62
63
64
65setMethod("whiches.max", "RasterStackBrick",
66	function(x) {
67
68		whichesMax <- function(i) {
69			m <- which(i == max(i, na.rm=TRUE))
70			sum(m * 10^(rev(seq_along(m)) - 1))
71		}
72
73		r <- raster(x)
74		nl <- nlayers(x)
75		if (nl > 9) {
76			stop('you can use only use this function for an object with less than 10 layers')
77		}
78
79		if (canProcessInMemory(x)) {
80			x <- values(x)
81			d <- dim(x)
82			i <- .rowSums(is.na(x), d[1], d[2]) < nl
83			y <- rep(NA, nrow(x))
84			if (sum(i) > 0) {
85				y[i] <- apply(x[i,], 1, whichesMax)
86			}
87			return( setValues(r, y) )
88		} else {
89			tr <- blockSize(x)
90			x <- readStart(x)
91			out <- raster(x)
92			out <- writeStart(out, '')
93			for (i in 1:tr$n) {
94				v <- getValues(x, row=tr$row[i], nrows=tr$nrows[i])
95				d <- dim(v)
96				j <- .rowSums(is.na(v), d[1], d[2]) < nl
97				y <- rep(NA, nrow(v))
98				if (sum(j) > 0) {
99					y[j] <- apply(v[j,], 1, whichesMax)
100				}
101				out <- writeValues(out, y, tr$row[i])
102			}
103			out <- writeStop(out)
104			x <- readStop(x)
105			return(out)
106		}
107	}
108)
109
110