Hello,
I am attempting to calculate the mean value in a raster (3432 x 8640) for
each of 689 polygons (code at end of email). To that end, I have a few
questions:
First, the for loop runs fine for a while (2 hours, maybe) before
terminating itself and returning this error:
Error in unlist(lapply(rings, function(x, pts) pointsInPolygon(pts, x, :
negative length vectors are not allowed
Based on the output matrix, it seems to have gotten hung up on polygon 24.
Any ideas why this might be and how I can fix it? In reading the shapefile I
both delete null objects and force rings.
Second, this process (if it works) will take forever. I've tried an approach
where I covert the polygons to a raster, but this seems to result in data
loss and produce a lot of null values. I've also read on this list a
discussion about just using those points in the bounding box. However, it
was entirely unclear to me how to implement that recommendation.
Any help would be greatly appreciated.
Cheers,
Evann
R Code:
# Projection
prj <- CRS("+proj=longlat +ellps=WGS84")
# Shapefile
polys <- readShapePoly("GeoEPR/geoepr-20100212-1248.shp", proj4string=prj,
verbose=TRUE, delete_null_obj=TRUE, force_ring=TRUE)
summary(polys)
# Raster
popimg <- readGDAL("PopDens/glds05ag.bil")
proj4string(popimg) <- prj
summary(popimg)
# Values
y <- slot(popimg,"data")
head(y)
# Total polygons
n <- length(slot(polys,"polygons"))
n
# Empty matrix
output1 <- as.data.frame(matrix(NA, nrow=n, ncol=4))
names(output1) <- c("cowgroup", "start", "end", "popdens")
head(output1)
# Beastly for loop
for(i in 1:n) {
output1[i,1] <- slot(polys, "data")[i,1]
output1[i,2] <- slot(polys, "data")[i,2]
output1[i,3] <- slot(polys, "data")[i,3]
overlap <- overlay(popimg,polys[i,])
output1[i,4] <- mean(y[overlap==1 & !is.na(overlap),1], na.rm=TRUE)
}
============================================
Evann Smith
Doctoral Student
Department of Government, Harvard University
============================================
[[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
