[R-sig-Geo] Local Moran'I for raster & focalFilter

Wed Apr 13 19:46:42 CEST 2011

Dear Robert,

Thanks for your reply. Now, the result is the same as with localmoran in spdep package. But I checked this with Getis & Ord (1996), seems the formula, as you implement in the script, is not the same as what has been explained in the Getis & Ord. First difference is that, Xi is excluded in calculation of the global mean and variance. The second difference is that, in the line "lz <- mean(z[-i], na.rm=na.rm)", 'mean' should be replaced with 'sum'. I re-implemented this function. The example we used here came from that reference (p.266). The calculation for the cell 136 has also been provided in the book. The result that this new function returns is the same as with the book, now. But the strange thing to me is that it is not the same as spdep. Any comment?

Babak

Getis,A. and Ord, J. K. 1996 Local spatial statistics: an overview. In P. Longley and M. Batty (eds) Spatial
analysis: modelling in a GIS environment (Cambridge: Geoinformation International), 261-277.

#----------------------------------------------
focalMoran<-function(ras,filter) {
	globsum <- cellStats(ras, sum)
	glob.sqsum <- cellStats(ras^2, sum)
	n<-ncell(tmp)
	lomo <- function (x, na.rm=TRUE, sigmaX=globsum, sigmaX2=glob.sqsum,N=n) {
		i <- trunc(length(x)/2)+1
	     	xbar=(sigmaX-x[i])/(N-1)
		s2<-((sigmaX2-x[i]^2)-((sigmaX -x[i])^2/(N-1)))/(N-1)
		z <- x - xbar
		lz <- sum(z[-i], na.rm=na.rm)
	    	return((z[i]/s2) * lz)
		}
	ras2 <- raster::expand(ras, extent(ras)+2*res(ras))
	ras2 <- focalFilter(ras2,filter,lomo) 
	ras2 <- crop(ras2, ras)
	return(ras2) 
}
#--------
locmor<-function(ras,d1,d2) { # function to calculate local Moran using spdep 
        require(raster) 
        require(spdep) 
        df<-data.frame(cbind(cells=1:ncell(ras),xyFromCell(ras,1:ncell(ras)),values=ras[1:ncell(ras)])) 
        df<-na.omit(df)
        df.nb <- dnearneigh(as.matrix(df[,2:3]), d1, d2) 
        df.listw <- nb2listw(df.nb) #turns neighbourhood object into a weighted list 
        lmor <- localmoran(df$values, df.listw) 
        ras[df$cells]<-as.data.frame(lmor)[[1]] 
        return(ras) 
} 

#####--------------- 

tmp <- c(55, 56, 54, 58, 65, 75, 82, 77, 74, 74, 69, 61, 62, 71, 73, 63, 62,
63, 64, 59, 85, 88, 95, 106, 110, 99, 89, 82, 79, 84, 97, 79, 55, 55, 56,
60, 91, 95, 86, 98, 115, 105, 110, 107, 101, 89, 85, 68, 55, 54, 53, 61, 82,
102, 88, 93, 96, 94, 110, 114, 109, 103, 92, 68, 59, 58, 60, 64, 88, 99, 82,
81, 71, 80, 89, 89, 89, 102, 104, 75, 63, 57, 58, 58, 77, 92, 82, 71, 59,
90, 105, 92, 79, 98, 110, 83, 62, 55, 56, 56, 80, 90, 99, 88, 64, 91, 112,
94, 76, 91, 100, 85, 62, 59, 55, 61, 99, 97, 93, 80, 65, 87, 107, 80, 59,
60, 67, 66, 65, 62, 68, 72, 102, 94, 90, 83, 74, 81, 96, 69, 52, 50, 51, 54,
62, 62, 86, 85, 55, 59, 64, 72, 75, 70, 70, 62, 66, 61, 55, 57, 52, 59, 61,
56, 41, 40, 43, 44, 46, 48, 50, 52, 68, 69, 60, 61, 42, 43, 44, 43, 42, 41,
42, 44, 43, 43, 44, 47, 58, 59, 55, 61, 44, 41, 39, 42, 44, 43, 42, 42, 42,
43, 43, 49, 56, 53, 53, 61, 43, 42, 40, 42, 42, 42, 41, 42, 42, 43, 42, 53,
66, 61, 51, 62, 40, 42, 41, 40, 43, 49, 46, 42, 42, 43, 43, 49, 59, 62, 53,
62, 40, 41, 42, 43, 49, 54, 47, 44, 42, 44, 43, 46, 52, 56, 56, 61) 
tmp <- matrix(tmp, ncol=16, nrow=16, byrow=T) 
tmp <- raster(tmp) 

filter<-c(NA,1,NA,1,1,1,NA,1,NA)
filter<-matrix(f,nc=3,nr=3,byrow=T)

tmp.focal<-focalMoran(tmp,filter) 
tmp.spdep<-locmor(tmp,d1=0,d2=0.0625)

tmp.focal[136] # the calculation for this cell has been provided in Getis and Ord(1996, p.268)
tmp.spdep[136]

par(mfrow=c(1,2))
plot(tmp.focal)	
plot(tmp.spdep)

