[R-sig-Geo] Local Moran'I for raster & focalFilter
Robert J. Hijmans
r.hijmans at gmail.com
Wed Apr 13 21:34:09 CEST 2011
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.
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)
>
>
>
More information about the R-sig-Geo
mailing list