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

[Babak Naimi]

Dear All,

I'm trying to calculate Local Moran's I for each grid cell in a raster given a distance (d1,d2) within which the neighbour cells are identified. I wrote a function which uses localmoran of spdep package and returns a raster of Moran's I measures. But in case of large dataset, it is a time consuming function. 

I am wondering if there is any efficient solution to calculating local moran's I for raster data. Maybe, using focalFilter is a solution. Therefore, 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 Naimi
PhD Candidate, Natural Resources Department,
Faculty of Geo-Information Science and Earth Observations (ITC), University of Twente, 
P.O. Box 217, 7500 AE Enschede, The Netherlands
Tel: +31 53 4874212


########

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) 
}


locMor.focal<-function(mtx) { # function to calculate local Moran using focalFilter
  require(raster)
  mtx.v<-as.vector(mtx)
  mtx.v.c<-mtx.v[c(trunc(length(mtx.v)/2)+1)] #value of the center of matrix (Xi)
  mtx.v.Wc<-mtx.v[-c(trunc(length(mtx.v)/2)+1)] # values of the matrix excluding center value (vector of Xj)
  mtx.v<-mtx.v[which(mtx.v != 0)]
  mtx.v.Wc<-mtx.v.Wc[which(mtx.v.Wc != 0)]
  
  z<-function(v) { #to calculate deviation from the mean for a value in the matrix
    return(v-mean(mtx.v,na.rm=T))
    }
  z.sum<-sum(sapply(mtx.v.Wc,z))
  mI<- (z(mtx.v.c)/var(mtx.v,na.rm=T))*z.sum
  return(mI)} 

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

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)

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

plot(tmp.spdep)
plot(tmp.focal)