Nähe in Raum und Zeit

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
... ...         ...       ...

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.

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 .