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

12

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

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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?
Mike Dolan Fliss
źródło
Doceń przekierowanie - czy powinienem usunąć stąd i opublikować tam ponownie? Jaki jest najlepszy ruch? Dzięki.
Mike Dolan Fliss

Odpowiedzi:

5

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.

Edzer Pebesma
źródło
Bardzo pomocny, dziękuję. W szczególności chcę uczcić Twoją ofertę rozwiązania, które działa przed i po naprawie. Czy miałbyś coś przeciwko opracowaniu: (1) Na czym polega błąd, w który uderzam tutaj-rgeos :: over, zwraca przestrzenną geografię wielokątów, a nie przestrzenną ramkę danych poli? Czy niektóre funkcje po prostu nie zwracają ramek danych ...? (2) Jak to ogólnie powinno działać z agregacją i więcej? Jestem trochę zdezorientowany co do ich zamierzonych różnic i przypadków użycia. Naprawdę doceniam wasze ważenie, dziękuję. I sidenote: jakieś sugestie dotyczące zrozumienia cyklu wydawania CRAN?
Mike Dolan Fliss
Ponadto, jeśli chodzi o pierwotne pytanie: muszę policzyć liczbę ekspozycji, ale też naprawdę muszę je zsumować - na przykład liczbę świń w każdej ekspozycji. Liczenie nakładek to początek ... ale wygląda na to, że potrzebuję rozwiązania polegającego na ściągnięciu najnowszych rgeo, tak? Nie da się tego zrobić bez agregacji funkcjonalnej (nie tylko liczenia)?
Mike Dolan Fliss
(1) rgeos :: over do podpisu SpatialPolygons,SpatialPolygonsDataFramepowinien zwrócić a data.frame, ale zwraca wektor indeksu identyczny z tym, kiedy ybyłby SpatialPolygons. sp::aggregaterobi to, co robisz, w bardziej przyjazny dla użytkownika sposób, zwracając Spatialobiekt zamiast data.frame. Pakiety CRAN są utrzymywane przez wolontariuszy.
Edzer Pebesma
OK, dzięki Edzer. Brzmi to tak, jakby agregacja opierała się na rgeos, więc aby uzyskać tę funkcjonalność przed cyklem wydawania CRAN (gdziekolwiek to jest), będę musiał dowiedzieć się, jak pobrać najnowsze rgeo i z tego zrezygnować. Dziękuję Ci. I dziękuję za całą pracę nad pakietem !!
Mike Dolan Fliss
Ponadto, Edzer, wielkie dzięki za notatkę na temat R-sis-geo. Nie byłem pewien, gdzie jest lepsze miejsce do publikowania, więc cieszę się, że wątek teraz tutaj wskazuje.
Mike Dolan Fliss,
1

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)
}
Mike Dolan Fliss
źródło
1
Dodałem alternatywę, aby naprawić rgeos 0.3-8, do mojej odpowiedzi.
Edzer Pebesma