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