Load Packages

pkgTest <- function(x)
{
  if (!require(x,character.only = TRUE))
  {
    install.packages(x,dep=TRUE)
    if(!require(x,character.only = TRUE)) stop("Package not found")
  }
}
`%notin%` <- Negate(`%in%`)
# Plotting and Data Wrangling
pkgTest("tidyverse")
pkgTest("RColorBrewer")
pkgTest("ggpubr")
pkgTest("lubridate")
pkgTest("rnaturalearth")
pkgTest("maptools")
pkgTest("ggmap")
pkgTest("metR")
pkgTest("kableExtra")
pkgTest("rnaturalearthdata")
pkgTest("sf")
pkgTest("sp")
pkgTest("moveHMM")
pkgTest("nlme")
pkgTest("stats")
pkgTest("AICcmodavg")
pkgTest('MetricsWeighted')
pkgTest("MASS")
pkgTest("ggspatial")

# load a custom function for extracting the Confidence intervals for stationary probabilities
source('./functions/moveHMMstationaryCI.R')

# Function to find the Pooled Standard Deviation of a Sample of Samples
# i.e. ILI by dive has a stdev, but a group of dives will have a pooled stdev
# Input: 
#   stdevs = standard deviations of sample
#   numSamples = sample sizes for each observation
pooledstd <- function(stdevs,numSamples) {
  psd <- sqrt( sum((numSamples-1)*(stdevs^2))/sum(numSamples-1) )
  return(psd)
} 

# Function to convert between Logit and Probability
logit2prob <- function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}
# Function to Convert Probability to Logit
prob_logit <- function(x) {
  out <- log(x / (1 - x))
  return(out)
}
# Function to perform a KS Test and returns a dataframe with the result
ksResults <- function(data1, data2, alt, Distribution, ID) {
  resKS <- ks.test(x=data1,y=data2,alternative=alt)
  res <- data.frame(ID = ID,
                    Distribution = Distribution,
                    test = resKS$alternative,
                    pvalue = as.numeric(resKS$p.value),
                    dvalue = as.numeric(resKS$statistic),
                    sampleSizeX = as.numeric(length(data1)),
                    meanX = mean(data1,na.rm=TRUE),
                    sdX= sd(data1,na.rm = TRUE),
                    sampleSizeY = as.numeric(length(data2)),
                    meanY = mean(data2,na.rm=TRUE),
                    sdY= sd(data2,na.rm = TRUE))
  return(res)
}

Load and Process Data

## Load Pre-processed Data
load("./dataProcessed/diveDataset.RData") # Load the processed dive dataframe
load("./dataProcessed/datasetFTLE.RData") # Load the processed dive+FTLE dataframe
load("./dataProcessed/datasetSimsFTLE.RData") # Load the Null Model data
load("./dataProcessed/datasetBackgroundFTLE.RData") # Load the Background Points data
load("./rmdFiles/world.RData") # load world SF basemap
load("./rmdFiles/bf.RData") # load Bathymetric Features
load("./rmdFiles/hfRadarStations.RData")
load("./rmdFiles/HFRadarCoverageMCP.RData") 
# HF Radar Coverage Area and Station Locations
hfRadarStns$name <- "Radar Station"

## Load CAR1 models
# Real Data Model
if (file.exists("./Models/feed_glmmPQL_FTLE_CAR1_Hrs.rds")) {
  feed_glmmPQL_FTLE_CAR1_Hrs <- readRDS("./Models/feed_glmmPQL_FTLE_CAR1_Hrs.rds")
} else { 
  cat(paste0("Could not find ./Models/feed_glmmPQL_FTLE_CAR1_Hrs.rds, will need to recreate below."))
}

