La proximité dans l'espace et le temps

J'ai quelques données ponctuelles qui représentent les emplacements lat-lat quotidiens d'un animal, avec un horodatage associé.

Je voudrais identifier tous les points où STATIONARY = TRUE. Un point est considéré comme stationnaire si un tampon de 100 km autour de lui chevauche 5 (disons) 5 points adjacents dans le temps supplémentaires . Donc, si le jour 10 est mon point d'intérêt, je veux demander si 5 jours temporellement adjacents sont dans un tampon de 100 km de ce point. Si les jours 5,6,7,8 & 9; OU jours 11, 12, 13, 14 et 15; OU les jours 8,9,11,12,13 (etc.) sont dans le tampon, puis STATIONARY = TRUE. Si, cependant, les jours 5, 7, 9, 11 et 13 se trouvent dans le tampon, mais pas les jours alternés (pairs) entre les deux, STATIONARY = FALSE

Je pense qu'une sorte de tampon de fenêtre mobile fournira la solution, mais je ne sais pas comment l'implémenter.

J'ai essayé de résoudre ce problème dans ArcGIS et R, mais je n'ai pas eu d'ondes cérébrales jusqu'à présent. C'est la solution la plus proche que j'ai trouvée, mais elle ne correspond pas tout à fait, je ne pense pas: Identification de points consécutifs dans un tampon spécifié

Voici quelques données fictives, qui se rapprochent de ma structure de données (bien qu'en réalité, j'ai deux emplacements par jour (midi et minuit) avec certains emplacements manquants - mais je m'inquiéterai à ce sujet plus tard)

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

Décomposons cela en morceaux simples. Ce faisant, tout le travail est accompli en seulement une demi-douzaine de lignes de code facilement testées.

Tout d'abord, vous devrez calculer les distances. Parce que les données sont en coordonnées géographiques, voici une fonction pour calculer les distances sur une donnée sphérique (en utilisant la formule 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)))
}

Remplacez-le par votre implémentation préférée si vous le souhaitez (comme celle qui utilise une donnée ellipsoïdale).

Ensuite, nous devrons calculer les distances entre chaque "point de base" (dont la stabilité est vérifiée) et son voisinage temporel. Il s'agit simplement de postuler distau quartier:

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

Troisièmement - c'est l'idée clé - les points stationnaires sont trouvés en détectant des quartiers de 11 points ayant au moins cinq dans une rangée dont les distances sont suffisamment petites. Implémentons cela un peu plus généralement en déterminant la longueur de la plus longue sous-séquence de valeurs vraies dans un tableau logique de valeurs booléennes:

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

(Nous trouvons les emplacements des fausses valeurs, dans l'ordre, et calculons leurs différences: ce sont les longueurs des sous-séquences de valeurs non fausses. La plus grande de ces longueurs est renvoyée.)

Quatrièmement, nous appliquons max.subsequenceà détecter des points stationnaires.

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

Ce sont tous les outils dont nous avons besoin.

À titre d'exemple, créons des données intéressantes ayant quelques blocs de points stationnaires. Je vais faire une promenade au hasard près de l'équateur.

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)

Les tableaux lonet latcontiennent les coordonnées, en degrés, des npoints en séquence. L'application de nos outils est simple après la première conversion en radians:

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

L'argument p[max(1,i-5):min(n,i+5), ]dit de regarder aussi loin que 5 pas de temps ou aussi loin que 5 pas de temps à partir du point de base p[i,]. Y compris k=5dit de rechercher une séquence de 5 ou plus d'affilée qui sont à moins de 100 km du point de base. (La valeur de 100 km a été définie par défaut dans is.stationarymais vous pouvez la remplacer ici.)

La sortie p.stationaryest un vecteur logique indiquant la stationnarité: nous avons ce pour quoi nous sommes venus. Cependant, pour vérifier la procédure, il est préférable de tracer les données et ces résultats plutôt que d'inspecter des tableaux de valeurs. Sur l'intrigue suivante, je montre l'itinéraire et les points. Chaque dixième point est étiqueté afin que vous puissiez estimer combien peuvent se chevaucher au sein des mottes stationnaires. Les points stationnaires sont redessinés en rouge uni pour les mettre en valeur et entourés de leurs tampons 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")

Pour d'autres approches (basées sur des statistiques) pour trouver des points stationnaires dans les données suivies, y compris le code de travail, veuillez visiter /mathematica/2711/clustering-of-space-time-data .