---
title: "anipaths: spline-based interpolation"
author: "Henry Scharf"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
    toc: true
vignette: >
  %\VignetteIndexEntry{anipaths}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
```{r, include = F}
options(rmarkdown.html_vignette.check_title = FALSE)
knitr::opts_chunk$set(eval = F)
```

## Animating animal trajectories

The package `anipaths` contains a collection of telemetry observations for turkey vultures originally analyzed in:

> Dodge S, Bohrer G, Bildstein K, Davidson SC, Weinzierl R, Mechard MJ, Barber D, Kays R, Brandes D, Han J (2014) Environmental drivers of variability in the movement ecology of turkey vultures (Cathartes aura) in North and South America. Philosophical Transactions of the Royal Society B 20130195.

To animate the locations, we first need to create a time stamp variable of class `numeric` or `POSIX`. One advantage to using the `POSIX` class is that we can specify the gaps in the interpolation (`delta.t`) using convenient character strings like `"hour"` or `"week"`.

```{r, fig.show='hold', fig.cap="test"}
library(anipaths)
vultures$POSIX <- as.POSIXct(vultures$timestamp, tz = "UTC")
vultures_paths <- vultures[format(vultures$POSIX, "%Y") == 2009, ] ## limit attention to 2009
animate_paths(paths = vultures_paths, 
              delta.t = "day",
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier")
```

### Using background maps

Often it may be useful to add a map of the relevant study area. There are lots of ways to incorporate a map. The simplest way is to set `background` to `TRUE`, in which case `anipaths` will do the best it can to select a map based on the data. In the next example, we've changed the time step to help the animations load a little faster.
```{r map_logical}
library(ggmap)
animate_paths(paths = vultures_paths, 
              delta.t = 2 * 24 * 60 * 60, ## number of seconds in two days
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier", 
              background = TRUE, img.name = "background_TRUE")
```

You can also give a long/lat `location`, `zoom` level (3-21; see `?ggmap::get_map()`), and `maptype` (`satellite`, `terrain`, `hybrid`) to be passed to `ggmap::get_map()`, and `anipaths` will make a background for you. 

As shown below, the value of `delta.t` can be specified as a numeric, which will be interpreted in whatever units used by `as.numeric(paths['Time.name'])` (in our case, this is seconds). 
```{r map_google}
vultures_paths <- vultures[format(vultures$POSIX, "%Y") == 2009:2010, ]
background <- list(center = c(-90, 10), 
                   zoom = 3, 
                   maptype = "satellite")
animate_paths(paths = vultures_paths, 
              delta.t = 3 * 24 * 60 * 60, ## number of seconds in three days
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier", 
              background = background, img.name = "background_google")
```

You can also supply your own background image. The projection and units should match the data.
```{r map_simple}
background <- geodata::world(path = ".")
sf::st_crs(background)$proj4string ## matches default projection in animate_paths()
animate_paths(paths = vultures_paths, 
              delta.t = "week",
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier", 
              background = background, img.name = "background_user")
system("rm -r gadm") ## remove geodata map from machine
```

The function `animate_paths()` must project the data to match Google's map tiles, so if your data aren't in long/lat format, make sure you update the `paths.proj` variable with the correct string.

### Adding individual-level covariate information

If you have a covariate of interest for each individual, `animate_paths()` can display that information as a colored ring around each individual. We don't have any natural individual-level covariates available for the vultures, so we'll make one up for the purposed of demonstration. The following code assigns each individual a random interval for each of three behaviors: *exploratory*, *directed*, and *stationary*.

```{r make_covariates}
behaviors <- c("exploratory", "directed", "stationary")
set.seed(1)
vultures_paths$behavior <- 
  unlist(sapply(unique(vultures_paths$individual.local.identifier), function(id){
    v_id <- vultures_paths[vultures_paths$individual.local.identifier == id, ]
    switches <- c(0, sort(sample(1:nrow(v_id), 2)), nrow(v_id))
    rep(behaviors[sample(1:3, 3)], diff(switches))
  }))
```

The covariates have now been appended to the `paths` data frame. We can let `animate_paths()` know we would like it to display this information by setting the argument `covariate` to match the name of the appropriate column in `paths` (i.e., `behavior`). The default colors are a gray scale. Any collection of colors can be provided (e.g., `covariate.colors = viridis::viridis(3)`) which will be turned into a palette to match the support of the covariate. 

```{r plot_covariates}
delta.t <- "day"
background <- geodata::world(path = ".")
sf::st_crs(background)$proj4string ## matches default projection in animate_paths()
animate_paths(paths = vultures_paths, 
              delta.t = delta.t,
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              covariate = "behavior", covariate.colors = viridis::viridis(3),
              ID.name = "individual.local.identifier",
              background = background, img.name = "covariates")
system("rm -r gadm") ## remove geodata map from machine
```

*Note*: if a covariate is represented numerically but should be treated as a factor, changing the class using `as.factor()` should prevent `animate_paths()` from interpolating to "intermediate" values (e.g., if male and female are coded as 0 and 1, defining the column of the data frame as a factor will constrain individuals to only take on values of 0 and 1).

