Punkte mit R zu Gittern zusammenfassen

14

Ich habe eine Frage zur räumlichen Aggregation in R. Ich versuche, einen Punktdatensatz zu einem Raster zusammenzufassen. Ich bin mir jedoch nicht sicher, wie ich das machen soll, da ich wenig Erfahrung mit solchen Sachen habe. Ich hatte gehofft, dass irgendjemand von Ihnen eine nützliche Anleitung / eine mögliche Lösung hat.

Mein Standpunkt ist ein Datensatz mit georeferenzierten Daten zu Konfliktereignissen in Afrika (siehe www.acleddata.com). Die Punkte sind mit Breiten- / Längenkoordinaten georeferenziert und enthalten Daten zu Ereignistyp und -zeit. Ich möchte diese Punkte zu einem 1x1-Grad-Raster zusammenfassen.

Daher sollte eine Gitterzelle die Informationen der Datenpunkte enthalten, wenn ein Ereignis in dieser Gitterzelle aufgetreten ist. Das Endprodukt sollte ein Datenrahmen sein oder etwas, das ich in eine CSV-Datei exportieren kann, da die Daten in einem Panel-Datensatz für statistische Analysen verwendet werden sollen.

Bisher habe ich die Daten und das Shapefile mit dem folgenden Code geladen und geplottet. Ich glaube, dass ich die Over- Funktion aus dem SP- Paket verwenden sollte, um zu aggregieren, aber ich weiß nicht wie. Hoffe, einer von euch kann helfen.

Den Code, den ich bisher verwendet habe, finden Sie hier mit dem entsprechenden visuellen Ergebnis dort .

Vorschläge dazu in QGIS sind ebenfalls willkommen.

PferdedasJahr
quelle
Dies ist eine schnelle und einfache Operation, die nur ein wenig Arithmetik erfordert. Aber in welchem ​​Format soll die Ausgabe erfolgen? "CSV" legt nur nahe, dass es sich um eine relationale Tabelle handeln sollte, dies stellt jedoch ein Problem dar: Wenn Sie aggregieren, entspricht jede Zelle möglicherweise einer unterschiedlichen Anzahl von Punkten. Normalerweise wählen Sie eine von zwei Optionen: Sie geben entweder einen Datensatz pro Punkt aus (einschließlich der ID der enthaltenen Zelle) oder Sie geben einen Datensatz pro Zelle aus und fügen statistische Zusammenfassungen der darin enthaltenen Punkte hinzu. Welche brauchst du
Whuber
1
Entschuldigung, das habe ich nicht angegeben. Was ich brauche, ist ein Datensatz pro Zelle . Ich verwende die CSV-Datei zu Paneldaten in machen zell Jahr Format.
Pferd des Jahres

Antworten:

12

Die heruntergeladenen Daten enthalten einige offensichtliche Positionsfehler. Als Erstes beschränken Sie die Koordinaten auf sinnvolle Werte:

data.df <- read.csv("f:/temp/All_Africa_1997-2011.csv", header=TRUE, sep=",",row.names=NULL)
data.df <- subset(data.df, subset=(LONGITUDE >= -180 & LATITUDE >= -90))

Das Berechnen von Gitterzellenkoordinaten und -identifikatoren ist lediglich eine Frage des Abschneidens der Dezimalstellen von den Breiten- und Längengradenwerten. (Bei beliebigen Rastern zentrieren und skalieren Sie diese zunächst auf die Einheitsgröße, kürzen die Dezimalstellen und skalieren und zentrieren sie dann wieder in ihre ursprüngliche Position zurück, wie im Code jiunten gezeigt.) Wir können diese Koordinaten zu eindeutigen Bezeichnern kombinieren. Hängen Sie sie an den Eingabedatenrahmen an, und schreiben Sie den erweiterten Datenrahmen als CSV-Datei. Es wird einen Datensatz pro Punkt geben:

ji <- function(xy, origin=c(0,0), cellsize=c(1,1)) {
  t(apply(xy, 1, function(z) cellsize/2+origin+cellsize*(floor((z - origin)/cellsize))))
}
JI <- ji(cbind(data.df$LONGITUDE, data.df$LATITUDE))
data.df$X <- JI[, 1]
data.df$Y <- JI[, 2]
data.df$Cell <- paste(data.df$X, data.df$Y)

Möglicherweise möchten Sie stattdessen eine Ausgabe, die Ereignisse in jeder Rasterzelle zusammenfasst. Um dies zu veranschaulichen, berechnen wir die Anzahl pro Zelle und geben diese aus, einen Datensatz pro Zelle:

counts <- by(data.df, data.df$Cell, function(d) c(d$X[1], d$Y[1], nrow(d)))
counts.m <- matrix(unlist(counts), nrow=3)
rownames(counts.m) <- c("X", "Y", "Count")
write.csv(as.data.frame(t(counts.m)), "f:/temp/grid.csv")

Ändern Sie für andere Zusammenfassungen das functionArgument in der Berechnung von counts. (Verwenden Sie alternativ eine Tabellenkalkulations- oder Datenbanksoftware, um die erste Ausgabedatei nach Zellenkennung zusammenzufassen.)

Lassen Sie uns zur Kontrolle die Zählungen mithilfe der Rastermitten zuordnen , um die Kartensymbole zu lokalisieren. (Die Punkte im Mittelmeer, in Europa und im Atlantik haben verdächtige Standorte: Ich vermute, dass viele davon auf eine Vermischung von Breiten- und Längengraden bei der Dateneingabe zurückzuführen sind.)

count.max <- max(counts.m["Count",])
colors = sapply(counts.m["Count",], function(n) hsv(sqrt(n/count.max), .7, .7, .5))
plot(counts.m["X",] + 1/2, counts.m["Y",] + 1/2, cex=sqrt(counts.m["Count",]/100),
     pch = 19, col=colors,
     xlab="Longitude of cell center", ylab="Latitude of cell center",
     main="Event counts within one-degree grid cells")

Afrika karte

Dieser Workflow ist jetzt

  • Gründlich dokumentiert (durch den RCode selbst),

  • Reproduzierbar (durch erneutes Ausführen dieses Codes),

  • Erweiterbar (durch offensichtliche Änderung des Codes) und

  • Ziemlich schnell (die gesamte Operation dauert weniger als 10 Sekunden, um diese 53052 Beobachtungen zu verarbeiten).

whuber
quelle
Code ist perfekt reproduzierbar. Ich habe jedoch noch eine weitere Frage. Wie hänge ich anstelle einer Zusammenfassung die Informationen aus der Eingabedatendatei an die Zelle im erstellten Raster an?
Horseoftheyear
1
Das ist nicht möglich , mit einer Leistung zu tun Tabelle , da die vollständigen Informationen für Zellen mit variabler Länge haben. Der richtige Weg, dies aufzuzeichnen, ist mit der ersten von mir gezeigten Ausgabeform: ein Datensatz pro Punkt mit einem Zellidentifizierungsattribut. Eines dieser beiden Formate - die Pro-Punkt- und Pro-Zellen-Tabelle - wird von dem von Ihnen verwendeten Statistikprogramm erwartet.
whuber
1
Ach ok Ich verstehe was du meinst. Es muss nur ein Raster für alle Zellen erstellt und zusammengeführt werden. Danke für die Hilfe.
Horseoftheyear
3

Nun, was Sie wollen, ist ein sogenannter "Spatial Join", der zwei Shapefiles miteinander vergleicht und die Summe (Zählnummer) der resultierenden Attributtabelle zuordnet. Wenn Sie nach "Spatial Join in R" suchen, finden Sie auch hier auf GIS.Stackexchange zahlreiche Beispiele. Ich googelte schnell und fand zum Beispiel diesen Code auf einer Mailingliste.

Wenn Sie einen räumlichen Attribut-Join in QGIS erzielen möchten, gehen Sie wie folgt vor:

  • Speichern Sie Ihre Shapes als .shp-Dateien (Befehl writeOGR aus dem rgdal-Paket)
  • Laden Sie sie in QGIS. Erstellen Sie Ihr Vektorraster über das MMQGIS-Plugin (Erstellen -> Rasterebene erstellen) mit entsprechender Skalierung neu.
  • Verwenden Sie das Werkzeug "Attribute verbinden" aus dem Menü Vektor -> Datenverwaltung. Wählen Sie ein Attribut Ihres Punkt-Layers aus (dies kann eine einfache Spalte sein, die TRUE (1) oder FALSE (0) für verschiedene Konfliktereignisse darstellt).
  • Wählen Sie Ihr Raster und summieren Sie alle Vorkommen und führen Sie aus. Danach würde ich auch Ihr Raster mit einer Form des afrikanischen Kontinents beschneiden.

Wenn der Join irgendwie fehlschlägt (funktioniert bei mir nicht immer), bleiben Sie bei SEXTANTE und suchen Sie nach der SAGA-Toolbox, die auch sehr gute Join-Funktionen bietet.

Brachvogel
quelle
Obwohl dies eine Lösung ist, ist sie besonders komplex und ineffizient, da das Zusammenfassen von Punkten zu einem Raster nur eine Frage von wenigen einfachen arithmetischen Operationen ist, die Rsich bei auszeichnen. Die Verwendung von Shapefiles, rgdalQGIS und Sextante ist ein bisschen wie die Empfehlung, dass jemand eine moderne automatisierte Industrieanlage anmietet, um zwei Boards zusammenzunageln :-).
Whuber
Ich werde diesen Ansatz dieses Wochenende versuchen. In naher Zukunft möchte ich möglicherweise verschiedene Formdateien miteinander kombinieren, damit dies nützlich sein kann. Danke für den Input und die Anregungen.
Horseoftheyear
@whuber: Das stimmt, aber wenn Sie Ihre Ausgabe verteilen und vielleicht stylen möchten, ist ein Shapefile die naheliegende Wahl. Trotzdem schönes R-Beispiel!
Brachvogel
Ich habe es endlich versucht. Das Problem bei diesem Ansatz ist jedoch, dass alle Beobachtungen zu einem Polygon zusammengefasst werden. Ich möchte zwar im Idealfall die Informationen über verschiedene Ereignisse im Laufe der Zeit aufbewahren. Aber es könnte sein, dass ich etwas falsch gemacht habe.
Horseoftheyear