# R/satellite.model.R In tripEstimation: Metropolis sampler and supporting functions for estimating animal movement from archival tags and satellite fixes

```"satellite.model" <-
function(day, X,
proposal.x,proposal.z,
fix.release=TRUE,fix.recapture=TRUE,
start.x,start.z,
posn.sigma=1,
behav.dist="gamma", behav.mean, behav.sd,
proj.string = "+proj=longlat", distfunction = NULL,
coordtransform = NULL

) {

if (length(day) != nrow(X)) stop("length of times should match number of locations")

##
## Utility functions
##

if (!is.null(distfunction)) {
dist <- distfunction
} else {
## (Great Circle) Distance in kilometres
dist <- function(a,b) {
r <- cos(pi/180*a[,2])*cos(pi/180*b[,2])*
cos(pi/180*(b[,1]-a[,1]))+sin(pi/180*a[,2])*sin(pi/180*b[,2])
6378.137*acos(pmin.int(r, 1))
}
}
if (!is.null(coordtransform)) { transf <- coordtransform } else {
# no transformation
transf <- function(x) x
}
if (!any(grep("longlat", proj.string))) {
##    require(rgdal) removed by namespace addition 2011-10-09 MDS
stop("projections not supported, to do this provide reprojection and distance functions via coordstransform and distfunction")
##    if (!any(grep("km", proj.string))) {
##      warning("Distances will be calculated in the units of the coordinate system, but kilometres are assumed for speed.")
##      print(CRSargs(CRS(proj.string)))
##    }
##    dist <- function(a, b) {
##      sqrt(rowSums((a - b)^2)) }

##    transf <-  function(x, inv = FALSE) {
##      rubbish <- capture.output(res <- project(x, proj.string, inv))
##      res
##    }
}
##
## Positional contribution to the log posterior
##
logp.position <- function(x) {
x[,1:2] <- transf(x[,1:2], inv = TRUE)
rowSums(dnorm(X, x, posn.sigma, log=TRUE))
}

##
## Behavioural contribution to the log posterior
##

## Times between observations
dt <- diff(unclass(day)/3600)

if(behav.dist=="gamma") {
alpha <- behav.mean^2/behav.sd^2
beta <- behav.mean/behav.sd^2
logp.behavioural <- function(k,x1,z,x2) {
## Average speed from x1 to z to x2
spd <- pmax.int(dist(x1,z)+dist(z,x2), 1e-06)/dt[k]
dgamma(spd, alpha, beta, log = TRUE)
}
} else {
log.sigma <- sqrt(log(1+behav.mean^2/behav.sd^2))
log.mu <- log(behav.mean)-log.sigma^2/2
logp.behavioural <- function(k, x1, z, x2) {
## Average speed from x1 to z to x2
spd <- (dist(x1, z) + dist(z, x2))/dt[k]
dnorm(spd, log.mu, log.sigma, log = TRUE)
}
}

##
## Locations to be held fixed
##
n <- length(day)
fixed.x <- rep(FALSE,n)
fixed.x[1] <- fix.release
fixed.x[n] <- fix.recapture

## Return a list with all model components
list(## Number of locations
n = n,
## The function for generating proposal x's
proposal.x=proposal.x,
## The function for generating proposal z's
proposal.z=proposal.z,
## The mask for the x's
## The mask for the z's
## Positional contribution to the log posterior
logp.position=logp.position,
## Behavioural contribution to the log posterior
logp.behavioural=logp.behavioural,
## Locations to be held fixed
fixed.x=fixed.x,
## Suggested starting points
start.x=start.x,
start.z=start.z)
}
```

## Try the tripEstimation package in your browser

Any scripts or data that you put into this service are public.

tripEstimation documentation built on May 2, 2019, 4:59 p.m.