Nähe in Raum und Zeit

10

Ich habe einige Punktdaten, die die täglichen Lat-Lon-Standorte eines Tieres mit einem zugehörigen Zeitstempel darstellen.

Ich möchte alle Punkte identifizieren, an denen STATIONARY = TRUE ist. Ein Punkt gilt als stationär, wenn ein 100 km langer Puffer um ihn herum zusätzliche (z. B.) 5 zeitlich benachbarte Punkte überlappt . Wenn also Tag 10 mein Punkt von Interesse ist, möchte ich fragen, ob 5 zeitlich benachbarte Tage innerhalb eines Puffers von 100 km von diesem Punkt liegen. Wenn die Tage 5,6,7,8 & 9; ODER Tage 11, 12, 13, 14 & 15; ODER Tage 8,9,11,12,13 (usw.) befinden sich im Puffer, dann STATIONARY = TRUE. Wenn sich jedoch die Tage 5, 7, 9, 11 und 13 innerhalb des Puffers befinden, jedoch nicht die alternativen (geraden) Tage dazwischen, dann ist STATIONARY = FALSE

Ich denke, eine Art beweglicher Fensterpuffer wird die Lösung bieten, aber ich weiß nicht, wie ich das implementieren soll.

Ich habe versucht, dieses Problem sowohl in ArcGIS als auch in R in den Griff zu bekommen, hatte aber bisher keine Gehirnwellen. Dies ist das, was ich einer Lösung am nächsten komme, aber es passt nicht ganz, glaube ich nicht: Identifizierung aufeinanderfolgender Punkte innerhalb eines bestimmten Puffers

Hier sind einige Dummy-Daten, die sich meiner Datenstruktur annähern (obwohl ich in Wirklichkeit zweimal täglich Standorte habe (Mittag und Mitternacht), wobei einige Standorte fehlen - aber darüber werde ich mich später Gedanken machen).

x<-seq(0,15,length.out=20)
y<-seq(10,-10,length.out=20)
t<-seq(as.POSIXct('2013-07-01'), length.out = 20, by = "days")
data<-data.frame(cbind(x,y,t=as.data.frame.POSIXct(t)))


            x           y          t
1   0.0000000  10.0000000 2013-07-01
2   0.7894737   8.9473684 2013-07-02
3   1.5789474   7.8947368 2013-07-03
4   2.3684211   6.8421053 2013-07-04
5   3.1578947   5.7894737 2013-07-05
6   3.9473684   4.7368421 2013-07-06
7   4.7368421   3.6842105 2013-07-07
... ...         ...       ...
Tom Finch
quelle
1
Frage? Angenommen, alle 10 Punkte befinden sich im Puffer und Sie haben eine Datumstrennung (ab Tag 1) von 1-3-4-12-13-20-21-22-29-30. Dann sagen Sie, dass Sie nur daran interessiert sind, Punkte auszuwählen das sind in den Tagen 1,2,3,4 & 12?
Hornbydd
Nein, ich würde mich nur für die Tage 1 bis 4 interessieren. Wenn das Tier den Puffer "verlässt" und dann am 12. Tag (oder 6. Tag) zurückkehrt, würde dies diese stationäre Periode "aufheben" - dh das Tier muss sich am Tag 1-2-3-4-5 für im Puffer befinden der Punkt in der Mitte des zu zählenden Puffers. Sinn ergeben? Ich bin mir nicht sicher ..
Tom Finch
1
Nur um zu überprüfen, ob der Punkt von Interesse Tag 7 war, dann wären Sie an Punkten interessiert, die für die Tage 7,8,9,10 und 11 innerhalb von 100 km liegen?
Hornbydd
Punkt 7 würde als stationärer Punkt ausgewählt, wenn die Tage 8, 9, 10, 11 und 12 100 km zurücklegen würden. Oder Tage 5,6,8,9,10. Ein Punkt wird also ausgewählt, wenn sich alle anderen 5 zeitlich benachbarten Punkte (die 5 vorherigen Tage, die 5 folgenden Tage oder einige Tage auf beiden Seiten) alle im Puffer befinden. Ich denke, das bewegliche Fenster ist die beste Art, es zu konzipieren. Für jeden "Brennpunkt" kann jeder Punkt vergessen werden, der länger als 5 Tage in der Vergangenheit / Zukunft liegt. Ich kann meine ursprüngliche Frage aktualisieren, da ich sie jetzt etwas besser verstehe ...
Tom Finch
Wie ist das Format der Daten? Haben Sie beispielsweise jede Zeit / jeden Ort als Vektorpunkt in einem Shapefile und eine Attributtabelle, in der die Zeit gespeichert ist? Oder werden jedes Mal / jeder Ort separat in verschiedenen Shapefiles gespeichert? Sind die Daten nicht einmal in einem Geodatenformat und einfach in einer Excel-Datei? Das zu wissen würde uns helfen zu antworten.

Antworten:

12

Lassen Sie uns dies in einfache Teile zerlegen. Auf diese Weise wird die gesamte Arbeit in nur einem halben Dutzend Zeilen leicht zu testenden Codes erledigt.

Zunächst müssen Sie Entfernungen berechnen. Da die Daten in geografischen Koordinaten angegeben sind, gibt es hier eine Funktion zum Berechnen von Entfernungen auf einem sphärischen Datum (unter Verwendung der Haversine-Formel):

