Moran's I : list of neighbour and generating residual maps.

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20090513/77898ed2/attachment.pl>
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20090514/7f63e2a9/attachment.pl>

Dear all,

Sorry for I forward this message to the list, but I got one error during
yesterday when send it, and I am not sorry if it worked fine.
Yes, the line with the overlay() got corrupted, and still is.

Try:

require(rgdal)
require(adehabitat)
(file1 <-  paste(system.file(package = "adehabitat"),
                "ascfiles/elevation.asc", sep = "/"))
el <- readGDAL(file1)
image(el, axes=T)
el.pix <- as(el, "SpatialPixelsDataFrame")
el.pol <- as.SpatialPolygons.SpatialPixels(el.pix)
rn <- sapply(slot(el.pol, "polygons"), function(x) slot(x, "ID"))
df <- data.frame(rn=rn, row.names=rn)
el.pol.spdf<-SpatialPolygonsDataFrame(el.pol, data=df)
el.pol.spdf$pixvalue <- el.pix$band1
spplot(el.pix, col.regions=heat.colors(20))
x11()
spplot(el.pol.spdf, "pixvalue", col.regions=heat.colors(20))

moves the values without overlaying - since the SpatialPixels and 
SpatialPolygons objects match one-to-one, this is OK.
From here you need a neighbour object, for example:
gsz <- gridparameters(el)[,2]
d <- ceiling(sqrt(sum(gsz^2)))
library(spdep)
nb_queen <- dnearneigh(coordinates(el.pol.spdf), 0, d)
nb_queen

before calculating Moran's I:

moran.test(el.pol.spdf$pixvalue, nb2listw(nb_queen, style="C"))

There will typically be a lot of autocorrelation if the raster cell 
resolution is much finer than the "size" of the entities.

Hope this helps,

Roger
Bests

milton

---------- Forwarded message ----------
From: milton ruser <milton.ruser at gmail.com>
Date: Wed, May 13, 2009 at 6:55 PM
Subject: Moran's I : list of neighbour and generating residual maps.
To: r-sig-geo at stat.math.ethz.ch

Dear GeoR Gurus,

I have some raster maps, with continuous values, and I would like
to (1) calculate Local Moran's I ;  (2) generate a residual map considering
different lags. I am trying to follow the Bivand and colleague's book,
but I don't know how to generate the neighbour list for each lag.

On the page 268 of book one can find a example of how to
plot the moran I index like moran.plot(NY8$Cases, listw=nb2listw (NY_b,
style="C")).

On the code I read a raster map, vectorize it (each pixel is one polygon),
and
the pixel value are stored as attribute for the polygon.

So, if somebody can help to solve the TWO questions listed above
I would be very thanks.

Cheers

milton
====
require(rgdal)
require(adehabitat)
(file1 <-  paste(system.file(package = "adehabitat"),
              "ascfiles/elevation.asc", sep = "/"))
el <- readGDAL(file1)
image(el, axes=T)
el.pix <- as(el, "SpatialPixelsDataFrame")
el.pol <- as.SpatialPolygons.SpatialPixels(el.pix)
str(el.pol, max.level=2)
el.pol at polygons[[1]]
plot(el.pol)

rn <- sapply(slot(el.pol, "polygons"), function(x) slot(x, "ID"))
df <- data.frame(rn=rn, row.names=rn)

el.pol.spdf<-SpatialPolygonsDataFrame(el.pol, data=df)
el.pol.spdf at data$pixvalue<-overlay(el<el.pol.spdf at data$pixvalue%3C-overlay(el>,
el.pix)@data[,1]
image(el, axes=T)
plot(el.pol, add=T)
head(el.pol.spdf at data)

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
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