10

I have some point data which represent daily lat-lon locations of an animal, with an associated timestamp.

I would like to identify all points where STATIONARY = TRUE. A point qualifies as being stationary if a 100km buffer around it overlaps an additional (say) 5 temporally adjacent points. So if day 10 is my point of interest, I want to ask whether 5 temporally adjacent days are within a 100km buffer of this point. If days 5,6,7,8 & 9; OR days 11,12,13,14 & 15; OR days 8,9,11,12,13 (etc.) are within the buffer, then STATIONARY = TRUE. If, however, days 5,7,9,11 & 13 are within the buffer, but not the alternate (even) days in between, then STATIONARY = FALSE

I think some kind of moving window buffer will provide the solution, but I don't know how to implement this.

I have been trying to get my head around this problem in both ArcGIS and R, but have had no brain waves so far. This is the closest I've got to a solution, but it doesn't quite fit, I don't think: Identification of consecutive points within a specified buffer

Here is some dummy data, which approximates my data structure (though in reality I have twice daily locations (midday and midnight) with some locations missing - but I'll worry about that later)

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
... ...         ...       ...
  • 1
    Question? Assuming all 10 points are within the buffer and you have a date separation (starting from day 1) of 1-3-4-12-13-20-21-22-29-30 then are you saying you are only interested selecting points that are in days 1,2,3,4 & 12? – Hornbydd Sep 16 '13 at 21:37
  • No, I'd only be interested in days 1-4. If the animal 'leaves' the buffer then returns on day 12 (or day 6), then that would 'cancel' that stationary period - i.e. the animal has to be in the buffer on day 1-2-3-4-5 for the point at the center of the buffer to be counted. Make sense? I'm not sure myself.. – Tom Finch Sep 16 '13 at 21:52
  • 1
    Just to check, if the point of interest was day 7 then you would be interested points that fall within 100Km for days 7,8,9,10 & 11? – Hornbydd Sep 16 '13 at 22:01
  • Point 7 would be selected as a stationary point if days 8,9,10, 11 & 12 were withing 100km. Or days 5,6,8,9,10. So any one point gets selected if any other 5 temporally adjacent points (the 5 previous days, the 5 subsequent days, or a few of days either side) are all within the buffer. I think the moving window is the best way of conceptualising it. For each 'focal' point any point more than 5 days into the past/future can be forgotten. I may update my original question, as I now understand it a bit more... – Tom Finch Sep 16 '13 at 22:11
  • What's the format of the data? For example, do you have each time/location as a vector point in a shapefile and an attribute table that stores the time? Or are each time/location stored separately in different shapefiles? Are the data not even in a geospatial format and simply in an Excel file? Knowing this would help us answer. – user21951 Sep 16 '13 at 22:27
12

Let's break this down into simple pieces. By doing so, all the work is accomplished in just a half dozen lines of easily tested code.

First, you will need to compute distances. Because the data are in geographic coordinates, here is a function to compute distances on a spherical datum (using the Haversine formula):

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

Replace this with your favorite implementation if you wish (such as one using an ellipsoidal datum).

Next, we will need to compute distances between each "base point" (being checked for staionarity) and its temporal neighborhood. That's simply a matter of applying dist to the neighborhood:

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

Third--this is the key idea--stationary points are found by detecting neighborhoods of 11 points having at least five in a row whose distances are sufficiently small. Let us implement this a little more generally by determining the length of the longest subsequence of true values within a logical array of boolean values:

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

(We find the locations of the false values, in order, and compute their differences: these are the lengths of subsequences of non-false values. The largest such length is returned.)

Fourth, we apply max.subsequence to detect stationary points.

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

Those are all the tools we need.


As an example, let's create some interesting data having a few clumps of stationary points. I will take a random walk near the Equator.

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)

The arrays lon and lat contain the coordinates, in degrees, of n points in sequence. Applying our tools is straightforward after first converting into 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))

The argument p[max(1,i-5):min(n,i+5), ] says to look as far back as 5 time steps or as far forward as 5 time steps from the base point p[i,]. Including k=5 says to look for a sequence of 5 or more in a row that are within 100 km of the base point. (The value of 100 km was set as the default in is.stationary but you could override it here.)

The output p.stationary is a logical vector indicating stationarity: we have what we came for. However, to check the procedure it's best to plot the data and these results rather than inspecting arrays of values. On the following plot I show the route and the points. Every tenth point is labeled so you can estimate how many might overlap within the stationary clumps. Stationary points are redrawn in solid red to highlight them and surrounded by their 100-km buffers.

Figure

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

For other (statistically-based) approaches to finding stationary points in tracked data, including working code, please visit https://mathematica.stackexchange.com/questions/2711/clustering-of-space-time-data.

  • wow, thanks! looking forward to getting my head around this. thanks again for your time and effort – Tom Finch Sep 18 '13 at 9:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.