Usando R para calcular el área de múltiples polígonos en un mapa que se interseca con otro polígono superpuesto

22

Tengo un shapefile descargado de la Encuesta de Artillería que proporciona límites de distrito electoral (división) para un condado del Reino Unido. He utilizado con éxito R para cargar el shapefile y he trazado varios mapas utilizando ggplot2 como se describe en esta pregunta . Todo está funcionando bastante bien.

Ahora me gustaría crear un nuevo polígono de forma arbitraria, agregarlo al mapa y luego calcular la población que vive en el área que se encuentra debajo de la forma, que podría cubrir o cubrir parcialmente varias divisiones. Tengo la población para cada división electoral y puedo hacer el supuesto simplificador de que la población en cada distrito está distribuida uniformemente. Eso sugiere los siguientes pasos.

1) Superponer una nueva forma en el mapa que cubre parcialmente varias divisiones electorales. Digamos que hay 3 divisiones, por el bien de la discusión. Se vería algo como esto. [Editar: excepto que en la imagen debajo de la forma se extienden 5 divisiones en lugar de 3]

2)Calculeelporcentajedeláreadecadaunadeestas3divisionesqueseintersecanconelpolígonosuperpuesto.

3)Calculelapoblaciónobteniendoelporcentajedeláreadecadadivisióncubiertaporlaformasuperpuestaymultiplicándolaporlapoblacióndecadadivisión.

Creoqueprobablementepuedaaveriguarcómocrearelpolígonoysuperponerloenelmapa,esdecir,agregarloalmarcodedatosexistenteusandolarespuestaútila this y otras preguntas. Lo que me preocupa es la tarea de calcular el porcentaje de cada división que está cubierta por la forma superpuesta. Las columnas lat y long en el marco de datos son esas extrañas cifras de OrdDance Survey OpenData (Eastings y Northings o algo así).

Entonces, mi primera pregunta es: ¿Cómo buscaría el área (o un subconjunto del área) de los polígonos que definen las fronteras de una división electoral utilizando estos datos? Porque incluso una el subconjunto significativo de este marco de datos es grande. He usado dput para crear un archivo de 500k ( que puede copiarse y pegado o descargado desde aquí ) en lugar de publicarlo en esta pregunta. El mapa que forma la base de la imagen anterior se creó con lo siguiente:

require(ggplot2)
ggplot(smalldf, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "grey50", size = 1, aes(fill = smalldf$bin))

Mi segunda pregunta es: ¿estoy usando las herramientas adecuadas? Actualmente estoy usando readShapePoly del paquete maptools para leer el shapefile. Luego uso fortify para crear un marco de datos de aproximadamente 130k líneas, adecuado para usar en ggplot . Tal vez debería estar usando un paquete diferente si hay uno con herramientas útiles para tales procesos?

    
pregunta SlowLearner 07.11.2012 - 12:02

2 respuestas

16

La respuesta de Spacedman y las sugerencias anteriores fueron útiles, pero en sí mismas no constituyen una respuesta completa. Después de un trabajo de detective por mi parte, me acerqué a una respuesta, aunque todavía no he logrado obtener gIntersection de la forma que deseo (consulte la pregunta original más arriba). Sin embargo, he logrado obtener mi nuevo polígono en SpatialPolygonsDataFrame.

ACTUALIZACIÓN 2012-11-11: Parece que he encontrado una solución viable (ver más abajo). La clave fue envolver los polígonos en una llamada SpatialPolygons cuando se usa gIntersection del paquete rgeos . La salida se ve así:

[1] "Haverfordwest: Portfield ED (poly 2) area = 1202564.3, intersect = 143019.3, intersect % = 11.9%"
[1] "Haverfordwest: Prendergast ED (poly 3) area = 1766933.7, intersect = 100870.4, intersect % = 5.7%"
[1] "Haverfordwest: Castle ED (poly 4) area = 683977.7, intersect = 338606.7, intersect % = 49.5%"
[1] "Haverfordwest: Garth ED (poly 5) area = 1861675.1, intersect = 417503.7, intersect % = 22.4%"

Insertar el polígono fue más difícil de lo que pensé porque, sorprendentemente, no parece haber un ejemplo fácil de seguir de insertar una nueva forma en un shapefile existente derivado de Topnance Survey. He reproducido mis pasos aquí con la esperanza de que sea útil para otra persona. El resultado es un mapa como este.

Si/cuandoresuelvoelproblemadelaintersección,editaréestarespuestayagregarélospasosfinales,amenosque,porsupuesto,alguienmeganeyproporcioneunarespuestacompleta.Mientrastanto,todosloscomentarios/consejossobremisoluciónsonbienvenidos.

Elcódigosigue.

