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