# Time-Shifted Models
if (file.exists("./Models/feed_glmmPQL_p24_CAR1_Hrs.rds")) {
  feed_glmmPQL_p24_CAR1_Hrs <- readRDS("./Models/feed_glmmPQL_p24_CAR1_Hrs.rds")
} else { 
  cat(paste0("Could not find ./Models/feed_glmmPQL_p24_CAR1_Hrs.rds, will need to recreate below."))
}
if (file.exists("./Models/feed_glmmPQL_p48_CAR1_Hrs.rds")) {
  feed_glmmPQL_p48_CAR1_Hrs <- readRDS("./Models/feed_glmmPQL_p48_CAR1_Hrs.rds")
} else { 
  cat(paste0("Could not find ./Models/feed_glmmPQL_p48_CAR1_Hrs.rds, will need to recreate below."))
}
if (file.exists("./Models/feed_glmmPQL_p96_CAR1_Hrs.rds")) {
  feed_glmmPQL_p96_CAR1_Hrs <- readRDS("./Models/feed_glmmPQL_p96_CAR1_Hrs.rds")
} else { 
  cat(paste0("Could not find ./Models/feed_glmmPQL_p96_CAR1_Hrs.rds, will need to recreate below."))
}
if (file.exists("./Models/feed_glmmPQL_p192_CAR1_Hrs.rds")) {
  feed_glmmPQL_p192_CAR1_Hrs <- readRDS("./Models/feed_glmmPQL_p192_CAR1_Hrs.rds")
} else { 
  cat(paste0("Could not find ./Models/feed_glmmPQL_p192_CAR1_Hrs.rds, will need to recreate below."))
}

# Simulated Tracks Model
if (file.exists("./Models/feed_glmmPQL_CRW_CAR1_Hrs.rds")) {
  feed_glmmPQL_CRW_CAR1_Hrs <- readRDS("./Models/feed_glmmPQL_CRW_CAR1_Hrs.rds")
} else { 
  cat(paste0("Could not find ./Models/feed_glmmPQL_CRW_CAR1_Hrs.rds, will need to recreate below."))
}

# Create a dataframe with no missing FTLE vales and columns for grouping and plotting
data_Gamm_ALL <- datasetFTLE %>% 
  # Only include Dives with FTLE data
  dplyr::filter(is.nan(ftle48)==0, is.na(ftle48)==0) %>% 
  group_by(depid) %>% 
  # Add a time Variable for the AR1 correlation structure 
  mutate(time = startI-min(startI),
         timeHrs = time/3600,
  # Add a binary Foraging field for Logistic Regression 
         Foraging = if_else(LungeCount > 0,1,0),
  # Add a hourly feeding rate for use in HMM
         total_Duration = dive_Duration + surf_Duration,
         feedingRateHR = LungeCount/total_Duration*60) %>% 
  ungroup()

# Create a feeding and non-feeding DF for plotting
bwFeeding <- data_Gamm_ALL %>% dplyr::filter(LungeCount>0) %>% 
  mutate(Distribution = "Feeding") %>% 
   group_by(depid) %>% 
  mutate(weight = 1/(n()*length(unique(data_Gamm_ALL$depid)))) %>% 
  ungroup()
bwNoFeeding <- data_Gamm_ALL %>% dplyr::filter(LungeCount==0) %>% mutate(Distribution = "Not Feeding")
# Create an overall feeding dataframe for calculating feeding rates and feeding rate percentiles
bwFeedingOverall <- diveDataset %>% 
  dplyr::filter(LungeCount >0) %>%  
  mutate(Distribution = "Feeding") %>% 
  group_by(depid) %>% 
  mutate(weight = 1/(n()*length(unique(data_Gamm_ALL$depid))),
         # Add a hourly feeding rate for use in HMM
         total_Duration = dive_Duration + surf_Duration,
         feedingRateHR = LungeCount/total_Duration*60) %>% 
  ungroup()

## Process the background FTLE data
data_GammBackground <- datasetBackgroundFTLE %>% 
  dplyr::filter(depid %in% unique(data_Gamm_ALL$depid), 
                Integration == "ftle48") %>%
  rename(ftle48 = FTLE) %>% 
  left_join(dplyr::select(data_Gamm_ALL, depid, DiveNum,Foraging,LungeCount,dttzStart)) %>% 
  dplyr::filter(!is.na(dttzStart)) %>%  # only include records in data_Gamm
  mutate(method = Distribution) %>%
  arrange(depid,dttzStart) %>% 
  group_by(depid) %>% 
  mutate(time = as.numeric(dttzStart)-min(as.numeric(dttzStart)),
         timeHrs = time/3600) %>% 
  ungroup() %>% 
  group_by(time) %>% 
  mutate(FTLE_ID = paste0(depid, seq(1,n()),sep="-B")) %>% 
  ungroup()
  
