1# Author: Robert J. Hijmans 2# Date : April 2011 3# Version 1.0 4# Licence GPL v3 5 6 7..moran <- function(x, directions=8) { 8 stopifnot(directions %in% c(4,8)) 9 # not memory safe 10 adj <- adjacent(x, 1:ncell(x), target=1:ncell(x), directions=8, pairs=TRUE) 11 z <- x - cellStats(x, mean) 12 wZiZj <- stats::na.omit(z[adj[,1]] * z[adj[,2]]) 13 z2 <- cellStats(z*z, sum) 14 NS0 <- (ncell(z)-cellStats(z, 'countNA')) / length(wZiZj) 15 mI <- NS0 * sum(wZiZj) / z2 16 return(mI) 17} 18 19 20 21Moran <- function(x, w=matrix(c(1,1,1,1,0,1,1,1,1),3,3) ) { 22 23 z <- x - cellStats(x, mean) 24 wZiZj <- focal(z, w=w, fun='sum', na.rm=TRUE, pad=TRUE) 25 wZiZj <- overlay(wZiZj, z, fun=function(x,y){ x * y }) 26 wZiZj <- cellStats(wZiZj, sum) 27 z2 <- cellStats(z*z, sum) 28 n <- ncell(z) - cellStats(z, 'countNA') 29 # weights 30 if (sum(! unique(w) %in% 0:1) > 0) { 31 zz <- calc(z, fun=function(x) ifelse(is.na(x), NA ,1)) 32 W <- focal( zz, w=w, fun='sum', na.rm=TRUE, pad=TRUE) 33 } else { 34 w2 <- w 35 w2[w2==0] <- NA 36 W <- focal( z, w=w2, fun=function(x, ...){ as.double(sum(!is.na(x))) }, pad=TRUE) 37 } 38 NS0 <- n / cellStats(W, sum) 39 mI <- NS0 * wZiZj / z2 40 return(mI) 41} 42 43 44MoranLocal <- function(x, w=matrix(c(1,1,1,1,0,1,1,1,1),3,3)) { 45 46 z <- x - cellStats(x, mean) 47 #weights 48 #w <- .getFilter(w) 49 if (sum(! unique(w) %in% 0:1) > 0) { 50 zz <- calc(z, fun=function(x) ifelse(is.na(x), NA ,1)) 51 W <- focal( zz, w=w, na.rm=TRUE, pad=TRUE) 52 } else { 53 w2 <- w 54 w2[w2==0] <- NA 55 W <- focal( z, w=w2, fun=function(x, ...){ sum(!is.na(x)) }, na.rm=TRUE, pad=TRUE) 56 } 57 lz <- focal(z, w=w, na.rm=TRUE, pad=TRUE) / W 58 59 n <- ncell(x) - cellStats(x, 'countNA') 60 s2 <- cellStats(x, 'sd')^2 61 # adjust variance denominator from n-1 to n 62 #s2 <- (s2 * (n-1)) / n 63 64 (z / s2) * lz 65} 66 67 68