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

[Robert Hijmans]

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



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Local-Moran-I-for-raster-focalFilter-tp6250323p6266969.html
Sent from the R-sig-geo mailing list archive at Nabble.com.