1# Author: Robert J. Hijmans 2# Date : April 2011 3# Version 1.0 4# Licence GPL v3 5 6 7 8.getFilter <- function(w, warn=TRUE) { 9 if (!is.matrix(w)) { 10 w <- .checkngb(w) 11 w <- matrix(1, nrow=w[1], ncol=(w[2])) 12 w[ceiling(dim(w)[1]/2), ceiling(dim(w)[2]/2)] <- 0 13 } else { 14 if (w[ceiling(dim(w)[1]/2), ceiling(dim(w)[2]/2)] != 0) { 15 if (warn) { 16 warning('central cell of weights matrix (filter) was set to zero') 17 } 18 w[ceiling(dim(w)[1]/2), ceiling(dim(w)[2]/2)] <- 0 19 } 20 stopifnot(all(w >= 0)) 21 } 22 if (min(dim(w) %% 2)==0) { 23 stop('dimensions of weights matrix (filter) must be uneven') 24 } 25 w 26} 27 28 29 30Geary <- function(x, w=matrix(c(1,1,1,1,0,1,1,1,1),3,3)) { 31 32 w <- .getFilter(w, warn=FALSE) 33 34 i <- trunc(length(w)/2)+1 35 n <- ncell(x) - cellStats(x, 'countNA') 36 37 fun <- function(x,...) sum(w*(x-x[i])^2, ...) 38 w2 <- w 39 w2[] <- 1 40 Eij <- cellStats(focal(x, w=w2, fun=fun, na.rm=TRUE, pad=TRUE), sum) 41 42 if (sum(! unique(w) %in% 0:1) > 0) { 43 xx <- calc(x, fun=function(x) ifelse(is.na(x), NA ,1)) 44 W <- focal(xx, w=w, na.rm=TRUE, pad=TRUE ) 45 } else { 46 w[w==0] <- NA 47 W <- focal(x, w=w, fun=function(x, ...){ sum(!is.na(x)) }, pad=TRUE ) 48 } 49 z <- 2 * cellStats(W, sum) * cellStats((x - cellStats(x, mean))^2, sum) 50 51 (n-1)*Eij/z 52} 53 54 55 56 57GearyLocal <- function(x, w=matrix(c(1,1,1,1,0,1,1,1,1),3,3)) { 58 w <- .getFilter(w) 59 i <- trunc(length(w)/2)+1 60 fun <- function(x,...) sum(w*(x-x[i])^2, ...) 61 w2 <- w 62 w2[] <- 1 63 Eij <- focal(x, w=w2, fun=fun, na.rm=TRUE, pad=TRUE) 64 65 s2 <- cellStats(x, 'sd')^2 66 if (ncell(x) < 1000000) { n <- ncell(x) - cellStats(x, 'countNA' ) 67 } else { n <- ncell(x) } 68 69 #s2 <- (s2 * (n-1)) / n 70 Eij / s2 71} 72 73