Example analysis of albatross trajectories

Animating

Using the tools from gganimate and ggspatial we can also add more context to the map and animate it. Here annotate_map_tile is used to add the OpenStreetMap background map and annotate_scale to add a scale bar. A variety of different animations are possible, here we show birds originating from different study_site’s. These kind of animations help to gain insights into the differences between subgroups of the data set. Besides tagging location, also different years, life stages or sexes are clear candidates to compare.

require(gganimate)
require(ggspatial)
animation_site <- ggplot() +
  annotation_map_tile(zoom = 5, progress = "none") +
  geom_sf(
    data = mt_track_lines(data),
    mapping = aes(group = individual_local_identifier),
    color = "black"
  ) +
  transition_states(study_site, state_length = 2) +
  enter_fade() +
  exit_fade() +
  ease_aes("cubic-in-out") +
  labs(title = "{closest_state}") +
  annotation_scale()
animation_site

We can also take this one step further by animating the map view. To do that we use a local equal area projection and define for each state a manual zoom.

animation_site +
  coord_sf(
    crs = sf::st_crs("+proj=aeqd +lon_0=-83 +lat_0=-6 +units=km")
  ) +
  view_zoom_manual(
    pause_length = 2,
    xmin = c(-850, -900, -950),
    ymin = c(-350, -800, -350),
    ymax = c(610, 700, 610),
    xmax = c(500, 600, 500)
  )

It is also possible to animate the movement of individuals. To make the animation smoother, we linearly interpolate the position of individuals. Here we interpolate to a location every 60 minutes. We avoid interpolating longer time lags (in this case larger than 3 hours) to avoid inferring movement between location to long apart, this causes some individuals to temporarily disappear from the animation. To ensure that the data set is fully regular we omit the existing locations. For visual purposes we only show a short period of a few days, however, this can be adjusted to what is most appropriate for the intended visualization. For each timestamp we render one frame. Here we color by individual but other attributes like the speed could also be used.

data_interpolated <- mt_interpolate(
  data[!sf::st_is_empty(data), ],
  time = seq(
    as.POSIXct("2008-7-27"),
    as.POSIXct("2008-8-1"), "60 mins"
  ),
  max_time_lag = units::as_units(3, "hours"),
  omit = TRUE
)
animation <- ggplot() +
  annotation_map_tile(zoom = 5, progress = "none") +
  annotation_scale() +
  theme_linedraw() +
  geom_sf(
    data = data_interpolated, size = 3,
    aes(color = individual_local_identifier)
  ) +
  transition_manual(timestamp) +
  labs(
    title = "Galapagos Albatrosses",
    subtitle = "Time: {current_frame}",
    color = "Individual"
  )
animate(animation,
  nframes = length(unique(data_interpolated$timestamp))
)

The previous animation was made using a relatively simple approach. However, with some additional changes it is possible to add a tail to the movement of individuals and use a dark color scheme. To visualize the tail we use shadow_wake which does not work in combination with transition_manual, therefore we switch here to the usage of transition_time. To keep the continuity of movement, this case we do in interpolate over longer time spans.

date_range <- as.POSIXct(c("2008-7-29", "2008-8-1"))
ts <- mt_time(data)
times <- sort(unique((c(date_range, ts[ts < max(date_range) & ts > min(date_range)]))))
data_interpolated <-
  mt_interpolate(
    data[!sf::st_is_empty(data), ],
    times,
    omit = TRUE
  )
label_df <- data.frame(
  timestamp = date_range,
  display_time = lubridate::with_tz(date_range, "America/Lima")
)
animation <- ggplot() +
  annotation_map_tile("cartodark", zoom = 5, progress = "none") +
  annotation_scale(bar_cols = c("gray80", "gray40"), text_col = "gray80") +
  geom_sf(data = mt_track_lines(data), color = "grey40") +
  theme_linedraw() +
  geom_sf(
    data = data_interpolated, size = 3,
    aes(color = individual_local_identifier)
  ) +
  scale_color_brewer(palette = "Set1") +
  guides(color = "none") +
  xlab("") +
  ylab("") +
  geom_text(
    data = label_df,
    aes(label = display_time, x = -10100000, y = -1370000),
    color = "grey80", size = 3, hjust = 0
  ) +
  transition_time(timestamp) +
  shadow_wake(0.2, exclude_layer = 6)

