¿Cómo espacial polígono %sobre% polígono trabajo para cuando la agregación de los valores en r?

Question

Estoy trabajando en un epidemiología ambiental proyecto en el que he punto de exposiciones (~ 2,000 industrial granjas de cerdos - Iho). Estos Iho spray en los campos cercanos, pero las heces de las gotas de agua y el olfato pueden viajar millas. Así que estos punto de exposiciones obtener 3mi de memoria, y quiero saber el número de exposiciones de la OHI (de varios tipos - suma de la cantidad de estiércol, el número de cerdos, lo que sea, más simple, sólo el número de la superposición de la exposición buffers) por NC bloques censales (~200,000). La exclusión de los bloques censales (azul) son: (1) nada en el top 5 de las ciudades más pobladas y (2) de los condados que no frontera del condado con un IHO (nota: esto fue hecho con el gRelate función y DE-9IM códigos - muy resbaladiza!). Vea a continuación la imagen visual de

El último paso es agregar el búfer de exposición a la representación de cada bloque censal. Aquí es donde estoy perplejo.

He tenido buenos momentos con el %sobre% funciones de la sp paquete hasta ahora, pero entiende que a partir de más de viñeta con la que la poli poli poli-línea se implementan en rgeos. La viñeta sólo cubre la línea de poli y la auto-referencia poli, y no con la agregación, así que estoy un poco confundido sobre cuáles son mis opciones de poli-poly con la función de agregación, como la suma o la media.

Para un caso de prueba, considere el siguiente, algo detallado fragmento de trabajar con el mundo de las fronteras del país de archivo. Este debe ser capaz de copiarse a cabo y se ejecuta como es, ya que estoy usando una semilla aleatoria de los puntos y ya estoy de descargar y descomprimir el mundo de archivo en el código.

En primer lugar, vamos a crear 100 puntos, a continuación, utilizar la función más con el fn argumento para añadir el elemento en el marco de datos. Hay un montón de puntos, pero echa un vistazo a Australia: 3 puntos, número 3 como una etiqueta. Tan lejos, tan bueno.

Ahora podemos transformar las geometrías por lo que podemos crear buffers, transformar la espalda, y el mapa de los búferes. (Incluido en el mapa anterior, ya que estoy limitado a dos enlaces.) Queremos saber cuántos de los búferes de superposición de cada país - en el caso de Australia, por el ojo, que es de 4. Yo no puedo por la vida de mí la figura de lo que está pasando, aunque para conseguir que con la función más. Ver mi desorden de un intento, en la final de líneas de código.

EDIT: tenga en cuenta que un comentarista en r-sis-geo menciona la función de agregado - también se hace referencia en el intercambio de la pila pregunta 63577 - por lo que un trabajo de torno / de flujo podría ser a través de esa función, pero no entiendo por qué yo tendría que ir a agregado para polypoly cuando más parece que la 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?

Answer 1

Gracias por la pregunta clara y reproducible ejemplo.

Su comprensión es correcta, y esto se reduce a un error en rgeos::otra vez, que fue fijado hace un mes, pero no ha hecho en un CRAN lanzamiento. El siguiente es un trabajo en torno a si usted está interesado sólo en el número de intersecciones:

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

Estoy usando NROW aquí en lugar de length por lo que se trabaja con el mal rgeos (0.3-8, de CRAN), así como la corrección (0.3-10, r-forge). La anterior sugerencia de uso

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

también cuenta el número de intersecciones, pero sólo con el fijo rgeos versión instalada. Su ventaja, además de una más intuitiva nombre, es que se devuelve directamente un Spatial de un objeto, con la geometría de world.map.

Para obtener rgeos 0.3-8 de trabajo, agregar

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

a la secuencia de comandos, antes de utilizar over.

Answer 2

Destapé una rápida (y mal codificado)-sustituto en el ínterin que crea el marco de datos que necesito, ya que mi pregunta no es respondida por el anterior recuento única solución o "el trabajo fuera de la nueva rgeos", que no estoy lo suficientemente capacitado para entender cómo hacerlo.

Esta función es claramente (1) incompleta (observe cómo puedo ignorar el fn argumento) y (2) ineficiente, ya que estoy entrando en él sin R poderoso conjunto de manipulaciones / sapply... (claramente yo soy procedentes de otras lenguas, sin que el poder) pero, sinceramente, yo todavía estoy confundido lo que la estructura de la sobre el retorno de la función (lista de listas de...? Y en blanco listas si NA?). Para lo que vale (ediciones de bienvenida), esta función hace el trabajo que necesita realizar, con éxito, e imita la acción de las otras funciones.

Ediciones de bienvenida:

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)
}