Verwenden von R zum Berechnen der Fläche mehrerer Polygone auf einer Karte, die sich mit einem anderen überlagerten Polygon schneiden

Ich habe ein Shapefile von der Ordnance Survey heruntergeladen, das die Grenzen der Wahlbezirke (Abteilungen) für eine Grafschaft des Vereinigten Königreichs angibt. Ich habe erfolgreich R verwendet, um das Shapefile zu laden, und verschiedene Maps mit ggplot2der in dieser Frage beschriebenen Methode gezeichnet . Es funktioniert alles ziemlich gut.

Jetzt möchte ich ein neues Polygon mit einer beliebigen Form erstellen, es der Karte hinzufügen und dann die Bevölkerung in dem Bereich berechnen, der unter der Form liegt und möglicherweise mehrere Bereiche abdeckt oder teilweise abdeckt. Ich habe die Bevölkerung für jede Wahlabteilung und kann vereinfachend davon ausgehen, dass die Bevölkerung in jeder Gemeinde gleichmäßig verteilt ist. Das legt die folgenden Schritte nahe.

1) Überlagern Sie eine neue Form auf der Karte, die teilweise mehrere Wahlkreise abdeckt. Nehmen wir an, es gibt 3 Abteilungen, um das Argument zu verdeutlichen. Es würde ungefähr so ​​aussehen. [Bearbeiten: außer dass im Bild unter der Form 5 Unterteilungen statt 3 überspannt]

2) Berechnen Sie den Prozentsatz der Fläche jeder dieser drei Unterteilungen, die sich mit dem überlagerten Polygon schneidet.

3) Schätzen Sie die Bevölkerung, indem Sie den Prozentsatz der Fläche jeder Abteilung ermitteln, die von der überlagerten Form bedeckt ist, und diesen Prozentsatz mit der Bevölkerung jeder Abteilung multiplizieren.

Ich denke, ich kann wahrscheinlich herausfinden, wie das Polygon erstellt und auf der Karte überlagert wird, dh es wird dem vorhandenen Datenrahmen hinzugefügt, wobei die nützliche Antwort auf diese und andere Fragen verwendet wird. Das bisschen, das mich beunruhigt, ist die Aufgabe, den Prozentsatz jeder Abteilung zu berechnen, die von der überlagerten Form abgedeckt wird. Die Spalten latund longim Datenrahmen sind die seltsamen OpenData-Zahlen von Ordnance Survey (Eastings und Northings oder so).

Meine erste Frage lautet also: Wie würde ich den Bereich (oder eine Teilmenge des Bereichs) der Polygone, die die Grenzen einer Wahldivision definieren, anhand dieser Daten ermitteln? Da selbst eine bedeutende Teilmenge dieses dputDatenrahmens groß ist, habe ich eine 500-KB-Datei erstellt ( die hier kopiert und eingefügt oder heruntergeladen werden kann ), anstatt sie in dieser Frage zu veröffentlichen. Die Karte, die die Basis für das obige Bild bildet, wurde wie folgt erstellt:

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

Meine zweite Frage ist: verwende ich die richtigen Werkzeuge? Derzeit benutze ich readShapePolyaus dem maptoolsPaket, um das Shapefile zu lesen. Anschließend fortifyerstelle ich einen Datenrahmen mit ca. 130.000 Zeilen, der für die Verwendung in geeignet ist ggplot. Vielleicht sollte ich ein anderes Paket verwenden, wenn es eines mit nützlichen Werkzeugen für solche Prozesse gibt?

Die obigen Antworten und Hinweise von Spacedman waren nützlich, stellen jedoch an sich keine vollständige Antwort dar. Nach einigen Detektivarbeiten von meiner Seite bin ich einer Antwort näher gekommen, obwohl ich es noch nicht geschafft habe, so zu werden, gIntersectionwie ich es will (siehe ursprüngliche Frage oben). Aber ich habe es geschafft , meine neue Polygon in die SpatialPolygonsDataFrame zu bekommen.

UPDATE 11.11.2012: Ich habe anscheinend eine brauchbare Lösung gefunden (siehe unten). Der Schlüssel bestand darin, die Polygone in einen SpatialPolygonsAufruf einzuschließen, wenn sie gIntersectionaus dem rgeosPaket verwendet wurden. Die Ausgabe sieht folgendermaßen aus:

[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%"

Das Einfügen des Polygons war schwieriger als gedacht, da es überraschenderweise kein leicht verständliches Beispiel für das Einfügen einer neuen Form in ein vorhandenes Shapefile gibt, das von Ordnance Survey abgeleitet wurde. Ich habe meine Schritte hier in der Hoffnung reproduziert, dass es jemand anderem nützlich sein wird. Das Ergebnis ist eine Karte wie diese.

Wenn ich das Kreuzungsproblem löse, bearbeite ich diese Antwort und füge die letzten Schritte hinzu, es sei denn, jemand schlägt mich und gibt eine vollständige Antwort. In der Zwischenzeit sind alle Kommentare / Ratschläge zu meiner Lösung willkommen.

Code folgt.

require(sp) # the classes and methods that make up spatial ops in R
require(maptools) # tools for reading and manipulating spatial objects
require(mapdata) # includes good vector maps of world political boundaries.
require(rgeos)
require(rgdal)
require(gpclib)
require(ggplot2)
require(scales)
gpclibPermit()

## Download the Ordnance Survey Boundary-Line data (large!) from this URL:
## https://www.ordnancesurvey.co.uk/opendatadownload/products.html
## then extract all the files to a local folder.
## Read the electoral division (ward) boundaries from the shapefile
shp1 <- 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
newpoly@polygons[[1]]@Polygons[[1]]@coords <- my.coords

## Change the name as well
newpoly@data$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(shp4@polygons)
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(shp4@polygons[3]), SpatialPolygons(shp4@polygons[6])))

## Calculate and print out intersecting area as % total area for each polygon
areas.list <- sapply(all.polys, function(x) {
    my.area <- shp4@polygons[[x]]@Polygons[[1]]@area # the OS data contains area
    intersected.area <- gArea(gIntersection(SpatialPolygons(shp4@polygons[x]), SpatialPolygons(shp4@polygons[overlaid.poly])))
    print(paste(shp4@data$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
      })

Verwenden Sie nicht readShapePoly - es ignoriert die Projektionsspezifikation. Verwenden Sie readOGR aus dem sp-Paket.

Überprüfen Sie für geografische Operationen wie Ihre Polygonüberlagerung das Rgeos-Paket.

Buchstäblich das Letzte, was Sie tun sollten, ist mit Festung und Verschwörung zu spielen. Bewahren Sie Ihre Daten in sp-class-Objekten auf, zeichnen Sie sie mit Basisgrafiken und lassen Sie den ggplot-Zucker bis zum Ende eines Projekts, und Sie benötigen einige hübsche Plots.