Odległość w przestrzeni i czasie

10

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
... ...         ...       ...
Tom Finch
źródło
1
Pytanie? Zakładając, że wszystkie 10 punktów znajduje się w buforze, a przedział dat (od 1 dnia) wynosi 1-3-4-12-13-20-20-21-22-29-30, to czy mówisz, że jesteś zainteresowany jedynie wyborem punktów które są w dniach 1,2,3,4 i 12?
Hornbydd,
Nie, interesują mnie tylko dni 1-4. Jeśli zwierzę „opuści” bufor, a następnie powróci w dniu 12 (lub dniu 6), wówczas „anuluje” ten okres stacjonarny - tj. Zwierzę musi znajdować się w buforze w dniach 1-2-3-4-5 przez punkt na środku bufora, który należy policzyć. Ma sens? Sam nie jestem pewien ...
Tom Finch,
1
Wystarczy sprawdzić, czy punktem zainteresowania był dzień 7, to czy interesują Cię punkty, które mieszczą się w odległości 100 km od dni 7,8,9,10 i 11?
Hornbydd,
Punkt 7 zostałby wybrany jako punkt stacjonarny, gdyby dni 8,9,10, 11 i 12 trwały 100 km. Lub dni 5,6,8,9,10. Tak więc wybierany jest dowolny punkt, jeśli wszystkie pozostałe 5 tymczasowo sąsiadujących punktów (5 poprzednich dni, 5 kolejnych dni lub kilka dni po obu stronach) znajdują się w buforze. Myślę, że ruchome okno jest najlepszym sposobem na jego konceptualizację. Dla każdego punktu centralnego można zapomnieć o dowolnym punkcie dłuższym niż 5 dni w przeszłość / przyszłość. Mogę zaktualizować moje pierwotne pytanie, ponieważ teraz rozumiem je nieco bardziej ...
Tom Finch,
Jaki jest format danych? Na przykład, czy masz za każdym razem / lokalizację jako punkt wektorowy w pliku kształtu i tabelę atrybutów, która przechowuje czas? A może za każdym razem / miejsce jest przechowywane osobno w różnych plikach kształtu? Czy dane nie są nawet w formacie geoprzestrzennym i po prostu w pliku Excel? Wiedza na ten temat pomoże nam odpowiedzieć.

Odpowiedzi:

12

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.

Postać

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 .

Whuber
źródło
wow, dzięki! nie mogę się doczekać, aż obejdę tę sprawę. jeszcze raz dziękuję za poświęcony czas i wysiłek
Tom Finch,