Proximidad en espacio y tiempo.

Tengo algunos datos de puntos que representan las ubicaciones diarias de un animal, con una marca de tiempo asociada.

Me gustaría identificar todos los puntos donde STATIONARY = TRUE. Un punto califica como estacionario si un búfer de 100 km a su alrededor se superpone (digamos) 5 puntos temporalmente adyacentes adicionales . Entonces, si el día 10 es mi punto de interés, quiero preguntar si 5 días temporalmente adyacentes están dentro de un búfer de 100 km de este punto. Si los días 5,6,7,8 y 9; O días 11,12,13,14 y 15; O los días 8, 9, 11, 12, 13 (etc.) están dentro del búfer, luego ESTACIONARIO = VERDADERO. Sin embargo, si los días 5,7,9,11 y 13 están dentro del búfer, pero no los días alternos (pares) intermedios, entonces ESTACIONARIO = FALSO

Creo que algún tipo de búfer de ventana móvil proporcionará la solución, pero no sé cómo implementar esto.

He estado tratando de entender este problema tanto en ArcGIS como en R, pero hasta ahora no he tenido ondas cerebrales. Esto es lo más cercano que tengo a una solución, pero no encaja, no creo: identificación de puntos consecutivos dentro de un búfer especificado

Aquí hay algunos datos ficticios, que se aproximan a mi estructura de datos (aunque en realidad tengo dos ubicaciones diarias (mediodía y medianoche) con algunas ubicaciones faltantes, pero me preocuparé de eso más 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 a dividir esto en partes simples. Al hacerlo, todo el trabajo se realiza en solo media docena de líneas de código fácilmente probado.

Primero, necesitará calcular distancias. Debido a que los datos están en coordenadas geográficas, aquí hay una función para calcular distancias en un dato esférico (usando la fórmula de 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)))
}

Reemplace esto con su implementación favorita si lo desea (como una que usa un dato elipsoidal).

A continuación, necesitaremos calcular las distancias entre cada "punto base" (que se verifica para determinar la estaionaridad) y su vecindad temporal. Eso es simplemente una cuestión de aplicar distal vecindario:

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

Tercero, esta es la idea clave: los puntos estacionarios se encuentran al detectar vecindarios de 11 puntos que tienen al menos cinco en una fila cuyas distancias son lo suficientemente pequeñas. Implementemos esto un poco más generalmente determinando la longitud de la subsecuencia más larga de valores verdaderos dentro de una 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)))

(Encontramos las ubicaciones de los valores falsos , en orden, y calculamos sus diferencias: estas son las longitudes de subsecuencias de valores no falsos. Se devuelve la longitud más grande).

Cuarto, aplicamos max.subsequencepara detectar puntos estacionarios.

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

Esas son todas las herramientas que necesitamos.

Como ejemplo, creemos algunos datos interesantes que tengan algunos grupos de puntos estacionarios. Haré una caminata aleatoria cerca del ecuador.

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)

Las matrices lony latcontienen las coordenadas, en grados, de npuntos en secuencia. Aplicar nuestras herramientas es sencillo después de la primera conversión a radianes:

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

El argumento p[max(1,i-5):min(n,i+5), ]dice mirar hacia atrás hasta 5 pasos de tiempo o hacia adelante hasta 5 pasos de tiempo desde el punto base p[i,]. Incluyendo k=5dice buscar una secuencia de 5 o más en una fila que estén dentro de los 100 km del punto base. (El valor de 100 km se configuró como predeterminado, is.stationarypero puede anularlo aquí).

La salida p.stationaryes un vector lógico que indica estacionariedad: tenemos lo que buscamos. Sin embargo, para verificar el procedimiento, es mejor trazar los datos y estos resultados en lugar de inspeccionar los conjuntos de valores. En la siguiente gráfica muestro la ruta y los puntos. Cada décimo punto está etiquetado para que pueda estimar cuántos podrían superponerse dentro de los grupos estacionarios. Los puntos estacionarios se vuelven a dibujar en rojo sólido para resaltarlos y rodeados por sus zonas de amortiguamiento 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 otros enfoques (basados ​​en estadísticas) para encontrar puntos estacionarios en datos rastreados, incluido el código de trabajo, visite /mathematica/2711/clustering-of-space-time-data .