SCM

SCM Repository

[raster] Annotation of /pkg/raster/R/cellStats.R
ViewVC logotype

Annotation of /pkg/raster/R/cellStats.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3417 - (view) (download)

1 : rhijmans 2239 # Author: Robert J. Hijmans
2 :     # Date : March 2009 / April 2012
3 :     # Version 1.0
4 : rhijmans 316 # Licence GPL v3
5 :    
6 : rhijmans 1105
7 : rhijmans 2239
8 :     .csTextFun <- function(fun) {
9 : rhijmans 3417 if (class(fun)[1] != 'character') {
10 : rhijmans 2239 if (is.primitive(fun)) {
11 :     test <- try(deparse(fun)[[1]], silent=TRUE)
12 :     if (test == '.Primitive(\"sum\")') { fun <- 'sum'
13 :     } else if (test == '.Primitive(\"min\")') { fun <- 'min'
14 :     } else if (test == '.Primitive(\"max\")') { fun <- 'max'
15 :     }
16 :     } else {
17 :     f <- paste(deparse(fun), collapse = "\n")
18 :     if (f == paste(deparse(mean), collapse = "\n")) {
19 :     fun <- 'mean'
20 : rhijmans 3208 } else if (f == paste(deparse(stats::sd), collapse = "\n")) {
21 : rhijmans 2239 fun <- 'sd'
22 :     } else if (f == paste(deparse(range), collapse = "\n")) {
23 :     fun <- 'range'
24 :     }
25 :     }
26 :     }
27 :     return(fun)
28 :     }
29 :    
30 :    
31 : rhijmans 2221
32 :     if (!isGeneric("cellStats")) {
33 :     setGeneric("cellStats", function(x, stat, ...)
34 :     standardGeneric("cellStats"))
35 :     }
36 : rhijmans 1510
37 : rhijmans 2221
38 : rhijmans 2239 setMethod('cellStats', signature(x='RasterStackBrick'),
39 :     function(x, stat='mean', na.rm=TRUE, asSample=TRUE, ...) {
40 : rhijmans 2221
41 :     stopifnot(hasValues(x))
42 : rhijmans 1150
43 : rhijmans 2221 makeMat <- FALSE
44 :     if (nlayers(x) == 1) {
45 :     makeMat <- TRUE
46 :     #return( cellStats(raster(x, values=TRUE, stat=stat, ...) )
47 :     }
48 : rhijmans 1609
49 : rhijmans 2239 stat <- .csTextFun(stat)
50 : rhijmans 2221
51 : rhijmans 1510 if (!inMemory(x)) {
52 :     if (canProcessInMemory(x)) {
53 :     x <- readAll(x)
54 :     }
55 :     }
56 :     if (inMemory(x) ) {
57 : rhijmans 1193 x <- getValues(x)
58 :     if (makeMat) {
59 :     x <- matrix(x, ncol=1)
60 : rhijmans 1105 }
61 : rhijmans 1150
62 : rhijmans 1193 if (class(stat) == 'character') {
63 : rhijmans 1150 if (stat == "mean" ) {
64 : rhijmans 2221 return( colMeans(x, na.rm=na.rm) )
65 :    
66 : rhijmans 1150 } else if (stat == "sum" ) {
67 : rhijmans 2221 return( colSums(x, na.rm=na.rm) )
68 : rhijmans 2881
69 :     } else if (stat == "min" ) {
70 : rhijmans 2907 v <- .colMin(x, na.rm=na.rm)
71 :     names(v) <- names(x)
72 :     return(v)
73 : rhijmans 2881
74 :     } else if (stat == "max" ) {
75 : rhijmans 2907 v <- .colMax(x, na.rm=na.rm)
76 :     names(v) <- names(x)
77 :     return(v)
78 : rhijmans 2221
79 :     } else if (stat == 'countNA') {
80 : rhijmans 3413 warning ("'countNA' is deprecated. Use 'freq(x, value=NA)' instead")
81 : rhijmans 2239 return( colSums(is.na(x)) )
82 :    
83 :     } else if (stat == 'sd') {
84 :    
85 : rhijmans 3208 st <- apply(x, 2, stats::sd, na.rm=na.rm)
86 : rhijmans 2239 if (! asSample) {
87 :     if (na.rm) {
88 :     n <- colSums(! is.na(x))
89 :     } else {
90 :     n <- nrow(x)
91 :     }
92 : rhijmans 3340 st <- sqrt(st^2 * ((n-1)/n))
93 : rhijmans 2239 }
94 : rhijmans 3340
95 : rhijmans 2239 return(st)
96 : rhijmans 2483
97 :     } else if (stat == 'rms') {
98 :     if (na.rm) {
99 :     n <- colSums(! is.na(x))
100 :     } else {
101 :     n <- nrow(x)
102 :     }
103 :     if (asSample) {
104 :     n <- n-1
105 :     }
106 : rhijmans 2872 # st <- apply(x, 2, function(x) sqrt(sum(x^2)/n))
107 :     return( sqrt( apply(x, 2, function(x) sum(x^2))/n ) )
108 :    
109 : rhijmans 2483
110 :     } else if (stat == 'skew') {
111 :     if (na.rm) {
112 :     n <- colSums(! is.na(x))
113 :     } else {
114 :     n <- nrow(x)
115 :     }
116 : rhijmans 2872 if (asSample) {
117 : rhijmans 3208 sdx <- apply(x, 2, stats::sd, na.rm=na.rm)
118 : rhijmans 2872 } else {
119 :     sdx <- apply(x, 2, function(x) sqrt(sum((x-mean(x, na.rm=na.rm))^2, na.rm=na.rm)/n))
120 :     }
121 :     return( colSums(t(t(x) - colMeans(x, na.rm=na.rm))^3, na.rm=na.rm) / (n * sdx^3) )
122 : rhijmans 2239 }
123 : rhijmans 2881 } # else
124 :    
125 :     return(apply(x, 2, stat, na.rm=na.rm, ...))
126 : rhijmans 1149 }
127 : rhijmans 1193
128 : rhijmans 2221 if (class(stat) != 'character') {
129 :     stop('cannot use this function for large files')
130 : rhijmans 316 }
131 : rhijmans 2221
132 :     st <- NULL
133 :     counts <- FALSE
134 :     if (stat == 'sum') {
135 :     fun <- sum
136 :     st <- 0
137 :     } else if (stat == 'min') {
138 : rhijmans 2881 st <- Inf
139 : rhijmans 2221 } else if (stat == 'max') {
140 : rhijmans 2881 st <- -Inf
141 : rhijmans 2221 } else if (stat == 'range') {
142 :     fun <- range
143 :     } else if (stat == 'countNA') {
144 : rhijmans 2483 warning ("'countNA' is depracted. Use freq(x, 'value=NA') instead")
145 : rhijmans 2221 st <- 0
146 :     counts <- TRUE
147 :     } else if (stat == 'skew') {
148 : rhijmans 2239
149 :     zmean <- cellStats(x, 'mean')
150 :     cnt <- 0
151 : rhijmans 2872 d3 <- 0
152 : rhijmans 2239 sumsq <- 0
153 : rhijmans 2221 counts <- TRUE
154 : rhijmans 2239
155 : rhijmans 2483 } else if (stat == 'mean' | stat == 'sd' | stat == 'rms') {
156 : rhijmans 2221 st <- 0
157 :     sumsq <- 0
158 :     cnt <- 0
159 :     counts <- TRUE
160 : rhijmans 2881
161 : rhijmans 2221 } else {
162 : rhijmans 2907 stop("invalid 'stat'. Should be 'sum', 'min', 'max', 'sd', 'mean', 'rms', or 'skew'")
163 : rhijmans 1193 }
164 : rhijmans 2221
165 : rhijmans 1105
166 : rhijmans 2221 tr <- blockSize(x)
167 : rhijmans 2562 pb <- pbCreate(tr$n, label='cellStats', ...)
168 : rhijmans 1193
169 : rhijmans 2221 for (i in 1:tr$n) {
170 :     d <- getValues(x, row=tr$row[i], nrows=tr$nrows[i])
171 :     if (makeMat) {
172 :     d <- matrix(d, ncol=1)
173 :     }
174 :     if (counts) {
175 : rhijmans 2902 if (na.rm & stat != 'countNA') {
176 : rhijmans 2906 nas <- colSums( is.na(d) )
177 : rhijmans 2221 if (min(nas) == nrow(d)) {
178 :     next
179 :     }
180 :     cells <- nrow(d) - nas
181 : rhijmans 2906 } else {
182 :     if (stat == 'countNA') {
183 :     nas <- colSums( is.na(d) )
184 :     } else {
185 :     cells <- nrow(d)
186 :     }
187 : rhijmans 2221 }
188 :     }
189 :    
190 :     if (stat=='mean') {
191 :     st <- colSums(d, na.rm=na.rm) + st
192 :     cnt <- cnt + cells
193 :    
194 :     } else if (stat=='sum') {
195 :     st <- colSums(d, na.rm=na.rm) + st
196 : rhijmans 1149
197 : rhijmans 2221 } else if (stat == 'sd') {
198 :     st <- colSums(d, na.rm=na.rm) + st
199 :     cnt <- cnt + cells
200 : rhijmans 2239 sumsq <- colSums(d^2, na.rm=na.rm) + sumsq
201 : rhijmans 1149
202 : rhijmans 2221 } else if (stat=='countNA') {
203 :     st <- st + nas
204 :    
205 : rhijmans 2902 } else if (stat=='rms') {
206 :    
207 :     sumsq <- colSums(d^2, na.rm=TRUE) + sumsq
208 :     cnt <- cnt + cells
209 :    
210 : rhijmans 2221 } else if (stat=='skew') {
211 : rhijmans 2872
212 :     d <- t( t(d) - zmean )
213 :     sumsq <- colSums(d^2, na.rm=TRUE) + sumsq
214 :     d3 <- colSums(d^3, na.rm=TRUE) + d3
215 : rhijmans 2239 cnt <- cnt + cells
216 :    
217 : rhijmans 2881 } else if (stat=='min') {
218 :     tmp <- .colMin(d, na.rm=na.rm)
219 :     st <- pmin(st, tmp, na.rm=na.rm)
220 :    
221 :     } else if (stat=='max') {
222 :     tmp <- .colMax(d, na.rm=na.rm)
223 :     st <- pmax(st, tmp, na.rm=na.rm)
224 :    
225 : rhijmans 2221 } else {
226 : rhijmans 2881 # range
227 : rhijmans 2221 st <- apply(rbind(d, st), 2, fun, na.rm=na.rm)
228 :     }
229 : rhijmans 1105
230 : rhijmans 2221 pbStep(pb, i)
231 : rhijmans 1193 }
232 : rhijmans 1105
233 : rhijmans 2221
234 :     if (stat == 'sd') {
235 :     meansq <- (st/cnt)^2
236 : rhijmans 2952 st <- sqrt(( (sumsq / cnt) - meansq ) * (cnt/(cnt-1)))
237 :     if (!asSample) {
238 : rhijmans 3340 #st <- sqrt( st^2 * (cnt / (cnt-1)))
239 :     st <- sqrt( st^2 * ((cnt-1) / cnt))
240 :    
241 :    
242 : rhijmans 2483 }
243 : rhijmans 2221 } else if (stat == 'mean') {
244 :     st <- st / cnt
245 : rhijmans 2483 } else if (stat == 'rms') {
246 :     if (asSample) {
247 :     st <- sqrt(sumsq/(cnt-1))
248 :     } else {
249 :     st <- sqrt(sumsq/cnt)
250 :     }
251 :    
252 : rhijmans 2221 } else if (stat == 'skew') {
253 : rhijmans 2872
254 : rhijmans 2239 if (asSample) {
255 : rhijmans 2872 stsd <- sqrt(sumsq/(cnt-1))^3
256 :     } else {
257 :     stsd <- sqrt(sumsq/cnt)^3
258 : rhijmans 2239 }
259 : rhijmans 2872 st <- d3 / (cnt*stsd)
260 : rhijmans 2907 } else if (stat %in% c('min', 'max')) {
261 :     names(st) <- names(x)
262 : rhijmans 2221 }
263 : rhijmans 1105
264 : rhijmans 2221 pbClose(pb)
265 :     return(st)
266 : rhijmans 316 }
267 : rhijmans 2221 )
268 : rhijmans 427
269 : rhijmans 2239
270 :    
271 :    
272 :    
273 :    
274 :     setMethod('cellStats', signature(x='RasterLayer'),
275 :     function(x, stat='mean', na.rm=TRUE, asSample=TRUE, ...) {
276 :    
277 :     stopifnot(hasValues(x))
278 :     stat <- .csTextFun(stat)
279 :    
280 : rhijmans 2906 if (! inMemory(x) ) {
281 : rhijmans 2239 if (canProcessInMemory(x)) {
282 :     x <- readAll(x)
283 :     }
284 :     }
285 :     if (inMemory(x) ) {
286 :     x <- getValues(x)
287 :    
288 :     if (class(stat) == 'character') {
289 :     if (stat == "mean" ) {
290 :     return( mean(x, na.rm=na.rm) )
291 :     } else if (stat == "sum" ) {
292 : rhijmans 3110 return( sum(as.double(x), na.rm=na.rm ) )
293 : rhijmans 2239 } else if (stat == 'countNA') {
294 :     return( sum(is.na(x)) )
295 :     } else if (stat == "range" ) {
296 :     return( range(x, na.rm=na.rm) )
297 :     } else if (stat == "min" ) {
298 :     return( min(x, na.rm=na.rm) )
299 :     } else if (stat == "max" ) {
300 :     return( max(x, na.rm=na.rm) )
301 :     } else if (stat == "sd" ) {
302 : rhijmans 3208 st <- stats::sd(x, na.rm=na.rm)
303 : rhijmans 2239 if (! asSample) {
304 :     if (na.rm) {
305 : rhijmans 3208 n <- length(stats::na.omit(x))
306 : rhijmans 2239 } else {
307 :     n <- length(x)
308 :     }
309 : rhijmans 3340 #st <- sqrt(st^2 * (n/(n-1)))
310 :     st <- sqrt(st^2 * ((n-1)/n))
311 :    
312 : rhijmans 2239 }
313 :     return(st)
314 : rhijmans 2902 } else if (stat == 'rms') {
315 :     if (na.rm) {
316 :     n <- sum(! is.na(x))
317 :     } else {
318 :     n <- length(x)
319 :     }
320 :     if (asSample) {
321 :     n <- n-1
322 :     }
323 :     # st <- apply(x, 2, function(x) sqrt(sum(x^2)/n))
324 : rhijmans 2907 return( sqrt( sum(x^2)/n ) )
325 : rhijmans 2902
326 :    
327 : rhijmans 2239 } else if (stat == "skew" ) {
328 :     if (na.rm) {
329 : rhijmans 3208 x <- stats::na.omit(x)
330 : rhijmans 2239 }
331 : rhijmans 2872 if (asSample) {
332 : rhijmans 3208 sdx <- stats::sd(x)
333 : rhijmans 2872 } else {
334 :     sdx <- sqrt(sum((x-mean(x))^2)/(length(x)))
335 :     }
336 :     return( sum( (x - mean(x))^3 ) / (length(x) * sdx^3) )
337 : rhijmans 2239 }
338 :     } else {
339 :     return( stat(x, na.rm=na.rm) )
340 :     }
341 :     }
342 :    
343 :    
344 :     if (class(stat) != 'character') {
345 :     stop('cannot use this function for large files')
346 :     }
347 :    
348 :     st <- NULL
349 :     counts <- FALSE
350 :     if (stat == 'sum') {
351 :     fun <- sum
352 :     st <- 0
353 :     } else if (stat == 'min') {
354 :     fun <- min
355 :     } else if (stat == 'max') {
356 :     fun <- max
357 :     } else if (stat == 'range') {
358 :     fun <- range
359 :     } else if (stat == 'countNA') {
360 :     st <- 0
361 :     counts <- TRUE
362 : rhijmans 2902
363 : rhijmans 2239 } else if (stat == 'skew') {
364 :     zmean <- cellStats(x, 'mean')
365 :     cnt <- 0
366 :     sumsq <- 0
367 : rhijmans 2872 d3 <- 0
368 : rhijmans 2239 counts <- TRUE
369 :    
370 : rhijmans 2902 } else if (stat == 'mean' | stat == 'sd' | stat == 'rms') {
371 : rhijmans 2239 st <- 0
372 :     sumsq <- 0
373 :     cnt <- 0
374 :     counts <- TRUE
375 :     } else {
376 :     stop("invalid 'stat'. Should be sum, min, max, sd, mean, or 'countNA'")
377 :     }
378 :    
379 :    
380 :     tr <- blockSize(x)
381 : rhijmans 2562 pb <- pbCreate(tr$n, label='cellStats', ...)
382 : rhijmans 2239
383 :     for (i in 1:tr$n) {
384 :     d <- getValues(x, row=tr$row[i], nrows=tr$nrows[i])
385 :     if (counts) {
386 : rhijmans 2902 if (na.rm & stat != 'countNA') {
387 : rhijmans 2906 nas <- sum(is.na(d) )
388 : rhijmans 2239 if (nas == length(d)) { # only NAs
389 :     next
390 :     }
391 :     cells <- length(d) - nas
392 : rhijmans 2906 } else {
393 :     if (stat == 'countNA') {
394 :     nas <- sum(is.na(d) )
395 :     } else {
396 :     cells <- length(d)
397 :     }
398 : rhijmans 2239 }
399 :     }
400 :    
401 :     if (stat=='mean') {
402 :     st <- sum(d, na.rm=na.rm) + st
403 :     cnt <- cnt + cells
404 :    
405 :     } else if (stat=='sum') {
406 : rhijmans 3110 st <- sum(as.double(d), na.rm=na.rm) + st
407 : rhijmans 2239
408 :     } else if (stat == 'sd') {
409 :     st <- sum(d, na.rm=na.rm) + st
410 :     cnt <- cnt + cells
411 :     sumsq <- sum( d^2 , na.rm=na.rm) + sumsq
412 :    
413 :     } else if (stat=='countNA') {
414 :     st <- st + nas
415 :    
416 :     } else if (stat=='skew') {
417 :    
418 : rhijmans 2872 d <- (d - zmean)
419 : rhijmans 2902 sumsq <- sum(d^2, na.rm=na.rm) + sumsq
420 :     d3 <- sum(d^3, na.rm=na.rm) + d3
421 : rhijmans 2239 cnt <- cnt + cells
422 : rhijmans 2902
423 :     } else if (stat=='rms') {
424 :     sumsq <- sum( d^2, na.rm=na.rm) + sumsq
425 :     cnt <- cnt + cells
426 : rhijmans 2239
427 : rhijmans 2902 } else {
428 :     st <- fun(d, st, na.rm=na.rm)
429 : rhijmans 2239 }
430 :    
431 :     pbStep(pb, i)
432 :     }
433 :     pbClose(pb)
434 :    
435 :     if (stat == 'sd') {
436 :     meansq <- (st/cnt)^2
437 : rhijmans 2952 st <- sqrt(( (sumsq / cnt) - meansq ) * (cnt/(cnt-1)))
438 :     if (!asSample) {
439 : rhijmans 3340 #st <- sqrt( st^2 * (cnt / (cnt-1)))
440 :     st <- sqrt( st^2 * ((cnt-1) / cnt))
441 :    
442 : rhijmans 2239 }
443 :     } else if (stat == 'mean') {
444 :     st <- st / cnt
445 : rhijmans 2902
446 :     } else if (stat == 'rms') {
447 :     if (asSample) {
448 :     st <- sqrt(sumsq/(cnt-1))
449 :     } else {
450 :     st <- sqrt(sumsq/cnt)
451 :     }
452 :    
453 : rhijmans 2239 } else if (stat == 'skew') {
454 :     if (asSample) {
455 : rhijmans 2872 stsd <- sqrt(sumsq/(cnt-1))^3
456 :     } else {
457 :     stsd <- sqrt(sumsq/cnt)^3
458 : rhijmans 2239 }
459 : rhijmans 2872 st <- d3 / (cnt*stsd)
460 :     }
461 : rhijmans 2239 return(st)
462 :     }
463 :     )
464 :    

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge