Wie funktioniert das räumliche Polygon% über% Polygon, wenn Werte in r aggregiert werden?

12

Ich arbeite an einem Umweltepidemiologieprojekt, bei dem ich Punktexpositionen habe (~ 2.000 industrielle Schweineoperationen - IHOs). Diese IHOs sprühen auf Felder in der Nähe, aber die Fäkalien, Wassertropfen und Gerüche können kilometerweit wandern. Diese Punktbelichtungen erhalten also 3-mi-Puffer, und ich möchte die Anzahl der IHO-Belichtungen (verschiedener Arten - Summe der Mistmenge, Anzahl der Schweine, was auch immer; am einfachsten, nur die Anzahl der überlappenden Belichtungspuffer) pro NC-Zensusblöcke wissen (~ 200.000). Ausschlusszählungsblöcke (blau) sind (1) alles in den Top 5 der bevölkerungsreichsten Städte und (2) Landkreise, die nicht an einen Landkreis mit IHO grenzen (Anmerkung: Dies wurde mit der Funktion gRelate und den Codes DE-9IM durchgeführt). sehr glatt!). Siehe Bild unten für eine visuelle Darstellung

Bildbeschreibung hier eingeben

Der letzte Schritt besteht darin, die gepufferte Expositionsdarstellung für jeden Zensusblock zu aggregieren. Hier bin ich ratlos.

Ich hatte bisher gute Zeiten mit den% over% -Funktionen im sp-Paket, verstehe aber aus der over-Vignette, dass poly-poly und poly-line over in rgeos implementiert sind. Die Vignette deckt nur Linien-Poly und selbstreferenzierendes Poly ab und nicht die Aggregation. Daher bin ich etwas verwirrt darüber, welche Optionen für Poly-Poly mit Funktionsaggregation wie Summe oder Mittelwert zur Verfügung stehen.

Betrachten Sie für einen Testfall das folgende, etwas ausführliche Snippet, das mit der Datei mit den Weltlandgrenzen arbeitet. Dies sollte kopiert und ausgeführt werden können, da ich einen zufälligen Startwert für die Punkte verwende und die Weltdatei im Code herunterlade und entpacke.

Zuerst erstellen wir 100 Punkte und verwenden dann die over-Funktion mit dem fn-Argument, um das Element im Datenrahmen zu addieren. Hier gibt es viele Punkte, aber werfen Sie einen Blick auf Australien: 3 Punkte, Nummer 3 als Label. So weit, ist es gut.

Bildbeschreibung hier eingeben

Jetzt transformieren wir Geometrien, damit wir Puffer erstellen, zurücktransformieren und diese Puffer zuordnen können. (In der vorherigen Karte enthalten, da ich auf zwei Links beschränkt bin.) Wir möchten wissen, wie viele Puffer sich in jedem Land überlappen - im Falle Australiens sind das 4. Ich kann für mein Leben nicht abschätzen, was los ist aber um das mit der over funktion hinzubekommen. In den letzten Codezeilen sehen Sie, wie ich einen Versuch durcheinanderbringe.

BEARBEITEN: Beachten Sie, dass ein Kommentator von r-sis-geo die Aggregatfunktion erwähnte - auf die auch in der Stapelaustauschfrage 63577 verwiesen wird -, so dass ein Workaround / Flow möglicherweise durch diese Funktion erfolgt, aber ich verstehe nicht, warum ich gehen muss Für Polypolyse zu aggregieren, wenn diese Funktionalität für andere räumliche Objekte zu gelten scheint.

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
quelle
Schätzen Sie die Weiterleitung - sollte ich von hier löschen und dort neu posten? Was ist der beste Zug? Vielen Dank.
Mike Dolan Fliss

Antworten:

5

Danke für die klare Frage und das reproduzierbare Beispiel.

Ihr Verständnis ist korrekt, und dies führt zu einem Fehler in rgeos :: over, der vor einem Monat behoben wurde, aber noch nicht zu einer CRAN-Veröffentlichung geführt hat. Das Folgende ist eine Problemumgehung, wenn Sie nur an der Anzahl der Kreuzungen interessiert sind:

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

Ich benutze NROWhier stattdessen, lengthdamit es mit den falschen Rgeos (0.3-8, von CRAN) sowie den korrigierten (0.3-10, von R-Forge) funktioniert. Der frühere Vorschlag der Verwendung

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

zählt auch die Anzahl der Kreuzungen, jedoch nur mit der installierten festen rgeos-Version. Neben einem intuitiveren Namen hat dies den Vorteil, dass ein SpatialObjekt mit der Geometrie von direkt zurückgegeben wird world.map.

Fügen Sie hinzu, um Rgeos 0.3-8 zum Laufen zu bringen

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

zu Ihrem Skript, bevor Sie verwenden over.

Edzer Pebesma
quelle
Sehr hilfreich, danke. Ich möchte besonders Ihr Angebot einer Lösung würdigen, die vor und nach dem Fix funktioniert. Würde es Ihnen etwas ausmachen, Folgendes näher zu erläutern: (1) Was ist der Fehler, den ich bei here-rgeos :: over treffe, ist die Rückgabe einer räumlichen Polygongeographie, nicht eines räumlichen Polydatenrahmens? Geben einige Funktionen nicht einfach Datenrahmen zurück ...? (2) Wie soll das generell mit Aggregat und mehr funktionieren? Ich bin ein bisschen verwirrt über die beabsichtigten Unterschiede und Anwendungsfälle. Schätzen Sie wirklich Ihr Abwägen, danke. Und Nebenbemerkung: Irgendwelche Vorschläge zum Verständnis des CRAN-Freigabezyklus?
Mike Dolan Fliss
Auch in Bezug auf die ursprüngliche Frage: Ich muss die Anzahl der Expositionen zählen, aber ich muss sie auch wirklich summieren - Dinge wie die Anzahl der Schweine in jeder Exposition. Das Zählen von Überlappungen ist ein Anfang ... aber es hört sich so an, als ob ich die neuesten Rgeos einsammeln müsste, ja? Keine Möglichkeit, diese funktionale Aggregation (nicht nur Zählen) ohne sie durchzuführen?
Mike Dolan Fliss
(1) rgeos :: over für Signatur SpatialPolygons,SpatialPolygonsDataFramesollte a zurückgeben data.frame, gibt jedoch einen Indexvektor zurück, der mit dem Zeitpunkt identisch ist, zu dem ydies geschehen wäre SpatialPolygons. sp::aggregatemacht das, was Sie damit machen, benutzerfreundlicher und gibt das SpatialObjekt anstelle des data.frame. CRAN-Pakete werden von Freiwilligen gepflegt.
Edzer Pebesma
Okay, danke Edzer. Es hört sich so an, als ob Aggregat sich auf Rgeos verlässt. Um diese Funktionalität vor dem CRAN-Veröffentlichungszyklus (wann immer dies der Fall ist) zu erhalten, muss herausgefunden werden, wie die neuesten Rgeos heruntergeladen und daraus gearbeitet werden können. Vielen Dank. Und vielen Dank für all Ihre Arbeit an dem Paket!
Mike Dolan Fliss
Edzer, vielen Dank für den Hinweis zu R-sis-geo. War mir nicht sicher, wo der bessere Ort zum Posten war, also bin ich froh, dass der Thread jetzt hier zeigt.
Mike Dolan Fliss
1

Ich habe in der Zwischenzeit einen schnellen (und schlecht codierten) Über-Ersetzer entwickelt, der den Datenrahmen erstellt, den ich benötige, da meine Frage durch die oben genannte Nur-Zählen-Lösung oder "Abarbeiten der neuen Probleme", die ich habe, nicht ganz beantwortet wird bin nicht so geschickt, um zu verstehen, wie man das macht.

Diese Funktion ist eindeutig (1) unvollständig (beachte, wie ich das fn-Argument ignoriere) und (2) ineffizient, da ich ohne die mächtigen Array-Manipulationen von R / sapply darauf komme ... (eindeutig komme ich aus anderen Sprachen ohne diese Macht), aber ehrlich gesagt, ich bin immer noch verwirrt, was die Struktur der Überfunktion zurückgibt (Liste der Listen ...? Und leere Listen, wenn NA?). Für das, was es wert ist (Änderungen erwünscht), erledigt diese Funktion die Arbeit, die ich tun muss, erfolgreich und ahmt die Aktion der anderen Funktionen nach.

Bearbeitungen erwünscht:

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
quelle
1
Ich habe meiner Antwort eine Alternative hinzugefügt, um RGEOS 0.3-8 zu beheben.
Edzer Pebesma