#------------------------------------------
#-----------------------------------------

Date: Tue, 12 Apr 2011 15:14:05 -0700 (PDT)
From: Robert Hijmans <r.hijmans at gmail.com>
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Local Moran'I for raster & focalFilter
Message-ID: <1302646445119-6266969.post at n2.nabble.com>
Content-Type: text/plain; charset=us-ascii

> Dear All, 
> 
> I tried to write a function that can be used in focalFilter of raster
> package, but the output
> result is different from what was calculated using spdep for a sample
> dataset. 
> It might be my mistake in the formula that is used in the function.
> Following you will
> find both functions as well as the resultss for a sample raster. 
>
> I would appreciate if anyone can help. 
> 
> Best regards, 
> Babak 

Babak, 
I think the mistake you made is that the local moran uses the _global_ mean
and the _global_ variance, whereas you used the _local_ mean and the _local_
variance. I have changed your function and now the result with focalFilter
is the same as with spdep.
Robert

locmor<-function(ras,d1=0,d2=6000) { # function to calculate local Moran
using spdep 
        require(raster) 
        require(spdep) 

df<-data.frame(cbind(cells=1:ncell(ras),xyFromCell(ras,1:ncell(ras)),values=ras[1:ncell(ras)])) 
        df.sub<-subset(df,!is.na(values)) 
        df.nb <- dnearneigh(as.matrix(df.sub[,2:3]), d1, d2) 
        df.listw <- nb2listw(df.nb) #turns neighbourhood object into a
weighted list 
        lmor <- localmoran(df.sub$values, df.listw) 
        ras[df.sub$cells]<-as.data.frame(lmor)[[1]] 
        return(ras) 
} 

#--------------- 

tmp <- c(55, 56, 54, 58, 65, 75, 82, 77, 74, 74, 69, 61, 62, 71, 73, 63, 62,
63, 64, 59, 85, 88, 95, 106, 110, 99, 89, 82, 79, 84, 97, 79, 55, 55, 56,
60, 91, 95, 86, 98, 115, 105, 110, 107, 101, 89, 85, 68, 55, 54, 53, 61, 82,
102, 88, 93, 96, 94, 110, 114, 109, 103, 92, 68, 59, 58, 60, 64, 88, 99, 82,
81, 71, 80, 89, 89, 89, 102, 104, 75, 63, 57, 58, 58, 77, 92, 82, 71, 59,
90, 105, 92, 79, 98, 110, 83, 62, 55, 56, 56, 80, 90, 99, 88, 64, 91, 112,
94, 76, 91, 100, 85, 62, 59, 55, 61, 99, 97, 93, 80, 65, 87, 107, 80, 59,
60, 67, 66, 65, 62, 68, 72, 102, 94, 90, 83, 74, 81, 96, 69, 52, 50, 51, 54,
62, 62, 86, 85, 55, 59, 64, 72, 75, 70, 70, 62, 66, 61, 55, 57, 52, 59, 61,
56, 41, 40, 43, 44, 46, 48, 50, 52, 68, 69, 60, 61, 42, 43, 44, 43, 42, 41,
42, 44, 43, 43, 44, 47, 58, 59, 55, 61, 44, 41, 39, 42, 44, 43, 42, 42, 42,
43, 43, 49, 56, 53, 53, 61, 43, 42, 40, 42, 42, 42, 41, 42, 42, 43, 42, 53,
66, 61, 51, 62, 40, 42, 41, 40, 43, 49, 46, 42, 42, 43, 43, 49, 59, 62, 53,
62, 40, 41, 42, 43, 49, 54, 47, 44, 42, 44, 43, 46, 52, 56, 56, 61) 
tmp <- matrix(tmp, ncol=16, nrow=16, byrow=T) 
tmp <- raster(tmp) 

tmp.spdep<-locmor(tmp,d1=0,d2=0.08838835) 

# new function for use with focalFilter
# first compute global mean and global variance
globmean <- cellStats(tmp, mean)
globvar <-  cellStats(tmp, sd)^2
# adjust variance denominator from n-1 to n
globvar <- (globvar * (ncell(tmp)-1)) / ncell(tmp)

lomo <- function (x, na.rm=TRUE, xbar=globmean, s2=globvar) {
	i <- trunc(length(x)/2)+1
    z <- x - globmean
	lz <- mean(z[-i], na.rm=na.rm)
    (z[i]/s2) * lz
}

filter<-matrix(1, nrow=3, ncol=3) 
tmp.focal<-focalFilter(tmp,filter,lomo) 

par(mfrow=c(1,2))
plot(tmp.spdep) 
plot(tmp.focal) 

# note that you could, but should not, do this:
tmp.focal<-focal(tmp,fun=lomo) 
# because  the trick to get the focal cell
#     i <- trunc(length(x)/2)+1
# does not work for the border cells. 
# To get the border cells:

tmp2 <- raster::expand(tmp, extent(tmp)+2*res(tmp))
tmpfoc2 <- focalFilter(tmp2,filter,lomo) 
tmpfoc2 <- crop(tmpfoc2, tmp)

par(mfrow=c(1,2))
plot(tmp.spdep) 
plot(tmpfoc2) 