## Process the Simulated Track Data
data_GammSim <- datasetSimsFTLE %>% 
  # dplyr::filter(FTLE_ID %in% unique(data_Gamm_ALL$depid)) %>% 
  mutate(method = if_else(method=="BrownianBridge","CRW",method),
         Foraging = if_else(LungeCount > 0,1,0)) %>% 
  dplyr::filter(is.nan(ftle48)==0, is.na(ftle48)==0) %>% 
  arrange(method,depid,dttzStart) %>% 
  group_by(depid) %>% 
  # Add a time Variable for the AR1 correlation structure
  mutate(time = as.numeric(dttzStart)-min(as.numeric(dttzStart)),
         timeHrs = time/3600) %>% 
  ungroup()

data_GammCRW <- data_GammSim %>% 
  dplyr::filter(method=="CRW") %>% 
  group_by(depid) %>% 
  mutate(SimNum = as.numeric(unlist(str_split(depid, pattern = "_SimmCRW-"))[2])) %>% 
  ungroup()

data_GammPlus24 <- data_GammSim %>% 
  dplyr::filter(method=="_plus24hrs") %>%
  mutate(method = "plus24hrs") %>% 
  ungroup()
data_GammPlus48 <- data_GammSim %>% 
  dplyr::filter(method=="_plus48hrs") %>%
  mutate(method = "plus48hrs") %>% 
  ungroup()
data_GammPlus96 <- data_GammSim %>% 
  dplyr::filter(method=="_plus96hrs") %>%
  mutate(method = "plus96hrs") %>% 
  ungroup()
data_GammPlus192 <- data_GammSim %>% 
  dplyr::filter(method=="_plus192hrs") %>%
  mutate(method = "plus192hrs") %>% 
  ungroup()

# Combine all data for Hypothesis testing
data_H_Tests <- data_Gamm_ALL %>% 
  dplyr::select(depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long) %>% 
  mutate(method = "Real Data",
         deployment = depid) %>% 
# add CRW
rbind(.,dplyr::select(data_GammCRW,depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long,method,deployment=FTLE_ID)) %>% 
# add plus24
rbind(.,dplyr::select(data_GammPlus24,depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long,method,deployment=FTLE_ID)) %>% 
  # add plus48
rbind(.,dplyr::select(data_GammPlus48,depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long,method,deployment=FTLE_ID)) %>% 
  # add plus96
rbind(.,dplyr::select(data_GammPlus96,depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long,method,deployment=FTLE_ID)) %>% 
  # add plus192
rbind(.,dplyr::select(data_GammPlus192,depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long,method,deployment=FTLE_ID)) %>% 
# add Background Points
rbind(.,dplyr::select(data_GammBackground,depid,Foraging,LungeCount,ftle48,time,timeHrs,dttzStart,Lat,Long,method,deployment=FTLE_ID))

# saveRDS(data_H_Tests, file = "./rmdFiles/data_H_Tests.rds")

Overview

Marine predators face the challenge of reliably finding prey that is patchily distributed in space and time. Predators make decisions at multiple hierarchically organized spatial and temporal scales, yet we have limited understanding of how habitat selection at multiple scales translates into foraging performance. There is mounting evidence that submesoscale (i.e., <100 km) ocean processes drive the formation of dense prey patches that should hypothetically provide feeding hot spots and increase predator foraging success. Here we integrate environmental remote-sensing with high-resolution (≥32 Hz) animal-borne biologging data to evaluate submesoscale surface current features in relation to the dive-scale foraging performance of blue whales in the California Current System.

Methods

Blue Whale Data

