In popovs/OPPtools: This package contains tools and templates to analyze seabird tracking data

knitr::opts_chunk$set(echo = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      dpi = 150)

library(OPPtools)
library(dplyr)
library(track2KBA)

report_references <- 'Beal, M., Oppel, S., Handley, J., Pearmain, E. J., Morera-Pujol, V., Carneiro, A. P. B., Davies, T. E., Phillips, R. A., Taylor, P. R., Miller, M. G. R., Franco, A. M. A., Catry,
I., Patrício, A. R., Regalla, A., Staniland, I., Boyd, C., Catry, P., and Dias, M. P. (2021). track2KBA: An R package for identifying important sites for biodiversity from tracking data. Methods in Ecology and Evolution, 12, 2372–2378.
\n\n
Calenge, C. (2006) The package adehabitat for the R software: a tool for the analysis of space and habitat use by animals. Ecological Modelling, 197, 516-519
\n\n
Environment and Climate Change Canada. 2022. British Columbia Seabird Colony Inventory: Digital dataset. Canadian Wildlife Service, Environment and Climate Change Canada, Pacific Region.
\n\n
Gaston, A. J. and S. B. Dechesne (2020). Rhinoceros Auklet (Cerorhinca monocerata), version 1.0. In Birds of the World (A. F. Poole and F. B. Gill, Editors). Cornell Lab of Ornithology, Ithaca, NY, USA. 
\n\n
Johnson, D. S., London, J. M., Lea, M.-A. and Durban, J. W. (2008) Continuous-time correlated random walk model for animal telemetry data. Ecology, 89: 1208-1215. 
\n\n
Lascelles, B. G., Taylor, P. R., Miller, M. G. R., Dias, M. P., Oppel, S., Torres,
L., Hedd, A., Corre, M. L., Phillips, R. A., Shaffer, S. A., Weimerskirch, H., & Small, C. (2016). Applying global criteria to tracking data to define important areas for marine conservation. Diversity and Distributions, 22, 422–431.
\n\n
R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.'

dat <- opp_download_data(study = params$movebank_id)

\newpage

Data Summary

r if(min(dat$year) == max(dat$year)) {paste("In", min(dat$year))} else {paste("Between", min(dat$year), "and", max(dat$year))}, r length(unique(dat$individual_id)) individual r params$spps at the r unique(dat$study_site) colony were tagged and tracked with r as.character(unique(dat$sensor_type)) tags. These r length(unique(dat$individual_id)) individuals generated r nrow(dat) detections between r as.Date(min(dat$timestamp)) and r as.Date(max(dat$timestamp)). This report provides a detailed summary of the tracking data available for this population, as well as a description of how data were processed prior to analysis to identify high-use areas where birds may be at risk from marine spills.

