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 .