Between 2016 and 2020, we deployed 23 high-resolution digital tags on blue whales in the California Current System (CCS). To assess the relationship between blue whale feeding and sub-mesoscale environmental features, we selected a subset of deployments that sampled data for > 24hours, collected GPS positions for > 66% of dives, and overlapped with the sampling footprint of the coastal HF Radar network. Several individuals initiated a southbound migration during the deployment and due to the shifts in behavior documented in Oestreich et al. 2020 were excluded from this study.

The resulting dataset contained:

d0 <- data.frame('Summary of Deployments' = c("Number of Individuals", "Number of Years Represented",
                        "Total Number of Dives", "Total Number of Dives with GPS Locations",
                        "Total Number of Dives with Corresponding HF Radar Data"),
                 Total = c(length(unique(data_Gamm_ALL$depid)), 
                           length(unique(lubridate::year(data_Gamm_ALL$dttzStart))),
                           length(diveDataset$depid),
                           length(datasetFTLE$depid[!is.na(datasetFTLE$Lat)]),
                           length(data_Gamm_ALL$depid) ))
x<-kable(d0,format='pandoc', align = 'lr', row.names = FALSE, col.names = gsub("[.]", " ", names(d0))) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
 cat(x[3:7], sep="\n")
## Number of Individuals                                        10
## Number of Years Represented                                   5
## Total Number of Dives                                      9605
## Total Number of Dives with GPS Locations                   8162
## Total Number of Dives with Corresponding HF Radar Data     7791
 rm(x,d0)

Deployment Overview

## Note: We use a pooled Standard Deviation for the overall summary across individuals.
ftleSummary <- datasetFTLE %>% 
  group_by(depid) %>% 
  summarize(Region=first(region),
            'Animal ID' = first(animalID),
            divesWithGPS = n(), 
            divesWithFTLE = sum(!is.na(ftle48)))

diveFeedingRates <- diveDataset %>%
  dplyr::filter(LungeCount >0) %>% 
  group_by(depid) %>%
  mutate(# Add a hourly feeding rate for use in HMM
         total_Duration = dive_Duration + surf_Duration,
         feedingRateHR = LungeCount/total_Duration*60) %>% 
  summarize('Mean\nDive-By-Dive\nFeeding Rate' = round(mean(feedingRateHR),1),
            'sd\nDive-By-Dive\nFeeding Rate'= round(sd(feedingRateHR),1))

overallSummary <- diveDataset %>% 
  mutate(Year = lubridate::year(dttzStart),
         yday= yday(dttzStart), 
         hr = hour(dttzStart),
         Foraging = if_else(LungeCount > 0,1,0)) %>% 
  group_by(depid) %>% 
  mutate(totalDives = n(),
         DeploymentLength = round(difftime(last(dtStart),first(dtStart),units = "days"),1),
         percentFeeding = sum(Foraging)/n()) %>% 
  ungroup() %>% 
  group_by(depid,yday) %>% 
  summarize(DeploymentLength = first(DeploymentLength),
            totalDives = first(totalDives),
            percentFeeding = first(percentFeeding),
            dailyLunges = sum(LungeCount),
            firstHour = first(hr),
            lastHour = last(hr),
            Year = first(Year)) %>% 
  slice(-1,-n()) %>%
  group_by(depid) %>% 
  summarize(Year = first(Year),
            DeploymentLength = first(DeploymentLength),
            totalDives = first(totalDives),
            percentFeeding = round(first(percentFeeding),2),
            'Mean\nDaily\nLunges' = round(mean(dailyLunges),0),
            'sd\nDaily\nLunges' = round(sd(dailyLunges),0),
            'Mean\nOverall Daily\nFeeding Rate' =  round(mean(dailyLunges/24),1),
            'sd\nOverall Daily\nFeeding Rate' = round(sd(dailyLunges/24),1)) %>% 
  left_join(ftleSummary,by="depid") %>% 
  left_join(diveFeedingRates,by="depid") %>% 
  dplyr::select("Deployment\nID"=depid, Region, `Animal ID`, Year, 
                'Deployment\nLength\n(days)' = DeploymentLength,
                'Total #\nDives' = totalDives,'Total #\nDives\nwith GPS'=divesWithGPS,
                'Total #\nDives\nwith FTLE'=divesWithFTLE,
                '%\nDives\nFeeding'= percentFeeding,everything())

