[R-sig-Geo] Spatial overlay vs. extract and aggregation methods

[Juta Kawalerowicz]

Hi,

I have a spatial polygon data frame (spdf1) and try to aggregate
information from it over polygons in a larger spatial polygon data
frame(spdf2). I wanted to run some simple test to see that my code is sound
(on nc.sids from maptools), but run into a problem which I don't
understand. As far as I can see, there are 2 ways to do this:

1) rasterise the smaller spdf1, extract and calculate over larger spdf2

library(maptools)
library(raster)
nc.sids <- readShapeSpatial(system.file("shapes/sids.shp",
package="maptools")[1],
                            IDvar="FIPSNO", proj4string=CRS("+proj=longlat
+ellps=clrk66"))

#plot(nc.sids)
r<-raster(ncol=180, nrow=180)
extent(r)<-extent(nc.sids)
rp<-rasterize(nc.sids, r, 'BIR74')
plot(rp)
plot(nc.sids, add=TRUE)

#this will take time on larger files...
v <- extract(rp, nc.sids, weights=TRUE)
res<-sapply(v, function(x) if (!is.null(x)) {sum(apply(x, 1, prod),
na.rm=TRUE) / sum(x[,2])} else NA )
cor(res, nc.sids$BIR74)
#seems to be working ok.

2) use over function from sp package

overlay<-over(nc.sids, nc.sids)
#this should return the same indices of nc.sids which fall inside nc.sids -
basically should be falling within itself for each polygon, right? But this
is not what I get.

Any hints would be greatly appreciated.

Cheers,
Juta

	[[alternative HTML version deleted]]