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

Robert J. Hijmans r.hijmans at gmail.com
Thu Apr 14 19:19:48 CEST 2011

```Hi Babak,

I would propose the below for the local Moran with a RasterLayer
because it deals with NA values in a memory safe manner (and will be a
bit slower because of that)

library(raster)

MoranLocal <- function(x) {
n <- ncell(x) - cellStats(x, 'countNA')
s2 <-  cellStats(x, sd)^2
# adjust variance denominator from n-1 to n
s2 <- (s2 * (n-1)) / n
z  <- x - cellStats(x, mean)
#weights
w  <- focal( z, fun=function(x, ...){ sum(!is.na(x))-1 } )
lz <- (focal(z, fun='sum', na.rm=TRUE) - z) / w
(z / s2) * lz
}

I looked up the Getis and Ord paper you mention (I happen to have the
book), and they do indeed subtract the focal cell from the global mean
& variance. That seems intuitive. But I do not think it was is
suggested in Anselin's LISA paper (Formula 12)  (and implemented in
spdep).

Also, it should be true that  globalMoran == mean(localMoran)

I think the below produces the global Moran (for the queen case
neighboring grid cells)

Moran <- function(x) {
z <- x - cellStats(x, mean)
wZiZj <- focal(z, fun='sum', na.rm=TRUE)
wZiZj <- overlay(wZiZj, z, fun=function(x,y){ (x-y) * y })
wZiZj <- cellStats(wZiZj, sum)
z2 <- cellStats(z*z, sum)
n <- ncell(z) - cellStats(z, 'countNA')
# weights
w <- focal( z, fun=function(x, ...){ max(0, sum(!is.na(x))-1) } )
NS0 <- n / cellStats(w, sum)
mI <- NS0 * wZiZj / z2
return(mI)
}

With your 'tmp' data, I get:

> Moran(tmp)
[1] 0.8526906

> cellStats(MoranLocal(tmp), mean)
[1] 0.8192615

Close, but not quite the same.

I'll put these functions in the raster package.

Robert

On Thu, Apr 14, 2011 at 7:34 AM, Babak Naimi <naimi at itc.nl> wrote:
> Dear Robert and Roger,
>
> Thank both of you. After all, I came up with the following function (focalMoran) that efficiently works. I tested it for a raster layer with 122345 not-NA cells, and compared with locmor (based on localmoran in spdep package). The results are almost identical (cor = 0.99987). I think the only difference between two function is that in the focalMoran, the center pixel is omitted in calculation of global mean and variance. However, focalMoran was so quicker than locmor. Following you will find the difference in computation as well as the functions.
>
> Many thanks,
> Babak
>
>
> b1<-raster(----) # a raster with dimension 355, 617, 219035  (nrow, ncol, ncell) and resolution of 843.8446 m
> filter<-matrix(1,nc=9,nr=9)
>
> system.time(
> b.foc<-focalMoran(b1,filter)
> )
>   user  system elapsed
>  44.63    0.45   48.91
>
> system.time(
> b.spdep<-locmor(b1,d1=0,d2=4863.92)
> )
>   user  system elapsed
> 5738.56    6.02 5832.11
>
> ###########################
>
> focalMoran<-function(ras,filter) {
>        globsum <- cellStats(ras, sum)
>        glob.sqsum <- cellStats(ras^2, sum)
>        n<-length(na.omit(ras[]))
>        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 <- mean(z[-i], na.rm=na.rm)
>                return((z[i]/s2) * lz)
>                }
>        ras2 <- raster::expand(ras, extent(ras)+dim(filter)[1]*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
>        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)
> }
>
>
>
>
> -----Original Message-----
> From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
> Sent: Wednesday, April 13, 2011 11:15 PM
> To: Robert J. Hijmans
> Cc: Babak Naimi; r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] Local Moran'I for raster & focalFilter
>
> On Wed, 13 Apr 2011, Robert J. Hijmans wrote:
>
>> Babak,
>>
>> I looked at the original paper,
>> http://onlinelibrary.wiley.com/doi/10.1111/j.1538-4632.1995.tb00338.x/abstract,
>> and the spdep manual, and other places, but they all seemed to be
>> saying somewhat different on the local moran statistic. I followed the
>> logic of the spdep::localmoran function, assuming that this is
>> "correct" and I should have phrased my answer more carefully: "based
>> onthe spdep::localmoran function, mean and variance should be
>> global...".  I do note that the two results are different but that
>> their correlation is very high R2 = 0.97; so perhaps they are both
>> right. Perhaps Roger can weigh in on this.
>
> I don't have access to my library now. My understanding is that you are
> short-circuiting the construction of the spatial weights (or the n
> separate star matrices, one for each i, in Michael Tiefeldorf's terms). In
> general, Moran's I is not often used with non-zero values on the diagonal
> of the weights matrix. Non-zero values are seen in the local Getis-Ord
> Gi_star, but not in local Gi. My guess (without references) would be that
> the weight for i itself should be zero (the centre cell of the raster).
> The difference between mean() and sum() is the difference between
> row-standardised (summing to 1) and binary (all unity) weights. The
> function as it is doesn't permit other weighting schemes (such as IDW in
> either row-standardised or binary form).
>
> Hope this helps,
>
> Roger
>
>>
>> To make it run, in your code I replaced this:
>> filter<-c(NA,1,NA,1,1,1,NA,1,NA)
>> filter<-matrix(f,nc=3,nr=3,byrow=T)
>>
>> with:
>> filter<-matrix(1 ,nc=3,nr=3)
>>
>> Robert
>>
>>
>> On Wed, Apr 13, 2011 at 10:46 AM, Babak Naimi <naimi at itc.nl> wrote:
>>> 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)
>>>
>>>
>>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of