Odległość w przestrzeni i czasie

Mam pewne dane punktowe, które reprezentują codzienne lokalizacje zwierzęcia w całej długości, wraz z powiązanym znacznikiem czasu.

Chciałbym zidentyfikować wszystkie punkty, w których STACJONARNY = PRAWDA. Punkt kwalifikuje się jako nieruchomy, jeśli wokół niego bufor o długości 100 km zachodzi na dodatkowe (powiedzmy) 5 tymczasowo sąsiadujących punktów. Więc jeśli moim celem jest dzień 10 , chcę zapytać, czy 5 tymczasowo sąsiadujących dni znajduje się w odległości 100 km od tego punktu. Jeśli dni 5,6,7,8 i 9; LUB dni 11,12,13,14 i 15; LUB dni 8,9,11,12,13 (itp.) Znajdują się w buforze, a następnie STATIONARY = TRUE. Jeśli jednak dni 5, 7, 9, 9, 13 i 13 znajdują się w buforze, ale nie są to dni alternatywne (parzyste) pomiędzy nimi, wówczas STACJONARNE = FAŁSZ

Myślę, że jakiś bufor ruchomego okna zapewni rozwiązanie, ale nie wiem jak to zaimplementować.

Próbowałem obejść ten problem w ArcGIS i R, ale jak dotąd nie miałem fal mózgowych. Jest to najbliższe mi rozwiązanie, ale nie do końca pasuje, nie sądzę: Identyfikacja kolejnych punktów w określonym buforze

Oto niektóre fikcyjne dane, które zbliżają się do mojej struktury danych (chociaż w rzeczywistości mam lokalizacje dwa razy dziennie (w południe i północ), a niektórych lokalizacji brakuje - ale będę się tym martwić później)

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

Podzielmy to na proste kawałki. Dzięki temu cała praca jest wykonywana w zaledwie pół tuzina linii łatwo testowanego kodu.

Najpierw musisz obliczyć odległości. Ponieważ dane są we współrzędnych geograficznych, tutaj jest funkcja do obliczania odległości na kulistym układzie odniesienia (przy użyciu wzoru Haversine):

#
# 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)))
}

Jeśli chcesz, zamień to na swoją ulubioną implementację (na przykład wykorzystującą elipsoidalną bazę danych).

Następnie będziemy musieli obliczyć odległości między każdym „punktem bazowym” (sprawdzanym pod kątem statyczności) a jego czasowym sąsiedztwem. To po prostu kwestia zastosowania distdo okolicy:

#
# 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, ...))

Po trzecie - jest to kluczowa idea - punkty stacjonarne można znaleźć, wykrywając sąsiedztwa 11 punktów mających co najmniej pięć z rzędu, których odległości są wystarczająco małe. Zaimplementujmy to nieco bardziej ogólnie, określając długość najdłuższego podsekwencji prawdziwych wartości w logicznej tablicy wartości boolowskich:

#
# 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)))

(Znajdujemy lokalizacje fałszywych wartości w kolejności i obliczamy ich różnice: są to długości podsekwencji wartości niefałszowych. Zwracana jest największa taka długość.)

Po czwarte, stosujemy się max.subsequencedo wykrywania punktów stacjonarnych.

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

To są wszystkie potrzebne narzędzia.

Jako przykład, stwórzmy interesujące dane z kilkoma skupiskami punktów stacjonarnych. Pójdę przypadkowym spacerem w pobliżu równika.

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)

Tablice loni latzawierają współrzędne npunktów w stopniach . Zastosowanie naszych narzędzi jest proste po pierwszym przekształceniu w radiany:

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

Argument p[max(1,i-5):min(n,i+5), ]mówi, aby spojrzeć wstecz na 5 kroków czasowych lub nawet na 5 kroków czasowych od punktu bazowego p[i,]. W tym k=5mówi, aby szukać sekwencji 5 lub więcej w rzędzie, które znajdują się w odległości 100 km od punktu bazowego. (Wartość 100 km została ustawiona jako domyślna w, is.stationaryale można ją tutaj zastąpić.)

Wynik p.stationaryjest logicznym wektorem wskazującym na stacjonarność: mamy to, po co przyszliśmy. Jednak, aby sprawdzić procedurę, najlepiej wykreślić dane i te wyniki, zamiast sprawdzać tablice wartości. Na poniższej działce pokazuję trasę i punkty. Co dziesiąty punkt jest oznaczony, dzięki czemu można oszacować, ile z nich może zachodzić na siebie w obrębie stacjonarnych grup. Punkty stacjonarne są przerysowane na stałe na czerwono, aby je podświetlić, i otoczone 100-kilometrowymi buforami.

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")

W przypadku innych (opartych na statystyce) podejść do znajdowania punktów stacjonarnych w śledzonych danych, w tym działającego kodu, odwiedź /mathematica/2711/clustering-of-space-time-data .