Dear all
Replies below
On 14 April 2011 19:19, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> 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).
According to my reading of Anselin (1995) wii should be excluded. See
page 98 in the text between equations 8 and 9 "...,and by convention
wii = 0".
Roger's post yesterday also suggested it should be excluded. But what
is the situation with 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.
>
This is a good idea. There is also interest in using other local
statistics (e.g., Getis) in remote sensing. This is what I was on my
way to calculating when I found the problem with focalFilter. Raster
is MUCH faster than using spdep for images - because it is not
necessary to recompute the neighbours for every pixel - as Babak has
shown.
> Robert
>
>