Find Marine Animal Dives Below the Water’s Surface and Summarize Them

Usage

In this vignete we will learn how to use the function find_dives to gain insight about a depth profile, calculate the mean dive duration and mean dive depth, mark the beginnings and endings of dives, as well as summarize dives using the function dive_stats().

MATLAB
D = find_dives(p,sampling_rate,mindepth,surface,findall)
R
D <- find_dives(p, mindepth, sampling_rate = NULL, surface = 1, findall = 0)

The order of inputs is sligtly different in R and Matlab

Inputs
  • p: is s sensor list or vector of depth/altitude in meters.
  • sampling_rate: the sampling rate of the sensor data in Hz (samples per second). This is only needed if p and t are not a sensor lists. The depth and temperature must both have the same sampling rate (use decdc() if needed to achieve this).
  • mindepth: the threshold in meters at which to recognize a dive or flight. Dives shallow or flights lower than mindepth will be ignored.
  • surface: (optional) the threshold in meters at which the animal is presumed to have reached the surface. Default value is 1. A smaller value can be used if the dive/altitude data are very accurate and you need to detect shallow dives/flights.
  • findall: (optional) when 1 forces the algorithm to include incomplete dives or flights at the start and end of the record. Default is 0 which only recognizes complete dives/flights.
Outputs
  • D: is a structure array with size equal to the number of dives/flights found. The fields of D are:
    • start time in seconds of the start of each dive/flight
    • end time in seconds of the start of each dive/flight
    • max maximum depth/altitude reached in each dive/flight
    • tmax time in seconds at which the animal reaches the max depth/altitude
Notes
  • If there are n dives/flights beyond mindepth in p, then D will be a structure containing n-element vectors, e.g., to access the start time of the kth dive/flight use D.start(k).
Theoretical Description

Diving and flying behaviour of animals is frequently studied in ecology to better understand foraging strategies, diving physiology, biomechanics… The main studied characteristics are mean and max dive/flight duration, mean and max dive/flight depth/attitude, surface interval (inter-dive interval or inter-flight interval), along with other parameters such as speed, and distance covered.

Technical Description

In this vignette, you will learn how to calculate the mean and max dive/flight duration and mean and max dive/flight depth/attitude, surface interval, mark the beginnings and endings of dives.

Verification

You could simply verify that find_dives works by simply plotting the the original depth, p and select by hand with ginput() the beginnings and endings of dives/flights, the deepest and highest point of each dive and flight, respectively. ginput() is a useful function in Matlab and Octave for measuring data on a plot – there isn’t a great equivalent in R, where interactive plots are not really commonly used. (How would this work for a flying animal? how could you check that is working? for a seabird that rests in the water would be the same but what for a flying animal that rests at different altitudes)

Caveats

For the find_dives function to work perfectly, you will need to remove any off-animal tag data before running it on the data.


Load and visualize the test dataset

Load the test dataset zc11_267a.nc which belongs to the data recorded from a suction cup tag attached to the back of a humpback whale [Note]. This dataset is built into the tagtools package, so you can access it using system.file.

MATLAB
zc_file_path = xxx
zc = load_nc(zc_file_path)
fs = zc.P.sampling_rate
R
zc_file_path<- system.file("extdata", "XXX", package = 'tagtools', mustWork = TRUE)
zc <- load_nc(zc_file_path)
 
fs <- zc$P$sampling_rate

Then use plott() to inspect the depth/(altiude in case of a flying animal):

MATLAB
figure
plott(zc.P)
ylabel('Depth m')
xlabel('Time hours')
R
plott(zc$P)
MATLAB Output
R Output

What do you notice?

There are some very deep dives in this profile! This is perhaps to be expected for Cuvier’s beaked whales, with their reputation for mysteriously extreme excursions.

As with all raw depth data, there are some problems with this dive profile. To know how to remove depth outliers, correct pressrue offsets and temperature effects, as well as how to remove periods of data when the tag is not on the animal, please see the Depth checks and corrections tutorial.

Here, we are going to crop the data before further analysis, because there is a period at the start of the recording when the tag was not yet deployed on the whale. To verify that the cropping has been done correct, plot the depth data after cropping it.

MATLAB
zc.Pc = crop(zc.P)
plott(zc.Pc)
ylabel('Depth m')
xlabel('Time hours')
R
zcCr = crop(zc$P)# if you want to review history/other fields
str(zcCr, max.level = 1)
plott(zcCr)
Expand to show figure …
MATLAB Output
R Output

What minimum depth threshold do you think you would use to detect dives this animal’s dives? Consider how you would justify your choice.

No less than fifty meters seems like a reasonable choice to detect this animal’s dives. After all, most of the dives this animal performs are more than 100 meters in depth.

Exercise: calculate the mean duration of dives deeper than 50m

Our goal with these data is to calculate the mean duration of dives deeper than 50 m. If you can think of a way to do this already, please go ahead and try! You can then compare your answer to the step-by-step procedure below.

You could measure each dive by hand on the depth plot (ginput is a useful function in Matlab and Octave for measuring data on a plot – there isn’t a great equivalent in R, where interactive plots are not really commonly used). But there is a toolbox function for this called find_dives(). See the help on this function to find out what it does and what options it has.

To find dives deeper than 50 m in your compensated dive data, type:

MATLAB
help find_dives

mindepth=50
surface=3
d = find_dives(zc.Pc,mindepth,surface)
R
? find_dives
d <- find_dives(Pc$p,mindepth=50,surface=3)

