[R-sig-Geo] Local Moran'I for raster & focalFilter
Babak Naimi
naimi at itc.nl
Thu Apr 14 16:34:59 CEST 2011
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
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list