Nothing
##' Filter satellite locations
##'
##' Apply a three stage algorithm to eliminate erroneous locations, based
##' on established procedures.
##'
##' These functions implement the location filtering procedure outlined in
##' Austin et al. (2003). \code{grpSpeedFilter} and \code{rmsDistFilter}
##' can be used to perform only the first stage or the second and third
##' stages of the algorithm on their own, respectively. Alternatively, the
##' three filters can be run in a single call using \code{austFilter}.
##'
##' The first stage of the filter is an iterative process which tests every
##' point, except the first and last (\var{w}/2) - 1 (where \var{w} is the
##' window size) points, for travel velocity relative to the
##' preceeding/following (\var{w}/2) - 1 points. If all \var{w} - 1 speeds
##' are greater than the specified threshold, the point is marked as
##' failing the first stage. In this case, the next point is tested,
##' removing the failing point from the set of test points.
##'
##' The second stage runs McConnell et al. (1992) algorithm, which tests
##' all the points that passed the first stage, in the same manner as
##' above. The root mean square of all \var{w} - 1 speeds is calculated,
##' and if it is greater than the specified threshold, the point is marked
##' as failing the second stage (see Warning section below).
##'
##' The third stage is run simultaneously with the second stage, but if the
##' mean distance of all \var{w} - 1 pairs of points is greater than the
##' specified threshold, then the point is marked as failing the third
##' stage.
##'
##' The speed and distance threshold should be obtained separately (see
##' \code{\link{distSpeed}}).
##'
##' @param time \code{POSIXct} object with dates and times for each point.
##' @param lon numeric vectors of longitudes, in decimal degrees.
##' @param lat numeric vector of latitudes, in decimal degrees.
##' @param id A factor grouping points in different categories
##' (e.g. individuals).
##' @param speed.thr numeric scalar: speed threshold (m/s) above which
##' filter tests should fail any given point.
##' @param dist.thr numeric scalar: distance threshold (km) above which the
##' last filter test should fail any given point.
##' @param window integer: the size of the moving window over which tests
##' should be carried out.
##' @param ... Arguments ultimately passed to \code{\link{distSpeed}}.
##' @section Warning:
##' This function applies McConnell et al.'s filter as described in Freitas
##' et al. (2008). According to the original description of the algorithm
##' in McConnell et al. (1992), the filter makes a single pass through all
##' locations. Austin et al. (2003) and other authors may have used the
##' filter this way. However, as Freitas et al. (2008) noted, this causes
##' locations adjacent to those flagged as failing to fail also, thereby
##' rejecting too many locations. In diveMove, the algorithm was modified
##' to reject only the \dQuote{peaks} in each series of consecutive
##' locations having root mean square speed higher than threshold.
##' @return
##' \code{rmsDistFilter} and \code{austFilter} return a matrix with 2 or 3
##' columns, respectively, of logical vectors with values TRUE for points
##' that passed each stage. For the latter, positions that fail the first
##' stage fail the other stages too. The second and third columns returned
##' by \code{austFilter}, as well as those returned by \code{rmsDistFilter}
##' are independent of one another; i.e. positions that fail stage 2 do not
##' necessarily fail stage 3.
##' @author Sebastian Luque \email{spluque@@gmail.com} and Andy Liaw.
##' @references
##' McConnell BJ, Chambers C, Fedak MA. 1992. Foraging ecology of southern
##' elephant seals in relation to bathymetry and productivity of the
##' Southern Ocean. \emph{Antarctic Science} 4:393-398.
##'
##' Austin D, McMillan JI, Bowen D. 2003. A three-stage algorithm for
##' filtering erroneous Argos satellite locations. \emph{Marine Mammal
##' Science} 19: 371-383.
##'
##' Freitas C, Lydersen, C, Fedak MA, Kovacs KM. 2008. A simple new
##' algorithm to filter marine mammal ARGOS locations. Marine Mammal
##' Science DOI: 10.1111/j.1748-7692.2007.00180.x
##' @keywords manip iteration
##' @seealso \code{\link{distSpeed}}
##' @examples
##' ## Using the Example from '?readLocs':
##' utils::example("readLocs", package="diveMove",
##' ask=FALSE, echo=FALSE)
##' ringy <- subset(locs, id == "ringy" & !is.na(lon) & !is.na(lat))
##'
##' ## Examples below use default Meeus algorithm for computing distances.
##' ## See ?distSpeed for specifying other methods.
##' ## Austin et al.'s group filter alone
##' grp <- grpSpeedFilter(ringy[, 3:5], speed.thr=1.1)
##'
##' ## McConnell et al.'s filter (root mean square test), and distance test
##' ## alone
##' rms <- rmsDistFilter(ringy[, 3:5], speed.thr=1.1, dist.thr=300)
##'
##' ## Show resulting tracks
##' n <- nrow(ringy)
##' plot.nofilter <- function(main) {
##' plot(lat ~ lon, ringy, type="n", main=main)
##' with(ringy, segments(lon[-n], lat[-n], lon[-1], lat[-1]))
##' }
##' layout(matrix(1:4, ncol=2, byrow=TRUE))
##' plot.nofilter(main="Unfiltered Track")
##' plot.nofilter(main="Group Filter")
##' n1 <- length(which(grp))
##' with(ringy[grp, ], segments(lon[-n1], lat[-n1], lon[-1], lat[-1],
##' col="blue"))
##' plot.nofilter(main="Root Mean Square Filter")
##' n2 <- length(which(rms[, 1]))
##' with(ringy[rms[, 1], ], segments(lon[-n2], lat[-n2], lon[-1], lat[-1],
##' col="red"))
##' plot.nofilter(main="Distance Filter")
##' n3 <- length(which(rms[, 2]))
##' with(ringy[rms[, 2], ], segments(lon[-n3], lat[-n3], lon[-1], lat[-1],
##' col="green"))
##'
##' ## All three tests (Austin et al. procedure)
##' austin <- with(ringy, austFilter(time, lon, lat, speed.thr=1.1,
##' dist.thr=300))
##' layout(matrix(1:4, ncol=2, byrow=TRUE))
##' plot.nofilter(main="Unfiltered Track")
##' plot.nofilter(main="Stage 1")
##' n1 <- length(which(austin[, 1]))
##' with(ringy[austin[, 1], ], segments(lon[-n1], lat[-n1], lon[-1], lat[-1],
##' col="blue"))
##' plot.nofilter(main="Stage 2")
##' n2 <- length(which(austin[, 2]))
##' with(ringy[austin[, 2], ], segments(lon[-n2], lat[-n2], lon[-1], lat[-1],
##' col="red"))
##' plot.nofilter(main="Stage 3")
##' n3 <- length(which(austin[, 3]))
##' with(ringy[austin[, 3], ], segments(lon[-n3], lat[-n3], lon[-1], lat[-1],
##' col="green"))
"austFilter" <- function(time, lon, lat, id=gl(1, 1, length(time)),
speed.thr, dist.thr, window=5, ...)
{
## FIRST STAGE ********************************************************
locs <- data.frame(time, lon, lat)
## Do first stage over each seal's data, returns vector as long as locs
first <- unlist(by(locs, id, grpSpeedFilter, speed.thr, window, ...),
use.names=FALSE)
## SECOND AND THIRD STAGES ********************************************
good <- which(first) # native subscripts that passed
last <- do.call(rbind, by(locs[good, ], id[good], rmsDistFilter,
speed.thr, window, dist.thr, ...))
filter123 <- cbind(firstPass=first,
secondPass=first, # 2nd and 3rd start the same as 1st
thirdPass=first)
filter123[good, 2:3] <- last
filter123
}
##' @describeIn austFilter Do stage one on 3-column matrix \code{x}
##' @param x 3-column matrix with column 1: \code{POSIXct} vector; column
##' 2: numeric longitude vector; column 3: numeric latitude vector.
##' @return \code{grpSpeedFilter} logical vector indicating those lines
##' that passed the test.
"grpSpeedFilter" <- function(x, speed.thr, window=5, ...)
{
## Value: Do stage one on matrix x (assuming it's a single unit),
## return a logical; whether each observation in x passed the test
## --------------------------------------------------------------------
## Arguments: x=matrix with cols: POSIXct, lon, lat; speed.thr=speed
## threshold (m/s), window=size of window to test; ...=arguments passed
## to distSpeed(), namely only 'method' for now.
## --------------------------------------------------------------------
## Author: Sebastian Luque
## --------------------------------------------------------------------
if (!window %% 2) stop ("window size must be an odd number")
if (nrow(x) < window) stop ("there are fewer rows than window size")
tpos <- window %/% 2 # test subscript - 1
testfun <- function(k) { # k=matrix with group to test
mid <- tpos + 1 # subscript of pt to test
ref <- c(-seq(tpos), seq(tpos)) # subscripts of pts to test against
mid.rep <- rep(mid, length(ref))
speeds <- distSpeed(k[mid.rep, ], k[mid + ref, ], ...)[, 3]
all(speeds > speed.thr, na.rm=TRUE) # TRUE if all speeds > thr
}
pass <- !logical(nrow(x)) # all pass at start up
## define all test rows and subscript for forward movement of window
testrows <- seq(1 + tpos, nrow(x) - tpos); i <- 1
for (j in testrows) {
test <- testfun(x[c(i:(i + tpos - 1), j:(j + tpos)), ])
if (test) {
pass[j] <- FALSE
} else {
i <- i + 1
}
}
pass
}
##' @describeIn austFilter Apply McConnell et al's filter and Austin et
##' al's last stage
"rmsDistFilter" <- function(x, speed.thr, window=5, dist.thr, ...)
{
## Value: Run McConnell et al's filter and Austin et al's last stage,
## return 2-col matrix of logicals; whether each observation passed
## each test.
## --------------------------------------------------------------------
## Arguments: x=matrix with cols: POSIXct, lon, lat; speed.thr=speed
## threshold (m/s), window=size of window to test; dist.thr=distance
## threshold (km); ...=arguments passed to distSpeed(), namely only
## 'method' for now.
## --------------------------------------------------------------------
## Author: Sebastian Luque
## --------------------------------------------------------------------
if (!window %% 2) stop ("window size must be an odd number")
wdw.errmess <- "there are fewer rows than window size"
if (nrow(x) < window) stop (wdw.errmess)
tpos <- window %/% 2 # test subscript - 1
ref <- c(-seq(tpos), seq(tpos)) # reference points for test
travel.fun <- function(k) { # k=subscripts
xmid <- k[1] # 1st is the middle
xmid.rep <- rep(xmid, length(k) - 1)
others <- k[-1] # 1st against all other positions
tr <- distSpeed(x[xmid.rep, ], x[others, ], ...)
tr[, c(1, 3)]
}
testrows <- seq(tpos + 1, nrow(x) - tpos) # subscripts of locs to test
rmsSwitch <- distPass <- !logical(length(testrows)) # to begin looping
while (any(rmsSwitch)) { # stop when switch is all FALSE
switchidx <- which(rmsSwitch)
if ((length(switchidx) + (tpos * 2)) < window) stop (wdw.errmess)
testrows.new <- testrows[rmsSwitch]
idx <- seq_along(testrows.new) # index the above
testidx.mtx <- testidx <- c(idx, sapply(ref, "+", idx))
testidx.inner <- testidx >= 1 & testidx <= length(idx)
testidx.mtx[testidx.inner] <- testrows.new[testidx[testidx.inner]]
testidx.low <- testidx < 1
testidx.mtx[testidx.low] <- testrows[1] - (1 - testidx[testidx.low])
testidx.high <- testidx > length(idx)
testidx.mtx[testidx.high] <- testrows[length(testrows)] +
(testidx[testidx.high] - length(idx))
testidx.mtx <- matrix(testidx.mtx, nrow=length(idx))
travel <- apply(testidx.mtx, 1, travel.fun)
if (dim(travel)[1] > 2) {
dist.refs <- seq(length(ref))
speed.refs <- seq(length(ref) + 1, nrow(travel))
} else {
dist.refs <- 1
speed.refs <- 2
}
dists <- as.matrix(travel[dist.refs, ])
speeds <- as.matrix(travel[speed.refs, ])
## root mean square value
rms <- apply(speeds, 2, function(k) { # do this for every test group
sqrt(sum(k ^ 2, na.rm=TRUE) / length(k))
})
## prelim passing locs in current set
rmsPass <- rms <= speed.thr
## find series of adj prelim failing locs in current set
rmsPass.rle <- rle(rmsPass)
bad <- which(!rmsPass.rle$values & rmsPass.rle$lengths > 1)
if (length(bad) < 1) break # stop looping if no adjc failing locs
beg <- rev(length(rmsPass) + 1 - cumsum(rev(rmsPass.rle$lengths)))
end <- cumsum(rmsPass.rle$lengths)
for (j in bad) {
bads <- seq(beg[j], end[j])
Vi <- rms[bads]
maxVi.idx <- which(Vi == max(Vi))
peaks <- bads[maxVi.idx]
ok <- bads[!bads %in% peaks]
rmsPass[ok] <- TRUE
}
rmsSwitch[switchidx[!rmsPass]] <- FALSE # set failing locs in orig set
## Distance filter
distt <- apply(dists, 2, function(k) sum(k, na.rm=TRUE) / length(k))
distFail <- distt > dist.thr
distPass[switchidx[distFail]] <- FALSE
}
rmsPass <- rmsSwitch
## Top and bottom ends pass as these can't be tested
untested <- matrix(!logical(tpos * 2), ncol=2)
rbind(untested, cbind(rmsPass, distPass), untested)
}
## TEST ZONE --------------------------------------------------------------
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.