require(sp)#theclassesandmethodsthatmakeupspatialopsinRrequire(maptools)#toolsforreadingandmanipulatingspatialobjectsrequire(mapdata)#includesgoodvectormapsofworldpoliticalboundaries.require(rgeos)require(rgdal)require(gpclib)require(ggplot2)require(scales)gpclibPermit()##DownloadtheOrdnanceSurveyBoundary-Linedata(large!)fromthisURL:##https://www.ordnancesurvey.co.uk/opendatadownload/products.html##thenextractallthefilestoalocalfolder.##Readtheelectoraldivision(ward)boundariesfromtheshapefileshp1<-readOGR("C:/test", layer = "unitary_electoral_division_region")
## First subset down to the electoral divisions for the county of Pembrokeshire...
shp2 <- shp1[shp1$FILE_NAME == "SIR BENFRO - PEMBROKESHIRE" | shp1$FILE_NAME == "SIR_BENFRO_-_PEMBROKESHIRE", ]
## ... then the electoral divisions for the town of Haverfordwest (this could be done in one step)
shp3 <- shp2[grep("haverford", shp2$NAME, ignore.case = TRUE),]

## Create a matrix holding the long/lat coordinates of the desired new shape;
## one coordinate pair per line makes it easier to visualise the coordinates
my.coord.pairs <- c(
                    194500,215500,
                    194500,216500,
                    195500,216500,
                    195500,215500,
                    194500,215500)

my.rows <- length(my.coord.pairs)/2
my.coords <- matrix(my.coord.pairs, nrow = my.rows, ncol = 2, byrow = TRUE)

## The Ordnance Survey-derived SpatialPolygonsDataFrame is rather complex, so
## rather than creating a new one from scratch, copy one row and use this as a
## template for the new polygon. This wouldn't be ideal for complex/multiple new
## polygons but for just one simple polygon it seems to work
newpoly <- shp3[1,]

## Replace the coords of the template polygon with our own coordinates
[email protected][[1]]@Polygons[[1]]@coords <- my.coords

## Change the name as well
[email protected]$NAME <- "zzMyPoly" # polygons seem to be plotted in alphabetical
                                 # order so make sure it is plotted last

## The IDs must not be identical otherwise the spRbind call will not work
## so use the spCHFIDs to assign new IDs; it looks like anything sensible will do
newpoly2 <- spChFIDs(newpoly, paste("newid", 1:nrow(newpoly), sep = ""))

## Now we should be able to insert the new polygon into the existing SpatialPolygonsDataFrame
shp4 <- spRbind(shp3, newpoly2)

## We want a visual check of the map with the new polygon but
## ggplot requires a data frame, so use the fortify() function
mydf <- fortify(shp4, region = "NAME")

## Make a distinction between the underlying shapes and the new polygon
## so that we can manually set the colours
mydf$filltype <- ifelse(mydf$id == 'zzMyPoly', "colour1", "colour2")

## Now plot
ggplot(mydf, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "black", size = 1, aes(fill = mydf$filltype)) +
    scale_fill_manual("Test", values = c(alpha("Red", 0.4), "white"), labels = c("a", "b"))

## Visual check, successful, so back to the original problem of finding intersections
overlaid.poly <- 6 # This is the index of the polygon we added
num.of.polys <- length([email protected])
all.polys <- 1:num.of.polys
all.polys <- all.polys[-overlaid.poly] # Remove the overlaid polygon - no point in comparing to self
all.polys <- all.polys[-1] ## In this case the visual check we did shows that the
                           ## first polygon doesn't intersect overlaid poly, so remove

## Display example intersection for a visual check - note use of SpatialPolygons()
plot(gIntersection(SpatialPolygons([email protected][3]), SpatialPolygons([email protected][6])))

## Calculate and print out intersecting area as % total area for each polygon
areas.list <- sapply(all.polys, function(x) {
    my.area <- [email protected][[x]]@Polygons[[1]]@area # the OS data contains area
    intersected.area <- gArea(gIntersection(SpatialPolygons([email protected][x]), SpatialPolygons([email protected][overlaid.poly])))
    print(paste([email protected]$NAME[x], " (poly ", x, ") area = ", round(my.area, 1), ", intersect = ", round(intersected.area, 1), ", intersect % = ", sprintf("%1.1f%%", 100*intersected.area/my.area), sep = ""))
    return(intersected.area) # return the intersected area for future use
      })
    
respondido por el SlowLearner 09.11.2012 - 09:33
15

No utilice readShapePoly: ignora la especificación de la proyección. Utilice readOGR del paquete sp.

Para operaciones geográficas como la superposición de polígonos, consulte el paquete rgeos.

Literalmente, lo último que debes hacer es jugar con fortify y ggplot. Mantenga sus datos en objetos de clase sp, trátelos con gráficos básicos y deje el ggplot sugar hasta el final de un proyecto y necesita algunos gráficos bonitos.

    
respondido por el Spacedman 07.11.2012 - 15:51

Lea otras preguntas en las etiquetas