Близость в пространстве и времени

У меня есть некоторые точечные данные, которые представляют дневные местоположения животного в течение долгого времени, с соответствующей отметкой времени.

Я хотел бы определить все точки, где STATIONARY = TRUE. Точка считается стационарной, если 100-километровый буфер вокруг нее перекрывает дополнительные (скажем, 5) смежные во времени точки. Таким образом, если интересует 10-й день , я хочу спросить, находятся ли 5 ​​соседних во времени дней в 100-километровом буфере этой точки. Если дни 5,6,7,8 и 9; ИЛИ дни 11, 12, 13, 14 и 15; ИЛИ 8,9,11,12,13 (и т. Д.) Находятся в буфере, тогда STATIONARY = TRUE. Однако, если дни 5,7,9,11 и 13 находятся в пределах буфера, но не чередуются (даже) дни между ними, то STATIONARY = FALSE

Я думаю, что некоторый буфер движущегося окна предоставит решение, но я не знаю, как реализовать это.

Я пытался разобраться с этой проблемой в ArcGIS и R, но пока у меня не было мозговых волн. Это самое близкое мне решение, но оно не совсем подходит, я не думаю: идентификация последовательных точек в указанном буфере

Вот некоторые фиктивные данные, которые приблизительно соответствуют моей структуре данных (хотя на самом деле у меня есть два раза в день местоположения (полдень и полночь) с отсутствующими некоторыми местоположениями - но я буду беспокоиться об этом позже)

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

Давайте разберем это на простые части. Таким образом, вся работа выполняется всего за полдюжины строк легко тестируемого кода.

Во-первых, вам нужно будет рассчитать расстояния. Поскольку данные представлены в географических координатах, здесь есть функция для вычисления расстояний на сферическом уровне (с использованием формулы Хаверсайна):

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

Замените это вашей любимой реализацией, если хотите (например, с использованием эллипсоидальных данных).

Далее нам нужно вычислить расстояния между каждой «базовой точкой» (проверяемой на стационарность) и ее временной окрестностью. Это просто вопрос применения distк соседству:

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

В-третьих, это ключевая идея - стационарные точки находятся путем обнаружения окрестностей из 11 точек, имеющих не менее пяти рядов, расстояния которых достаточно малы. Давайте реализуем это немного более широко, определяя длину самой длинной подпоследовательности истинных значений в логическом массиве логических значений:

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

(Мы находим местоположения ложных значений по порядку и вычисляем их различия: это длины подпоследовательностей не ложных значений. Возвращается наибольшая такая длина.)

В-четвертых, мы применяем max.subsequenceдля обнаружения стационарных точек.

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

Это все инструменты, которые нам нужны.

В качестве примера давайте создадим некоторые интересные данные, имеющие несколько скоплений стационарных точек. Я совершу случайную прогулку возле экватора.

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)

Массивы lonи latсодержат координаты, в градусах, nточек в последовательности. Применение наших инструментов просто после первого преобразования в радианы:

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

Аргумент p[max(1,i-5):min(n,i+5), ]говорит, что нужно смотреть назад на 5 временных шагов или на 5 временных шагов от базовой точки p[i,]. Включая k=5говорит, чтобы искать последовательность из 5 или более подряд, которые находятся в пределах 100 км от базовой точки. (Значение 100 км было установлено по умолчанию, is.stationaryно вы можете переопределить его здесь.)

На выходе p.stationaryполучается логический вектор, указывающий на стационарность: у нас есть то, для чего мы пришли. Однако для проверки процедуры лучше составить график данных и этих результатов, а не проверять массивы значений. На следующем участке я показываю маршрут и точки. Каждая десятая точка помечена, чтобы вы могли оценить, сколько из них может перекрываться в стационарных скоплениях. Стационарные точки перерисовываются сплошным красным, чтобы выделить их, и окружены их 100-километровыми буферами.

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

Для других (основанных на статистике) подходов к поиску стационарных точек в отслеживаемых данных, включая рабочий код, посетите страницу /mathematica/2711/clustering-of-space-time-data .