Proximidade no espaço e no tempo

Eu tenho alguns dados pontuais que representam locais lat-lon diários de um animal, com um timestamp associado.

Gostaria de identificar todos os pontos em que STATIONARY = TRUE. Um ponto se qualifica como estacionário se um buffer de 100 km ao seu redor se sobrepuser a (digamos) 5 pontos temporariamente adjacentes adicionais . Portanto, se o dia 10 for o meu ponto de interesse, quero perguntar se 5 dias temporalmente adjacentes estão dentro de um buffer de 100 km desse ponto. Se dias 5,6,7,8 e 9; Dias OR 11,12,13,14 e 15; Os dias OR 8,9,11,12,13 (etc.) estão dentro do buffer e, em seguida, STATIONARY = TRUE. Se, no entanto, os dias 5,7,9,11 e 13 estiverem dentro do buffer, mas não os dias alternativos (pares) intermediários, então STATIONARY = FALSE

Acho que algum tipo de buffer de janela em movimento fornecerá a solução, mas não sei como implementar isso.

Eu tenho tentado resolver esse problema no ArcGIS e no R, mas até agora não tive ondas cerebrais. Este é o mais próximo que eu tenho de uma solução, mas não se encaixa, não acho: Identificação de pontos consecutivos dentro de um buffer especificado

Aqui estão alguns dados fictícios, que se aproximam da minha estrutura de dados (embora na realidade eu tenha locais duas vezes ao dia (meio-dia e meia-noite) com alguns locais ausentes - mas vou me preocupar com isso mais tarde)

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

Vamos dividir isso em pedaços simples. Ao fazer isso, todo o trabalho é realizado em apenas meia dúzia de linhas de código facilmente testado.

Primeiro, você precisará calcular as distâncias. Como os dados estão em coordenadas geográficas, aqui está uma função para calcular distâncias em um dado esférico (usando a fórmula 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)))
}

Substitua por sua implementação favorita, se desejar (como uma que usa um dado elipsoidal).

Em seguida, precisaremos calcular as distâncias entre cada "ponto base" (que está sendo verificado quanto à estagnação) e sua vizinhança temporal. Isso é simplesmente uma questão de se aplicar distao bairro:

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

Terceiro - essa é a idéia principal - pontos estacionários são encontrados detectando bairros de 11 pontos com pelo menos cinco em uma linha cujas distâncias são suficientemente pequenas. Vamos implementar isso um pouco mais geralmente, determinando o comprimento da subsequência mais longa de valores verdadeiros dentro de uma matriz lógica de valores booleanos:

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

(Nós encontramos os locais dos valores falsos , em ordem, e calculamos suas diferenças: estes são os comprimentos das subsequências de valores não falsos. O maior tamanho é retornado.)

Quarto, aplicamos max.subsequencepara detectar pontos estacionários.

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

Essas são todas as ferramentas que precisamos.

Como exemplo, vamos criar alguns dados interessantes com alguns grupos de pontos estacionários. Vou dar um passeio aleatório perto do Equador.

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)

As matrizes lone latcontêm as coordenadas, em graus, dos npontos em sequência. A aplicação de nossas ferramentas é simples após a primeira conversão em radianos:

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

O argumento p[max(1,i-5):min(n,i+5), ]diz para olhar para trás em 5 etapas de tempo ou até 5 etapas de tempo a partir do ponto base p[i,]. A inclusão k=5diz para procurar uma sequência de 5 ou mais em uma fileira que esteja a 100 km do ponto base. (O valor de 100 km foi definido como padrão, is.stationarymas você pode substituí-lo aqui.)

A saída p.stationaryé um vetor lógico indicando estacionariedade: nós temos o que buscamos. No entanto, para verificar o procedimento, é melhor plotar os dados e esses resultados em vez de inspecionar matrizes de valores. No gráfico a seguir, mostro a rota e os pontos. Cada décimo ponto é rotulado para que você possa estimar quantos podem se sobrepor dentro dos grupos estacionários. Os pontos estacionários são redesenhados em vermelho sólido para destacá-los e cercados por seus amortecedores de 100 km.

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

Para outras abordagens (baseadas em estatística) para encontrar pontos estacionários em dados rastreados, incluindo código de trabalho, visite /mathematica/2711/clustering-of-space-time-data .