In FlukeAndFeather/breath: Analyze whale breathing behavior

knitr::opts_chunk$set(echo = TRUE)
library(cowplot)
library(DT)
library(glue)
library(ncdf4)
library(RcppRoll)
library(shiny)
library(tidyverse)

Trying a new approach to breath identification that isolates surface intervals first.

Load data

Read a medium-term blue whale tag NetCDF file and clean the data up a bit.

# load data
nc_path <- "~/Documents/GitHub/exploration/breath-rate/data/Bm181021-TDR11_prh8.nc"
prh_nc <- nc_open(nc_path)

# read variables
dn <- ncvar_get(prh_nc, "DN")
dt <- breath::dn_to_posix(dn, "US/Pacific")
secs <- as.numeric(dt - dt[1], unit = "secs")
p <- ncvar_get(prh_nc, "P")
pitch <- ncvar_get(prh_nc, "pitch")
speed <- ncvar_get(prh_nc, "speedJJ")
fs <- ncatt_get(prh_nc, "P", "sampling_rate")$value

# Create features (smooth over 2s and calculate rate of change)
feat <- tibble(dt, secs, p, pitch, speed) %>%
  mutate(p2 = roll_mean(p, n = fs * 2, fill = NA),
         pitch2 = roll_mean(pitch, n = fs * 2, fill = NA),
         speed2 = roll_mean(speed, n = fs * 2, fill = NA),
         ddt_p = (lead(p2) - p2) * fs,
         ddt_p2 = roll_mean(ddt_p, n = fs * 2, fill = NA),
         ddt_pitch = (lead(pitch2) - pitch2) * fs,
         ddt_pitch2 = roll_mean(ddt_pitch, n = fs * 2, fill = NA))

Isolate surface intervals

Use run-length encoding to identify surface intervals (<10m for >2s).

# find regions of depth<10m for >2s
depth_rle <- rle(feat$p2 < 10)
depth_rle$surface <- depth_rle$lengths > 2 * fs & depth_rle$values == TRUE
depth_rle$start <- cumsum(lag(depth_rle$lengths, default = 1))
surfaces <- tibble(
  lengths = depth_rle$lengths,
  values = depth_rle$values,
  is_surface = lengths > 10 * fs & values == TRUE,
  surface_idx = ifelse(is_surface, cumsum(is_surface), NA),
  surface_start = cumsum(lag(lengths, default = 1)),
  surface_end = lead(surface_start, default = nrow(feat) + 1) - 1
) %>%
  filter(is_surface) %>%
  select(surface_idx:surface_end)

Identify breaths

Breaths are classified as local depth minima shallower than a meter. I tried various approaches with pitch and rate of pitch change, but they resulted in too many false negatives. Simplest approach seems to work best.

# Find local minima in window
find_peaks <- function(x, width, prominence, idx_off = 0) {
  if (length(x) < width)
    return(NULL)
  peak_idx <- which.max(x)

  # Skip edges
  if (peak_idx %in% c(1, length(x)))
    return(NULL)

  peak_val <- x[peak_idx]
  left_bound <- last(which(x[1:(peak_idx - 1)] > peak_val))
  if (is.na(left_bound))
    left_bound <- 1
  right_bound <- peak_idx + first(which(x[(peak_idx + 1):length(x)] > peak_val))
  if (is.na(right_bound))
    right_bound <- length(x)
  left_low <- min(x[left_bound:(peak_idx - 1)])
  right_low <- min(x[(peak_idx + 1):right_bound])
  low_val <- max(left_low, right_low)
  peak <- c(peak_idx + idx_off, peak_val - low_val)

  left <- right <- NULL

  # recurse on the left
  new_left <- peak_idx - floor(width / 2)
  if (new_left >= width)
    left <- find_peaks(x[1:new_left],
                       width, prominence,
                       idx_off + 0)

  # recurse on the right
  new_right <- (peak_idx + floor(width / 2))
  if (length(x) - new_right + 1 >= width)
    right <- find_peaks(x[new_right:length(x)],
                        width, prominence,
                        idx_off + new_right - 1)

  if (peak[2] >= prominence)
    rbind(peak, left, right)
  else
    rbind(left, right)
}

find_breaths <- function(surface_idx, surface_start, surface_end, period) {
  df <- slice(feat, surface_start:surface_end)

  # find local depth minima in 2 x period second window at surface
  peaks <- find_peaks(-df$p2, width = 2 * period * fs, prominence = 0.5)
  if (is.null(peaks))
    return(NULL)
  peaks_p <- df$p2[peaks[, 1]]
  breath_idx <- peaks[, 1][peaks_p < 0.5]
  slice(df, breath_idx) %>%
    mutate(surface_idx = surface_idx,
           row_idx = surface_start + breath_idx - 1)
}

# Find breaths (assumes period >=5, v. conservative for blue whale)
breaths <- pmap_dfr(surfaces, find_breaths, period = 5)

Visualize breaths

Plot depth, pitch, and speed to validate breaths in surface intervals

# Plot depth, pitch, and speed with breaths for a surface interval
plot_surface <- function(surface_idx) { 
  i <- surfaces$surface_start[surface_idx]:surfaces$surface_end[surface_idx]
  df <- slice(feat, i)
  b <- filter(breaths, surface_idx == { surface_idx })
  suppressWarnings(
    plot_grid(
      ggdraw() + 
        draw_label(glue("Surface {surface_idx}"), 
                   fontface = 'bold', 
                   x = 0, 
                   hjust = 0) +
        theme(plot.margin = margin(0, 0, 0, 7)),
      ggplot(df, aes(dt, p2)) +
        geom_line(size = 0.2, color = "blue") +
        geom_vline(aes(xintercept = dt), b, linetype = 3) +
        scale_x_datetime(date_labels = "%b-%d %H:%M") +
        scale_y_reverse("Depth (m)") +
        theme_minimal() +
        theme(axis.text.x = element_blank(),
              axis.title.x = element_blank()),
      ggplot(df, aes(dt, pitch2 * 180 / pi)) +
        geom_line(size = 0.2, color = "red") +
        geom_vline(aes(xintercept = dt), b, linetype = 3) +
        scale_x_datetime(date_labels = "%b-%d %H:%M") +
        scale_y_continuous("Pitch (deg)", limits = c(-15, 15)) +
        theme_minimal() +
        theme(axis.text.x = element_blank(),
              axis.title.x = element_blank()),
      ggplot(df, aes(dt, speed2)) +
        geom_line(size = 0.2, color = "purple") +
        geom_vline(aes(xintercept = dt), b, linetype = 3) +
        scale_x_datetime(date_labels = "%b-%d %H:%M") +
        scale_y_continuous("Jiggle speed (m/s)") +
        theme_minimal() +
        theme(axis.title.x = element_blank()),
      ncol = 1,
      align = "v"
    )
  )
}

# Spot check five surface intervals
set.seed(934)
surf_idx <- sort(sample(unique(breaths$surface_idx), 5))
walk(surf_idx, ~ print(plot_surface(.x)))

Shiny app

renderPlot({
  hist(faithful$eruptions, probability = TRUE,
       breaks = as.numeric(input$n_breaks),
       xlab = "Duration (minutes)",
       main = "Geyser eruption duration")

  dens <- density(faithful$eruptions, adjust = input$bw_adjust)
  lines(dens, col = "blue")
})