Usage
For the full theoretical description, refer to dead-reckoned track – technical explanation.
Inputs
A:
is a sensor data structure or annx3
acceleration matrix with columns[ax ay az]
, or acceleration sensor list (in R).A
can be in any consistent unit, e.g., g or m/s^2.M
: is a sensor data structure or annx3
magnetometer matrix with columns [mx,my,mz] or magnetometer sensor list (in R). M can be in any consistent unit (e.g., in uT or Gauss) A and M must have the same size (and so are both measured at the same sampling rate).s
: is the forward speed of the animal in m/s. s can be a single number meaning that the animal is assumed to travel at a constant speed. s can also be a vector with the same number of rows as A and M, e.g., generated by ocdr().sampling_rate:
The sampling rate of the sensor data in Hz (samples per second). This input will be ignored if A and/or M are sensor structures or sensor lists (in R), in which case the sampling rate will be extracted from them-
fc:
(optional) The cut-off frequency in Hz of a low-pass filter to apply toA
and M before computing bodypointing angle. The filter cut-off frequency is in Hz. The filter length is 4*sampling_rate/fc. Filtering adds no group delay. If fc is empty or not given, the default value of 0.2 Hz (i.e., a 5 second time constant) is used.
Outputs
T:
is the estimated track in a local level frame. The track is defined as meters of northward and eastward movement (termed ‘northing’ and ‘easting’, i.e, T=[northing,easting]) relative to the animal’s position at the start of the measurements (which is defined as [0,0]). The track sampling rate is the same as for the input data and so each row of T defines the track coordinates at times 0,1/fs,2/fs,… relative to the start time of the measurements.pe:
is the estimated depth or altitude predicted from the speed and pitch angle. This can be compared against the measured depth/altitude to assess errors in the dead-reckoned track. Note that even if pe matches the observed depth, this does not guarantee that the track is accurate.
Theoretical Description
Dead reckoning is the process of deducing a track from periodic estimates of heading and speed. For the full theoretical description, refer to dead-reckoned track – technical explanation.
Technical Description
In this vignette, you will learn how to estimate a simple dead-reckoned track (pseudo-track) based on speed and body pointing angle using the ptrack()
function. For the full technical description, refer to dead-reckoned track – technical explanation.
Caveats
The big disadvantage of dead-reckoning is that errors accumulate over time, resulting in inaccurate positions unless occasional GPS fixes are available to correct the track. For the full caveats description, refer to dead-reckoned track – technical explanation.
This practical covers how to compute a dead-reckoned track from accelerometer and magnetometer data. If you have GPS data and would like to merge the GPS and dead-reckoned tracks to produce a high-resolution track with lower error and try visualizing the track coloured by another parameter, this vignette will also show you how to do it.
Load and Visualize the Test Data Set
Load the test dataset testset7
, which belongs to the data recorded from a suction cup tag attached to the back of a humpback whale [Note]This dataset has already been converted from a source file that was offloaded from the tag into a NetCDF file. In doing so, some metadata was gleaned from the file and added to the data. Other metadata was added by hand.. This dataset is built into the tagtools package, so you can access it using system.file
.
hw_file_path = xxx
hw = load_nc(hw_file_path)
fs = hw.P.sampling_rate
hw_file_path<- system.file("extdata", "testset7.nc", package = 'tagtools', mustWork = TRUE)
hw <- load_nc(hw_file_path)
fs <- hw$P$sampling_rate
Then, check the contents of the data structures to answer the following questions.
What is the sampling rate of the accelerometer data?
1 Hz.
What processing steps have already been applied to the magnetometer data?
sens_struct
, do_cal
, tag2animal
, crop_to
, decdc(5)
… So, it has been converted to a sensor structure, calibrated, put into the animal frame, cropped to a certain segment of data, and decimated by a factor of five (as a result, the sampling rate is now the same as the acceleration data, 1 Hz).
What is in the 3 columns in the POS (GPS position) data?
The three columns are the time, the latitude, and the longitude.
Which frame is the accelerometer data in? Which frame is the magnetometer data in?
The accelerometer data is oriented in the animal frame, and so is the magnetometer data.
Plot the pressure and accelerometer data with plott()
to get a sense for what the animal might be doing in this data segment. Note that the code example below assumes you have called the data set hw
:
plott(hw.P, hw.Aa)
subplot(211)
ylabel('Depth')
subplot(212)
ylabel('Acc')
legend({'x-axis','y-axis','z-axis'})
plott(X = list(Depth = hw$P, Acc = hw$Aa))
Expand to show figure …