### Output formats

The default format for displaying the animation is an `index.html` file. This is great for broad functionality, and the generated images may be used as the basis of a variety of other animation formats (e.g., the "animate" package in LaTeX).  If you have `ffmpeg` installed on your system, you can also set `method = "mp4"` to create a stand-alone video that is easy to share with others. For more information, see https://www.ffmpeg.org/. 

## Checking the interpolation

As a way to check that `anipaths` is producing reasonable interpolations of the telemetry observations, a generic `plot()` function is provided that takes an argument of class `paths_animation` produced by calling `animate_paths()` with `return.paths = TRUE`. See `example(plot.paths_animation)`. Adjusting the parameters of the interpolation can be done by modifying the `s_args` argument.

The following code demonstrates the use of `s_args` by intentionally over-smoothing the trajectories.
```{r}
interp <- animate_paths(paths = vultures_paths, 
                        delta.t = "day",
                        coord = c("location.long", "location.lat"),
                        Time.name = "POSIX",
                        ID.name = "individual.local.identifier",
                        s_args = rep(list(list(k = 10)), 10),
                        return.paths = T)
plot(interp, i = 2)
```

This piece of code shows the number of observations for each individual vulture, in the same order they are handled by `animate_paths()` (i.e., the same order as `unique(vultures_paths$individual.local.identifier)`). The wide variability in duration and number of observations required individual-specific numbers of knots in the GAM interpolation. In this case, the number of observations divided by 4 is used as the number of knots for each individual, with a maximum of 306 knots because each frame in the animation corresponds to a single day between March 1 and December 31.
```{r}
obs_counts <- merge(data.frame(individual.local.identifier = unique(vultures_paths$individual.local.identifier)), 
                    aggregate(timestamp ~ individual.local.identifier, data = vultures_paths, FUN = length),
                    by = "individual.local.identifier", sort = F)
obs_counts
s_args <- lapply(obs_counts$timestamp, function(x) c(k = floor(min(x / 4, 306))))
interp <- animate_paths(paths = vultures_paths, 
                        delta.t = "day",
                        coord = c("location.long", "location.lat"),
                        Time.name = "POSIX",
                        ID.name = "individual.local.identifier",
                        s_args = s_args,
                        return.paths = T)
plot(interp, i = 2)
```

As of version 0.10.2, the default value for the number of knots also works well for this example.

```{r}
interp <- animate_paths(paths = vultures_paths, 
                        delta.t = "day",
                        coord = c("location.long", "location.lat"),
                        Time.name = "POSIX",
                        ID.name = "individual.local.identifier",
                        return.paths = T, verbose = T)
plot(interp, i = 2)
```

## Additional features

### Showing complete trajectories

Set `whole.path = TRUE`. It is not required, but probably clearer visually, if we also set `tail.length = 0`.

```{r whole_path_anim}
animate_paths(paths = vultures_paths, 
              delta.t = delta.t,
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier", 
              background = background, 
              whole.path = TRUE, tail.length = 0, img.name = "whole_traj")
```

### Dimming selected individuals

```{r dim_anim}
vultures_paths <- vultures[format(vultures$POSIX, "%Y") == 2009 & 
                             vultures$location.lat > 32, ]
animate_paths(paths = vultures_paths, 
              delta.t = delta.t,
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier", 
              background = background, dimmed = c(1, 3, 5), img.name = "dim")
```

### Showing relational information (beta)

It is also possible to visualize dynamic pair-wise relational information with `anipaths`. This information needs to be coded as an array of adjacency matrices (binary and weighted edges allowed, no support yet for directed information). If the network information is on a different time scale than the position information, it will be interpolated using smoothing splines to match.

We don't have existing relationships among the vultures to display, so we'll make some up for the purposes of demonstration. 

```{r make_network}
vultures_paths <- vultures[format(vultures$POSIX, "%Y") == 2009, ]
set.seed(1)
n_indiv <- length(unique(vultures_paths$individual.local.identifier))
change_pts <- 5
network.times <- seq(min(vultures_paths$POSIX), max(vultures_paths$POSIX), l = change_pts + 2)
network <- array(NA, dim = c(n_indiv, n_indiv, length(network.times)))
for(time_i in 2:length(network.times)){
  network_mat <- matrix(sample(1:0, n_indiv^2, prob = c(0.1, 0.9), replace = T), 
                        n_indiv, n_indiv)
  network_mat[lower.tri(network_mat)] <- t(network_mat)[lower.tri(network_mat)]
  diag(network_mat) <- 1
  network[, , (time_i - 1):time_i] <- network_mat
}
```

```{r network_anim}
delta.t <- 3600*24*2
animate_paths(paths = vultures_paths, 
              delta.t = delta.t,
              coord = c("location.long", "location.lat"),
              Time.name = "POSIX",
              ID.name = "individual.local.identifier",
              background = background, network = network, network.times = network.times,
              img.name = "network")
```

```{r remove files, eval = F}
system("rm -r js; rm -r css; rm -r images; rm index.html")
```