d should return a data frame with the start, end, and maximum depth of 29 dives. Note that we have added surface equals to 3 meters in the find_dives function. This is because the surface period for the first 4 dives is not accurate, it is shifted 2 meters, and so if not surface is not defined, d would have returned a data frame of 25 dives. How can you get the mean dive duration & mean (maximum) dive depth from this structure? Code below provides the answer.

MATLAB
total_dive_duration = d.end-d.start
mean_dive_duration = mean(total_dive_duration)
mean_dive_duration % since this is using data from find_dives, the mean duration is in seconds. If you want it in minutes then type
mean(total_dive_duration)/60

% and mean depth
total_dive_depth = d.max
mean_dive_depth = mean(total_dive_depth)
mean_dive_depth
R
total_dive_duration <- matrix(0)
for(n in 1:nrow(d)) {
  total_dive_duration <- total_dive_duration + d[n,2] - d[n,1]
  }
mean_dive_duration <- total_dive_duration/nrow(d)
mean_dive_duration
# since this is using data from find_dives, the mean duration is in seconds
# and mean depth
total_dive_depth <- matrix(0)
for(n in 1:nrow(d)) {
  total_dive_depth <- total_dive_depth + d[n,3]
}
mean_dive_depth <- total_dive_depth/nrow(d)
mean_dive_depth

When you have got the mean dive depth, try plotting the start and end of the dives on the depth plot:

MATLAB
plott(zc.Pc)
hold on 

figure 
plott(zc.Pc)
hold on 
plot(d.start/3600, 0,'g^')
plot(d.end/3600, 0,'r^')
R
plott(X=list(Pcmf=Pcmf$p), r=TRUE)
points(d$start/(3600*24),rep(0,nrow(d)),col='green', pch=19)
points(d$end/(3600*24),rep(0,nrow(d)),col='red', pch=17)
MATLAB Output
R Output

Again, this plot might be rather small. As a result it might be tricky to make sense of the markers that designate the starts and ends of dives. You should be able to click Zoom to view the plot in a larger window.

Note: if you cropped the time such that the units of the x-axis are not in days, you will have to adjust the multipliers in the points code accordingly. In the example above, the start and end times returned by find_dives are in seconds so we needed to divide them by 3600 to match the unit (days) automatically selected for time by plott

Use dive_stats

Now use dive_stats() to compute summary statistics for all the dives you detected.

In addition to the maximum excursion and duration, dive_stats divides each excursion into three phases: “to” (descent for dives, ascent for flights), “from” (ascent for dives, descent for flights), and “destination“. The “destination” (bottom for dives and top for flights) phase of the excursion is identified using a “proportion of maximum depth/altitude” method, whereby for example the bottom phase of a dive lasts from the first to the last time the depth exceeds a stated proportion of the maximum depth. Average vertical velocity is computed for the to and from phases using a simple method: total depth/altitude change divided by total time. If an angular data variable is also supplied (for example, pitch, roll or heading), then the circular mean and variance are also computed for each dive phase and the dive as a whole. Could anyone add more info in this regard here… why is this important or useful for?. As I do not understand this, I have not add the optional part regarding calculating dive_stats including pitch, as is done in https://animaltags.github.io/tagtools_r/articles/dive-stats.html.

MATLAB
dive_cues = [d.start,d.end];% start and end of dive in seconds since the start of the recording
ds = dive_stats(zc.Pc, dive_cues) 
ds
R
ds <- dive_stats(ZPCr, dive_cues=dt[,c('start', 'end'),]) 
str(ds)

Have a look at the dive stats and perhaps make a plot of some or all of them.

Do you notice anything interesting?

If you plot duration (ds.dur) vs. the maximum depth (ds.maxz) for each dive you would be albe to see how the diving behaviour of beaked whales is bimodal, with longest dives being also the deepest.

The destination phase from the output above has been calculated from the first to the last time the depth of each excursion exceeds 85% of the maximum depth, as a default. You could change this proportion and compare your results with already published results from the same or differnt species that uses this method to identify and clasify the destination or bottom phase.

Lets set the proportion to 70% and compare the dive_stats results with the ones obtained above:

MATLAB
dive_cues = [d.start,d.end];% start and end of dive in seconds since the start of the recording

prop=0.70
ds_70 = dive_stats(zc.Pc, dive_cues,[],[],0.70) 
ds_70 
R
ds_70 <- dive_stats(ZPCr, dive_cues=dt[,c('start', 'end'),],[],[],prop=0.70) 
str(ds)

Lets plot the dive profile coloring the start and end of the bottom phase for both scenarios to compare them visually.

MATLAB
figure
plott(zc.Pc)
hold on
start_085=(ds.dest_st+d.start)/3600
end_085=(ds.dest_et+d.start)/3600
start_070=(ds_70.dest_st+d.start)/3600
end_070=(ds_70.dest_et+d.start)/3600
plot(start_085, 0,'go')%start of destination phase with 0.85 proportion in green circles. 
hold on 
plot(end_085,0,'g*')%end of destination phase with 0.85 proportion in green asterix.
hold on 
plot(start_070,0,'ro')%start of destination phase with 0.85 proportion in red circles. 
hold on 
plot(end_070,0,'r*')%end of destination phase with 0.85 proportion in red asterix. 
ylabel('Depth, m')
xlabel('Time, hours')
R

Expand to show figure …
MATLAB Output
R Output

Further Resources

to continue working through these practicals, consider rotation-test or fine-scale-tracking ( these two pages are not done yet!) .