Then, plot the GPS positions:
figure
plot(hw.POS.data(:,3),hw.POS.data(:,2))
xlabel('Latitude')
ylabel('Longitude')
plot(hw$POS$data[,3],hw$POS$data[,2])
Expand to show figure …

Dead-reckoning
What is Dead-reckoning?
Hearken back to your earliest trigonometry course. Chances are, you used angles and distances to calculate other distances (on a right triangle, especially). Then, in calculus, you’ve probably used angles and rates of change to calculate distance or other rates of change. Dead-reckoning is essentially this—using just a heading and a forward speed in order to plot a forward path. The technique is fairly commonplace in sailing and aviation, and can be made accurate with the help of known position data. If you want more background, consult this video for an example with some simple trigonometry.
However, per Wikipedia (which simply says it very well), dead-reckoning “is subject to significant errors of approximation. For precise positional information, both speed and direction must be accurately known at all times during travel. Most notably, dead reckoning does not account for directional drift during travel through a fluid medium. These errors tend to compound themselves over greater distances”. We will see how this is true over the course of this vignette—water is seldom entirely still, and there is quite a current present in the data we will soon investigate.
Why use Dead-reckoning?
The GPS plot (plotted above) shows a mix of intensive and extensive movements, but the constraint of only getting positions when the animal is at the surface means we cannot infer much about the movement behaviour within dives. Dead-reckoning from accelerometer and magnetometer data is the only way to estimate movement within dives without requiring external tracking infrastructure.
Estimating Animal Speed
Dead-reckoning uses the accelerometer and magnetometer to calculate the direction of travel. In what frame do these data need to be?
They both need to be in the animal frame. Thankfully, they both already are!
An estimate of the forward speed is also required. We don’t have a speed sensor—doing this would actually be quite cumbersome, as compared to an acceleration sensor. Nevertheless, we can compute the vertical speed (i.e., the differential of the depth) during descents and ascents, which might be a good starting guess (for a detailed explanation on how depth_rate() works, refer to Estimate vertical speed and forward speed tutorial):
v = depth_rate(hw.P)
figure
plott(hw.P,v,fs)
%add some labels
subplot(211)
ylabel('Depth m')
subplot(212)
ylabel('Vertical speed, m/s')
v <- depth_rate(hw$P)
plott(X = list('Depth m)' = hw$P, 'Vertical speed, m/s'= v),fsx = fs)
Expand to show figure …

Zoom in to individual dives on this plot to get a rough idea of the descent and ascent speed of the whale. Here we are zooming to see the first dive in detail.
figure
plott(hw.P,fs,v,fs)
xlim([0 6])% zoom in the first 6 minutes
%add some labels
subplot(211)
ylabel('Depth m')
subplot(212)
ylabel('Vertical speed, m/s')
#Set interactive = TRUE and zoom in to individual dives on this plot to get a rough idea of the descent and ascent speed of the whale.
plott(X = list('Depth m)' = hw$P, 'Vertical speed, m/s'= v),fsx = fs, interactive = TRUE)
#Or, if the interactive figure drives you a bit crazy, just set the axis limits:
plott(X = list('Depth m)' = hw$P, 'Vertical speed, m/s'= v), fsx = fs, xl = c(0,0.1)
Expand to show figure …

Which is roughly the estimated speed during descents and ascents?
2 m/s
Computing a Dead-Reckoned Track
Call your speed estimate speed spd
and use it in the following line to compute the dead-reckoned track:
figure
spd=2
DR = ptrack(hw.Aa, hw.Ma,spd)
plot(DR(:,2),DR(:,1))
xlabel('Easting, m')
ylabel('Northing, m')
spd <- #???? # fill in your estimate here, we used 2 for the example
DR <- ptrack(A = hw$Aa, M = hw$Ma, s = spd)
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track')
Expand to show figure …

A dead-reckoned track is a series of distances north (or south if negative) and east (or west if negative) of the starting point which is position (0, 0) in the plot. The first two columns of DR are these ‘northing’ and ‘easting’ values. The DR track is defined in a Local Level Frame (LLF) as opposed to latitude and longitude. An LLF is like a map—a region that is small enough that we can assume the earth is flat—centered on the starting point.
How does the spatial resolution of the dead-reckoned track look compared to the surface GPS positions?
At first glance, the shape is quite similar, so this might not concern us yet. Also, the spatial resolution looks much, much better! So, making this track seems like a good idea.
Adding GPS positions
To plot the GPS positions on the same plot, we need to first convert them into the same LLF. The first GPS position is only 0.8s into the dataset (this is hw.POS.data(1,1)
) so we can say that the dead-reckoned track starts from this point. To convert latitude and longitude into the LLF use:
figure
POSLLF = lalo2llf(hw.POS.data(:,2:3));
plot(DR(:,2),DR(:,1))
hold on
plot(POSLLF(:,2),POSLLF(:,1),'r.-')
xlabel('Easting, m')
ylabel('Northing, m')
POSLLF <- lalo2llf(trk = hw$POS$data[,c(2:3)])
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track',
yl = c(-1000, 5000))
lines(POSLLF$easting, POSLLF$northing, type = 'b', col = 'red', pch = 20)
Expand to show figure …


