---
title: "Trajectory analysis of Turkey vultures"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Trajectory analysis of Turkey vultures}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, fig.width = 7, fig.height = 7,
  comment = "#>"
)
if (Sys.info()["user"] != "bart") {
  if (Sys.getenv("MBPWD") != "") {
    options(keyring_backend = "env")
    move2::movebank_store_credentials("move2_user", Sys.getenv("MBPWD"))
  } else {
    knitr::opts_chunk$set(eval = FALSE)
  }
}
Sys.setlocale("LC_TIME", "en_US.utf8")
```

## Retrieve data

```{r setup}
library(move2)
vulture_data <-
  movebank_download_study("Turkey vultures in North and South America")
vulture_data
```

In this case some tracks have very long time gaps, to prevent distant points to
be connected by a line we split those tracks by creating a new id.

```{r map_of_vulture_track, fig.alt="A map showing the trajectories of vultures using coast lines as a reference"}
library(dplyr, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(rnaturalearth, quietly = TRUE)
library(units, quietly = TRUE)
vulture_lines <- vulture_data %>%
  mutate_track_data(name = individual_local_identifier) %>%
  mutate(
    large_gaps = !(mt_time_lags(.) < set_units(1500, "h") |
      is.na(mt_time_lags(.))),
    track_sub_id = cumsum(lag(large_gaps, default = FALSE)),
    new_track_id = paste(mt_track_id(.), track_sub_id)
  ) %>%
  mt_set_track_id("new_track_id") %>%
  mt_track_lines()
ggplot() +
  geom_sf(data = ne_coastline(returnclass = "sf", 50)) +
  theme_linedraw() +
  geom_sf(
    data = vulture_lines,
    aes(color = name)
  ) +
  coord_sf(
    crs = sf::st_crs("+proj=aeqd +lon_0=-83 +lat_0=8 +units=km"),
    xlim = c(-3500, 3800), ylim = c(-4980, 4900)
  )
```

## Categorize seasons

```{r categorize_seasons}
library(magrittr, quietly = TRUE)
library(lubridate, quietly = TRUE)
vulture_data <- vulture_data %>% mutate(
  month = month(mt_time(), label = TRUE, abbr = FALSE),
  season = recode_factor(month,
    January = "Wintering", February = "Wintering",
    March = "Wintering", April = "North migration",
    May = "North migration", June = "Breeding",
    July = "Breeding", August = "Breeding",
    September = "South migration", October = "South migration",
    November = "Wintering", December = "Wintering"
  ),
  # Here we change season to NA if the next location is either from
  # a different track or season
  season = if_else(season == lead(season, 1) &
    mt_track_id() == lead(mt_track_id(), 1),
  season, NA
  )
)
```

Annotate speed and azimuth to the trajectory.

```{r annotate_track}
vulture_data <- vulture_data %>% mutate(azimuth = mt_azimuth(.), speed = mt_speed(.))
```
## Seasonal distribution per individual
```{r azimuth_distribution, fig.height=6, fig.alt="A plot of the azimuths of the tracks split out per season and individual"}
library(circular, quietly = TRUE)
vulture_azimuth_distributions <- vulture_data %>%
  filter(speed > set_units(2, "m/s"), !is.na(season)) %>%
  group_by(season, track_id = mt_track_id()) %>%
  filter(n() > 50) %>%
  summarise(azimuth_distribution = list(density(
    as.circular(
      drop_units(set_units(
        azimuth,
        "degrees"
      )),
      units = "degrees",
      modulo = "asis",
      zero = 0,
      template = "geographic", rotation = "clock", type = "angles"
    ),
    bw = 180, kernel = "vonmises"
  )))

# Load purrr for map function
library(purrr, quietly = TRUE)
# Load tidy r for unnest function
library(tidyr, quietly = TRUE)
vulture_azimuth_distributions %>%
  mutate(
    x = map(azimuth_distribution, ~ as.numeric(.$x)),
    y = map(azimuth_distribution, ~ .$y)
  ) %>%
  select(-azimuth_distribution) %>%
  unnest(c(x, y)) %>%
  ggplot() +
  geom_path(aes(x = as.numeric(x), y = y, color = season)) +
  coord_polar(start = pi / 2) +
  theme_linedraw() +
  facet_wrap(~ factor(track_id)) +
  scale_x_continuous(
    name = NULL, breaks = (-2:1) * 90,
    labels = c("S", "W", "N", "E")
  ) +
  scale_y_continuous(name = NULL, limits = c(-0.8, 1.4), expand = c(0L, 0L)) +
  labs(color = "Season")
```

## Individual trajectory

```{r speed_direction, fig.height=5, fig.alt="A plot exploring the relation between direction and speed for one track"}
leo <- vulture_data |>
  filter_track_data(individual_local_identifier == "Leo") |>
  mutate(speed_categorical = cut(speed, breaks = c(2, 5, 10, 15, 35)))
leo |> ggplot(aes(x = azimuth, y = speed)) +
  geom_point() +
  scale_x_units(unit = "degrees", breaks = -2:2 * 90, expand = c(0L, 0L)) +
  theme_linedraw()
```
```{r seasonal_speed_distribution, fig.alt="Plot that visualizes the speed and directions per season for one individual"}
leo |>
  filter(speed > set_units(2L, "m/s"), !is.na(season)) |>
  ggplot() +
  coord_polar(start = pi) +
  geom_histogram(
    aes(
      x = set_units(azimuth, "degrees"),
      fill = speed_categorical
    ),
    breaks = set_units(seq(-180L, 180L, by = 10L), "degrees"),
    position = position_stack(reverse = TRUE)
  ) +
  scale_x_units(
    name = NULL,
    limits = set_units(c(-180L, 180L), "degrees"),
    breaks = (-2L:2L) * 90L
  ) +
  facet_wrap(~season) +
  scale_fill_ordinal("Speed") +
  theme_linedraw()
```

Plot turn angle distribution.

```{r turnangle_plot, fig.height=4, fig.width=5, fig.alt="Plot of the turning angle distribution"}
pi_r <- set_units(pi, "rad")
leo %>%
  mutate(turnangle = mt_turnangle(.)) %>%
  filter(speed > set_units(2L, "m/s"), lag(speed, 1L) > set_units(2L, "m/s")) %>%
  ggplot() +
  geom_histogram(
    aes(
      x = turnangle,
      fill = speed_categorical
    ),
    position = position_stack(reverse = TRUE)
  ) +
  scale_fill_ordinal("Speed") +
  coord_polar(start = pi) +
  scale_x_units(limits = c(-pi_r, pi_r), name = NULL) +
  scale_y_continuous(limits = c(-500L, 650L), breaks = c(0L, 250L, 500L)) +
  theme_linedraw()
```

## Net displacement

Here we visualize the distance to the first location of each trajectory. Notice the distances from the first location directly have the correct units associated due to the projection that is included.

```{r nsd, fig.height=5, fig.alt="Plot of the net squared displacement overtime per individual"}
vulture_data <- vulture_data %>%
  group_by(mt_track_id()) %>%
  mutate(displacement = c(st_distance(
    !!!syms(attr(., "sf_column")),
    (!!!syms(attr(., "sf_column")))[row_number() == 1]
  )))

vulture_data %>% ggplot() +
  geom_line(aes(
    x = timestamp,
    y = set_units(displacement, "km"),
    color = individual_local_identifier
  )) +
  ylab("Distance from start") +
  theme_linedraw()
```

Using `dplyr` it is also possible to do more complex operations per individual (or other grouping variables for that matter). 
Here we do a spatial intersection for each time an individual crosses a circular buffer a 1000 kilometers away from its initial location. 
```{r buffer}
vulture_data %>%
  # Group by individual/track_iu
  group_by(mt_track_id()) %>%
  # Create a list with each individual being a element of the list
  group_split() %>%
  # Use purrr::map to do some operation on each track
  purrr::map(~ {
    # Create a circular buffer around the first location
    buffer <- sf::st_cast(sf::st_buffer(
      sf::st_geometry(head(.x, 1)),
      units::set_units(1000, "km")
    ), "LINESTRING")
    # Do a spatial interpolation to the locations the track crosses the buffer
    mt_interpolate(
      .x,
      sf = buffer, omit = TRUE
    )
  }) |>
  # Use purrr compact to omit individuals that never reach a 1000 kilometers from the first location and thus do not have an crossing
  purrr::compact() |>
  # Combine back into one move2
  mt_stack()
```

Using `dplyr::group_map()` we can calculate statistics for groups, such as here a yearly convex hull for each individual. 

```{r convex}
yearly_convex_hull <- vulture_data %>%
  group_by(year = factor(lubridate::year(timestamp)), track=mt_track_id()) %>%
  group_map(~ st_sf(.y, geom = sf::st_convex_hull(sf::st_union(.x)))) %>%
  bind_rows()
yearly_convex_hull |> ggplot() +
  geom_sf(data = ne_coastline(returnclass = "sf", 50)) +
  geom_sf(aes(fill = year), alpha = .1) +
  theme_linedraw() +
  coord_sf(
    crs = sf::st_crs("+proj=aeqd +lon_0=-83 +lat_0=8 +units=km"),
    xlim = c(-3500, 3800), ylim = c(-4980, 4900)
  )
```