The time period selected is 3 days and we generate one frame every 30 minutes. The detail argument is used to estimate additional intermediate locations so that the tail gets a smooth shape.

animate(animation,
  nframes = 3 * 24 * 2 + 1, detail = 5
)

Note that the timezone is set to Lima, this provides us the possibility to easily identify that most movement occurs during the day and that the birds are relatively stationary and floating with the current during the night.

Spatial analysis

Now lets take a deeper dive into the albatross data. We have seen that these birds breed on the Galapagos islands and forage on the coast of South America. For some of these birds only one foraging trip has been recorded while for others multiple trips are available. In this example we explore the tracks within four different stages of these trips (in the breeding area, the outbound and inbound flight and the foraging area). To do this we do spatial intersections with the defined regions, split the track into different sections and summarize each track.

First we download the data again, however this time we additionally include the accelerometer data, this gives us one trajectory per individual. Some individuals are marked in the deployment data as being ‘not used for analysis’, so we omit these individuals from the dataset.

require(units)
require(dplyr)
require(sf)
data <- movebank_download_study("Galapagos Albatrosses",
  sensor_type_id = c("gps", "acceleration")
)
data <- data %>%
  filter_track_data(deployment_comments != "not used in analysis")

The dataset now contains both accelerometer and gps data, these records are not necessarily recorded at the same time and are thus reported as separate rows in the data. The accelerometer data are not associated with locations. In order to associate the measurements to a specific region we do a linear interpolation for all missing locations (the default of mt_interpolate). The records as downloaded by default are not sorted by time and individual, therefore we first need to ensure this ordering is correct.

data <- data[order(mt_track_id(data), mt_time(data)), ]
data <- mt_interpolate(data)

Next, we want to associate records with the respective regions (breeding and foraging), we do this by a spatial intersection (st_join). For simplicity here the breeding region is defined as with a fixed distance from any of the deployment locations of the tracks. The foraging region is defined as the area within a 100 kilometer from the coast of South America. As a result a new column called region is generated where all records outside of either of the regions has a NA value.

library(rnaturalearth)
breeding_area <- st_buffer(mt_track_data(data)$deploy_on_location, as_units(25, "km")) %>%
  st_union()
foraging_area <- ne_countries(110,
  returnclass = "sf",
  continent = "South America"
) %>%
  st_union() %>%
  st_buffer(as_units(100, "km"))
regions <- st_sf(data.frame(
  region = c("Breeding", "Foraging"),
  polygon = c(breeding_area, foraging_area)
))
data <- st_join(data, regions)

Next, we need to find those trajectories that are either inbound or outbound. To identify these we look for series of NA records where the previous location was within the breeding region and the next location in the foraging region, or vice versa. This approach has the advantage that only full trips from one to the other region are defined as commutes. Shorter trips outside of the breeding or foraging region can be redefined to be either foraging or breeding so that these sections are not cut up.

data <- data %>%
  group_by(mt_track_id(.)) %>%
  mutate(
    region_change = paste(
      vctrs::vec_fill_missing(region),
      vctrs::vec_fill_missing(region, "up")
    ),
    region = case_match(region_change,
      "Foraging Breeding" ~ "Inbound",
      "Breeding Foraging" ~ "Outbound",
      "Breeding Breeding" ~ "Breeding",
      "Foraging Foraging" ~ "Foraging",
      .default = region
    )
  )

Next, we redefine tracks, up to now the tracks corresponded to all tracking data from one individual. Now we want to convert all continuous data within one region as being one track. Here rle is used to identify continues section within one region. We do omit locations that have no region associated, these occur at the start or end of trajectories

data <- data %>%
  mutate(
    sequence_number = with(rle(region), rep(seq_along(lengths), lengths)),
    track = paste(individual_local_identifier, region, sequence_number)
  ) %>%
  ungroup() %>%
  mutate_track_data(individual = droplevels(individual_local_identifier)) %>%
  mt_set_track_id("track") %>%
  filter(!is.na(region))