# Summary dataframe for the calculation of the pooled std dev of  
# dive feeding rates 
frSummary <- bwFeedingOverall %>% 
  group_by(depid) %>% 
  summarize(meanFR = mean(feedingRateHR),
            sdFR = sd(feedingRateHR),
            sampleSize = n())
# Summary dataframe for the calculation of the pooled std dev of 
# daily feeding rates 
dailyFRSummary <- diveDataset %>% 
  mutate(yday= yday(dttzStart)) %>% 
  group_by(depid,yday) %>% 
  summarize(dailyLunges = sum(LungeCount)) %>% 
  slice(-1,-n()) %>%
  group_by(depid) %>% 
  summarize(meanDailyLunges = round(mean(dailyLunges),0),
            sdDailyLunges = round(sd(dailyLunges),0),
            meanOverallDailyFeedingRate =  round(mean(dailyLunges/24),1),
            sdOverallDailyFeedingRate = round(sd(dailyLunges/24),2) ,
            numDays = n())

totals <- overallSummary %>% 
  summarize("Total",
            "",# Region
            "  ",# Animal ID
            length(unique(Year)),
            sum(`Deployment\nLength\n(days)`),
            sum(`Total #\nDives`),
            sum(`Total #\nDives\nwith GPS`),
            sum(`Total #\nDives\nwith FTLE`),
            round(mean(`%\nDives\nFeeding`),2),
            round(mean(`Mean\nDaily\nLunges`),0),
            # Daily Lunges Pooled standard dev
            sdDLunges=round(pooledstd(stdevs = dailyFRSummary$sdDailyLunges,numSamples =  dailyFRSummary$numDays),0),
            round(mean(`Mean\nOverall Daily\nFeeding Rate`),1),
           # Daily Feeding Rate Pooled standard dev
            sdDFR=round(pooledstd(stdevs = dailyFRSummary$sdOverallDailyFeedingRate,numSamples = dailyFRSummary$numDays),1),
            round(mean(`Mean\nDive-By-Dive\nFeeding Rate`),3),
            # Dive by Dive Feeding Rate Pooled standard dev
            sdDbDFR=round(pooledstd(stdevs = frSummary$sdFR,numSamples = frSummary$sampleSize),1)
            )
colnames(totals) <- colnames(overallSummary)
d2<-rbind(overallSummary,totals)
kable(d2, align = 'lcccccccccccc', row.names = FALSE) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Deployment ID Region Animal ID Year Deployment Length (days) Total # Dives Total # Dives with GPS Total # Dives with FTLE % Dives Feeding Mean Daily Lunges sd Daily Lunges Mean Overall Daily Feeding Rate sd Overall Daily Feeding Rate Mean Dive-By-Dive Feeding Rate sd Dive-By-Dive Feeding Rate
Bm160523-A20 Central CA CRC-3308 2016 4.1 853 770 770 0.44 140 45 5.8 1.9 24.10 8.5
Bm160716-A20 Southern CA CRC-3090 2016 3.9 462 335 335 0.58 365 50 15.2 2.1 29.20 7.5
Bm160918-A08 Southern CA CRC-3348 2016 4.1 907 791 670 0.65 407 65 16.9 2.7 25.50 9.1
Bm170622-TDR12 Southern CA CRC-0542 2017 17.7 2728 2588 2569 0.53 346 66 14.4 2.8 23.20 5.6
Bm170925-TDR12 Southern CA CRC-3427 2017 4.1 505 456 255 0.18 84 67 3.5 2.8 22.00 9.2
Bm170926-TDR14 Southern CA CRC-3107 2017 4.3 373 259 259 0.60 330 14 13.8 0.6 25.90 7.7
Bm181021-TDR14 Central CA CRC-3520 2018 9.1 1152 908 908 0.32 191 82 8.0 3.4 21.00 5.4
Bm190709-TDR14 Central CA CRC-1789 2019 8.0 731 635 633 0.40 104 101 4.4 4.2 17.60 6.7
Bm190710-TDR15 Northern CA CRC-3697 2019 13.5 1502 1062 1039 0.50 258 136 10.7 5.7 24.40 8.8
Bm201028-TDR17 Central CA CRC-3732 2020 2.6 392 358 353 0.39 331 190 13.8 7.9 24.90 8.4
Total 5 71.4 9605 8162 7791 0.46 256 94 10.7 3.9 23.78 7.4

