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