¿Cómo funciona el polígono espacial% sobre el polígono% cuando se agregan valores en r?

Estoy trabajando en un proyecto de epidemiología ambiental donde tengo exposiciones puntuales (~ 2,000 operaciones de cerdos industriales - OHI). Estas OHI rocían en los campos cercanos, pero las heces y el olor de las heces pueden viajar millas. Entonces, estas exposiciones puntuales obtienen amortiguadores de 3mi, y quiero saber la cantidad de exposiciones de la OHI (de varios tipos: suma de la cantidad de estiércol, cantidad de cerdos, lo que sea; lo más simple, solo la cantidad de amortiguadores de exposición superpuestos) por bloques de censo de Carolina del Norte (~ 200,000). Los bloques censales de exclusión (azul) son (1) cualquier cosa en las 5 ciudades más pobladas y (2) condados que no bordean un condado con una OHI (nota: eso se hizo con la función gRelate y los códigos DE-9IM - muy hábil!). Vea la imagen a continuación para un visual

El último paso es agregar la representación de exposición amortiguada a cada bloque del censo. Aquí es donde estoy perplejo.

He tenido buenos momentos con las funciones% over% en el paquete sp hasta ahora, pero entiendo por la viñeta over que poly-poly y poly-line over se implementan en rgeos. La viñeta solo cubre line-poly y poly de autorreferencia, y no con agregación, por lo que estoy un poco confundido sobre cuáles son mis opciones para poly-poly con agregación de funciones, como sum o mean.

Para un caso de prueba, considere el siguiente fragmento, algo detallado, que funciona con el archivo de fronteras de países del mundo. Esto debería poder copiarse y ejecutarse tal cual, ya que estoy usando una semilla aleatoria para los puntos y desde que estoy descargando y descomprimiendo el archivo mundial en código.

Primero, creamos 100 puntos, luego usamos la función over con el argumento fn para sumar el elemento en el marco de datos. Aquí hay muchos puntos, pero eche un vistazo a Australia: 3 puntos, número 3 como etiqueta. Hasta aquí todo bien.

Ahora transformamos geometrías para poder crear buffers, transformar de nuevo y mapear esos buffers. (Incluido en el mapa anterior, ya que estoy limitado a dos enlaces). Queremos saber cuántos búferes se superponen en cada país; en el caso de Australia, a simple vista, eso es 4. No puedo entender lo que está pasando. aunque para conseguir eso con la función over. Vea mi lío de un intento en las líneas finales de código.

EDITAR: tenga en cuenta que un comentarista en r-sis-geo mencionó la función agregada, también mencionada en la pregunta 63577 de intercambio de pila, por lo que una solución / flujo podría ser a través de esa función, pero no entiendo por qué tendría que ir agregar para polypoly cuando over parece tener esa funcionalidad para otros objetos espaciales.

require(maptools)
require(sp)
require(rgdal)
require(rgeos)

download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip", destfile="world.zip")
unzip("world.zip")
world.map = readOGR(dsn=".", "TM_WORLD_BORDERS_SIMPL-0.3", stringsAsFactors = F)
orig.world.map = world.map #hold the object, since I'm going to mess with it.

#Let's create 500 random lat/long points with a single value in the data frame: the number 1
set.seed(1)
n=100
lat.v = runif(n, -90, 90)
lon.v = runif(n, -180, 180)
coords.df = data.frame(lon.v, lat.v)
val.v = data.frame(rep(1,n))
names(val.v) = c("val")
names(coords.df) = c("lon", "lat")
points.spdf = SpatialPointsDataFrame(coords=coords.df, proj4string=CRS("+proj=longlat +datum=WGS84"), data=val.v)
points.spdf = spTransform(points.spdf, CRS(proj4string(world.map)))
plot(world.map, main="World map and points") #replot the map
plot(points.spdf, col="red", pch=20, cex=1, add=T) #...and add points.

#Let's use over with the point data
join.df = over(geometry(world.map), points.spdf,  fn=sum)
plot(world.map, main="World with sum of points, 750mi buffers") #Note - happens to be the count of points, but only b/c val=1.
plot(points.spdf, col="red", pch=20, cex=1, add=T) #...and add points.
world.map@data = data.frame(c(world.map@data, join.df))
#world.map@data = data.frame(c(world.map@data, over(world.map, points.spdf, fun="sum")))
invisible(text(getSpPPolygonsLabptSlots(world.map), labels=as.character(world.map$val), cex=1))
#Note I don't love making labels like above, and am open to better ways... plus I think it's deprecated/ing