Map of Deployments

hourlySum <- datasetFTLE %>% 
  mutate(yday = yday(dttzStart),
         hr = hour(dttzStart)) %>% 
  group_by(depid,yday, hr) %>% #View()
  summarize(LungeCount = sum(LungeCount),
            Lat = mean(Lat,na.rm=TRUE),
            Long = mean(Long,na.rm=TRUE))
p<-ggplot() +
  geom_sf(data = world,color=NA,fill=NA) +
  geom_polygon(data=fortify(HFRadarCoverageMCP),
               aes(x=long,y=lat, alpha = name),
               color = NA, fill = "blue") +
  # add 200m contour
  geom_contour(data = bf, 
               aes(x=x, y=y, z=z),
               breaks=c(-200),
               size=c(0.4),
               colour="darkgrey", show.legend = FALSE) +
  geom_text_contour(data = bf, aes(x=x, y=y,z = z),breaks=c(-200), 
                    show.legend = FALSE, size = 2.2, alpha = .6, nudge_y = -.002) +
  # Add whale Locations
  geom_path(data=data_Gamm_ALL, aes(x=Long,y=Lat,group=depid,color=depid),
            alpha = 0.95,size = .9, show.legend = TRUE) +
  geom_point(data=hourlySum[hourlySum$LungeCount>0,],
             aes(x=Long,y=Lat,group = depid,
                 size = LungeCount, fill = depid),
             alpha = 0.75, shape=21,show.legend = TRUE) +
  geom_sf(data = world) +
  geom_point(data=hfRadarStns,
             aes(st_coordinates(hfRadarStns$geometry)[,1],st_coordinates(hfRadarStns$geometry)[,2],
                 shape = name),size=2.5, fill = 'magenta') +
  scale_alpha_manual(values=c("Coverage Area"=.15))+
  scale_shape_manual(values = c("Radar Station" = 23))+
  coord_sf(xlim = c(min(data_Gamm_ALL$Long,na.rm = TRUE)-.05, max(data_Gamm_ALL$Long,na.rm = TRUE)+.05),
           ylim = c(min(data_Gamm_ALL$Lat,na.rm = TRUE)-.05, max(data_Gamm_ALL$Lat,na.rm = TRUE)+.05),
           expand = FALSE) +
  guides(size = guide_legend(order=4, override.aes = list(linetype = 0)),
         alpha = guide_legend(order=2, override.aes = list(linetype = 0,size=0,color=NA)),
         shape = guide_legend(order=1, override.aes = list(linetype = 0)))+
  labs(size = "Hourly Lunge Count",alpha=NULL, shape="High Frequency Radar", color="Deployment",fill="Deployment") +
  xlab("Longitude") + 
  ylab("Latitude") + 
  theme(axis.title = element_text(family="Times",face="bold", size=20),
        axis.text = element_text(family="Times", face="bold", size=18),
        panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), 
        panel.background = element_rect(fill = "aliceblue"),
        legend.text = element_text(face="bold", size=18),
        legend.title = element_text(face="bold", size=20),
        plot.title = element_text(face="bold", size=24,hjust = 0.5),
        legend.position= c(.8,.7),
        legend.key.size = unit(.5, "cm"),
        legend.box.background = element_rect(color="white", size=1,fill="white"),
        # legend.box.margin = margin(6, 6, 6, 6),
        legend.box = "vertical", legend.key = element_blank(),
        legend.margin = margin(0.25,0.15,0.25,0.15, unit="cm"),
        legend.spacing.y = unit(.1, "cm"))
p