[R-sig-Geo] Fwd: Moran's I : list of neighbour and generating residual maps.

[Roger Bivand]

On Thu, 14 May 2009, milton ruser wrote:

> 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