Używanie R do obliczania powierzchni wielu wielokątów na mapie, które przecinają się z innym nałożonym wielokątem

Mam plik kształtu pobrany z Ordnance Survey, który wyznacza granice okręgu wyborczego dla hrabstwa Zjednoczonego Królestwa. Z powodzeniem użyłem R, aby załadować plik kształtu i sporządziłem różne mapy, używając ggplot2opisu opisanego w tym pytaniu . Wszystko działa całkiem dobrze.

Teraz chciałbym stworzyć nowy wielokąt o dowolnym kształcie, dodać go do mapy, a następnie obliczyć populację zamieszkującą obszar leżący pod tym kształtem, który może obejmować lub częściowo obejmować wiele dywizji. Mam populację dla każdego działu wyborczego i mogę uprościć założenie, że populacja na każdym okręgu jest równomiernie rozmieszczona. Sugeruje to następujące kroki.

1) Nałóż nowy kształt na mapę, który częściowo obejmuje wiele dywizji wyborczych. Powiedzmy, że istnieją 3 dywizje, ze względu na argument. Wyglądałoby to mniej więcej tak. [Edycja: oprócz tego, że na obrazku poniżej kształt otacza 5 podziałów zamiast 3]

2) Oblicz procent powierzchni każdego z tych 3 działów, która przecina się z nałożonym wielokątem.

3) Oszacuj populację, uzyskując procent powierzchni każdego działu objętej nałożonym kształtem i mnożąc go przez populację każdego działu.

Myślę, że prawdopodobnie uda mi się ustalić, jak utworzyć wielokąt i nałożyć go na mapę, tj. Dodać go do istniejącej ramki danych, używając użytecznej odpowiedzi na to i inne pytania. Niepokoi mnie to zadanie ustalenia procentu każdego podziału objętego nałożonym kształtem. latI longkolumn w ramce danych są te dziwne postacie Ordnance Survey OpenData (Eastings i Northings lub coś).

Moje pierwsze pytanie brzmi: w jaki sposób mógłbym znaleźć obszar (lub jego podzbiór) wielokątów, które określają granice podziału wyborczego przy użyciu tych danych? Ponieważ nawet znaczący podzbiór tej ramki danych jest duży, dputstworzyłem plik 500k ( który można skopiować i wkleić lub pobrać stąd ) zamiast zamieszczać go w tym pytaniu. Mapa stanowiąca podstawę dla powyższego obrazu została utworzona w następujący sposób:

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

Moje drugie pytanie brzmi: czy używam odpowiednich narzędzi? Obecnie używam readShapePolyz maptoolspakietu do odczytu pliku shapefile. Następnie używam fortifydo utworzenia ramki danych o długości około 130 000 linii, odpowiedniej do użycia w ggplot. Może powinienem używać innego pakietu, jeśli jest taki z przydatnymi narzędziami do takich procesów?

Odpowiedź Spacedmana i powyższe wskazówki były przydatne, ale same w sobie nie stanowią pełnej odpowiedzi. Po kilku pracach detektywistycznych z mojej strony zbliżyłem się do odpowiedzi, chociaż nie udało mi się dostać gIntersectiontak, jak chcę (patrz oryginalne pytanie powyżej). Mimo to udało mi się przenieść mój nowy wielokąt do SpatialPolygonsDataFrame.

AKTUALIZACJA 11.11.2012: Wydaje mi się, że znalazłem realne rozwiązanie (patrz poniżej). Kluczem było zawinięcie wielokątów w SpatialPolygonswywołanie podczas korzystania gIntersectionz rgeospakietu. Dane wyjściowe wyglądają następująco:

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

Wstawienie wielokąta było trudniejsze niż myślałem, ponieważ, nieoczekiwanie, nie ma łatwego do naśladowania przykładu wstawiania nowego kształtu do istniejącego pliku kształtu pochodzącego z Ordnance Survey. Powtórzyłem tutaj swoje kroki w nadziei, że przyda się komuś innemu. Rezultatem jest mapa taka jak ta.

Jeśli / kiedy rozwiążę problem skrzyżowania, wyedytuję tę odpowiedź i dodam ostatnie kroki, chyba że, oczywiście, ktoś mnie pokona i dostarczy pełną odpowiedź. Tymczasem mile widziane są komentarze / porady dotyczące mojego rozwiązania.

Kod następuje.

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

Nie używaj readShapePoly - ignoruje specyfikację projekcji. Użyj readOGR z pakietu sp.

W przypadku operacji geograficznych, takich jak nakładka wielokąta, sprawdź pakiet rgeos.

Dosłownie ostatnią rzeczą, którą powinieneś zrobić, to grać z fortyfikacjami i ggplot. Przechowuj dane w obiektach klasy sp, drukuj je z grafiką bazową i pozostaw cukier ggplot do końca projektu, a potrzebujesz ładnych wykresów.