Raw tracking data are stored on Movebank (https://www.movebank.org/) in project number r as.character(params$movebank_id). All analysis were conducted with R (version 4.1.0, R Core Team 2021). Code used to generate this report is available at: https://github.com/popovs/OPPtools.

dat2 <- opp2KBA(dat)

# Raw trip data
opp_map_tracks(tracks = dat2$data,
               center = dat2$site,
               zoom = NULL,
               coast_scale = 10,
               viridis_option = "D",
               show_locs = T)

\newpage

opp_logger_dotplot(dat)

\newpage

Classifying trips off-colony

Estimating high-use areas with 95% utilization distributions (UD) requires defining periods when individuals are away from the colony and interpolating GPS locations at a consistent time-interval. We first identified trips where birds were away from the colony, and classified these trips based on how complete and consistent GPS fixes were acquired during these trips. Trips were identified as any continuous track at least r params$minDist km from the colony which lasted for a minimum of r params$minDur r ifelse(params$minDur == 1, paste("hour"), paste("hours")). Trips were considered incomplete if the trips started or ended more than r params$maxDist km from the colony. Trips were considered to be gappy if there were gaps between successive GPS fixes greater than r params$gapTime hr and the bird traveled farther than r params$gapDist km during that time.

Table 1 summarizes trip attributes for each type of trip. The appendix includes plots showing the trip classification from each individual deployment. r if (length(grep('storm-petrel', params$spp)) == 1) "Because of the low GPS sampling rate and long commuting distances of storm-petrels, some consecutive trips are potentially combined here. These merged trips were retained in order to include time spent close to the colony in estimated utilization distributions. For this reason, trip summary statistics in Tables 1 and 2 should be interpreted with caution." Because this analysis is primarily concerned with identifying at-sea areas where birds may be vulnerable to marine spills, no effort was made to distinguish between foraging trips and non-foraging trips.

trips <- opp_get_trips(dat2,
              innerBuff = params$minDist,
              returnBuff = params$maxDist,
              duration = params$minDur,
              gapTime = params$gapTime,
              gapDist = params$gapDist,
              gapLimit = params$gapLimit,
              showPlots = F
              )

ft <- trips@data %>%
  filter(Type != 'Non-trip') %>%
  group_by(ID, tripID, Type) %>%
  mutate(
    duration = as.numeric(difftime(max(DateTime), min(DateTime), units = 'hours'))
  ) %>%
  group_by(Type) %>%
  summarize(
    `Number of trips` = length(unique(tripID)),
    `Number of raw GPS fixes` = n(),
    `Median duration (hrs)` = round(median(duration), 1),
    `Median of max trip distance (km)` = round(median(max(ColDist))/1000, 1)
    ) %>%
  flextable::flextable(cwidth =1.25) %>% 
  flextable::align(align = 'center', part = 'all') %>%
  flextable::set_caption("Summary of trip attributes by trip classification.")

ft

Track interpolation

Locations collected during complete, incomplete, and gappy trips were interpolated at a r params$timestep interval using a continuous time correlated random walk model using the crawl package (version 2.2.1, Johnson et al. 2008) in R. Gaps in tracks during a trip, as defined above, were r ifelse(params$interpolateGaps == F, '**not**','') interpolated. Figure 2 shows maps comparing the raw and interpolated locations within trips for all deployments. Table 2 provides a summary of the number of raw and interpolated locations, trip duration, maximum trip distance and trip type for each identified foraging trip. Plots comparing the raw and interpolated data for each individual deployment are provided in the appendix.

interp <- ctcrw_interpolation(trips,
                    site = dat2$site,
                    type = c('Complete', 'Incomplete', 'Gappy'),
                    timestep = params$timestep,
                    interpolateGaps = params$interpolateGaps,
                    showPlots = F)

\newpage

# Raw trip data
p_raw <- opp_map_tracks(tracks = interp$data,
               center = dat2$site,
               zoom = NULL,
               coast_scale = 10,
               viridis_option = "D",
               show_locs = T)

# Interpolated data
p_int <- opp_map_tracks(tracks = interp$interp,
               center = dat2$site,
               zoom = NULL,
               coast_scale = 10,
               viridis_option = "D",
               show_locs = T)

print(ggpubr::ggarrange(p_raw, p_int,
                        ncol = 1,
                        labels = c("A", "B"),
                        align = "v")
      )

\newpage

my_trip_summary <- OPPtools::sum_trips(interp) %>% 
  dplyr::mutate(duration = signif(duration, 3),
                max_dist_km = signif(max_dist_km, 3),
                departure = as.Date(departure)) %>% 
  dplyr::select(tripID, raw_n_locs, interp_n_locs, departure, duration, max_dist_km, complete) %>% 
  as.data.frame() %>% 
  dplyr::rename(`Trip ID` = tripID,
                `GPS locations` = raw_n_locs,
                `Interpolated locations` = interp_n_locs,
                `Start date` = departure,
                `Maximum distance (km)` = max_dist_km,
                `Duration (hrs)` = duration,
                Type = complete) %>% 
  flextable::flextable(cwidth =c(1, 0.8, 0.8, 1, 0.8, 1, 0.8)) %>% 
  flextable::align(align = 'center', part = 'all') %>%
  flextable::set_caption("Summary of individual foraging trip attributes.")


my_trip_summary

# Report will end here if there are fewer than 20 tracks in the dataset
n_tracks <- length(unique(interp$interp$ID))

show.text <- n_tracks < 20
calc.kernels <- n_tracks >= 20

\newpage

Kernel density estimates

h <- opp_href(interp$interp)
s <- opp_step(interp$interp)
s <- ifelse(s < 1, 1, s)

Kernel density estimates are sensitive to the kernel smoother; therefore it is important to choose an appropriate smoother for the species and GPS fix rate. The appendix includes plots comparing kernel smoothers for each individual track. These plots were examined to identify any potential biases associated with using a particular smoother.

We considered three smoothers using a traditional kernel density estimate: href, href/2, and step. The href smoother is a data-driven method which is calculated from the variance in the X and Y directions in the data; in this case the median href was r signif(h, 2) km. Initial data exploration indicated that the href parameter could over-smooth tracks. We therefore also considered the href/2 smoother, which was r signif(h/2, 2) km. Finally, we tested a step smoother based on the median step length between interpolated locations across all tracks. r if(s == 1){"The step function resulted in a value of less than 1 km; therefore, step kernel density estimates were calculated with a smoother value of 1 km."} else {paste("The step smoother was", signif(s, 2), "km.")} All UDs were calculated using the the adehabitatHR package (version 0.4.19, Calenge 2006) in R.

href_ud <- opp_kernel(interp,
                      interpolated = TRUE,
                      smoother = "href",
                      extendGrid = params$gridExtend,
                      res = params$gridRes)

href0.5_ud <- opp_kernel(interp,
                      interpolated = TRUE,
                      smoother = (h/2),
                      extendGrid = params$gridExtend,
                      res = params$gridRes)

step_ud <- opp_kernel(interp,
                      interpolated = TRUE,
                      smoother = s,
                      extendGrid = params$gridExtend,
                      res = params$gridRes)

cat("# Low track sample size

The relatively small number of GPS tracks available for this breeding population is not sufficient for characterizing important at-sea areas (Beal et al 2021, Lascelles et al 2016).")

\newpage

References

r report_references

\newpage

Appendix - Diagnostic plots for individual tracks

This appendix includes plots showing the analysis steps described above for each individual track in the dataset for this population. On each page, the top plot (A) shows the trip classifications. Locations are plotted as distance from colony over time, with points coloured according to trip type. Complete trips start and end close to the colony, and do not contain any large gaps in GPS locations. Incomplete trips do not contain any large gaps in GPS locations during the trips, but either start or end far away from the colony. Gappy plots contain a prolonged period of missing locations where the bird also travelled a long distance. The middle plot (B) shows the interpolated GPS locations relative to raw GPS locations. Locations are plotted as distance from colony over time. Gaps in the interpolated (teal) tracks, where no dashed line is plotted, were times when the bird movement was classified as a 'non-trip'. Times where the interpolated track is shown with a dashed line, but no points are plotted, are large gaps in GPS tracking during a trip where no interpolation was performed over gaps in the GPS locations. The bottom plot (C) compares three different kernel density estimation (KDE) smoothers; each plot shows how the 50% and 95% utilization distribution changes depending on the smoothing parameter and method.

\newpage

trip_plots <- plot_trip_dist(data = trips, plotsPerPage = 1,
                             innerBuff = params$minDist,
                             returnBuff = params$maxDist,
                             showPlots = F)

interp_plots <- plot_interp_dist(data = interp,
                                 plotsPerPage = 1,
                                 showPlots = F)

plot_uds <- list(Step = step_ud,
                 Href = href_ud,
                 `Href/2` = href0.5_ud)

pp <- opp_map_indUD(plot_uds,
                    tracks = interp$interp,
                    center = dat2$site,
                    ud_levels = c(50, 95),
                    zoom = NULL,
                    coast_scale = 50,
                    viridis_option = "A")

for (i in 1:length(trip_plots)) {

  nn<- names(trip_plots)[i]

  if (nn %in% names(interp_plots) & nn %in% names(pp)) {

  p1 <- trip_plots[[nn]] +
    ggplot2::labs(title = paste0('Diagnostic plots for ', names(trip_plots),'.'))

  p2 <- interp_plots[[nn]]+
    ggplot2::labs(title = '')

  p3 <- pp[[nn]]+
    ggplot2::labs(title = '')

  ptop <- (ggpubr::ggarrange(p1, p2,
                           ncol = 1,
                           labels = c("A", "B"))
  )

  print(ggpubr::ggarrange(ptop, p3,
                          ncol = 1,nrow = 2,heights = c(1.5,1),
                          labels = c('', 'C'))
  )
  } else {
      p1 <- trip_plots[[nn]] +
    ggplot2::labs(title = paste0('Diagnostic plots for ', names(trip_plots),'.'),
                  caption = 'Interpolation and kernel density plots were not run because there were no off-colony trips.')
  }
}

popovs/OPPtools documentation built on July 8, 2023, 2:29 a.m.