I suggested by Tomislav Hengl, with very large datasets this might be 
better done in a dedicated GIS program. SAGA GIS is one good option 
which you can use from within R. Another, certainly when you work in 
Linux, is to do this in GRASS GIS using spgrass6. You would need to 
carry out the following steps (check each of the help function of each 
of the commands for additional options).

1) Open GRASS GIS
2) Run R on the command line within GRASS
3) library(spgrass6)
4) writeRAST6(your_raster_layer, "name_raster_in_GRASS") # writes your 
raster to the current GRASS mapset
5) writeVECT6(your_vector_layer, "name_of_vector_in_GRASS") # writes 
your vector to the current GRASS mapset
6) execGRASS(g.region, parameters=list(rast="name_raster_in_GRASS") # 
sets the region (resolution, extend) to that of your raster. Note that 
this is equivalent to the GRASS command g.region 
rast="name_raster_in_GRASS".
7) execGRASS(v.rast.stats, 
parameters=list(vector="name_of_vector_in_GRASS", 
raster="name_raster_in_GRASS", colprefix="prefix_column_name")) # This 
calculates zonal statistics, including the mean, and writes them to 
columns in your attribute table. The name of the column names are 
prefix_statistics
8) Last step, use readVECT6 to get your vectormap back in R.

Paulo


 > 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 at stat.math.ethz.ch
 > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
 >



	[[alternative HTML version deleted]]

