Menggunakan R untuk menghitung luas beberapa poligon pada peta yang berpotongan dengan poligon overlay lainnya

Saya memiliki shapefile yang diunduh dari Ordnance Survey yang memberikan batas-batas (divisi) bangsal pemilihan untuk sebuah daerah di Britania Raya. Saya telah berhasil menggunakan R untuk memuat shapefile dan merencanakan berbagai peta menggunakan ggplot2seperti yang dijelaskan dalam pertanyaan ini . Semuanya bekerja dengan baik.

Sekarang saya ingin membuat poligon baru dengan bentuk sewenang-wenang, menambahkannya ke peta, lalu menghitung populasi yang tinggal di daerah yang terletak di bawah bentuk, yang mungkin menutupi atau sebagian menutupi beberapa divisi. Saya memiliki populasi untuk setiap divisi pemilihan dan saya dapat membuat asumsi penyederhanaan bahwa populasi di setiap bangsal terdistribusi secara seragam. Itu menyarankan langkah-langkah berikut.

1) Overlay bentuk baru pada peta yang sebagian mencakup beberapa divisi pemilihan. Katakanlah ada 3 divisi, demi argumen. Akan terlihat seperti ini. [Sunting: kecuali bahwa pada gambar di bawah ini bentuknya mengangkang 5 divisi daripada 3]

2) Hitung persentase luas masing-masing dari 3 divisi yang bersinggungan dengan poligon yang dilapis.

3) Perkirakan populasi dengan mendapatkan persentase area masing-masing divisi yang dicakup oleh bentuk overlay dan mengalikannya dengan populasi masing-masing divisi.

Saya pikir saya mungkin bisa mengetahui cara membuat poligon dan melapisinya di peta yaitu menambahkannya ke bingkai data yang ada menggunakan jawaban yang berguna untuk ini dan pertanyaan lain. Bagian yang membuat saya khawatir adalah tugas menghitung persentase setiap divisi yang dicakup oleh bentuk overlay. The latdan longkolom dalam frame data adalah mereka yang aneh angka Ordnance Survey OpenData (Eastings dan Northings atau sesuatu).

Jadi pertanyaan pertama saya adalah: Bagaimana cara mencari area (atau subset area) poligon yang menentukan batas-batas divisi pemilu menggunakan data ini? Karena bahkan subset yang bermakna dari kerangka data ini besar, saya telah menggunakan dputuntuk membuat file 500k ( yang dapat disalin dan ditempelkan atau diunduh dari sini ) daripada mempostingnya dalam pertanyaan ini. Peta yang membentuk dasar untuk gambar di atas dibuat dengan yang berikut:

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

Pertanyaan kedua saya adalah: apakah saya menggunakan alat yang tepat? Saat ini saya menggunakan readShapePolydari maptoolspaket untuk membaca shapefile. Saya kemudian gunakan fortifyuntuk membuat bingkai data sekitar 130k garis, cocok untuk digunakan dalam ggplot. Mungkin saya harus menggunakan paket yang berbeda jika ada satu dengan alat yang berguna untuk proses seperti itu?

Jawaban dan petunjuk Spacedman di atas bermanfaat, tetapi tidak dengan sendirinya merupakan jawaban penuh. Setelah beberapa pekerjaan detektif pada bagian saya, saya sudah semakin dekat dengan jawaban walaupun saya belum berhasil mendapatkan gIntersectioncara yang saya inginkan (lihat pertanyaan asli di atas). Namun, saya telah berhasil memasukkan poligon baru saya ke SpatialPolygonsDataFrame.

UPDATE 2012-11-11: Saya sepertinya telah menemukan solusi yang bisa diterapkan (lihat di bawah). Kuncinya adalah untuk membungkus poligon dalam SpatialPolygonspanggilan saat menggunakan gIntersectiondari rgeospaket. Outputnya terlihat seperti ini:

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

Memasukkan poligon lebih sulit daripada yang saya pikir karena, anehnya, tampaknya tidak ada contoh yang mudah diikuti untuk memasukkan bentuk baru ke dalam formulir yang diturunkan oleh Survey Ordnance Survey. Saya telah mereproduksi langkah-langkah saya di sini dengan harapan akan bermanfaat bagi orang lain. Hasilnya adalah peta seperti ini.

Jika / ketika saya menyelesaikan masalah persimpangan, saya akan mengedit jawaban ini dan menambahkan langkah-langkah terakhir, kecuali, tentu saja, ada yang mengalahkan saya untuk itu dan memberikan jawaban lengkap. Sementara itu, komentar / saran tentang solusi saya sejauh ini semuanya diterima.

Kode berikut.

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

Jangan gunakan readShapePoly - ini mengabaikan spesifikasi proyeksi. Gunakan readOGR dari paket sp.

Untuk operasi geografis seperti overlay poligon Anda, lihat paket rgeos.

Secara harfiah hal terakhir yang harus Anda lakukan adalah bermain dengan fortifikasi dan ggplot. Simpan data Anda di objek kelas-sp, plot dengan grafis dasar, dan biarkan gula ggplot hingga akhir proyek dan Anda memerlukan beberapa plot yang cantik.