#Now buffer...
pointbuff.spdf = gBuffer(spTransform(points.spdf, CRS("+init=EPSG:3358")), width=c(750*1609.344), byid=T)
pointbuff.spdf = spTransform(pointbuff.spdf, world.map@proj4string)
plot(pointbuff.spdf, col=NA, border="pink", add=T)



#Now over with the buffer (poly %over% poly).  How do I do this?
world.map = orig.world.map
join.df = data.frame(unname(over(geometry(world.map), pointbuff.spdf, fn=sum, returnList = F)) ) #Seems I need to unname this...?
names(join.df) = c("val")
world.map@data = data.frame(c(world.map@data, join.df)) #If I don't mess with the join.df, world.map's df is a mess..
plot(world.map, main="World map, points, buffers...and a mess of wrong counts") #replot the map
plot(points.spdf, col="red", pch=20, cex=1, add=T) #...and add points.
plot(pointbuff.spdf, col=NA, border="pink", add=T)
invisible(text(getSpPPolygonsLabptSlots(world.map), labels=as.character(world.map$val), cex=1)) 
#^ But if I do strip it of labels, it seems to be misassigning the results?
# Australia should now show 4 instead of 3.  I'm obviously super confused, probably about the structure of over poly-poly returns.  Help?

Gracias por la pregunta clara y el ejemplo reproducible.

Su comprensión es correcta, y esto se reduce a un error en rgeos :: over, que se corrigió hace un mes pero aún no se ha convertido en una versión de CRAN. La siguiente es una solución alternativa si solo le interesa la cantidad de intersecciones:

world.map$val = sapply(over(geometry(world.map), pointbuff.spdf, returnList = TRUE), NROW)

Estoy usando NROWaquí en lugar de lengthpara que funcione con los rgeos incorrectos (0.3-8, de CRAN) así como con los corregidos (0.3-10, de r-forge). La sugerencia anterior de usar

a = aggregate(pointbuff.spdf, world.map, sum)

también cuenta el número de intersecciones, pero solo con la versión fija de rgeos instalada. Su ventaja, además de un nombre más intuitivo, es que devuelve directamente un Spatialobjeto, con la geometría de world.map.

Para hacer que rgeos 0.3-8 funcione, agregue

setMethod("over",
    signature(x = "SpatialPolygons", y = "SpatialPolygonsDataFrame"),
        rgeos:::overGeomGeomDF)

a su secuencia de comandos, antes de usar over.

Mientras tanto, preparé un reemplazo rápido (y mal codificado) que crea el marco de datos que necesito, ya que mi pregunta no está totalmente respondida por la solución de conteo anterior o "trabajar con los nuevos rgeos", que yo No soy lo suficientemente hábil como para entender cómo hacerlo.

Esta función es claramente (1) incompleta (observe cómo ignoro el argumento fn) y (2) ineficiente, ya que estoy llegando sin las poderosas manipulaciones de matriz / sapply de R ... (claramente vengo de otros idiomas sin ese poder) pero honestamente, todavía estoy confundido sobre lo que devuelve la estructura de la función over (¿lista de listas ...? ¿Y listas en blanco si NA?). Por lo que vale (ediciones bienvenidas), esta función hace el trabajo que necesito hacer, con éxito, e imita la acción de las otras funciones.

Ediciones bienvenidas:

overhelper <- function(pol, pol.df, fn=sum, verbose=F){
   if(verbose) {cat("Building over geometry...\n"); t=Sys.time(); t}
   geolist = over(geometry(pol), pol.df, returnList = T)
   if(verbose) {cat("Geometry done. Aggregating df. \n"); Sys.time()-t;t=Sys.time();t;}
   results = data.frame(matrix(0,nrow=length(pol), ncol=ncol(pol.df)))
   names(results) = names(pol.df)
   end = length(geolist)

   for (i in 1:end){
     if(verbose) cat(i, "...")
     results[i,] = sapply(pol.df@data[unlist(geolist[i]),], fn)
   }
   if(verbose) cat("Aggregation done! (", Sys.time()-t, ") \n Returning result vector.")
   return (results)
}