Jak działa% wielokąta przestrzennego w stosunku do% wielokąta podczas agregowania wartości wr?

Pracuję nad projektem epidemiologii środowiskowej, w którym mam narażenia punktowe (~ 2000 operacji wieprzowych - IHO). Te IHO rozpylają na pobliskie pola, ale krople wody i zapach kału mogą pokonywać kilometry. Te ekspozycje punktowe otrzymują 3 miliony buforów i chcę poznać liczbę ekspozycji IHO (różnego rodzaju - suma obornika, liczba świń, cokolwiek; najprościej, tylko liczba nakładających się buforów ekspozycji) na bloki spisu NC (~ 200 000). Bloki spisu wykluczającego (niebieskie) to (1) wszystko w 5 najbardziej zaludnionych miastach i (2) hrabstwach, które nie graniczą z hrabstwem z IHO (uwaga: zostało to zrobione z funkcją gRelate i kodami DE-9IM - bardzo zręczny!). Zobacz obraz poniżej

Ostatnim krokiem jest agregacja buforowanej reprezentacji ekspozycji dla każdego bloku spisu. Oto, gdzie jestem zakłopotany.

Do tej pory dobrze się bawiłem z funkcjami% over% w pakiecie sp, ale z winiety over rozumiem, że poly-poly i poly-line over są zaimplementowane w rgeos. Winieta obejmuje tylko line-poly i samo-odwołujące się poly, a nie agregację, więc jestem trochę zdezorientowany, jakie są moje opcje dla poly-poly z agregacją funkcji, takich jak suma lub średnia.

W przypadku testowym zapoznaj się z poniższym, nieco szczegółowym fragmentem, który działa z plikiem granic kraju na świecie. Powinno to być możliwe do skopiowania i uruchomienia w obecnej postaci, ponieważ używam losowego źródła dla punktów oraz pobieram i rozpakowuję plik świata w kodzie.

Najpierw tworzymy 100 punktów, a następnie używamy funkcji over z argumentem fn, aby dodać element do ramki danych. Jest tutaj wiele punktów, ale spójrz na Australię: 3 punkty, numer 3 jako etykieta. Jak na razie dobrze.

Teraz przekształcamy geometrie, abyśmy mogli tworzyć bufory, przekształcać je z powrotem i mapować te bufory. (Uwzględnione na poprzedniej mapie, ponieważ jestem ograniczony do dwóch linków.) Chcemy wiedzieć, ile buforów pokrywa się z każdym krajem - w przypadku Australii, na oko, to 4. Nie mogę z mojego życia zrozumieć, co się dzieje chociaż aby uzyskać to dzięki funkcji over. Zobacz mój bałagan związany z próbą w ostatnich wierszach kodu.

EDYCJA: Zauważ, że komentator r-sis-geo wspomniał o funkcji agregującej - również przywołanej w pytaniu 63577 dotyczącym wymiany stosów - więc obejście / przepływ może odbywać się za pośrednictwem tej funkcji, ale nie rozumiem, dlaczego muszę iść do agregowania dla polipowatości, gdy wydaje się, że ma tę funkcjonalność dla innych obiektów przestrzennych.

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?

Dzięki za jasne pytanie i powtarzalny przykład.

Twoje zrozumienie jest prawidłowe, a sprowadza się to do błędu w rgeos :: over, który został naprawiony miesiąc temu, ale nie został jeszcze udostępniony w wersji CRAN. Oto obejście, jeśli interesuje Cię tylko liczba skrzyżowań:

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

Używam NROWtutaj zamiast, lengthaby działał z niewłaściwymi rgeos (0,3-8, z CRAN), a także z poprawionymi (0,3-10, z r-forge). Wcześniejsza sugestia użycia

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

zlicza także liczbę skrzyżowań, ale tylko z zainstalowaną stałą wersją rgeos. Jego zaletą, oprócz bardziej intuicyjnej nazwy, jest to, że bezpośrednio zwraca Spatialobiekt o geometrii world.map.

Aby uruchomić rgeos 0.3-8, dodaj

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

do skryptu przed użyciem over.

W międzyczasie podniosłem szybki (i źle zakodowany) zamiennik, który tworzy ramkę danych, której potrzebuję, ponieważ na moje pytanie nie ma odpowiedzi w powyższym rozwiązaniu liczącym tylko „odpracować nowe rgeo”, które ja nie jestem wystarczająco wykwalifikowany, by zrozumieć, jak to zrobić.

Ta funkcja jest wyraźnie (1) niekompletna (zauważ, jak ignoruję argument fn) i (2) nieefektywna, ponieważ przychodzę do niej bez potężnych manipulacji tablicą R / sapply ... (najwyraźniej pochodzę z innych języków bez tej mocy), ale szczerze mówiąc, nadal jestem zdezorientowany, co zwraca struktura funkcji over (lista list ...? I puste listy, jeśli NA?). Jeśli chodzi o to, co jest warte (mile widziane edycje), ta funkcja wykonuje pracę, którą potrzebuję, z powodzeniem i naśladuje działanie innych funkcji.

Mile widziane zmiany:

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