[R-sig-Geo] Local Moran'I for raster & focalFilter
Robert J. Hijmans
r.hijmans at gmail.com
Thu Apr 14 21:40:28 CEST 2011
> 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".
>
Nick,
I believe that Wii is zero in these functions because it uses:
focal(z, fun='sum', na.rm=TRUE) - z)
The above first sums all focal values but then it subtracts the
central cell value. This alleviates us from the problem of finding the
central value, which is difficult around the borders of the raster
(for which the work-around is to use additional rows/columns with NA
values).
I look forward to seeing your implementation of the G statistics :)
Robert
Robert
On Thu, Apr 14, 2011 at 12:27 PM, Nick Hamm <nick at hamm.org> wrote:
> 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).
>
> 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
>>
>>
>
More information about the R-sig-Geo
mailing list