How well does the dead-reckoned track match up to the GPS positions?
They do not match up much at all. They touch at the beginning… and that’s about it. After that, the shapes are somewhat preserved, but the predicted track is way farther north, and quite a ways farther east, than the actual known positions.
The dead-reckoned track is always computed with respect to the water, not the ground, but we are plotting it here with respect to the ground. A more accurate track would take into account the water current.
Can you imagine what current direction would be needed to make the dead-reckoned track more closely match the GPS positions?
The current would need to be moving quite a bit south, and also some distance west, in order to compensate for the large errors we are seeing.
There are a number of ways to combine the GPS positions and the dead-reckoned track into a single track which has both the absolute accuracy of GPS and the high temporal resolution of dead-reckoning. A simple method is to force the dead-reckoned track to meet the GPS positions at each surfacing by adding a constant ‘current’ to each track point in the preceding dive. This is done by fit_tracks
():
[TRK,C] = fit_tracks(POSLLF,hw.POS.data(:,1),DR(:,1:2),fs);
POSLLF = lalo2llf(hw.POS.data(:,2:3));
plot(DR(:,2),DR(:,1))
hold on
plot(POSLLF(:,2),POSLLF(:,1),'r.-')
hold on
plot(TRK(:,2),TRK(:,1),'g') % add the fitted track to your plotc
xlabel('Easting, m')
ylabel('Northing, m')
TRK <- fit_tracks(P= POSLLF,T=hw$POS$data[,1],D = DR[,1:2],sampling_rate=fs)
# add to plot
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track',
yl = c(-1000, 5000))
lines(POSLLF$easting, POSLLF$northing, type = 'b', col = 'red', pch = 20)
lines(TRK$easting, TRK$northing, col = 'darkgreen')
Expand to show figure …


Interpreting the tracks
Zoom into the plot to see how the green fitted track interpolates the red GPS positions.
If the green track is to be believed, how effectively do the surface positions capture the movement of the animal?
They don’t capture the movement of the animal very well. So much happens between the positions that we’re able to collect at the surface, so we are missing out if surface positions is all we collect/show
We often want to plot a tag-derived variable as a function of position, and a coloured track plot is one way to do this. For example, it can be useful to see where the animal is diving along the track. To colour the track by depth, use col_line with P as the colour information:
figure % open a new figure window
col_line(TRK(:,2),TRK(:,1),hw.P.data)
hcb=colorbar
title(hcb,'Depth, m')
xlabel('Easting, m')
ylabel('Northing, m')
set(gcf,'color','white')
col_line(TRK[,2], TRK[,2], c = hw$P$data)
Expand to show figure …

Coming soon …
Zoom in to see what the scale is of the track tortuosity: Is there tortuousity within individual dives, or is the tortuousity occuring across dives?
While most of the tortuosity occurs across dives, there is some tortuosity within dives as well.
If you are done and want a challenge, try colouring the track by the absolute roll angle instead of the depth. Refer to the orientation estimation tutorial for a deep understanding on how a2pr() function works. You can use scatter3
in MATLAB or col_line3
in R to plot the 3-d positions (i.e., TRK(:,2)
, TRK(:,1)
, and P.data
), and then color the plot by absolute roll angle remembering to convert from radians to degrees.
figure
[pitch roll] = a2pr(hw.Aa)
roll_deg = roll/pi*180
scatter3(TRK(:,2), TRK(:,1), hw.P.data,[], roll_deg)
set( gca, 'ZDir', 'reverse' )
colorbar
hcb=colorbar
title(hcb,'roll, deg')
xlabel('Easting, m')
ylabel('Northing, m')
zlabel('Depth, m')
pitch_roll <- a2pr(hw$Aa)
roll_deg <- pitch_roll$r/pi*180
col_line3(x = TRK[,2], y =TRK[,1],
z = hw$P$data, c = roll_deg)
Expand to show figure …


Further Resources
For a more in-depth explanation of the usage and inner workings of the ptrack()
function, refer to the Technical Explanation of a dead-reckoned track.