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.ru...@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@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....@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.s...@data$pixvalue<-overlay(el<el.pol.s...@data$pixvalue%3c-overlay(el>,
el.pix)@data[,1]
image(el, axes=T)
plot(el.pol, add=T)
head(el.pol.s...@data)
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
