Trajectory analysis of Turkey vultures

Retrieve data

library(move2)
vulture_data <-
  movebank_download_study("Turkey vultures in North and South America")
vulture_data
#> A <move2> with `track_id_column` "individual_local_identifier" and `time_column`
#> "timestamp"
#> Containing 19 tracks lasting on average 1083 days in a
#> Simple feature collection with 215719 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -123.9255 ymin: -39.8214 xmax: -48.5788 ymax: 54.01817
#> Geodetic CRS:  WGS 84
#> # A tibble: 215,719 × 7
#>   sensor_type_id individual_local_ident…¹ manually_marked_outl…² timestamp          
#>          <int64> <fct>                    <lgl>                  <dttm>             
#> 1            653 Irma                     FALSE                  2004-09-06 17:00:00
#> 2            653 Irma                     FALSE                  2004-09-06 18:00:00
#> 3            653 Irma                     FALSE                  2004-09-06 19:00:00
#> 4            653 Irma                     FALSE                  2004-09-06 20:00:00
#> 5            653 Irma                     FALSE                  2004-09-07 00:00:00
#> # ℹ 215,714 more rows
#> # ℹ abbreviated names: ¹individual_local_identifier, ²manually_marked_outlier
#> # ℹ 3 more variables: event_id <int64>, visible <lgl>, geometry <POINT [°]>
#> First 5 track features:
#> # A tibble: 19 × 54
#>   deployment_id   tag_id individual_id animal_life_stage animal_mass attachment_type
#>         <int64>  <int64>       <int64> <fct>                     [g] <fct>          
#> 1      17225120 16883951      17002744 adult                    2372 harness        
#> 2      17225134 16883951      17002745 adult                    1951 harness        
#> 3      17225131 16883928      17002732 adult                    2012 harness        
#> 4      17225133 16883937      17002737 adult                    2108 harness        
#> 5      17225122 16883967      17002753 adult                      NA harness        
#> # ℹ 14 more rows
#> # ℹ 48 more variables: deployment_comments <chr>, deploy_off_timestamp <dttm>,
#> #   deploy_on_timestamp <dttm>, duty_cycle <chr>,
#> #   deployment_local_identifier <fct>, …

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.

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

A map showing the trajectories of vultures using coast lines as a reference

Categorize seasons

library(magrittr, quietly = TRUE)
#> 
#> Attaching package: 'magrittr'
#> The following object is masked from 'package:raster':
#> 
#>     extract
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.

vulture_data <- vulture_data %>% mutate(azimuth = mt_azimuth(.), speed = mt_speed(.))

Seasonal distribution per individual

library(circular, quietly = TRUE)
#> 
#> Attaching package: 'circular'
#> The following objects are masked from 'package:stats':
#> 
#>     sd, var
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"
  )))
#> `summarise()` has grouped output by 'season'. You can override using the `.groups`
#> argument.

# Load purrr for map function
library(purrr, quietly = TRUE)
#> 
#> Attaching package: 'purrr'
#> The following object is masked from 'package:magrittr':
#> 
#>     set_names
# Load tidy r for unnest function
library(tidyr, quietly = TRUE)
#> 
#> Attaching package: 'tidyr'
#> The following object is masked from 'package:magrittr':
#> 
#>     extract
#> The following object is masked from 'package:raster':
#> 
#>     extract
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")

A plot of the azimuths of the tracks split out per season and individual

Individual trajectory

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()
#> Warning: Removed 6683 rows containing missing values or values outside the scale range
#> (`geom_point()`).

A plot exploring the relation between direction and speed for one track

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 that visualizes the speed and directions per season for one individual

Plot turn 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()
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
#> Warning: Removed 8 rows containing missing values or values outside the scale range
#> (`geom_bar()`).

Plot of the turning angle distribution

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.

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

Plot of the net squared displacement overtime per individual

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.

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()
#> A <move2> with `track_id_column` "individual_local_identifier" and `time_column`
#> "timestamp"
#> Containing 15 tracks lasting on average 752 days in a
#> Simple feature collection with 166 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -120.2711 ymin: -28.16434 xmax: -62.0751 ymax: 45.63105
#> Geodetic CRS:  WGS 84
#> First 5 features:
#>   individual_local_identifier           timestamp                   geometry
#> 1                      Disney 2004-11-14 19:59:51 POINT (-81.45408 31.70607)
#> 2                      Disney 2004-11-25 13:38:30 POINT (-81.51575 31.86253)
#> 3                      Disney 2004-12-15 20:05:43 POINT (-81.27212 31.69371)
#> 4                      Disney 2005-01-12 17:30:58 POINT (-82.39753 32.12783)
#> 5                      Disney 2005-01-17 14:18:56 POINT (-82.37178 32.12628)
#> First 5 track features:
#> # A tibble: 15 × 54
#>   deployment_id   tag_id individual_id animal_life_stage animal_mass attachment_type
#>         <int64>  <int64>       <int64> <fct>                     [g] <fct>          
#> 1      17225133 16883937      17002737 adult                    2108 harness        
#> 2      17225132 16883941      17002739 adult                    1750 harness        
#> 3      17225135 16883943      17002740 adult                    1975 harness        
#> 4      17225136 16883945      17002741 adult                    1810 harness        
#> 5      17225130 16883947      17002742 adult                      NA harness        
#> # ℹ 10 more rows
#> # ℹ 48 more variables: deployment_comments <chr>, deploy_off_timestamp <dttm>,
#> #   deploy_on_timestamp <dttm>, duty_cycle <chr>,
#> #   deployment_local_identifier <fct>, …

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

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