Blue Whales and Lagrangian Features
Methods and Analysis
James A. Fahlbusch, Max F. Czapanskiy, John Calambokidis, David E. Cade, Briana Abrahms, Elliott L. Hazen, Jeremy A. Goldbogen
18 July, 2022
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
rm(hourlySum)Whale Data Processing
Tags used in this study were Wildlife Computers TDR10-F (n=7) and Acousonde acoustic tags (n=3). All tags sampled depth, accelerometry (>= 32Hz) and GPS. Tags were deployed from a 6-7m Rigid Hull Inflatable Boat using a 4m carbon fiber pole. Tags were attached to the animal with 3 or 4 stainless steel darts that were 6cm long with 1-2 rows of petals (Szesciorka et al. 2016, Calambokidis et al. 2019). Tags were recovered and the data downloaded after detaching from the animal.
Raw tag data was pre-processed in Matlab using custom scripts (following Cade et al 2016) to calculate the animal’s pitch, roll, and speed from the accelerometer sensor data. Individual feeding events (i.e. Lunges) were identified manually using stereotypical signatures in kinematic data. Raw GPS location data was filtered for unrealistic whale speeds (e.g. > 6 meters/second) using the ArgosFilter package (Version 0.62) in R (Version 3.5.1). All processed data was down-sampled to 1Hz for analysis.
Dives were identified as excursions to a depth of greater than 10 meters. Locations were interpolated for dives with a start time between 2 known positions less than 15 min from the dive’s start time. Only dives with a corresponding GPS location were included in this analysis.
Environmental Data
HF Radar Data
For each deployment, surface current data were downloaded from the Southern California Coastal Ocean Observing System (SCCOOS, http://www.sccoos.org/data/hfrnet/) at hourly resolution (cells 6 km on a side) for the deployment period \(\pm\) 7 days with a bounding box of \(\pm\) 1 degree around the deployment locations.
Data gaps in the raw HF radar surface current measurements were restored using the algorithms described in Ameli and Shadden (2019), selecting a concave hull (alpha shape radius of 10 km) with the detection and exclusion of land.
FTLE Calculation
Using the restored surface current data, we calculated the backward-in-time FTLE using TRACE (http://transport.me.berkeley.edu/trace/), a Lagrangian analysis tool that follows the methodology described in Shadden et al. (2005; 2009). We used an evenly-spaced grid of tracers having 10 times the spatial resolution of the HF radar data (sensu Shadden et al., 2009), and applied a free-slip boundary condition to tracers near land. Tracer advection used a bilinear spatial interpolation and an adaptive 4th order Runge-Kutta-Fehlberg integration method. At every hourly time-step, the trajectories of the evenly spaced grid of tracers were integrated for the preceding 48 hour period, and FTLE was calculated from the time-dependent movement of tracer trajectories.
- TRACE Parameters
- Initial Time: 1 day prior to Tag On Animal
- Number of Initial Times: Number of hours in the deployment + 48
- this ensures the data will envelop the deployment by $$1 day
- Initial Times Increment: 1 Hour
- Integration Duration: 48 Hours
- Reverse Direction of Time
- Number of Tracers: FTLE is calculated with a grid of tracers having
10x the resolution of the input data.
- e.g. the 6km HF radar data for a deployment may have a grid 39 X 49 cells (Long X Lat). The number of tracers for this deployment would be 390 X 490.
- TRACE Output
- The final resolution varies slightly by deployment, but is roughly 590m (i.e., 6km HF radar resolution/10).
Overlap with Environmental Data
The sampling unit for this analysis is Blue Whale Dives and are classified as either feeding or non-feeding. Dives also have an associated Feeding Rate, which is the # of lunges per dive minute + post-dive surface period.
FTLE is sampled at hourly intervals, with each FTLE value representing \(\pm\) 30 minutes from the sample timestamp. For each blue whale dive location, we extract the FTLE value from the corresponding FTLE layer.
Statistical Analysis
The movements of free-ranging predators can be described in terms of habitat utilization in which the predator uses a subset of the habitat available. The utilization is considered selective when a predator targets a certain set of features disproportionately to their availability in the environment (Johnson 1980). Within the hierarchical framework described by Johnson 1980, there have been many works describing the lower order selection types of blue whales using satellite tags (e.g. Abrahms et al 2018), which relates to the selection of a home range. Additionally, recent work by Cade et al 2021 has shed light on some of the highest order selection types, such as the selction for the densest prey within a patch (e.g. 4th order).
In this manusript, we are contributing to a knowledge gap of intermediate order selection processes. This includes the 3rd order, which pertains to the usage of habitat within a home range.
The pattern we explore quantitatively is shown in the animation below.