These used gps tracking devices were also equipped with accelerometers, here we do not go into the calibration of these measurements, however we can still calculate a proxy for energy expenditure, the dynamic body acceleration (DBA). The acceleration measurements, in this case collected in bursts, are stored in movebank as a text string which can be parsed to an expenditure for each burst as follows (notice that for these tags no calibration is available):

acc_to_dba <- function(x) {
  acc_mat <- matrix(as.numeric(unlist(strsplit(x, " "))), nrow = 2)
  mean(colSums(abs(acc_mat - rowMeans(acc_mat))))
}
data$dba <- unlist(lapply(data$eobs_accelerations_raw, acc_to_dba))

Now we can calculate summaries for each track, in this case we calculate the number of records in each track but also the start and end so the duration can be calculates. Furthermore we calculate summary statistics for the ground speed as measured by the gps and the DBA.

track_summary <- data %>%
  mt_track_lines(
    region = unique(region), n = dplyr::n(), start = min(timestamp),
    end = max(timestamp),
    across(
      all_of(c("ground_speed", "dba")),
      list(
        mean = function(x) mean(x, na.rm = TRUE),
        sd = function(x) sd(x, na.rm = TRUE)
      )
    )
  ) %>%
  mutate(duration = as_units(end - start))

We can first tabulate the number of tracks per region for each individual. Here we see that each individual has one track more in the breeding region then in the other regions. This makes sense as all individuals were equipped with transmitters in the breeding region and data is only retrieved once they returned.

table(track_summary$individual, track_summary$region)

Now we can plot the obtained summary statistics for the trajectories in each ‘region’. First we plot a map, to show where each track is, the inbound tracks are longer then the outbound tracks.

ggplot(track_summary) +
  geom_sf(data = ne_coastline(returnclass = "sf", 50)) +
  geom_sf(aes(color = region)) +
  theme_linedraw() +
  coord_sf(
    crs = st_crs("+proj=aeqd +lon_0=-83 +lat_0=-6 +units=km"),
    xlim = c(-1000, 600),
    ylim = c(-800, 700)
  ) +
  labs(color = "Region") +
  scale_color_brewer(type = "qual")

However the distance is covered in a shorter time as can be seen in the following graph. Here the units are changes to days to facilitate the interpretation.

ggplot(track_summary, aes(x = region, y = duration)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(aes(color = individual, group = individual),
    position = position_jitterdodge()
  ) +
  xlab("") +
  scale_y_units("Duration", unit = "days", trans = "log10") +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5)) +
  theme_linedraw() +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5)) +
  scale_color_brewer("Individual", type = "qual", palette = "Set1")

This corresponds with the fact that the ground speed measured by the gps sensor is on average higher (notice units get propagated from movebank).

ggplot(
  track_summary,
  aes(x = region, y = ground_speed_mean)
) +
  theme_linedraw() +
  geom_boxplot(outlier.shape = NA) +
  geom_point(aes(color = individual, group = individual),
    position = position_jitterdodge()
  ) +
  xlab("") +
  ylab("Mean ground speed") +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5)) +
  scale_color_brewer("Individual", type = "qual", palette = "Set1")

At the same time DBA is lower in these inbound flight corresponding to the tail winds these birds profit from. We also see that birds have a higher speed and DBA in the foraging region compared to the breeding region corresponding with their respective activities.

ggplot(
  track_summary,
  aes(x = region, y = dba_mean)
) +
  theme_linedraw() +
  geom_boxplot(outlier.shape = NA) +
  geom_point(aes(color = individual, group = individual),
    position = position_jitterdodge()
  ) +
  ylab("Mean DBA") +
  xlab("") +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5)) +
  scale_color_brewer("Individual", type = "qual", palette = "Set1")

This analysis provides an example on how movement data can be analyzed using move2 and how the retention of meta data is important to maintain consistency throughout the analysis. This enables tight integration with other packages facilitating visualization and spatial operations.