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