Bagaimana cara kerja poligon spasial% lebih dari% poligon ketika mengumpulkan nilai dalam r?

Saya sedang mengerjakan proyek epidemiologi lingkungan di mana saya memiliki paparan titik (~ 2.000 operasi babi industri - IHO). IHO ini disemprotkan pada ladang terdekat, tetapi tetesan air dan bau tinja dapat menempuh jarak bermil-mil. Jadi eksposur titik ini mendapatkan 3mi buffer, dan saya ingin tahu jumlah eksposur IHO (dari berbagai jenis - jumlah kotoran, jumlah babi, apa pun; paling sederhana, hanya jumlah buffer eksposur yang tumpang tindih) per blok sensus NC (~ 200.000). Blok sensus pengecualian (biru) adalah (1) apa pun di 5 kota terpadat teratas dan (2) kabupaten yang tidak membatasi kabupaten dengan IHO di dalamnya (catatan: yang dilakukan dengan fungsi gRelate dan kode DE-9IM - sangat apik!). Lihat gambar di bawah ini untuk visual

Langkah terakhir adalah untuk menggabungkan representasi keterpaparan yang disangga ke setiap blok sensus. Di sinilah saya bingung.

Saya sudah bersenang-senang dengan fungsi% over% dalam paket sp sejauh ini, tetapi pahami dari sketsa berlebihan bahwa poly-poly dan poly-line over diimplementasikan dalam rgeos. Sketsa hanya mencakup poli baris dan poli referensi sendiri, dan tidak dengan agregasi, jadi saya agak bingung tentang apa opsi saya untuk poli-poli dengan agregasi fungsi, seperti jumlah atau rata-rata.

Untuk kasus uji, pertimbangkan cuplikan di bawah ini, yang agak bertele-tele yang bekerja dengan file perbatasan negara dunia. Ini harus dapat disalin dan dijalankan apa adanya, karena saya menggunakan seed acak untuk poin dan karena saya mengunduh dan membuka ritsleting file dunia dalam kode.

Pertama, kita membuat 100 poin, lalu menggunakan fungsi over dengan argumen fn untuk menjumlahkan elemen dalam bingkai data. Ada banyak poin di sini, tetapi lihatlah Australia: 3 poin, nomor 3 sebagai label. Sejauh ini baik.

Sekarang kita mengubah geometri sehingga kita dapat membuat buffer, mengubah kembali, dan memetakan buffer itu. (Termasuk di peta sebelumnya, karena saya terbatas pada dua tautan.) Kami ingin tahu berapa banyak penyangga yang tumpang tindih dengan masing-masing negara - dalam kasus Australia, menurut mata, angka itu adalah 4. Saya tidak bisa seumur hidup mencari tahu apa yang terjadi meskipun untuk mendapatkan itu dengan fungsi over. Lihat kekacauan upaya saya di baris kode terakhir.

EDIT: Perhatikan bahwa komentator di r-sis-geo menyebutkan fungsi agregat - juga dirujuk pada pertanyaan pertukaran tumpukan 63577 - jadi pekerjaan di sekitar / aliran mungkin melalui fungsi itu, tapi saya tidak mengerti mengapa saya harus pergi untuk agregat untuk polipoly ketika over tampaknya memiliki fungsi itu untuk objek spasial lainnya.

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?

Terima kasih atas pertanyaan yang jelas dan contoh yang dapat direproduksi.

Pemahaman Anda benar, dan ini bermuara pada bug di rgeos :: over, yang diperbaiki sebulan lalu tetapi belum membuatnya menjadi rilis CRAN. Berikut ini adalah penyelesaian jika Anda hanya tertarik pada jumlah persimpangan:

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

Saya menggunakan di NROWsini bukannya lengthsehingga bekerja dengan rgeos yang salah (0,3-8, dari CRAN) serta yang diperbaiki (0,3-10, dari r-forge). Saran sebelumnya untuk menggunakan

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

juga menghitung jumlah persimpangan, tetapi hanya dengan versi rgeo tetap diinstal. Keuntungannya, selain nama yang lebih intuitif, adalah bahwa ia langsung mengembalikan Spatialobjek, dengan geometri world.map.

Untuk mengaktifkan rgeos 0.3-8, tambahkan

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

ke skrip Anda, sebelum Anda gunakan over.

Saya membuat over-replacer cepat (dan kode buruk) sementara itu yang menciptakan bingkai data yang saya butuhkan, karena pertanyaan saya tidak cukup dijawab oleh solusi penghitungan-hanya di atas atau "bekerja dari rgeo baru", yang saya Saya tidak cukup terampil untuk memahami bagaimana melakukannya.

Fungsi ini jelas (1) tidak lengkap (perhatikan bagaimana saya mengabaikan argumen fn) dan (2) tidak efisien, karena saya datang tanpa manipulasi array yang kuat R / sapply ... (jelas saya datang dari bahasa lain tanpa kekuatan itu) tetapi jujur, saya masih bingung apa struktur fungsi kembali kembali (daftar daftar ...? Dan daftar kosong jika NA?). Untuk apa nilainya (suntingan selamat datang), fungsi ini melakukan pekerjaan yang perlu saya lakukan, berhasil, dan meniru tindakan fungsi lainnya.

Suntingan selamat datang:

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