#
# Spherical distance.
# `x` and `y` are (long, lat) pairs *in radians*.
dist <- function(x, y, R=1) {
  d <- y - x
  a <- sin(d[2]/2)^2 + cos(x[2])*cos(y[2])*sin(d[1]/2)^2
  return (R * 2*atan2(sqrt(a), sqrt(1-a)))
}

Ersetzen Sie dies durch Ihre bevorzugte Implementierung, wenn Sie dies wünschen (z. B. eine mit einem ellipsoiden Datum).

Als nächstes müssen wir die Abstände zwischen jedem "Basispunkt" (der auf Staionarität geprüft wird) und seiner zeitlichen Nachbarschaft berechnen. Das ist einfach eine Frage der Bewerbung distfür die Nachbarschaft:

#
# Compute the distances between an array of locations and a base location `x`.
dist.array <- function(a, x, ...) apply(a, 1, function(y) dist(x, y, ...))

Drittens - dies ist die Schlüsselidee - werden stationäre Punkte gefunden, indem Nachbarschaften von 11 Punkten mit mindestens fünf in einer Reihe erkannt werden, deren Abstände ausreichend klein sind. Lassen Sie uns dies etwas allgemeiner implementieren, indem wir die Länge der längsten Teilsequenz wahrer Werte innerhalb eines logischen Arrays boolescher Werte bestimmen:

#
# Return the length of the longest sequence of true values in `x`.
max.subsequence <- function(x) max(diff(c(0, which(!x), length(x)+1)))

(Wir finden die Positionen der falschen Werte in der richtigen Reihenfolge und berechnen ihre Unterschiede: Dies sind die Längen der Teilsequenzen nicht falscher Werte. Die größte solche Länge wird zurückgegeben.)

Viertens wenden wir max.subsequencean, um stationäre Punkte zu erkennen.

#
# Determine whether a point `x` is "stationary" relative to a sequence of its
# neighbors `a`.  It is provided there is a sequence of at least `k`
# points in `a` within distance `radius` of `x`, where the earth's radius is
# set to `R`.
is.stationary <- function(x, a, k=floor(length(a)/2), radius=100, R=6378.137) 
  max.subsequence(dist.array(a, x, R) <= radius) >= k

Das sind alle Werkzeuge, die wir brauchen.


Lassen Sie uns als Beispiel einige interessante Daten mit einigen Klumpen stationärer Punkte erstellen. Ich werde einen zufälligen Spaziergang in der Nähe des Äquators machen.

set.seed(17)
n <- 67
theta <- 0:(n-1) / 50 - 1 + rnorm(n, sd=1/2)
rho <- rgamma(n, 2, scale=1/2) * (1 + cos(1:n / n * 6 * pi))
lon <- cumsum(cos(theta) * rho); lat <- cumsum(sin(theta) * rho)

Die Arrays lonund latenthalten die Koordinaten in nPunkten von aufeinanderfolgenden Punkten. Das Anwenden unserer Werkzeuge ist nach der ersten Umrechnung in Bogenmaß unkompliziert:

p <- cbind(lon, lat) * pi / 180 # Convert from degrees to radians
p.stationary <- sapply(1:n, function(i) 
  is.stationary(p[i,], p[max(1,i-5):min(n,i+5), ], k=5))

Das Argument p[max(1,i-5):min(n,i+5), ]besagt, dass bis zu 5 Zeitschritte oder bis zu 5 Zeitschritte vom Basispunkt zurückgeschaut werden soll p[i,]. Einschließlich k=5sagt, nach einer Folge von 5 oder mehr in einer Reihe zu suchen, die innerhalb von 100 km vom Basispunkt liegen. (Der Wert von 100 km wurde als Standard in festgelegt, is.stationaryaber Sie können ihn hier überschreiben.)

Die Ausgabe p.stationaryist ein logischer Vektor, der die Stationarität anzeigt: Wir haben das, wofür wir gekommen sind. Um das Verfahren zu überprüfen, ist es jedoch am besten, die Daten und diese Ergebnisse zu zeichnen, anstatt Arrays von Werten zu untersuchen. Auf dem folgenden Plot zeige ich die Route und die Punkte. Jeder zehnte Punkt ist beschriftet, damit Sie abschätzen können, wie viele sich innerhalb der stationären Klumpen überlappen könnten. Stationäre Punkte werden in durchgehendem Rot neu gezeichnet, um sie hervorzuheben, und von ihren 100 km langen Puffern umgeben.

Zahl

plot(p, type="l", asp=1, col="Gray", 
     xlab="Longitude (radians)", ylab="Latitude (radians)")
points(p)
points(p[p.stationary, ], pch=19, col="Red", cex=0.75)
i <- seq(1, n, by=10)
#
# Because we're near the Equator in this example, buffers will be nearly 
# circular: approximate them.
disk <- function(x, r, n=32) {
  theta <- 1:n / n * 2 * pi
  return (t(rbind(cos(theta), sin(theta))*r + x))
}
r <- 100 / 6378.137  # Buffer radius in radians
apply(p[p.stationary, ], 1, function(x) 
  invisible(polygon(disk(x, r), col="#ff000008", border="#00000040")))
text(p[i,], labels=paste(i), pos=3, offset=1.25, col="Gray")

Weitere (statistisch basierte) Ansätze zum Auffinden stationärer Punkte in verfolgten Daten, einschließlich Arbeitscode, finden Sie unter /mathematica/2711/clustering-of-space-time-data .

whuber
quelle
Wow, danke! Ich freue mich darauf, mich darum zu kümmern.
Tom Finch