Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You and Your Team.

Learn More →

Dead-reckoning animal movements in R: a reappraisal using Gundog.Tracks

Dead-reckoning animal movements in R: a reappraisal using Gundog.Tracks Background: Fine‑scale data on animal position are increasingly enabling us to understand the details of animal movement ecology and dead‑reckoning, a technique integrating motion sensor ‑ derived information on heading and speed, can be used to reconstruct fine ‑scale movement paths at sub ‑second resolution, irrespective of the environ‑ ment. On its own however, the dead‑reckoning process is prone to cumulative errors, so that position estimates quickly become uncoupled from true location. Periodic ground‑truthing with aligned location data (e.g., from global positioning technology) can correct for this drift between Verified Positions ( VPs). We present step ‑by‑step instruc‑ tions for implementing Verified Position Correction ( VPC) dead‑reckoning in R using the tilt ‑ compensated compass method, accompanied by the mathematical protocols underlying the code and improvements and extensions of this technique to reduce the trade‑ off between VPC rate and dead‑reckoning accuracy. These protocols are all built into a user‑friendly, fully annotated VPC dead‑reckoning R function; Gundog.Tracks, with multi‑functionality to reconstruct animal movement paths across terrestrial, aquatic, and aerial systems, provided within the Additional file 4 as well as online (GitHub). Results: The Gundog.Tracks function is demonstrated on three contrasting model species (the African lion Panthera leo, the Magellanic penguin Spheniscus magellanicus, and the Imperial cormorant Leucocarbo atriceps) moving on land, in water and in air. We show the effect of uncorrected errors in speed estimations, heading inaccuracies and infrequent VPC rate and demonstrate how these issues can be addressed. Conclusions: The function provided will allow anyone familiar with R to dead‑reckon animal tracks readily and accu‑ rately, as the key complex issues are dealt with by Gundog.Tracks. This will help the community to consider and imple‑ ment a valuable, but often overlooked method of reconstructing high‑resolution animal movement paths across diverse species and systems without requiring a bespoke application. Keywords: Animal behaviour, Animal movement, Global Positioning System, R (programming language), Track integration, Tri‑axial accelerometers, Tri‑axial magnetometers Background Reconstructing animal movement paths is an important *Correspondence: richard.m.g@hotmail.com 1 tool in ecology, providing insights into animal space-use, Swansea Lab for Animal Movement, Department of Biosciences, College of Science, Swansea University, Singleton Park, Swansea SA2 8PP, Wales, behaviour and habitat selection [1–3]. However, accurate UK estimation of paths at fine temporal scales has proved a Full list of author information is available at the end of the article © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Gunner et al. Anim Biotelemetry (2021) 9:23 Page 2 of 37 persistent challenge [4, 5]. Dead-reckoning is a method whilst the utility of transmission telemetry is restricted used to reconstruct animal movement paths, based on to periodic resurfacing events [58]. Moreover, speed can sequentially integrating the vector of travel from a pre- be more easily measured or approximated in water, with determined position using estimates of heading (also previous studies obtaining estimates via acoustic flow termed ‘bearing’ or ‘yaw’) and speed (and displacement noise [e.g., 59], passive sonar [e.g., 60], pitch and change about the vertical axis for 3-D movements), over an in depth [e.g., 11] and speed sensors [cf. 16, 61, 62]. The elapsed time interval [6–9]. In its most advanced form, it efficacy of such techniques diminishes within the aerial can provide positional data with sub-second resolution, environment, principally, due to the marked difference irrespective of the environment [e.g., 10, 11, 12] and it between water and air density [cf. 63] and the current therefore has huge potential for providing data that can speed and volatility of wind [cf. 64, 65]. Indeed, this may elucidate many fundamental behavioural and ecological explain why, in part (to our knowledge), only one study issues related to space-use. to date has dead-reckoned a volant species [66]. More The concept of dead-reckoning (also termed ‘track recently, dynamic body acceleration (DBA, see Wilson integration’) originated to aid nautical navigation [6, 13], et al. [67], for recent review) has been validated as a proxy though its utility to reconstruct uninterrupted fine-scale of speed for terrestrial animals [68, 69] although there are (in time and space) animal movement paths by inte- still very few studies that use the dead-reckoning method grating different sensors in animal-attached tags was in terrestrial animals [e.g., 9, 10, 26, 34, 70, 71]. suggested over three decades ago [14, 15]. Today, this We suggest that another reason that Verified Position typically involves the simultaneous deployment of tri- Correction (VPC) dead-reckoning has been little used axial accelerometers and magnetometers [e.g., 9, 10, 16, relates to the apparent difficulty and poor accessibility of 17–20], utilising the tilt-compensated compass method the analytical processes involved. With this in mind, the [21–24] (see “Glossary” for a definition of dead-reckon - primary aim of this paper is to provide potential users ing-related terminology used throughout). with a clear, concise roadmap for implementing dead- The utility of dead-reckoning depends on the accu - reckoning protocols. Specifically, we revisit the dead- racy of speed and heading estimates (see Table  1) and, reckoning methodology, from calibrating magnetometry due to the nature of vector integration, dead-reckoned data and deriving heading (tilt-compensated compass tracks accumulate errors (commonly termed ‘drift’) over method), to VPC dead-reckoning within both terrestrial time [15, 25, 26]. As a result, periodic ground-truthing and fluid media. We provide simplistic conceptual expla - by a secondary source is important for maintaining the nations and mathematical protocols and describe the accuracy of animal paths with all its underlying scales pitfalls within the procedure that can increase error. We and tortuosity of movement [9, 10, 27]. For this reason, also translate the relevant equations into complementary dead-reckoning data are normally enhanced by combin- R code [cf. 93, available at 94] throughout the text, with ing it with other methods for providing Verified Positions fully annotated scripts deposited in Additional files 2, 3, (VPs). These are primarily; direct observation [e.g., 28], 4, 5 and GitHub [available at 95]. light intensity-based geolocation [e.g., 29], VHF—[e.g., In addition to the above, we outline recent advances 30], acoustic—[e.g., 31] and GPS telemetry [e.g., 26]. to the VPC dead-reckoning technique. For use in ter- Other, less utilised, systems that may also have merit restrial environments, this includes implementing step at sites frequented by the tagged animals, include radio counts as a distance measure, by-passing dynamic body frequency identification (RFID) stations [cf. 32], camera acceleration (DBA) as a proxy of speed, and assessing the traps [cf. 33] and video footage, such as closed-circuit tel- value of ‘reverse dead-reckoning’ (useful when VPs are evision (CCTV) surveillance [e.g., 34]. Although all these concentrated to the latter end of an animal’s trajectory). systems are subject to a number of issues that can make For marine and aerial environments, we demonstrate their positional fixes temporally widely spaced [e.g., 4, 35, the value of integrating tidal-/air current data with dead- 36], inaccurate [e.g., 37, 38, 39] or impossible [e.g., 40, reckoned vectors (hereafter termed ‘current integration’) 41, 42], they can be critical in providing ground-truthed to reduce errors attributed to drift [cf. 46, 92]. Across all positions, even infrequently, with which to reset accumu- three media of travel (land, water and air), we show the lated drift [9, 26]. value of incorporating different speed coefficients accord - Of the above VP options, GPS-corrected dead-reck- ing to behaviour types. In addition, we provide examples oning is the most widely used and there is a marked bias of the various methods by which VP data can be under- towards marine studies [e.g., 10–12, 15–17, 19, 27, 43, sampled (relevant for high-res GPS datasets) prior to 44–56]. This is likely for logistical reasons, with many correcting dead-reckoned tracks and discuss the scales aquatic animals being larger (and thus can carry larger/ at which users should consider VP correction (which more devices) than their terrestrial counterparts [57], depend on the details of species-specific movement and Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 3 of 37 ‑ ‑ ff Table 1 Possible system errors that can affect the utility of animal deadreckoning within the ‘tilt compensated compass’ framework System error Reasons for error Underlying causes References Possible mitigation measures Derived heading Errors in deriving static acceleration (postural) During bouts of high centripetal (turning) accel [65, 72–75] Gyrointegrated data [cf. 76, 77] estimates eration Freefalling behaviour Using Euler angles (angle of rotation about each The orientation of the device with respect to [21, 79–82] Quaternion estimated heading [cf. 83, 84, 85] axis of a given coordinate system) the earth’s frame of reference (cf. Additional file 1: Text S2) can only be defined reasonably at angles less than perpendicular or less than a 180° inversion (dependent on pitch and roll equations used—cf. “VPC deadreckoning procedure in R” section) from their longitudinal and lateral axes of ‘normal’ posture (otherwise, unstable measures arise from the Gimbal lock singularity complex [cf. 78], whilst x, y and z values can become inversed and/or represent different ‘surge’, ‘sway ’ and ‘heave’ planes) Tag placement/dislodgment In line with the above—range for accurate angular [19, 86] Ensure tag orientation is noted during deployment (pitch and roll) measures are restricted in one or and retrieval operations (and subsequently used more dimensions in corrections) Heading will be biased according to the degree of displacement about the zaxis Variations in the strength and declination of Animals that undertake long journeys (regionally/ [87] Ensure at least one magnetic calibration procedure magnetic fields globally) is carried out (see Additional file 1: Text S3 for Environmental and manmade magnetic noise details) and apply magnetic declination offset to (iron distortions) heading values where required Derived speed Deviations of the DBA ~ speed relationship Load bearing [12, 65, 68, 69, 88–91] Iteratively modulate the gradient and/or intercept Moving over a deformable substrate/changeable within the DBA–speed linear regression according incline to environmental circumstance and mode of Changing gait movement [cf. 68] Moving within fluid media Bypass DBA (e.g., use speed/acoustic sensors, step/ Gliding/thermalling behaviour tail/wing beat frequency, vertical speed, etc. [e.g., 47]) Both External forces (e.g., current vectors in air and Decreases the signalto noise ratio of motion [11, 15, 46, 92] Smooth postural (and pre derivative data)/DBA water flow) sensor data estimates [cf. 67] Aects the relationship between an animal’s Incorporate current flow vectors within the dead (longitudinal axis) direction of travel from their reckoning procedure true vector of travel Some animals do not always move in the same direction as their anterior–posterior axis Gunner et al. Anim Biotelemetry (2021) 9:23 Page 4 of 37 User functionality length of data acquisition). We specifically demonstrate Gundog.Tracks is an all-encompassing dead-reckoning the above using our R-functions (Gundog.Tracks being function that can be used to dead-reckon animal paths the primary function for dead-reckoning), providing travelling terrestrially or through fluid media. Table  2 examples of its utility across various scenarios. Lastly, details all the function’s input requirements/options. we highlight the relevance of heading and distance cor- rection factors (derived from the VPC procedure), which Reverse dead‑reckoning can also be used to interrogate the animal–environ- Dead-reckoning backwards is useful when the start ment interaction and biases stemming from animal tag position is unknown, but the finishing coordinates are performance. known. For example, central-place foraging, diving ani- To illustrate our approach, we use three model species mals returning to land from the sea may not acquire a (the African lion Panthera leo, the Magellanic penguin satellite fix for an appreciable period of time following Spheniscus magellanicus and the Imperial cormorant submersion in water which can make determining the Leucocarbo atriceps) that cover almost two orders of size start position difficult. So, when VPs are skewed to the magnitude in body mass and that operate in markedly dif- latter part of the track, it may be beneficial to start the ferent environments and at different scales of movement. iterative dead-reckoning process from that end. This To make this review more broadly applicable to research- involves reversing the order of data to be dead-reckoned ers of varying dead-reckoning and R knowledge, we have and changing heading values by 180 degrees prior to departed from a traditional article format and instead, dead-reckoning. split this body of work into two distinct sections: firstly, we provide an overview of the critical Gundog.Tracks Integrating current vectors function and provide a brief review of the conceptual Wind or ocean currents can change the relationship workflow (“Implementation of Gundog.Tracks” section). between an animal’s (longitudinal axis) bearing and With respect to this, we discuss the relevant strengths speed of travel from their true vector of travel [46, 92]. and limitations of the current dead-reckoning framework This drift can be incorporated within movement paths by and the key considerations involved. Secondly, we detail advancing each iterated dead-reckoned vector according each ‘potential’ stage of the VPC dead-reckoning proce- to the direction and speed of the current at that point in dure with exemplar mathematical equations and R syntax space and time (cf. Fig. 1). (“VPC dead-reckoning procedure in R” section). DBA–speed derivation Implementation of Gundog.Tracks Given the approximate linear relationship between We expand on key concepts in Additional file  1 and DBA [sensu 67] and terrestrial animal speed provide complimentary R scripts (outlined below) in [speed = (DBA·m) + c], DBA estimates can be multi- Additional files 2, 3, 4 and 5. We also supply an exam- plied by a gradient, m (the multiplicative coefficient) ple data set of a Magellanic penguin walking out to sea and summed with an intercept, c (the constant) to in Additional file  6, which can be used to trial each of derive speed [10, 26]. These values are typically substi - the provided R scripts and perform the full dead-reck- tuted with results from DBA–speed linear regression oning process. Mathematical equations are referred to estimates, such as from treadmill tests or using GPS- as ‘Eqs.  1–33’ and R syntax as ‘R ’, where ‘ ’ is the refer- x x derived speed [26, 69, 97, 98]. The m -coefficient should ence number. To simplify concepts, we use base R syntax be selected such that (uncorrected) dead-reckoned (wherever possible) and typically use vectors to dem- tracks accord with the apparent straight-line distance onstrate points made, though ‘df$’ directly before the between VPs. Importantly, the DBA–speed relation- variable name indexes data retained within data frame ship may be a function of terrain-type (e.g., sand vs. columns (assuming data frame is called ‘df’). We note, concrete), animal state (e.g., weight variation) and however, that more efficient code implementations are mode of movement (e.g., running vs. climbing) [cf. 68]. possible (e.g., data.table [96] and lapply()) than pre- For instance, a condor gliding within a thermal would sented here, especially for large data, but here wanted have high speeds, despite having negligible DBA, while to make the code as readable as possible in this manu- an Ibex traversing across different substrate types and script, especially to persons not familiar with complex gradients would impart varying magnitudes of accel- coding. More efficient code will be implemented through eration that may scale non-linearly with a change in updated GitHub versions of the functions. See Additional stride gait. It may be of value, therefore, to iteratively file  1: Text S1 for our model species’ device setup and change the supplied m (and possibly c) values between capture protocol and the glossary for a definition of dead- VPs according to behaviour and environment. The user reckoning related terminology. Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 5 of 37 Table 2 Gundog.Tracks input fields and description of their role Funcon input Descripon Ref TS Timestamp - POSIXct object. No missing data (NA’s) permied - o o h Heading (0 to 360 ) - No missing data (NA’s) permied v - DBA (g or m/s ) or speed (m/s) - No missing data (NA’s) permied Elevaƒon / depth data (m) - No missing data (NA’s) permied 0 elv Pitch ( ) – Only supply if user wants radial distance modulated according to pitch (cf. "VPC dead-reckoning p NULL ). No Current speed (m/s) - Supplied as a single value or vector/column of changeable values. NA’s are replaced with cs NULL the most recent non-NA prior to it (observaƒons carried forward) o o Current heading (0 to 360 ) - Supplied as a single value or vector/column of changeable values. NA’s are replaced ch NULL with the most recent non-NA prior to it (observaƒons carried forward) Mulƒplicaƒve coefficient (gradient) - If speed (m/s) supplied for v, then m must be 1. m 1 Supplied as a single value or vector/column of changeable values Constant (y-intercept) – If speed supplied for v, then c must be 0. c 0 Supplied as a single value or vector/column of changeable values Marked Events – 0 denotes periods of sta�onary behaviour and 1 (or any integer number > 0) denotes periods of ME traversing movement. ME overrides ini�al speed input / DBA-derived speed (calculated within the func�on itself). 1 NA’s and character values are replaced with zero. lo Star�ng longitude coordinate to advance dead-reckon track from – Decimal format, e.g., 26.31989 0 Star�ng la�tude coordinate to advance dead-reckon track from – Decimal format e.g., -06.11995 0 la VP longitude coordinates – Decimal format. Missing reloca�on data expressed as either NA’s or 0’s. First (or last VP.lon NULL if reverse dead-reckoning) element/row allocated as lo within the func�on VP la�tude coordinates – Decimal format. Missing reloca�on data expressed as either NA’s or 0’s. First (or last if VP.lat NULL reverse dead-reckoning) element/row allocated as la within the func�on TRUE = Supplied VPs removed at �mes when ME = 0 (relevant for high-res VP datasets, when loca�on error is high VP.ME FALSE during rest). Note, this does not remove the element/row allocated as lo/la How the func�on under-samples VPs prior to correc�on (subsequent to the VP.ME subset, if set to TRUE) – “divide” = Fix kept every x (thresh) segments of supplied VPs, based on row number. The first and last fixes are always included “ me” (s) = F ix kept every x (thresh) accumulated seconds ( or t he n ext available fix a œer a period o f missing method NULL loca�onal data (≥ thresh). The first and last fixes are always included “distance” = Fix kept every x (thresh) p ropor�onal s egments of t he t otal a ccumulated distance ( m) b etween supplied VPs (using the stepping interval ‘dist.step’). The first and last fixes are always included “all” = Every supplied VP kept (irrespec�ve of thresh value) Threshold - Degree of V P under-sampling prior to d ead-reckon c orrec�on. The f requency of under sampling thresh 1 depends on the method selected The stepping interval used for calcula�ng distance b etween V Ps, both w ithin the VP s ummary d istance m etrics (see Table. 3) and within the ‘method = distance' VP under-sampling protocol prior to VPC. For example, dist.step dist.step 1 = 5 computes distance between every 5th VP (irrespec�ve of the �me difference between them) TRUE = VPC dead-reckoning is bounded by the first and last VP present bound FALSE = VPC dead-reckoning is unbounded by the last available VP. The last dead-reckoned track segment inherits TRUE the previous correc�on factors TRUE = ‘normal’ dead-reckoning procedure Outgoing FALSE = R everse dead-reckoning. N ote la a nd lo p osi�ons should n ow b e the finishing l ongitude a nd la�tude TRUE coordinates, respec�vely FALSE = No summary plots TRUE = R graphics window ini�alized: VPC = 4 summary plots / no VPC = 1 summary plot (cf. top leœ ) Top le) Uncorrected dead-reckoned track (blue) and VP track (red). If currents are supplied, the blue track has currents integrated and an addi�onal green track with no current integra�on is plo¤e d Plot Top right) VPC dead-reckoned track (blue) in rela�on to VP track (red) TRUE Bo�om le) Net error (m) between VPs and dead-reckoned posi�ons (un-corrected = red and corrected = black). If currents are supplied, then uncorrected w ith no c urrent i ntegra�on = green a nd uncorrected w ith current integra�on = red. Bo�om right) VPC corrected dead-reckoned track Ref. refers to the default value when no input is stated. Red shading represents required user inputs and green and orange shading reflect optional inputs (the latter change when using VPC). Note that if speed estimates (v) are directly inputted into the function then m and c (and possibly ME) defaults should not be changed. If either one of the VP.lon, VP.lat or method inputs is specified as NULL, then no VPC will occur Gunner et al. Anim Biotelemetry (2021) 9:23 Page 6 of 37 DBA is a weak proxy of speed for many marine animals because overall body tissue density changes with depth when air is associated so that speed may be invariant of the movement kinematics [cf. 103, 104]. DBA is also a weak proxy for flying animals that glide at constant veloc - ity, use thermals or bank [cf. 105]. One of the most com- mon methods for determining animal speed in water is via devices that estimate flow or resistance rate [16, 19, 64, 106]. These often have appreciable limitations, with currents, biofouling, blockage and turbulence affecting performance [64], and many of these issues are appli- cable to volant species, so that bird speed measures are typically restricted to GPS-derived estimates of ground Fig. 1 Schematic representation of a current flow vector (orange) speed [cf. 107]. In the absence of a reliable motion sen- (due to its speed and direction) being integrated to a given travel sor-derived speed proxy, previous reported approximated vector (blue). The x, y reflect the initial location of a dead‑reckoned speed estimates according to movement modes and/or track, x2 and y2 are the resultant location following the integration of topological whereabouts can be used [cf. 30]. For exam- a travel vector (prior to current integration) and xxx and yyy advance ple, for various diving animals such as penguins, a sim- these x 2 and y2 values a step further in the direction of the current flow vector. The dashed lines indicate the magnitude of the x and y ple depth threshold may prove effective to differentiate dimensions of travel (both pre‑ and post ‑ current flow integration) between various previous reported modal ‘surface-rest- and the green line reflects the actual travel vector ing’ and ‘underwater-commuting’ speeds [61, 108–110]. For volant species, whilst wingbeat frequency or ampli- tude does not scale reliably with air speed [cf. 111], the may also opt to supply a ‘marked events’ (ME) vector (a interplay between both can decipher various flight modes marked event is a term we use that refers to a number of (e.g., ‘cruising speed’ vs taking off/landing) [65, 112]. sequential (in time) data points within a dataset coded Furthermore, tail beat frequency has been shown to be by integer values) to ensure dead-reckoned tracks are a good predictor of swimming speed for various fish spe - not advanced with non-animal movement behaviours. cies [113–115]. For diving animals, a proxy for horizontal Within Gundog.Tracks, ME values of one or greater speed can be obtained based on animal pitch and rate of reflect progressive movement, and zero values code for change of depth [11, 116]. Specifically, the rate change of stationary behaviour—dead-reckoned tracks are not depth is divided by the tangent of the body pitch. advanced when ME = 0 (irrespective of the allocated In any case, when high-resolution VP data are available speed). For example, in its simplest form, ME could (e.g., 0.01–10  Hz GPS), for instance, during short-term be filled with binary 0  s and 1  s as governed by a DBA trial deployments, speed estimates can be compared threshold (labelling the ME vector 0 in sleep and resting alongside those derived between VPs and approximated behaviour). according to behaviour type (elucidated from, for exam- ple, accelerometer—[e.g., 117, 118], magnetometry— Pre‑determining speed [e.g., 105, 119], depth- [e.g., 120] or altitude—[e.g., 65] For terrestrial species (specifically bipeds and quadru - data), and uncorrected dead-reckoned tracks can be peds), the interplay between peak heave acceleration compared alongside VPs to determine where biases amplitude and periodicity may be a useful indicator for may occur visually. Furthermore, the correction factors the movement gait adopted [99], which may help decide obtained from the VPC process are viable comparators the m-coefficient in the DBA–speed relationship [69]. for detecting consistent under- or over-estimations of There may be times, however, when DBA is an unreliable speed and/or heading offsets (e.g., due to tag placement). proxy of terrestrial speed [cf. 68]. At this time, given that Essentially, when empirical speed evidence is unavail- the stride cycle can be easily detected by cyclic peaks in able, the user can ad hoc iteratively adjust allocated speed a given acceleration channel [e.g., 100, 101, 102], peak values or the underlying DBA–speed coefficients until periodicity (and amplitude) may be used as a proxy of uncorrected dead-reckoned track segments scale propor- distance moved by providing a distance per step estimate tionately to their aligned ground-truthed positions (pre- (assuming constant distance travelled between step gaits VPC). Within Gundog.Tracks, the user can modulate m, if only concerning step periodicity—cf. “VPC dead-reck- c and ME values to switch between predetermined speed oning procedure in R” section and Additional file  1: Text (m = 1, c = 0), DBA-derived speed (m > 0, c ≥ 0) and sta- S4). tionary behaviour (ME = 0). Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 7 of 37 VPC procedure with input modulated according to the animal in ques- Ground-truthing dead-reckoned tracks typically involves tion and data available (see Fig. 2). the linear drift correction method [cf. 26, 46], outlined The function outputs a data frame containing vari - in Constandache et al. [121] and Symington and Trigoni ous descriptive columns which, depending on the input, [122]. In essence, a shift vector aligns the starting dead- includes (but is not limited to): reckoned path segment with the VP at time point one, after which the difference between the VP and dead-reck - • The correction factors used oned path segment at time point two is calculated to pro- • Heading and radial distance estimates (both pre- and vide a correction vector that is applied linearly between post-current integration and/or VPC) time point one and time point two. Our method follows • Distance moved and speed estimates (both in 2-D the protocols outlined by Walker et  al. [9], whereby the and 3-D when elevation/depth data supplied) underlying correction coefficients (hereafter termed • Net error between dead-reckoned positions and VPs ‘factors’) for both heading and (radial) distance are cal- (both pre- and post-correction) culated—adjusting the length and heading at each dead- • Various VP summaries including notation of when reckoned path segment until the end points align to each VPs are present and which fixes were used in the cor - VP along the path. This process requires the trigonomet - rection process. ric ‘as the crow flies’ Haversine formulae [123–125] which allows one to translate a distance across the curvature of When specified, 2-D summary plots demonstrating the the Earth’s surface (detailed within “VPC dead-reckoning relationship between dead-reckoned positions and VPs procedure in R” section). The advantage of this method (both pre- and post-current integration and/or VPC) are is that, whilst correction factors are constant between provided (e.g., Fig.  3). Table  3 details all the function’s VPs, it does not assume that the dead-reckoned path available outputs (modulated according to input). Gun- deviates linearly over time from the true path because dog.Tracks uses the na.locf() function from the ‘zoo’ pack- (radial) distance is multiplied by the distance correction age [126] and the slice() function from the ‘dplyr’ package factor. This ensures that parts of track where the animal [127] (both are checked as dependencies and installed is determined to be stationary (e.g., ME = 0) are left unal- when required within this function). Output 2-D dis- tered. The function’s method of VPC, automatically han - tance/speed estimates are calculated with the Haversine dles NaN and Infinite (Inf ) values which can arise during formula. When depth/elevation data are supplied (and the derivation of the distance correction factors (when changes between sets of coordinates) 3-D distance/speed no dead-reckoned movement occurs between successive estimates are calculated with a variant of the Euclidean VPs—detailed within “VPC dead-reckoning procedure Formula—converting x, y, z from polar to Cartesian coor- in R” section). It is worth noting that even animals that dinates, and incorporating the Earth’s oblate spheroid [cf. travel in 3-D can be subject to the 2-D dead-reckoning World Geodetic System (WGS84)], via conversion from formulae and Haversine computation of distance correc- Geodetic- to Geocentric-latitude [cf. 128]. tion factors because we typically assume that both dead- The interplay between numerical precision in R, cor - reckoned- and VP positions are aligned in vertical space rection rate and net error can make more than one (assuming reliable pressure—[60]/depth [13] data) and round of adjustment necessary for dead-reckoning fixes attempt to control for the horizontal component of speed to accord exactly with ground-truthed locations (cf. [e.g., “VPC dead-reckoning procedure in R” section— Fig. 4a), particularly given that slight discrepancies accu- Eqs. (25, 27)] pre-correction. Although not covered here, mulate over time. Each iteration of the correction pro- we acknowledge that various state–space modelling tech- cess produces a tighter adherence between estimated and niques have also been developed to georeference dead- ground-truthed positions [cf. 9]. Typically, this does not reckoned tracks [e.g., 11, 47]. involve more than two rounds of VPC to achieve a maxi- mum net error of 0.01 m (the threshold used within Gun- Default inputs for calculations and outputs dog.Tracks) across a ca. (1 Hz) 2-week-long track. For an Gundog.Tracks default input takes the form: indication of processing times see Additional file  1: Text S6, Fig. S4; for example dead-reckoning a lion at 1  Hz Gundog.Tracks(TS, h, v, elv = 0, p = NULL, for 7 (continuous) days (with plot = TRUE, dist.step = 5, cs = NULL, ch = NULL, m = 1, c = 0, ME = 1, VP.ME = TRUE, method = “time” and thresh = 3600) lo = 0, la = 0, VP.lon = NULL, VP.lat = NULL, took 25 s to compute (on a MSI GP72 7RD Leopard lap- VP.ME = FALSE, method = NULL, thresh = 1, top with intel core i7 processor). Logically, the net error dist.step = 1, bound = TRUE, Outgo- between VPs and (corrected) dead-reckoned positions is ing = TRUE, plot = FALSE), positively correlated to the time between corrections (cf. Gunner et al. Anim Biotelemetry (2021) 9:23 Page 8 of 37 Fig. 2 Schematic of the conceptual workflow involved when dead‑reckoning using Gundog.Tracks—elaborated within “ VPC dead‑reckoning procedure in R” section. Note Gundog.Peaks (Additional file 3) is a peak finder function that locates peaks based on local signal maxima and Gundog.Compass (Additional file 2) is a function to correct iron distortions from tri‑axial magnetometry data and subsequently compute tilt‑ compensated heading. Both functions are elaborated within “VPC dead‑reckoning procedure in R ” section and in Additional file 1. The direction of workflow and key questions asked follow from green—(pre ‑processing and data alignment) to purple—(computing heading) sections, before splitting into blue—(air/water) and brown—(land) sections (computing speed) and culminating at the red section (final pre ‑ dead reckoning checks/data formats and post‑ dead‑reckoning checks/plots) in conjunction with the process of using Gundog.Tracks in R (yellow) Fig.  4b) [cf. 46], although the rate of net error ‘drop-off ’ according to the permissiveness of the environment, such is dependent on the accuracy of the initial (uncorrected) as high-density shrub or submersion in water [e.g., 38, dead-reckoned track (cf. Fig.  5), itself, modulated by the 129, 130]. GPS technology is arguably the most popular extent of system errors (Table  1) and initial user-defined and widely used method for determining estimates of track scaling. free-ranging animal movement [cf. 131, 132, 133]. This Within this process, people assume VPs to be perfect, is because inspection of data is less complex and time- however, across all VP determining methods, the rate consuming than some of the alternatives, whilst improve- and accuracy of data acquisition is highly moderated ments in design and battery longevity have enabled Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 9 of 37 Fig. 3 Dead‑reckoned (DR) movement path of lion as provided by Gundog.Tracks summary plots (within the initialised R graphics window). This is an approximate 2‑ week trajectory over an approximated total travel (DR) distance of > 142 km. (Pre‑filtered) GPS (red) was sampled at 1 Hz and derived heading and speed measurements were sub‑sampled to 1 Hz (initial acceleration/magnetometry data were recorded at 40 Hz). The VPC dead‑reckoned track (blue) was constructed using DBA–GPS‑ derived speed regression estimates and corrected approx. every 6 h. Note, for dead‑reckoning within fluid media, an additional green dead‑reckoned track with current integration and its associated distance estimates are also plotted (pre‑ correction) when wind/ocean currents are supplied (cf. Fig. 6). Accumulated 3‑D DR distance is shown when elevation/depth data are supplied GPS units to be attached to a plethora of animals (up to appreciably more depending upon the propagation of sig- almost four orders of magnitude in size and mass [cf. 11, nal quality and/or receiver reception capability [38, 138, 134]) and record at high frequencies (e.g., ≥ 1  Hz [131, 139]. As such, VP error becomes more relevant at smaller 135]). Consequently, GPS units are unparalleled for pro- scales of assessed movement and this is the reason why viding such detailed quantification of space-use out - VP distance-moved estimates can go from being typically side of the VPC dead-reckoning framework, and are the underestimated at low frequencies (due to linear interpo- most utilised VPC method within (including the case lation of tortuous movements) [26, 140, 141] to overesti- study datasets within this study). However, locational mated at high frequencies [97, 136] and result in highly accuracy (excepting precision error radius [cf. 136] and variable correction factors within the VPC dead-reckon- variable latency [cf. 137]) can vary by a few metres or be ing process [cf. 10]. Indeed, judicious selection of VPC Gunner et al. Anim Biotelemetry (2021) 9:23 Page 10 of 37 Table 3 Gundog.Tracks data frame output names and their parameters Funcon output Descripon Ref Row.number Row number Supplied mestamp - POSIXct object Timestamp DR.seconds Accumulated me (s) based on the supplied mestamp o o Heading Supplied heading (0 to 360 ) Marked.events Supplied Marked events (or replicated default) DBA.or.speed Supplied DBA (g or m/s ) or speed (m/s) Pitch Supplied pitch ( ) The calculated q-coefficient (prior to VPC) (see “VPC dead-reckoning procedure in R” secon—Eq. (26)) # Radial.distance Elevaon Supplied elevaon / depth (m) Elevaon.diff Rate change of supplied elevaon/depth (m/s) - (elevaon difference / me difference between rows) Supplied current speed (m/s) Current.strength o o Current.heading Supplied current heading (0 to 360 ) o o Heading.current.integrated # Updated heading (0 to 360 ) following addion of current vectors (prior to VPC) Radial.distance.current.integrated Updated q-coefficient following addion of current vectors (prior to VPC) # DR.longitude Dead-reckoned longitude coordinates – Decimal format (prior to VPC) DR.la tude Dead-reckoned latude coordinates – Decimal format (prior to VPC) Corrected dead-reckoned longitude coordinates – Decimal format (post VPC) DR.longitude.corr DR.la tude.corr Corrected dead-reckoned latude coordinates – Decimal format (post VPC) Dist.corr.factor Distance correcon factor (observaons carried forward) ↑ # o o Head.corr.factor ↑ # Heading correcon factor (0 to 360 ) (observaons carried forward) o o Heading.corr # Corrected heading (0 to 360 ) (post VPC) Corrected q-coefficient (post VPC) # Radial.distance.corr Distance (m) between uncorrected dead-reckoned posions and VPs (observaons carried forward), Distance.error.before.correc on ↑ subsequent to sub-sampling according to ME, if VP.ME = TRUE Distance (m) between corrected dead-reckoned posions and VPs (observaons carried forward), Distance.error.a er.correc on ↑ subsequent to sub-sampling according to ME, if VP.ME = TRUE DR.distance.2D Two-dimensional distance moved (m) between dead-reckoned fixes * DR.distance.3D Three-dimensional distance moved (m) between dead-reckoned fixes * DR.cumula ve.distance.2D Accumulated two-dimensional distance moved (m) between dead-reckoned fixes * DR.cumula ve.distance.3D Accumulated three-dimensional distance moved (m) between dead-reckoned fixes * Two-dimensional (straight-line) distance moved (m) from starng posion DR.distance.from.start.2D * Three-dimensional (straight-line) distance moved (m) from the starng posion * DR.distance.from.start.3D DR.speed.2D Horizontal speed (m/s) (DR.distance.2D / me difference between rows) * DR.speed.3D Total speed (m/s) (DR.distance.3D / me difference between rows) * Accumulated me (s) between supplied VPs (observaons carried forward) VP.seconds VP.longitude Supplied VP longitude values (observaons carried forward), sub-sampled according to ME, if VP.ME = TRUE ↑ VP.la tude Supplied VP latude values (observaons carried forward), sub-sampled according to ME, if VP.ME = TRUE ↑ Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 11 of 37 Table 3 (continued) Denotes when a fix was present (1) or absent (0), subsequent to sub-sampling according to ME, if VP.ME = VP.fix.present TRUE Denotes which VPs were used to correct (1) and which VPs were ignored (0) VP.used.to.correct Number.of.VPCs Increments by 1 each ƒme a VP was used to correct (observaƒons carried forward) ↑ Replicates the thresh value set (or default) or warns the user that addiƒonal VP under-sampling was VP.thresh required if ‘Inf’ values produced Two-dimensional distance moved (m) between VPs, subsequent to sub-sampling according to ME, if VP.ME VP.distance.2D = TRUE and using the stepping interval ‘dist.step’ Accumulated two-dimensional distance moved (m) between VPs, subsequent to sub-sampling according to VP.cumulave.distance.2D ME, if VP.ME = TRUE and using the stepping interval ‘dist.step’ The shading of Ref refers to when the outputs occur; red shading = always, purple shading = when pitch data are supplied, green shading = when elevation/depth data are supplied, blue = when current data are supplied and orange = when the user opts to undertake VPC. The symbol * demonstrates that the metrics will be derived from VPC tracks when correction is initialised. Note that, subsequent to reverse dead-reckoning, the data frame is reverted (back to original time order), though as a result, observations appear to be carried backwards in some instances, indicated by ↑. Due to the nature of reverse dead-reckoning (cf. Additional file 4), some input fields are shifted forward one row following the initial inversion of data. As such, fields indicated by #, are one row further forward in time (this is important when relating Head.corr.factor and Heading.corr output to the equivalent (uncorrected) Heading output. However, when currents are integrated, the Head. corr.factor and Heading.corr outputs refer to Heading.current.integrated and these are synchronised row-wise. All heading related data are rotated back 180 degrees following reverse dead-reckoning Fig. 4 Net error between (GPS‑ corrected) dead‑reckoned and GPS positions for a track from 5 African lions. a Maximum net error (m) between ground‑truthed GPS and time ‑matched dead‑reckoned positions after one iteration of correction, both as a function of GPS correction rate [one correction per 1—(red), 12—(green) and 24—(blue) hours] and underlying m‑ coefficient used to determine the DBA‑ derived speed. Data from 5 lions (individual denoted by symbol shape) over a period of 12 days. Note that the difference in error varies according to individual, initial speed estimate and the rate of correction. b Net error between dead‑reckoned positions and all available GPS fixes, according to correction rate (data from the same 5 lions), subsequent to the iterative procedure of GPS correction (maximum distance between GPS fix used in correction procedure and according dead‑reckoned position < 0.01 m). Boxes denote the median and 25–75% interquartile range with a blue ‘loess’ smooth line. Whiskers extend to 1.5 * interquartile range in both directions rate is critical in maximising dead-reckoned track accu- It was suggested by Bidder et  al. [10], that the next racy when relocation data are taken at fine spatial- and stage in this work is to derive a standardised set of rules temporal resolutions [26] (cf. Table  2—‘VP.ME’, ‘method’, to maximise the value of both GPS (though this applies to ‘thresh’ and ‘dist.step’ inputs to aid in modulating VPC any VP method) and dead-reckoned data in line with the rate). Likewise, the initial screening for location anoma- questions being asked. We argue that consistent trends lies, across all VP methods and sampling intervals, is in the magnitude and/or bias of correction factors can be important so as to prevent incorrect distortion of tracks. used as a diagnostic tool for elucidating: (i) VP inaccu- Put simply, the higher the quality of VP data input, the racy (e.g., possibly manifested by extremely high distance greater the robustness of the VPC dead-reckoning and heading correction factors), (ii) required alterations output. to the DBA–speed relationship [e.g., due to traversing Gunner et al. Anim Biotelemetry (2021) 9:23 Page 12 of 37 across different substrates (e.g., Fig.  5)] and (iii) drift due incorporated to some degree. Whilst not illustrated here, to current vectors [cf. 16, 46] (e.g., Fig. 6). advancing tracks without a VeDBA threshold would incur greater error still. Lastly, in this section, the dis- The case‑studies tance correction factor was consistently high (Fig. 5, inset An important question to address is how often to do VP b ) as the lion travelled along the Botswana fence bound- correction. This is obviously dependent upon the scales of ary, perhaps as a result of the animal walking on the com- movement elicited and the medium in/on which the ani- pact dirt road at this location (Fig. 5, inset a ), altering the mal in question navigates. Put simply, one should VP cor- VeDBA–speed relationship. Such patterns in correction rect as little as possible, but as much as is necessary and factors (whether consistent or highly variable) can high- we elaborate on this using our model species operating in light issues with the underlying track scaling. different media. Within Fig.  5, the 1 Hz GPS track (blue) Where animals move in water or air, obtaining accu- is plotted alongside two different dead-reckoned tracks; rate estimates of speed is more difficult without the use [(a) = uncorrected and (b) = corrected approx. every of speed sensors. Naturally, the resolution and accuracy 30  min (method = “time”)] from 12  days of data acquisi- of initial dead-reckoning track scaling (pre-VPC) reduces tion of one lion. There were two variations in the method when speed has to be approximated using constant values of scaling the dead-reckoned tracks; a track based on a according to behaviour type (a strategy used here). There Vectorial Dynamic Body Acceleration (VeDBA) threshold is a balance between initial dead-reckoning accuracy (red), and a track advanced based on periods of identified and required VPC. The lower the initial track accuracy, movement (purple). The m-coefficient and c-constant the more frequent it should be corrected, and additional values were determined from the VeDBA–GPS speed drift caused by external-force vectors compounds this relationship (Fig.  5, inset a ) and the Movement Verified issue. Within Fig.  6, we illustrate the value that current Filtering (MVF) protocol outlined by Gunner et  al. [97] correction, dependent on current information, brings to was used to depict movement and anomalous GPS fixes the VPC procedure if the derived track is to be superim- (green) and to compute reasonable GPS-derived speed posed on the environment. Here, one Magellanic pen- estimates. This case study demonstrates three impor - guin was dead-reckoned with and without tidal vector tant points. Firstly, on its own, dead-reckoning is subject integration (instantaneous tidal currents were deduced to substantial drift and so VPC is essential for resetting from a 3-D numerical model validated in the region this error. The more frequent a user corrects, the more [144], at hourly, 1  km grid nodes). Commuting speed accurate the dead-reckon track becomes (relative to VPs), was allocated 2.1 m/s [cf. 61, 145] and changed according though VP error can also be substantial, especially dur- to “VPC dead-reckoning procedure in R” section—R . ing rest behaviour (see Gunner et al. [97] for demonstra- Surface period ‘rest’ speed were allocated 0.416  m/s [cf. tion of this). For collared animals, heading measurements 108]. VP accuracy improved considerably both pre- and can become inaccurate at times of erratic collar roll (cf. post- VPC when currents were integrated which points Table 1) and conjointly, GPS performance is also reduced to the value of acquiring current data if possible, particu- when antenna position becomes compromised [e.g., 142]. larly if VPs are sparse. Notably the combination of dead- Secondly, and in conjunction to the above, irrespective reckoning and VP estimation of both movements relative of VPC rate, the initial allocation of speed is important. to the ground and fluid, may detail specific orientation Here, only dead-reckoning identified movement periods strategies used and thus can have value for assessing the resulted in greater accuracy than just advancing tracks ability of drift compensation in aquatic or volant animals based on a VeDBA threshold. This is because even sta - [46, 92] tionary behaviours can impart appreciable DBA [e.g., For all our case study animals, GPS units were set to 143] (beyond the threshold), and thus wrongly advance record at 1  Hz. With this temporal resolution (which tracks. The false patterns of tortuosity created from is not always possible anyway due to the high-power this, whilst scaled and possibly rotated with VPC (cf. requirements of the GPS), the value of dead-reckoning “VPC dead-reckoning procedure in R” section), remain would seem questionable. However, dead-reckoning can: (See figure on next page.) Fig. 5 Dead‑reckoned lion track in relation to GPS positions [a uncorrected and b corrected—approx. every 30 min (black circles)]. The start of the track (lo and la) is denoted with a black x. Three corresponding sections of each track are denoted with the same number and the finishing positions denoted with a circle (coloured according to its reference track). Note that the horizontal straight‑line sections (cf. yellow arrow) result from the lion following the Botswana boundary fence (which this individual eventually crossed). Mean net error between (corrected) dead‑reckoned positions and all available GPS fixes was higher for tracks resolved using a VeDBA threshold (0.11 g), than for tracks advanced only at times of depicted movement using the MVF method Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 13 of 37 Gunner et al. Anim Biotelemetry (2021) 9:23 Page 14 of 37 Fig. 6 One Magellanic penguin’s dead‑reckoned foraging trip at sea, lasting approximately 9 h (yellow arrow denotes the trajectory direction over time. Black track = GPS. Fifteen corrections (black circles) were made (method = “divide”). For comparison, the grey dotted track is the GPS‑ corrected dead‑reckoned track with current integration approx. every 1 min (where possible ‑method = “time”) (a). Note the difference of net error between dead‑reckoned positions and all available GPS fixes across the various tracks [insert = grey track] (b). Both uncorrected and corrected dead‑reckoned tracks had less error after current integration (black arrows vector every 5 min) and this was reflected in the direction and magnitude of heading correction factors required per unit time (c). Heading correction factors obtained from the track corrected approx. every 1 min; the colour of the scale bar indicates the extent of the heading correction factor required) Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 15 of 37 (i) work when GPS cannot—such as when an animal is error (relative to VPs). As such, and akin with VP under- underwater [e.g., 18] or in thick forest [cf. 2] and it can sampling, choice of under-sampling data to be dead-reck- (ii) by-pass the issues arising from GPS inaccuracies such oned may have implications to the resultant accuracy, as ‘jitter’ [cf. 97], allowing for more accurate and finer and this will be moderated according to the scales (and scale delineations of movement. This is illustrated in media) of movement elicited by the animal in question. Fig. 7, in which 12 outgoing (green) and incoming (blue) Beyond this, Fig.  7 also demonstrates the importance dead-reckoned trajectories from Magellanic penguins of initial track advancement, with three variants used, walking to and from their nest are plotted. Incoming including step counts instead of DBA. tracks were reverse-dead-reckoned (Outgoing = FALSE, Finally, obtaining accurate estimates of altitude or bound = FALSE), because the GPS did not always reg- depth allow users to plot and investigate scales of contin- ister fixes for minutes after birds left the water and uous movement in three dimensions and at times when because nest coordinates were known (Fig. 7, inset). This VP success rate fails completely (such as underwater). We explains why the blue tracks extend into the sea rather demonstrate this using the Imperial cormorant in Fig.  8. than encroach further inland when speed was overesti- After visual inspection of data, uncorrected tracks were mated. What is evident is that even ‘accurate’ GPS paths scaled according to the following speeds: periods of flying are coarsely resolved due to precision errors. Indeed, allocated 12 m/s, surface ‘rest’ periods allocated 0.1 m/s, even with little or no GPS error, this can greatly com- bottom phase of dives allocated 0.4 m/s and descent and promise movement estimates [cf. 136]. Conversely, the ascent speeds modulated according to “VPC dead-reck- precision of the dead-reckoned tracks is only limited by oning procedure in R” section—Eq.  (25). Note that ele- the amount of initial motion sensor data under-sampling vation was not resolved during flying periods (although (usually required in some capacity to make datasets more flying periods were dead-reckoned). Regardless of the manageable and less computationally expensive). Such current limitations, the VPC dead-reckoning proce- fine-scale estimates can therefore (with suitable VPC) dure represents a substantial advance for resolving, and allow users to define movement in space with unprec - thereby allowing investigation of, continuous, fine-scale, edented resolution. The benefit of this is that such reso - free-ranging 2- or 3-D space-use with all its underlying lution can resolve important metrics of movement, such scales of tortuosity and distances moved (e.g., Figs. 7 and as step duration [cf. 146] and the number and extent of 8). turns made [cf. 147]; useful parameters for investigat- ing navigation and foraging strategies according to envi- VPC dead‑reckoning procedure in R ronmental circumstance—though, such parameters are Preparing the three axes of rotation for derivation also useful without superimposing on the environment. of heading Moreover, even dead-reckoned tracks that are sparsely The tilt-compensated compass method is a well-known corrected or never corrected can detail important move- practice for deriving heading [e.g., 21, 22, 81]. Correct ment-specific behaviours [12], for example, circling coordinate system axis alignment and suitable calibra- behaviour [148]. tion of tri-axial magnetometry data [cf. 149] are crucial Ultimately, the higher the frequency at which dead- pre-processors, without which, heading estimates would reckoning is undertaken, the better the resolution and likely incorporate substantial error [cf. 21, 149]. The tilt- detail of reconstructed tracks. However, accuracy only compensated compass method described below (fol- improves up to a point because extrapolated travel vec- lowing the framework outlined by Pedley [21]), requires tors (heading and speed estimates) nearly always com- the aerospace (x-North, y-East, z-Down) (right-handed) prise some degree of error (no matter how small) and so, coordinate system, or ‘NED’ (cf. Additional file  1: Text with very high frequencies (> 1  Hz), more error is accu- S2, Fig. S1). We provide examples of axis alignment, out- mulated per unit time [cf. 16, 44]. In particular, when the line the importance of transforming between coordinate temporal resolution of dead-reckoning results in a spa- frames (relative to the Earths fixed frame) and recom - tial resolution dominated more by sensor noise than by mend a universal configuration calibration procedure to ‘actual’ movement of the animal in question, dead-reck- aid correct axis alignment within Additional file  1: Text oning accuracy will begin to decrease (at least pre-VPC). S2. The extent of this will depend on the size, speed and life - Multiple mathematically sophisticated algorithms style of the animal in question. For example, the benefits have been developed to correct distortions from each of dead-reckoning a lion at 40  Hz rather than 1  Hz are magnetometer channel’s output [e.g., 23, 149, 150, questionable (how often does a lion turn substantially 151, 152]. We provide an annotated R script—Gundog. within a second?), particularly given the additional com- Compass (Additional file 2) that corrects both soft and putation time (cf. Additional file  1: Text S6) and possible hard iron distortions from tri-axial magnetometry data Gunner et al. Anim Biotelemetry (2021) 9:23 Page 16 of 37 Fig. 7 Twelve outgoing (green) and incoming (blue) dead‑reckoned trajectories from Magellanic penguins walking to and from their nest. Three variants of track advancement were used: a A VeDBA threshold (0.1 g) and constant m‑ coefficient (1.4) (b), depicted movement periods using the LoCoD method to identify steps (cf. Wilson et al. 2018) and constant m‑ coefficient (1.4) and c depicted individual steps within depicted movement periods, from which a constant distance estimate (0.16 m) was multiplied by step frequency (x̄ no. steps/s) (full details within Additional file 1: Text S4) (c). Note that the accuracy with respect to the radial distance can be evaluated by examining the track stops in relation to the shore‑line. DBA‑ derived speed estimates were typically overestimated for incoming tracks, due to the birds being heavier (and thus impart greater DBA per stide cycle) after foraing. Tracks (from c) were GPS‑ corrected (d) (method = “ distance” , dist.step = 5, VP .ME = TRUE, thresh = between 8 and 15 (depending on track length) approx. every 50 m). A portion of the GPS‑ corrected dead‑reckoned tracks (bottom panel) are magnified (2 iterations) to show the difference in resolution of movement tortuosity, between GPS and dead‑reckoned tracks Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 17 of 37 Fig. 8 GPS‑ corrected dead‑reckoned tracks of Imperial cormorants foraging at sea: a 15 birds (blue = male, red = female). b Shows one of these tracks illustrated in 3‑D. Note gaps between dives are either associated with current drift, while the bird is resting at the sea surface, or periods of flight. c and d show the descent, bottom phase and ascent of a given dive in both 2‑D (c) and 3‑D, respectively and subsequently computes tilt-compensated heading Tilt‑compensated heading derivation (0° to 360°). Within this function, there are two main Device orientation is expressed in terms of a sequence of methods of correction to choose from, based on the Euler angle [roll (Φ), pitch (θ), yaw (Ψ)] rotations about mathematical protocols outlined by Vitali [153]—least- the x-, y- and z-axes, respectively, relative to the (iner- square error approximation (constructing an ellipsoid tial) Earths fixed frame of reference (e.g., Earth-Centre, rotation matrix) and Winer [154]—scale biases with Earth-Fixed (ECEF) system) [155]. Being a vector field simple orthogonal rescaling (avoiding matrices alto- sensor with two degrees of rotational freedom, acceler- gether). We expand on this user-defined function- ometers are insensitive to rotations about the gravity vec- ality, as well as outlining the causes of soft and hard tor and thus discerning heading requires the arctangent iron distortions and the initial calibration procedure of the ratio between the x- and y-orthogonal magnetom- required to correct such distortions within Additional eter measurements [156]. For the correct computation of file 1: Text S3. heading, these two channels need to be aligned parallel to Gunner et al. Anim Biotelemetry (2021) 9:23 Page 18 of 37 the earth’s surface. This is achieved by correcting any ori - animals, high frequency of body posture changes could entation (de-rotation) according to pitch and roll angles cause discrepancy between static acceleration data and (postural offsets) which can be deduced from accelera - magnetism data, which could consequently affect head - tion. These angles are typically approximated by deriving ing estimation [162]. Although this effect would not gravity-based (static) acceleration [see 72, 157] from each change general shapes of movement paths, we suggest channel by employing one of four approaches using: (i) that prior to the normalisation process (and magnetic a running mean [e.g., 72, 86], (ii) a Fast-Fourier transfor- calibration procedure), it may be of value to initially mation [e.g., 158], (iii) a high-pass filter [e.g., 159] or (iv) smooth out (see Eq. 1, R ) small deviations within mag- 1:4 a Kalman-filter [e.g., 160]. Here, we use a computation- netometry data, both to avoid this type of error and to ally simple running mean over 2 s [72] (Eq. 1): reduce the magnitude of anomalous spikes in magnetic inference. We used a smoothing window of 10 events for i+ the 40 Hz datasets used in this study. G = A x,y,z x,y,z (1)     j=i− NG G x x     NG � G = · y y G · G + G · G + G · G where w is an integer specifying the window size and x x y y z z NG G z z G and A represents the smoothed and raw compo- x,y,z x,y,z (2) nents of acceleration, respectively. In the absence of lin-     NM M x x ear (dynamic) acceleration [see 157, 161], values of G x,y,z 1     NM � M = · y y reflect the device orientation with respect to the earth’s M · M + M · M + M · M x x y y z z NM M z z reference frame (though see Table. 1), reading approx. (3) + 1  g when orientated directly towards the gravity vec- tor (down), − 1  g against the gravity vector (up) and 0  g (R5) NGx = Gx / sqrt(Gx^2 + Gy^2 + at perpendicular to it (horizontal). In R, the ‘zoo’ package Gz^2) [126] provides useful wrappers to apply arithmetic opera- (R6) NGy = Gy / sqrt(Gx^2 + Gy^2 + tions in a rolling fashion (R ). 1:4 Gz^2) (R7) NGz = Gz / sqrt(Gx^2 + Gy^2 + (R1) install.packages("zoo") ; Gz^2) library(zoo) (R8) NMx = Mx / sqrt(Mx^2 + My^2 + (R2) Gx = rollapply(Ax, width=w, Mz^2) FUN=mean, align="center", (R9) NMy = My / sqrt(Mx^2 + My^2 + fill="extend") Mz^2) (R3) Gy = rollapply(Ay, width=w, (R10) NMz = Mz / sqrt(Mx^2 + My^2 + FUN=mean, align="center", Mz^2) fill="extend") (R4) Gz = rollapply(Az, width=w, Depending on deployment position, the device-carried FUN=mean, align="center", NED coordinate frame (the x-, y-, z-axes) may not corre- fill="extend") spond with the animal’s body-carried NED frame. When this occurs, prior to deriving animal orientation, the nor- Here, w should be replaced with the window width malised gravity and magnetic vectors are required to be of choice (e.g., for 20  Hz data and a smoothing of 2  s corrected so that their measurements are expressed rela- required, replace w with 40). We use a centre-aligned tive to the body frame of the animal [45]. This requires index (compared to the rolling window of observations), three rotation sequences, using 3 by 3 rotation matrices with “extend” to indicate repetition of the leftmost or (Eqs. 4, 6) and involves two intermediate frames. The aer - rightmost non-NA value (though fill can equally be set as ospace sequence used here is as follows: NA, 0, etc.). Importantly, for correct trigonometric formulae out- 1. A right-handed rotation (C ), about the z-axis axis of put within the tilt-compensated compass method, the the device’s frame ( D ), through angle Ψ (Eq. 4), to get vectorial sum of static acceleration (G ) and calibrated x,y,z to the first intermediate frame (F ). magnetometry (M ) measurements across all three x,y,z 1 2. A right-handed rotation (C ) about the y-axis at F , spatial-dimensions must be normalised (to a unit vector) 1 through angle θ (Eq. 5), to get to the second interme- with a scaled magnitude (radius) of one (Eqs.  2, 3, R ). 5:10 diate frame ( F ). It was previously demonstrated that, for fast moving 2 Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 19 of 37 3. A right-handed rotation (C ) about the x-axis at F , The product of the conventionally used aerospace rota - through angle Φ (Eq.  6), to get to the animal’s body tion sequence outlined above (to get from the tag frame frame ( B). to the animal’s body frame) can be expressed as (Eq. 7).   (�,θ,�) (�) (θ) (�) cos (�) sin (�) 0 C = C · C · C (7) B/D B/F F /F F /D 2 2 1 1 (�)   − sin (�) cos (�) 0 C = (4) F /D 0 01 When matrix multiplied out, this yields (Eq.  8)—often referred to as a Direction Cosine Matrix (DCM). The   composition of this DCM varies according the (six pos- cos (θ) 0 − sin (θ) (θ) sible) orderings of the three rotation matrices (Eqs.  4,   01 0 C = (5) F /F 2 1 6) and the direction of intended rotation relative to the sin (θ) 0 cos (θ) direction of measured g within the NED system (see Additional file 1: Text S2).   cos (�) · cos (θ) sin (�) · cos (θ) − sin (θ) (�,θ ,�)   (8) C = cos (�) · sin (θ) · sin (�) − sin (�) · cos (�) sin (�) · sin (θ) · sin (�) + cos (�) · cos (�) cos (θ) · sin (�) B/D cos (�) · sin (θ) · cos (�) + sin (�) · sin (�) sin (�) · sin (θ) · cos (�) − cos (�) · sin (�) cos (θ) · cos (�) Note the left-handed rule of reading the vectorial nota-   (�,θ ,�) tion of ordered rotations, for example C means 10 0 B/D (�)   going from the device frame to the animal’s body frame, C = 0 cos (�) sin (�) (6) B/F by first rotating about the z-axis (though angle  ), fol- 0 − sin (�) cos (�) lowed by the y-axis (though angle θ) and then lastly the Note the right-handed rule of rotation; a positive x-axis (though angle  ). The device offset can be esti - reflects a clockwise rotation of the anterior–poste - mated from direct observation or deduced using photo- rior axis (relative to North), a positive θ reflects a nose- graphs or from the tag data itself. For example, assuming upward tilt of this axis and a positive  reflects a bank that ‘normal animal posture’ has no pitch and roll angle angle tilt to the right about this axis. Reversing the offset, then a tri-axial spherical plot of static acceleration direction of two axes causes a 180° inversion about the [164] would show a densely populated band of datapoints remaining axis and interchanging two axes (e.g., x with y) at the cross-sectional origin of 0  g about the x- and or reversing the direction of one or all three axes reverses y-axes, respectively, when the tag and body NED axes are the ‘handedness’ of rotation [right-handed—‘counter- in alignment. clockwise’ vs. left-handed—‘clockwise’ (when viewed NG and NM are pre-multiplied by the DCM to x,y,z x,y,z from the tip of the z-axis)]. Rotation matrices are orthog- compensate for offset. However, device offset is often onal (unitary), with every row and column being linearly parametrised by roll, pitch and/or yaw angles relative independent and normal to every other row and column. to the animal’s body frame and thus, the device actually The consequence of this is that the inverse of a rotation requires de-rotation (switching the ‘handedness’ of rota- matrix is its transpose [163] [which essentially reverses tion) according to these values. For example, a + 45° yaw the direction of rotation, and within (Eqs.  4, 6), this is offset requires an inverse rotation about the z-axis by achieved by negating the sign of the sines]. Importantly, − 45°, rather than a further + 45° rotation. This simply because rotation matrices are not symmetric, the order of involves taking the transpose of the DCM (Eq.  9), which matrix multiplication is important [45] (otherwise, Euler is the same as the transpose of each of the individual angles are without meaning for describing orientation). rotation matrices (Eqs. 10, 11).   cos � · cos θ cos � · sin θ · sin � − sin � · cos � cos � · sin θ · cos � + sin � · sin � ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (�,θ,�)   (9) sin � · cos θ sin � · sin θ · sin � + cos � · cos � sin � · sin θ · sin � − cos � · sin � C ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) B/D − sin (θ) cos (θ) · sin (�) cos (θ) · cos (�) Gunner et al. Anim Biotelemetry (2021) 9:23 Page 20 of 37     B D PitchSinAngle * RollSinAngle - YawSi- NGb NG x x T T T (�) (θ) (�) nAngle * RollCosAngle) + NMz *     NGb = C · C · C · NG (10) y y B/F F /F F /D 2 2 1 1 (YawCosAngle * PitchSinAngle * RollCo- NGb NG x x sAngle + YawSinAngle * RollSinAn- gle)     B D NMb NM (R21) NMby = NMx * YawSinAngle * Pitch- x x T T T (�) (θ) (�)     NMb = C · C · C · NM y y CosAngle + NMy * (YawSinAngle * B/F F /F F /D 2 2 1 1 NMb NM x x PitchSinAngle * RollSinAngle - YawSi- (11) nAngle * RollCosAngle) + NMz * (YawSinAngle * PitchSinAngle * RollSi- where T is the matrix transpose and resultant NGb x,y,z nAngle - YawCosAngle * RollSinAn- and NMb vectors are expressed in the animal’s body- x,y,z gle) carried NED frame. The input of these gravity and mag - (R22) NMbz = -NMx * PitchSinAngle + NMy netic vectors are supplied as 3 by 1 column matrices for * PitchCosAngle * RollSinAngle + true matrix multiplication, and when expanding out NMz * PitchCosAngle * RollCosAngle (Eq. 11), this results in (Eq. 12) [substituting NM with NG Here, Roll, Pitch and Yaw inputs denote the angular off - expands out (Eq. 10)].     B D NMb NM · cos (�) · cos (θ) + NM · (cos (�) · sin (θ) · sin (�) − sin (�) · cos (�)) + NM · (cos (�) · sin (θ) · cos (�) + sin (�) · sin (�)) x x y z      NMb   NM · sin (�) · cos (θ) + NM · (sin (�) · sin (θ) · sin (�) + cos (�) · cos (�)) + NM · (sin (�) · sin(θ) · sin(�) − cos(�) · sin(�))  y x y z NMb −NM · sin (θ) + NM · cos (θ) · sin (�) + NM · cos (θ) · cos(�) x x y z (12) In R then, the alignment of device to body axes for both set of the device, relative to the animal body frame. Note, gravity and magnetic vectors can be performed using the standard trigonometric functions operate in radians, not following procedure (R ). degrees. In base R, π = pi. Multiplying values by pi/180 11:22 coverts degrees into radians, whilst multiplying values (R11) RollSinAngle = sin(Roll * pi/180) by 180/pi does the reverse. This rotation correction pro - (R12) RollCosAngle = cos(Roll * pi/180) cedure is implemented within Gundog.Compass when (R13) PitchSinAngle = sin(Pitch * pitch, roll and/or yaw offsets are supplied (Additional pi/180) file 2). (R14) PitchCosAngle = cos(Pitch * Following the alignment of device and body axes, pitch pi/180) and roll of the animal are calculated from the DCM, and (R15) YawSinAngle = sin(Yaw * pi/180) because there are multiple variations in the order that (R16) YawCosAngle = cos(Yaw * pi/180) rotation sequences can be composed and applied, there (R17) NGbx = NGx * YawCosAngle * Pitch- are also different valid equations that output different CosAngle + NGy * (YawCosAngle * pitch and roll angle estimates, for equivalent static accel- PitchSinAngle * RollSinAngle - YawSi- eration input. The convention is to use formulae that nAngle * RollCosAngle) + NGz * have no dependence on yaw rotation and restrict either (YawCosAngle * PitchSinAngle * RollCo- the pitch or the roll angles within the range − 90° to + 90° sAngle + YawSinAngle * RollSinAn- (but not both), with the other axis of rotation able to lie gle) between − 180° and 180°, thereby eliminating duplicate (R18) NGby = NGx * YawSinAngle * Pitch- solutions at multiples of 360°. Multiplying (Eq.  8) by the CosAngle + NGy * (YawSinAngle * measured Earth’s gravitational field vector (+ 1  g when PitchSinAngle * RollSinAngle + YawCo- initially aligned downwards along the z-axis) simplifies sAngle * RollCosAngle) + NGz * down to (Eq.  13). The accelerometer output for this aer - (YawSinAngle * PitchSinAngle * RollSi- ospace rotation sequence is thus only dependent on the nAngle - YawCosAngle * RollSinAn- roll and pitch angles which can be solved (Eqs.  14, 15), gle) allowing roll angles the greater freedom [161]. This is (R19) NGbz = -NGx * PitchSinAngle + NGy relevant for studies using collar-mounted tags, whereby * PitchCosAngle * RollSinAngle + collar may roll > 90° in either direction from default NGz * PitchCosAngle * RollCosAngle orientation. (R20) NMbx = NMx * YawCosAngle * Pitch- CosAngle + NMy * (YawCosAngle * Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 21 of 37         0 cos (�) · cos (θ) sin (�) · cos (θ) − sin(θ) 0 − sin (θ)         NGb = · =  0   cos (�) · sin (θ) · sin (�) − sin (�) · cos (�) sin (�) · sin (θ) · sin (�) + cos (�) · cos (�) cos(θ) · sin(�)   0   cos (θ) · sin (�)  xyz 1 cos (�) · sin (θ) · cos (�) + sin (�) · sin (�) sin (�) · sin (θ) · cos (�) − cos (�) · sin (�) cos(θ) · cos(�) 1 cos (θ) · cos (�) (13)   � � � � −NGb 180   tan θ = � ⇒ θ = atan2 −NGb , NGb · NGb + NGb · NGb · (14) xyz x y y z z NGb + NGb 2 2 y z containing 1  s and − 1  s according to the direction of NGb 180 tan � = ⇒ � = atan2 NGb , NGb · xyz y z measured g from NGb (R ). NGb π 23 The magnetic vector of the device is then de-rotated (15) to the Earth frame (tilt-corrected) by pre-multiplying by The equation for roll (Eq.  15), however, has a region of the product of the inverse roll multiplied by inverse pitch instability at obtuse pitch angles (e.g., for NED systems, rotation matrix (Eq. 17), which when expanded out gives the x-axis points directly up or down, with respect to the (Eq. 18). Earth’s frame of reference). Whilst there is no ‘gold stand-     ard’ solution to this problem of singularity (using Euler NMbf cos (θ) sin (θ) · sin (�) sin (θ) · cos (�) angles), an attractive circumvention (detailed within     NMbf = 0 cos (�) − sin (�) [161]) is to modify (Eq. 15) and add a very small percent- NMbf − sin (θ) cos (θ) · sin (�) cos (θ) · cos (�)   age ( µ ) of the NGb reading into the denominator, pre- NMb venting it ever being zero and thus driving roll angles to   · NMb zero when pitch approaches −/+90° for stability (Eq. 16). NMb (17)     NMbf NMb · cos (θ) + NMb · sin (θ) · sin (�) + NMb · sin (θ) · cos (�) x x y z     NMbf = NM · cos (�) − NMb · sin (�)     y y z NMbf −NMb · sin (θ) + NMb · cos (θ) · sin (�) + NMb · cos (θ) · cos (�) z x y z (18) Here NMbf are the calibrated, normalised mag- x,y,z � = atan2 NGb , sign(NGb ) y z netometry data (expressed in the animal’s body-carried NED frame) after tilt-correction. Finally, yaw (ψ) (head- · (NGb · NGb + µ · NGb · NGb ) · z z x x ing—now defined by the compass convention, relative to (16) magnetic North) can be computed from the NMbf and where sign(NGb ) is allocated the value + 1 when NGb z z NMbf (Eq. 19) via; is non-negative and − 1, when NGb is negative (recovers directionality of NGb , subsequent to the square-root). ψ = atan2 −NMbf , NMbf · (19) y x Taken together then, in R, pitch and roll are computed according to, (R ) with outputs within the range of 24:25 We outline the R code for this procedure below (R ). 26:34 − 90° to + 90° for pitch and − 180° to + 180° for roll, and this is the formula we use in the tilt-compensated method (R26) RollSinAngle = sin(Roll * pi/180) outlined below (and within Additional file 2). (R27) RollCosAngle = cos(Roll * pi/180) (R28) PitchSinAngle = sin(Pitch * (R23) mu = 0.01 ; sign = ifelse(NGbz >= pi/180) 0, 1, -1) (R29) PitchCosAngle = cos(Pitch * (R24) Pitch = atan2(-NGbx, sqrt(NGby^2 pi/180) + NGbz^2)) * 180/pi (R30) NMbfx = NMbx * PitchCosAngle + (R25) Roll = atan2(NGby, sign * NMby * PitchSinAngle * RollSinAn- sqrt(NGbz^2 + mu * NGbx^2)) * gle + 180/pi NMbz * PitchSinAngle * RollCosAngle (R31) NMbfy = NMby * RollCosAngle – Here, prior to the derivation of pitch and roll, μ is allo- NMbz * RollSinAngle cated the value 0.01 and a vector termed ‘sign’ is created, Gunner et al. Anim Biotelemetry (2021) 9:23 Page 22 of 37 (R32) NMbfz = -NMbx * PitchSinAngle + acceleration (cf. Eq.  1, R ) from their raw equivalent 1:4 NMby * PitchCosAngle + RollSinAn- (R ). gle + NMbz * PitchCosAngle * RollCosAngle (R37) v = sqrt((Ax - Gx)^2 + (Ay - (R33) Yaw = atan2(-NMfby, NMfbx) * 180/ Gy)^2 + (Az - Gz)^2), pi (R34) Yaw = ifelse(Yaw < 0, Yaw + 360, where Ax, Ay, Az and Gx, Gy, Gz are the raw and static Yaw) (smoothed) values of each channel’s recorded acceleration. Note, yaw output from (R ) uses the scale − 180° to We recommend implementing a running mean (cf. + 180°. (R ) converts to the scale 0° to 360 (specifi - Eq.  1, R ) to raw VeDBA values to ensure that both 34 1:4 cally, 0° to 359°). This is also achieved by using a modulus acceleration and deceleration components of a stride (mod) operator (Eq.  20, R ), which in base R takes the cycle are incorporated together per unit time and to form %%. reduce the magnitude of small temporal spikes (likely not attributable to the scale of movement elicited [cf. 97]. ψ = mod(360 + ψ, 360) (20) Choice of smoothing window size is dependent on the scale of movement being investigated, though as a basic (R35) Yaw = (360 + Yaw) %% 360 rule, we suggest 1 to 2  s. For similar reasons, it is also worth post-smoothing raw pitch, roll and heading out- Magnetic declination is defined as the angle on the puts, although heading requires a circular mean (Eqs. 22, horizontal plane between magnetic north and true 23) [cf. 169]: north [165]. Prior to dead-reckoning, magnetic decli-   nation should be summed to heading values to convert n n � � � � � � 1 π 1 π from magnetic to true North [166]. There are many   θ = a tan 2 sin h · , cos h · p j j n 180 n 180 online sources to calculate the magnetic declination of j=i j=i an area [e.g., 167]. Notably, logical corrections may need (22) to be performed to ensure data does not exceed either circular direction after applying magnetic declination h = mod 360 + θ · , 360 p (23) (R ). where h and h are the unsmoothed and smoothed (R36) h = ifelse(h < 0, h + 360, h) ; h heading values, θ the arithmetic mean after converting = ifelse(h > 360, h - 360, h), degrees to cartesian coordinates and mod refers to the modulo operator. where h refers to the vector containing the heading In R, the above formula can be made into a function data. Should the user not correct for axis align- (R ), to be applied within the ‘rollapply’ wrapper (replac- ment between the device and animal body frame ing ‘FUN = mean’ with ‘FUN = Circ.Avg’) (cf. R ). 1:4 (cf. Eqs. 4–12, R ) then a reasonable post-cor- 11:22 rection for small discrepancies about the yaw axis (R38) Circ.Avg = function(x){ would be to subtract the difference to h values at this point. H.East = mean(sin(x * pi / 180)) Preparing speed estimates H.North = mean(cos(x * pi / 180)) The vectorial dynamic body acceleration (VeDBA) (Eq. 21) [cf. 67, 168] was our choice of DBA-based speed MH =(atan2(H.East, H.North)) * 180/pi proxy for terrestrial dead-reckoning purposes. This is given by: MH = (360 + MH) %% 360 2 2 2 v = D + D + D (21) x y z return(MH) where v represents VeDBA, D , D and D are the x y z } dynamic acceleration values from each axis, themselves Speed ( s ) can be estimated from VeDBA ( v ) via (Eq. 24). obtained by subtracting each axis’ static component of s = (v · m) + c (24) Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 23 of 37 where m is the multiplicative coefficient and c is a con- thus should only be calculated at times when the animal stant [10, 69]. Here, a user can define various bouts is travelling ‘ballistically’ (at considerable vertical speed). of movement from motion sensor data (e.g., via vari- (R41) s = ifelse(abs(p) >= 10, abs(RCD ous machine-learning approaches (for review see Far- / tan(p * pi/180)), s) rahi et al. [170]) or the Boolean-based LoCoD method [101]) and/or substrate condition (e.g., via GPS), to be In the above example (R ), nominal speed values are cross-referenced when allocating variants of the speed overwritten with the trigonometric formula output coefficients. As a simple example, in R, should walk- (Eq.  25) at times of ‘appreciable’ pitch (10°) [cf. 171], ing (coded for as 1) and running (coded for as 2) be where RCD is the rate change of depth and p is the pitch teased apart from all other (non-moving) data (coded (in radians). An upper limit should be imposed on speed for as 0) within a Marked Events vector (ME), then values derived in this way because values can become ME can be used to allocate various m (and if applica- highly inflated when the pitch angle is particularly acute. ble, c ) values using simple ‘ifelse’ statements (R ). 39:40 Converting speed to a distance coefficient (R39) m = ifelse(ME == 1, 1.5, Speed (s) estimates are multiplied by the time differ - ifelse(ME == 2, 3.5, 0)) ence between the values (TD ) to give a distance estimate (R40) c = ifelse(ME > 0, 0.1, 0) (units in metres) which, in turn, standardises coefficient comparisons across datasets sampled at different rates. These distance values are then divided by the approxi - Here, walking is given an arbitrary coefficient of 1.5 mate radius of the earth (R = 6,378,137 m) to give a radial and running, 3.5 with a value of 0.1 for their constants. distance coefficient ( q ) [se e 172] (Eq. 26)” All other ME values are given a 0 coefficient and 0 con - stant, which results in no speed at such times, regardless s · TD of DBA magnitude. q = (26) By‑passing DBA as a speed proxy Assuming that high-resolution depth data are not avail- Dividing the number of steps detected within a given able, but ‘absolute’ speed estimates have been obtained, rolling window length (cf. R ), by the window length then an alternative to Eq.  25, (in accordance with the 1:4 (s) gives an estimated step count per second. This can equal pitch assumption) is to derive horizontal distance be converted to speed by multiplying by a distance per estimates by multiplying the absolute distance by the step estimate (assuming constant distance travelled cosine of the pitch ( θ ) (converted from degrees to radi- between step gaits). We review this further in Addi- ans), which can equally be performed on the radial dis- tional file  1: Text S4, including a simple peak finder tance (Eq. 27): function—Gundog.Peaks (Additonal file 3) that locates peaks based on local signal maxima, using a given roll- q = q · cos θ · (27) ing window, with each candidate peak filtered accord- ing to whether it surpassed a threshold height (in In R, to determine accurate lengths of time between conjunction with other potential user-defined thresh- values, it is best to save date and time variables together olds). Note, this method can equally be applied to as POSIX class [173]. Creating timestamp (TS) objects non-terrestrial species, using flipper/tail beats instead, with POSIXct class enables greater control and manipula- where appropriate. tion of time data. This makes computing the rolling time For diving animals, a proxy for horizontal speed can difference ( TD ; units in seconds) between data points be obtained based on animal pitch and rate change in simple (R ): depth [47, 116]. Specifically, rate change of depth ( d ) (R42) TD = c(0, difftime(TS, lag(TS), (units in m/s) is divided by the tangent of pitch ( θ ) units = "secs")[-1]) (converted from degrees to radians) (Eq. 25): We detail how to create timestamp objects of POSIXct �d s = class within Additional file  1: Text S5, including format- (25) tan θ · ting with decimal seconds (important for infra-second datasets) and various codes useful for manipulating data Here, resultant speed values need to be made absolute to be dead-reckoned based on time. (positive). This calculation is only valid when the direction In R then, following the computation of TD , q is of movement is the same as the direction of the animal’s obtained via (R ). 43:44 longitudinal axis (equal pitch assumption) [cf. 47] and Gunner et al. Anim Biotelemetry (2021) 9:23 Page 24 of 37 (R43) s = (v * m) + c DR.lon[i] = DR.lon[i-1] + (R44) q = (s * TD) / 6378137 atan2(sin(h[i]) * sin(q[i]) * cos(DR.lat[i-1]), cos(q[i]) - sin(DR. Note, if a negative c intercept is used (e.g., to allow for lat[i-1]) * sin(DR.lat[i])) some body movement without translation), then any neg- ative speed values would need to be equated to zero as an additional step. Reverse dead‑reckoning As previously mentioned, the ME vector (progressive For this, firstly, the time difference is computed as usual movement coded by integer values greater than zero (R ) and the dimensions of each vector required in the (e.g., 1) and stationary behaviour coded by zero) can be dead-reckoning calculation are reversed. We bind all used to ensure q (essentially the distance moved) is zero relevant vectors into a data frame (df ) (R ), subsequent when ME reads zero, ensuring dead-reckoned tracks are to reversing data frame dimensions (R ); the last row not advanced at such times, regardless of the computed becomes the first row, second to last row becomes the speed (R ). second, etc. Note, this can equally be achieved by using the rev() function within base R, on each individual vec- (R45) q = ifelse(ME == 0, 0, q) tor. These reversed columns are now restored as vectors (R ) and shifted forward by one element (R ). This is 53 54 Derivation of coordinates required for correct alignment in time so that dead-reck- Once q and h are obtained, coordinates are advanced oning works in exactly the opposite manner to ‘forward’ using (Eqs. 28, 29); dead-reckoning. Lat = asin(sin Lat · cos q + cos Lat · sin q · cos h) i 0 0 (R50) TD = c(0, difftime(TS, lag(TS), (28) units = "secs")[-1]) Lon =Lon + atan2((sin h · sin q · cos Lat ), i 0 0 (R51) df = data.frame(TD, h, v, m, c, (cos q − sin Lat · sin Lat )) ME) 0 i (29) (R52) df = df[dim(df)[1]:1, ] where Lat , Lat and Lon , Lon are the previous and pre- 0 i 0 i (R53) TD = df[, ’TD’] ; h = df[, ’h’] ; sent latitude and longitude coordinates, respectively (in v = df[, ’v’] ; radians), h is the (present) heading (in radians) and q is m = df[, ’m’] ; c = df[, ’c’] ; ME = the (present) distance coefficient. df[, ’ME’] In R, the above can be performed iteratively within a (R54) TD = c(NA, TD[-length(TD)]) ; h = for-loop (iteration of code repeated per consecutive c(NA, h[-length(h)]) ; ith element of data; R ). Initialising the output latitude v = c(NA, v[-length(v)]) ; m = c(NA, (DR.lat) and longitude (DR.lon) variables to the required m[- length(m)]) ; length (e.g., as governed by the vector length of other c = c(NA, c[-length(c)]) ; ME = c(NA, input data (heading, speed, etc.) speeds up processing ME[-length(ME)]) time (R ). Within the trigonometric dead-reckoning for- The next step is to rotate heading 180° and correct for mulae, the starting latitude (la) and longitude (lo) coordi- its circular nature (R ). nates and heading ( h ) values must be supplied in radians (R ). The la and lo values are saved as the first elements (R55) h = h - 180 ; h = ifelse(h < 0, h of the DR.lat and DR.lon vectors to be advanced, respec- + 360, h) tively (R ). Lastly, is determined and DR.lon and DR.lat are advanced based on the dead-reckoning formula (cf. (R46) DR.lat = rep(NA, length(h)) ; R ), except in this instance, the first element of DR. 46:49 DR.lon = rep(NA, length(h)) lon and DR.lat needs to be supplied by the ‘known’ last lo (R47) la = la * pi/180 ; lo = lo * and la coordinates. pi/180 ; h = h * pi/180 (R48) DR.lat[1] = la DR.lon[1] = lo Integrating current vectors (R49) for(i in 2:length(DR.lat)) { In R, current vectors can be added according to (R ). 56:60 DR.lat[i] = asin(sin(DR.lat[i-1]) * Current speed (cs) is in m/s (ensure values are absolute) cos(q[i]) n+ cos(DR.lat[i-1]) * and current heading (ch) uses the scale 0° to 360°. Note sin(q[i]) * cos(h[i])) the use of ‘yy’ and ‘xx’ vectors, storing the previous DR. lat and DR.lon coordinates prior to implementing the Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 25 of 37 next ‘current drift’ vector per iteration. The current speed estimates from motion data—though they are essentially is also standardised according to the time period length the same). Lat − Lat Lon − Lon i 0 i 0 −1 2 2 (30) d = 2·R·sin sin + cos (Lat ) · cos (Lat ) · sin 0 i 2 2 where R is the Earth’s radius and d , the output in metres. sin (�Lon) · cos (Lat ), and Earth’s radius (analogous to the derivation of q ). i b = atan2 · cos (Lat ) · sin (�Lat) · cos (Lat ) · cos (�Lon) π 0 i When reverse dead-reckoning, it is important to ensure (31) that cs and ch are included in the steps outlined above (R ). where Lon represents Lon − Lon , Lat represents i 0 50:55 Lat − Lat and b output is in the scale − 180° to + 180°. To i 0 convert b to the conventional 0° to 360° scale, 360 should (R56) DR.lat = rep(NA, length(h)) ; be added to values < 0. DR.lon = rep(NA, length(h)) For each VP, the distance is divided by the dead-reck- (R57) xx <- rep(NA, length(cs)) ; yy <- oned distance providing a distance correction factor rep(NA, length(cs)) (ratio; Eq.  32). The heading correction factor is com - (R58) la = la * pi/180 ; lo = lo * puted by subtracting the dead-reckoned bearing from pi/180 ; the VP bearing (Eq.  33). To ensure that difference does h = h * pi/180 ; ch = ch * pi/180 not exceed 180° in either circular direction, 360 should (R59) DR.lat[1] = la DR.lon[1] = lo be added to values < − 180 and 360 subtracted from val- (R60) for(i in 2:length(DR.lat)) { ues > 180. A simple example of why this is relevant can be DR.lat[i] = asin(sin(DR.lat[i-1]) * illustrated by subtracting a dead-reckoned bearing value cos(q[i]) + cos(DR.lat[i-1]) * of 359° from a VP bearing value of 1°—post-correction, sin(q[i]) * cos(h[i])) the difference is + 2°. yy[i] = DR.lat[i] DR.lon[i] = DR.lon[i-1] + Distance VP atan2(sin(h[i]) * sin(q[i]) * Distance = corr.factor (32) Distance DR cos(DR.lat[i-1]), cos(q[i]) - sin(DR. lat[i-1]) * sin(DR.lat[i])) Heading = Bearing − Bearing (33) xx[i] = DR.lon[i] corr.factor VP DR DR.lat[i] = asin(sin(yy[i]) * All intermediate q values are multiplied by the distance cos((cs[i] * TD[i]) / 6378137) + correction factor and the heading correction factor is cos(yy[i]) * sin((cs[i] * TD[i]) / added to all intermediate h values (ensuring that h values 6378137) * cos(ch[i])) are in degrees). To ensure circular range is maintained DR.lon[i] = xx[i] + atan2(sin(ch[i]) * between 0° and 360°, 360 should be subtracted from val- sin((cs[i] * TD[i]) / ues > 360 and added to values < 0. 6378137) * cos(yy[i]), cos((cs[i] * Specifically, we follow the protocol illustrated within TD[i]) / 6378137) - sin(yy[i]) * Fig.  9 for intermediate values. Note the formulae to cal- sin(DR.lat[i])) culate both distance ( d ; Eq.  30) and bearing ( b ; Eq.  31) between two points, are also used to recalculate both the heading ( h ) and radial distance ( q ) between current-inte- VPC procedure grated dead-reckoned fixes (pre-VPC; cf. R ). Note that Specifically, this method entails calculating the difference the Haversine distance is required to be converted back of Haversine distance (net error) and bearing (from true to radial distance by dividing by R (R = 6,378,137). North) between consecutive VPs and the correspond- In R, the formulae to calculate the great-circle dis- ing time-matched dead-reckoned track positions. The tance and great circular bearing are saved within the trigonometric Haversine formulae (Eqs.  30, 31) are used disty (R ) and beary (R ) functions, respectively, where 61 62 to calculate the great-circle distance ( d ) and great circu- lon1, lat1, long2 and lat2 represent longitude and lati- lar bearing ( b ) between consecutive VPs and consecu- tude positions (decimal format) at t and t , ( t repre- i i+1 tive (time-matched) dead-reckoned positions (note we senting time). use the term ‘bearing’ to differentiate between heading Gunner et al. Anim Biotelemetry (2021) 9:23 Page 26 of 37 Fig. 9 Schematic diagram illustrating the order of fixes used when calculating the Distance and Heading (difference of both GPS corr.factor corr.factor and dead‑reckoned (DR) positions between arrow heads). Note the discrepancy with the order at which these correction factors are applied to intermediate DR positions (as denoted by colour shading). Known starting position denoted with asterisk (R61) disty = function(long1, lat1, d = 6378137 * c long2, lat2) { return(d) long1 = long1 * pi/180 ; long2 = long2 * pi/180 ; lat1 = lat1 * } (R62) beary = function(long1, lat1, pi/180 ; lat2 = lat2 * pi/180 long2, lat2) { long1 = long1 * pi/180 ; long2 = long2 * pi/180 ; lat1 = lat1 * pi/1 a = sin((lat2 - lat1) / 2) * sin((lat2 80 ; lat2 = lat2 * pi/180 - lat1) / 2) + cos(lat1) * a = sin(long2 - long1)*cos(lat2) b = cos(lat1) * sin(lat2) - sin(lat1) cos(lat2) * sin((long2 - long1) / 2) * * cos(lat2) * cos(long2 - long1) sin((long2 - long1) / 2) c = ((atan2(a, b) / pi)*180) return(c) c = 2 * atan2(sqrt(a), sqrt(1 - a)) } Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 27 of 37 Below, we outline an example of VPC in R and assume distances between VPs are stored within the column VP coordinates (decimal format) are aligned in the same termed VP.distance (R ). The VP.distance is divided by length vectors/columns as motion sensor-derived data, the DR.distance to provide the distance correction factor, e.g., heading, DBA/speed, etc., with the corresponding termed Dist.corr.factor (R ). Importantly here, an ifelse indexed (element-/row-wise) time. Typically, motion sen- statement is incorporated so that Dist.corr.factor defaults sor data are recorded at much higher frequency so that to zero at times when both VP.distance and DR.distance there are many dead-reckoned fixes between sequential are zero (otherwise dividing zero by zero in R produces VPs. As such, in the example below, we assume NAs are NaN’s). expressed in the VP longitude and latitude fields at times of missing locational data. This approach of synchronis - (R66) df.sub$DR.loni = c(df.sub[-1, ing VP—with motion sensor data also applies when inte- ’DR.longitude’], NA) grating current data; assuming ch and cs are element/ row-wise matched to the relevant VP grid node. df.sub$DR.lati = c(df.sub[-1, ’DR. Firstly, an indexing row number (Row.number) vector, latitude’], NA) the length of the data used in the dead-reckoning opera- tion (e.g., h ) is created (R ), which is relevant for merg- df.sub$VP.loni = c(df.sub[-1, ’VP. ing full-sized and under-sampled data frames together longitude’], NA) (seen later). Together, the row number, (uncorrected) dead-reckoned longitude and latitude coordinates, VP df.sub$VP.lati = c(df.sub[-1, ’VP. longitude and latitude coordinates, heading and the latitude’], NA) radial distance vectors are inputted column-wise into a (R67) df.sub$DR.distance= disty(df. ‘main’ data frame, termed ‘df ’ (R ; user-assigned column 64 sub$DR.longitude, names of each vector are within quotation marks). This df.sub$DR.latitude, df.sub$DR.loni, data frame is then filtered removing rows with missing df.sub$DR.lati) VP data and stored as df.sub (R ). This under-sampled 65 (R68) df.sub$VP.distance= disty(df. data frame thus, row-wise, contains the time-matched sub$VP.longitude, dead-reckoned and ground-truthed positions. The VPC df.sub$VP.latitude, df.sub$VP.loni, process is analogous for reverse dead-reckoned tracks— df.sub$VP.lati) although VP.lon and VP.lat must also be reversed [Row. (R69) df.sub$Dist.corr.factor = number remains in ascending order (not reversed)]. The ifelse(df.sub$VP.distance == 0 & first element of VP.lon and VP.lat must be the lo and la, df.sub$DR.distance == 0, 0, df.sub$VP. respectively (or for reverse dead-reckoning, the last ele- distance / df.sub$DR.distance) ment prior to reversing these vectors). Analogous to the distance correction, the bearings between consecutive dead-reckoned estimates are stored (R63) Row.number = rep(1:length(h)) within the column termed DR.head (R ) and the cor- (R64) df = data.frame(Row.number, ’DR. responding bearings between VPs are stored within the longitude’ = DR.lon, column termed VP.head (R ). Logical corrections are ’DR.latitude’ = DR.lat, ’VP.longitude’ performed to convert both to the 0° to 360° scale (R ), = VP.lon, DR.head is subtracted from VP.head providing the head- ’VP.latitude’ = VP.lat, h, q) ing correction factor, termed Head.corr.factor (R ) (R65) df.sub = df[!with(df, is.na(VP. and further logical corrections are performed to ensure longitude) | is.na(VP.latitude)) a minimum and maximum difference range between ,] − 180° to + 180 (R ). Both sets of dead-reckoned and VP coordinates are (R70) df.sub$DR.head = beary(df.sub$DR. shifted backwards one row within new columns termed; longitude, DR.loni, DR.lati, VP.loni, VP.lati (R ). Row-wise, these columns represent the consecutive fix at t with their i+1 df.sub$DR.latitude, df.sub$DR.loni, originals being t . This provides the correct format for df.sub$DR.lati) the inputs required within the disty (cf. R ) and beary 61 (R71) df.sub$VP.head = beary(df.sub$VP. (cf. R ) functions. The distances between consecu - 62 longitude, tive dead-reckoned estimates are stored within the col- df.sub$VP.latitude, df.sub$VP.loni, umn termed DR.distance (R ) and the corresponding 67 df.sub$VP.lati) Gunner et al. Anim Biotelemetry (2021) 9:23 Page 28 of 37 (R72) df.sub$DR.head = ifelse(df. coordinates, heading and radial distance each time) until sub$DR.head < 0, dead-reckoning fixes accord ‘exactly’ (Gundog.Tracks uses df.sub$DR.head + 360, df.sub$DR.head) a threshold of 0.01 m) with ground-truthed locations. An df.sub$VP.head = ifelse(df.sub$VP.head important pitfall of the correction process to consider is < 0, that dividing ‘any’ value (e.g., > 0) by 0 results in infinite df.sub$VP.head + 360, df.sub$VP.head) (Inf) values in R. This can arise during the correction (R73) df.sub$Head.corr.factor = process when there is a given distance between con- df.sub$VP.head - df.sub$DR.head secutive VPs, but no displacement between the accord- (R74) df.sub$Head.corr.factor = ing dead-reckoned positions. This can be a consequence ifelse(df.sub$Head.corr.factor < of ground-truthing too frequently (typically relevant to -180, high-res GPS studies), where positional noise is more (df.sub$Head.corr.factor + 360), apparent during rest periods [cf. 97] and/or wrongly df.sub$Head.corr.factor) assigned speed estimates/ME values. Gundog.Tracks df.sub$Head.corr.factor = ifelse(df. automatically resamples VPC rate where necessary to sub$Head.corr.factor > 180, avoid Inf values, essentially by changing the VPC rate to df.sub$Head.corr.factor - 360), avoid using successive VPs at times of no dead-reckoned df.sub$Head.corr.factor) track advancement. Lastly, Gundog.Tracks outputs mes- Only the relevant columns; Row.number, Dist.corr.fac- sages to the user’s console, detailing up to six stages of tor and Head.corr.factor are preserved (R ) and merged dead-reckoning progression, which includes reporting back into the main data frame (df ) based on the matching the maximum distance (units in metres) between dead- row numbers (R ). Both Dist.corr.factor and Head.corr. reckoned- and ground-truthed positions (used within factor express NA’s between VPs. These are replaced with the VPC procedure) at each iteration of correction and the most recent non-NA (observations carried forwards; whether automatic VPC resampling due to Inf values R ). Dist.corr.factor and Head.corr.factor values are occurred. shifted forward by one row (R ) for correct alignment purposes with respect to h and q values to be adjusted Conclusion (cf. Fig.  9; R ). A logical correction is performed to We have provided a comprehensive, fully integrated 79:80 ensure that a 0° to 360° circular scale is maintained after application of the dead-reckoning procedure within the the heading correction (R ). Note, the na.locf() function framework of the programming language, R, from pre- is required from the ‘zoo’ package, to replace NA values processing raw tri-axial accelerometery and magnetom- with the last non-NA value. etry data to VPC dead-reckoning. We have highlighted important considerations to increase the accuracy of the (R75) df.sub = df.sub[, c(’Row.number’, analytical procedure and to avoid misinterpretation of ’Dist.corr.factor’, ’Head.corr. error. We have also supplied extensive Additional files factor’)] 1, 2, 3, 4 and 5 and supporting functions to aid the pro- (R76) df = merge(df, df.sub, by = "Row. cess of deriving fine-scale movement paths, including number", all = TRUE) the protocols to correct magnetometry data and derive (R77) df$Dist.corr.factor = (tilt-compensated) heading. Importantly, we have dem- na.locf(df$Dist.corr.factor) onstrated the value of Gundog.Tracks; a multi-functional df$Head.corr.factor = na.locf(df$Head. and user-friendly tool to derive animal movement paths corr.factor) across all media of travel, with detailed input flexibil - (R78) df$Dist.corr.factor = c(NA, ity and output summaries. We suggest the next phase in df$Dist.corr.factor[-nrow(df)]) advancing the utility of animal dead-reckoning includes df$Head.corr.factor = c(NA, df$Head. looking for ‘track signatures’ that may signify a particu- corr.factor[-nrow(df)]) lar behaviour or reference a particular ‘ground-truthed’ (R79) q = (df$q * df$Dist.corr.factor) location. Lastly to advance the utility of Gundog.Tracks, (R80) h = (df$h + df$Head.corr.factor) we aim to optimise future iterations of the online code to (R81) h = ifelse(h > 360, h - 360, h) ; speed up computation time on larger datasets (e.g., sub- h = ifelse(h < 0, h + 360, h) second data collected over many months). Supplementary Information These updated coefficients are substituted into the The online version contains supplementary material available at https:// doi. dead-reckoning formula (cf. R ) and this process is 46:49 org/ 10. 1186/ s40317‑ 021‑ 00245‑z. repeated iteratively (using the updated dead-reckoned Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 29 of 37 Declarations Additional file 1. Methods expanded. Text S1. Device set up and cap ‑ ture protocol. Text S2. The importance of having the correct coordinate Ethics approval and consent to participate system and axis alignment. Text S3. Magnetometer calibration, rotation We thank the Conservation Agency from the Chubut Province, Argentina, correction and deriving yaw (heading)—Gundog.Compass() explained. for the permits to work at Punta León and Península Valdés protected areas Text S4. Step counts as a distance estimate—Gundog.Peaks() explained. (Disp No. 047/19‑SsCyAP). All penguin and cormorant handling procedures Text S5. Time Data in R (POSIXct). Text S6. VPC dead‑reckoning—Gundog. were reviewed and approved by the Dirección de Fauna y Flora Silvestre y el Tracks()—explained. Ministerio de Turismo y Áreas Protegidas de la Provincia de Chubut (permits Additional file 2. Gundog.Compass() (.R file). to work at San Lorenzo and Punta León, No. 060/19‑DFyFS‑MP and No. 047‑ SsCy/19). Ethical approval was also given by Animal Welfare Ethical Review Additional file 3. Gundog.Peaks() (.R file). Body (AWERB), approval number: SU‑Ethics‑Student ‑260919/1894, reference: Additional file 4. Gundog.Tracks() (.R file). IP‑1819‑30. Conditions and approvals for lion fieldwork were granted by the Animals Scientific Procedures Act (ASPA) at Queens University of Belfast (QUB‑ Additional file 5. Step by step guide of using Gundog.Tracks (.R file) (to BS‑AREC‑18‑006) and Pretoria University (NAS061‑19), permit authorisation use in conjunction with below). was given by South African National Parks (Permit Number SCAM 1550). Additional file 6. Raw sensor and GPS data frame (.txt) of a penguin walk ‑ ing out to sea from its nest. Consent for publication Not applicable. Acknowledgements Competing interests We thank South African National and the Department of Wildlife and National The authors declare no conflict of interest. Parks, Botswana, for allowing our research in the Kgalagadi Transfrontier Park. We are grateful to support and kind assistance of the staff and Rangers at the Kgalagadi National Park who were involved with this work, especially Steven Glossary Smith, Christa von Elling, Wayne Oppel and Corera Links. Pertaining to the field Acceleration The first der ivative (rate of change) work carried out in Argentina, we express our gratitude to Andrea Benvenuti, of an object’s velocity with respect to Fabian Gabelli, Monserrat Del Caño, La Chola, Miguel, Estancia El Pedral and Estancia San Lorenzo for assistance in various aspects of the research. We also time. Units are expressed as metres per thank the Instituto de Biología de Organismos Marinos (IBIOMAR‑ CONICET ) for second squared (m/s ) or in G-forces logistical support. HME is funded by an Irish Research Council Government of Ireland postgraduate scholarship. (g). A single G-force on Earth (though this does vary slightly with elevation) Authors’ contributions is 9.81 m/s . Tri-axial accelerometers The authors declare no conflict of interest and are in agreement to submit to Animal Biotelemetry. RMG conceived the study and RMG and RW wrote the measure acceleration in three orthogo- initial draft. PH constructed tag housings for all model species used. Data col‑ nal planes (surge—‘anterior–posterior’, lection for the lions was led by SF, DG, PV, LVS and AB and assisted by CJT, MFB, sway—‘medio-lateral’ and heave— DMS, SB, MVR, PH and RMG. Data collection for the penguins was led by FQ and data collection for the cormorants was led by AGL, with assistance from ‘dorsal–ventral’). Under non-moving KY, TY, and RMG. LB, MDH, RPW and RMG conceptualised the key considera‑ conditions, relative to gravity, the tions underlying the R code procedures and associated case‑studies, and RMG wrote the Gundog scripts and conducted the analysis of the case‑studies. device tilt (pitch and roll) can be cal- MHT supplied data for the Instantaneous tidal currents of the San Lorenzo culated directly from raw accelerom- region. All authors contributed to manuscript revision and LB, HJW, HME, AGL etery values since they are composed and LVS contributed to the testing and revision of the R syntax. All authors read and approved the final manuscript. entirely of the static force (gravity). Under linear acceleration, ‘moving’ Funding This research contributes to the CAASE project funded by King Abdullah forces applied to the device (e.g., due University of Science and Technology (KAUST ) under the KAUST Sensor to the animal moving) are superim- Initiative. Fieldwork in the Kgalagadi Transfrontier Park was supported in part posed to static readings and as such by a Department for Economy Global Challenges Research Fund grant to MS. Fieldwork within the Chubut Province was supported in part by the National measured animal acceleration is typi- Agency for Scientific and Technological Promotion of Argentina (PICT 2017‑ cally comprised of both a static and 1996 and PICT 2018‑1480), and the Grants‑in‑Aid for Scientific Research from the Japan Society for the Promotion of Science (16K18617). dynamic component. Barometric Pressure with the Earth’s atmosphere, Availability of data and materials pressure that is a measure of force per unit We provide a step‑by step example R script and an example data file of a Magellanic penguin walking out to sea (with time, raw acceleration and mag‑ area, often expressed as standard netometry data, marked events and aligned GPS positions) to demonstrate atmosphere (symbol: atm), defined as some of the key concepts outlined within “VPC dead‑reckoning procedure 101,325 Pa (1013.25 mbar; 1 Pa = 1 N/ in R” section when dead‑reckoning using Gundog.Tracks (including the initial calibration of magnetometry data with Gundog.Compass). All scripts, and the m ). The Earths mean sea-level atmos - example penguin data file have been uploaded to GitHub [95] and will be pheric pressure is approx. 1 atm. Baro- made available if the manuscript is accepted for publication. Online scripts will be continually updated and any queries, suggestions and/or reported bugs metric pressure decreases with eleva- should be emailed to the corresponding author. tion and increases with depth. Centripetal Inertial force caused by circular acceleration motion because an object is always Gunner et al. Anim Biotelemetry (2021) 9:23 Page 30 of 37 accelerating when either its direction (ECEF) system the Earth (this is often simplified to or magnitude (speed) changes, and in ‘Earth frame of reference’ or ‘Earths circular motion, the direction changes fixed frame’ in text). Its origin is fixed instantaneously. This can cause the at the Earth’s centre (the x-axis points animal to ‘pull g’, such as at times of towards the intersection of the Earth’s banking and cornering very fast. Greenwich Meridian and equatorial Coordinate In 3-D space, this is a set of three vectors plane, the y-axis pointing 90 degrees frame (x-, y-, z-axes) of unit length, perpen- East of the x-axis and the z pointing dicular (orthogonal) to each other. north, along the Earth’s rotation axis). Current flow (or ‘external’ current flow vectors). Note, this is different to the Earth- vectors The heading and speed of tidal-/ Centred Inertial (ECI) system, which air-currents. is non-rotating (and the x-axis instead Current Adding current flow vec tors to dead- always points towards the vernal integration reckoned travel vectors. equinox). De-rotation W ithin the tilt-compensated compass Equal pitch The animal moves in the same direction framework, this is the conversion of the assumption and angle as its anterior–posterior axis magnetic vector values through multi- (relative to North and the gravity vec- plying by the transpose (inverse) of the tor, respectively). pitch and then roll rotation matrices. Georeference W ithin the dead-reckoning frame- Distance The 2-D Haversine distance ratio work this is another term used for correction between successive Verified positions carrying out VPC-dead-reckoning, factor (VPs) (used in the VP-correction pro- or drift-correction or GPS-corrected cedure) and corresponding dead-reck- dead-reckoning. oned positions. This is multiplied to all Gimbal lock This is the loss of a degree of freedom intermediate (between VPs) radial dis- in 3-D, when two axes become paral- tance (q) values. lel to each other (locked in the same Drift The accumulation of spatial errors attitude, reflecting the same rotation). relative to a Verified Position, arising For example if the anterior–posterior from integrating incorrect dimensions axis (‘surge’ or ‘forward-back’—x-axis of travel. for NED coordinate frames) points in Dynamic The dynamic component of acceleration, the plane of the gravity vector (pitched body which is typically induced by the limb 90 degrees up or down), then the dor- acceleration and/or spine kinematics of the animal sal–ventral (‘heave’ or ‘up-down— (DBA) (and thus the attached accelerometer). z-axis for NED coordinate frames) and Generally, more mechanical work (via the medio-lateral axis (‘sway’ or ‘side- muscular contraction), corresponds to-side’—y-axis for NED coordinate to higher metabolic rate and greater frames) become parallel to each other, magnitudes in accelerometery read- and changes about the yaw can no ings (dependent on tag deployment longer be compensated for [changes site). Typically, dynamic values from in the roll (or ‘bank’) is equivalent to each multi-axial channel are integrated changing the heading]. into an overall metric, such as ‘Overall GPS-derived The Haversine distance calculated Dynamic Body Acceleration’ [ODBA = speed b etween successive GPS coordinates, │DBAx│ + │DBAy│ + │DBAz│] divided by the time taken between or ‘Vectorial Dynamic Body Accelera- locations. 2 2 tion’ [VeDBA = (DBAx + DBAy + D Ground- Empirical evidence (often information 2 0.5 BAz ) ]. Such derivatives have been truthing obtained by direct observation), as demonstrated as useful proxies for opposed to inference for validating movement-based power. something under investigation. Within Earth-Centre, This defines a non-inertial reference the dead-reckoning framework, VPs Earth-Fixed coordinate frame that rotates with such as GPS locations are used to Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 31 of 37 periodically ground-truth an animal’s system frame, the origin affixed as the devices position. centre of gravity and its axes oriented Haversine Computes the great-circular distance along the geodetic directions defined formula (units in metres) between two loca- by the Earth surface (the x- and y-axis tions (using their longitude and lati- pointing true north and East, respec- tude coordinates) on a sphere, apply- tively, parallel to the geoid surface ing trigonometry to map a triangle to and the z-axis pointing downwards the surface of a unit sphere. This for - towards the Earth’s surface). mula is only an approximation because Pitch The angle of device’s anterior–poste- the Earth is not a perfect sphere, with rior inclination or declination, relative numerical errors also arising at the to the horizontal plane of the Earth’s antipodal regions. surface. Pitch is often expressed as an Heading The difference of heading (or ‘bearing’) Euler angle, which describe the atti- correction factor f rom true North between consecutive tude and rotations of a device via a VPs (used in the VPC procedure) and given Euler angle sequence (yaw, pitch temporally aligned dead-reckoned posi- and roll) of rotations (using rotation tions. This is summed to all intermedi- matrices). Pitch can be derived from ate (between VPs) heading (h) values. the static component of accelera- Inf Results from numerical calculations tion. Assuming an NED system, pitch which are mathematically infinite (e.g., defines the degree of rotation about in R, dividing any value by zero results the y-axis. in Inf ). Radial Within the dead-reckoning framework, Linear drift At each path segment, the dead-reck- distance (q) this refers to a progression distance correction oned path is shifted to the position of accounting for the approximate cur- method the first VP encounter using a shift vature of the Earth (longlat projection vector. A correction vector then adds approximating the geoid to a sphere; the difference between the VP and radius (R) = 6,378,137 m). estimated dead-reckoned end points Right-handed The direction in which the linearly over this path segment period. coordinates/ fingers curl when pointing the right Multiplicative rotations thumb along the positive direction (m)-coefficient Within the dead-reckoning frame- (+ 1 g) of the z-axis (e.g., down for work, this refers to the gradient of a NED coordinates), reflect the direction linear regression, e.g., [speed =​ (Ve of rotation to be applied about each DBA· m) + c], where m is the multi- axis (for a given Euler angle sequence), plicative factor of VeDBA and c is the with the index finger representing the subsequently summed constant value x-axis and the middle finger repre- (reflecting the y-intercept). senting the y-axis, respectively, when NaN Non-n umeric (un-defined) values (e.g., splayed out at right angles to the in R, diving zero by zero results in thumb. NaN). Roll The angle of rotation about the device’s Net error Here, net error reflects the 2-D Haver - anterior–posterior axis. Roll is often sine distance (units in metres) between expressed as an Euler angle, which VPs and temporally aligned dead-reck- describe the attitude and rotations of a oned positions. device via a given Euler angle sequence Non-movement Behaviour performed while stationary, (yaw, pitch and roll) of rotations (using behaviours whereby the animal may be moving, rotation matrices). Roll is thus derived e.g., feeding on the spot, but there is after rotating by yaw and pitch and can no locomotion (not moving to a differ - be derived from the static component ent position in 3-D space). of acceleration. Assuming an NED sys- North–East– Often used in flight mechanics, tem, roll defines the degree of rotation Down (NED) this defines a non-inertial 3-D coordinate about the x-axis (also termed ‘bank Gunner et al. Anim Biotelemetry (2021) 9:23 Page 32 of 37 angle’). integration together. Assuming Cartesian coor- Static body The static component of dinates, vector addition is performed acceleration acceleration, due to gravity. Apart by adding the corresponding com- (SBA) from being used to calculate the angle ponents of the vectors together. E.g., of device tilt, increased inertial (cen- [A + B = (a + b , a + b , . . . , a + b )]. 1 1 2 2 n n tripetal) acceleration, e.g., when the Vertical speed Di stance travelled vertically up (at alti- animal ‘pulls g’, can be captured more tude) or down (at depth) divided by fully with static measures (rather than the time period between values. DBA estimates), and analogous to World Geodetic The typical model of the System 1984 VeDBA, the computation of the Vec- (WGS-84) E arth’s shape (standard for maps and torial Static Body Acceleration [VeSB satellite navigation), defining a coor - 2 2 2 0.5 A = (SBAx + SBAy + SBAz ) ] has dinate system that accounts for the been considered as a proxy of power. oblate spheroid. Tilt-compensated The compass heading compass Yaw The orientation of the device, gen- method (estimated using the arctangent ratio erally, with respect to true North between two orthogonal components (assuming any required magnetic dec- of the magnetic vector) is only accu- lination offset has been applied). Yaw, rate if the magnetometer outputs [typi- also termed ‘heading’ or ‘bearing’, is cally x, y channels—assuming the NED often expressed as an Euler angle, coordinate system is used (Additional which describe the attitude and rota- file  1: Text S2)] are taken when the tions of a device via a given Euler angle compass is level. Assuming the accel- sequence (yaw, pitch and roll) of rota- erometer-magnetometer approach, tions (using rotation matrices). Yaw static acceleration measures are used can be derived from the static com- to calculate the angles between the ponent of acceleration. Assuming an tag’s gravity (and thus magnetic) vec- NED system, yaw defines the degree of tor and the Earths frame of reference rotation about the z-axis. Yaw requires (e.g., Earth-Centered, Earth-Fixed the tilt-compensated compass method (ECEF) coordinate system). These to compute. angles are typically expressed as pitch Author details and roll Euler angles which are used to Swansea Lab for Animal Movement, Department of Biosciences, College of Science, Swansea University, Singleton Park, Swansea SA2 8PP, Wales, UK. compensate for variations in the mag- School of Biological Sciences, Queen’s University Belfast, Belfast, 19 Chlorine netometer output due to device tilt. Gardens, Belfast BT9 5DL, Northern Ireland, UK. Department of Agriculture, The tilt-compensated compass method Land Reform and Rural Development, Government of South Africa, Preto‑ ria 001, South Africa. Department of Migration, Max Planck Institute of Animal covers the procedures of adjusting the Behavior, 78315 Radolfzell, Germany. Department of Veterinary Tropical Dis‑ coordinate frame of the device to cor- eases, Faculty of Veterinary Science, University of Pretoria, Onderstepoort 0110, South Africa. School of Biology and Environmental Science, University College respond with a level inclination and Dublin, Belfield, Dublin, Ireland. Instituto de Biología de Organismos Marinos subsequently compute the compass (IBIOMAR), CONICET, Boulevard Brown 2915, U9120ACD Puerto Madryn, heading from the adjusted magnetom- Chubut, Argentina. Departamento de Ecología, Genética y Evolución & Insti‑ tuto de Ecología, Genética y Evolución de Buenos Aires (IEGEBA), CONICET, etry values. Pabellón II Ciudad Universitaria, C1428EGA Buenos Aires, Argentina. Centre Tortuosity The straight-line distance between for Biomathematics, College of Science, Swansea University, Swansea SA2 8PP, the start and end positions of a given UK. Graduate School of Environmental Studies, Nagoya University, Furo‑cho, Chikusa‑ku, Nagoya, Japan. Organization for the Strategic Coordination path segment, divided by the sum of of Research and Intellectual Properties, Meiji University, Nakano, Tokyo, Japan. the consecutive intermediate indi- Savanna and Grassland Research Unit, South African National Parks, Scientific Services Skukuza, Kruger National Park, Skukuza 1350, South Africa. Veteri‑ vidual distance steps that constituted nary Wildlife Services, South African National Parks, 97 Memorial Road, Old the total path segment’s length. Values Testing Grounds, Kimberley 8301, South Africa. Mammal Research Institute, closer to 0 (or conversely values closer Department of Zoology and Entomology, University of Pretoria, Pretoria 002, South Africa. Instituto Andino Patagónico de Tecnologías Biológicas y to 1 if subtracting the resultant ‘tortu- Geoambientales, Grupo GEA, IPATEC‑UNCO ‑ CONICET, San Carlos de Bariloche, osity’ value from 1) reflect more twists Río Negro, Argentina. Red Sea Research Centre, King Abdullah University of Science and Technology, Thuwal 23955, Saudi Arabia. Center for Zoo and turns in the movement path. and Wild Animal Health, Copenhagen Zoo, Roskildevej 38, 2000 Frederiksberg, Vector Adding vectors (of travel) Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 33 of 37 Denmark. Department of Zoology and Entomology, University of Fort Hare 20. Wilson RP, Shepard E, Liebsch N. Prying into the intimate details of Alice Campus, Ring Road, Alice 5700, South Africa. animal lives: use of a daily diary on animals. Endanger Species Res. 2008;4(1–2):123–37. Received: 8 March 2021 Accepted: 26 May 2021 21. Pedley M. eCompass‑build and calibrate a tilt ‑ compensating electronic compass. Circuit Cellar Mag Comput Appl. 2012;265:1–6. 22. Li Z, Li X, Wang Y. A calibration method for magnetic sensors and accelerometer in tilt‑ compensated digital compass. In: 2009 9th inter‑ national conference on electronic measurement & instruments, 16–19 Aug 2009; 2009. p. 2‑868‑862‑871. References 23. Ozyagcilar T. Implementing a tilt‑ compensated eCompass using 1. Browning E, Bolton M, Owen E, Shoji A, Guilford T, Freeman R. Predict‑ accelerometer and magnetometer sensors. Freescale semiconductor. ing animal behaviour using deep learning: GPS data alone accurately Application Note, AN4248; 2012. predict diving in seabirds. Methods Ecol Evol. 2018;9(3):681–92. 24. Gheorghe MV, Bodea MC. Calibration optimization study for tilt‑ 2. McClune DW. Joining the dots: reconstructing 3D environments compensated compasses. IEEE Trans Instrum Meas. 2018;67(6):1486–94. and movement paths using animal‑borne devices. Anim Biotelem. 25. Liu Y, Battaile BC, Trites AW, Zidek JV. Bias correction and uncertainty 2018;6(1):5. characterization of dead‑reckoned paths of marine mammals. Anim 3. Schlägel UE, Signer J, Herde A, Eden S, Jeltsch F, Eccard JA, Dammhahn Biotelem. 2015;3(1):51. M. Estimating interactions between individuals from concurrent animal 26. Dewhirst OP, Evans HK, Roskilly K, Harvey RJ, Hubel TY, Wilson AM. movements. Methods Ecol Evol. 2019;10(8):1234–45. Improving the accuracy of estimates of animal path and travel 4. Cagnacci F, Boitani L, Powell RA, Boyce MS. Animal ecology meets GPS‑ distance using GPS drift‑ corrected dead reckoning. Ecol Evol. based radiotelemetry: a perfect storm of opportunities and challenges. 2016;6(17):6210–22. Philos Trans R Soc B Biol Sci. 2010;365(1550):2157–62. 27. Mitani Y, Sato K, Ito S, Cameron MF, Siniff DB, Naito Y. A method for 5. Williams HJ, Taylor LA, Benhamou S, Bijleveld AI, Clay TA, de Grissac reconstructing three‑ dimensional dive profiles of marine mammals S, Demšar U, English HM, Franconi N, Gómez‑Laich A. Optimizing using geomagnetic intensity data: results from two lactating Weddell the use of biologgers for movement ecology research. J Anim Ecol. seals. Polar Biol. 2003;26(5):311–7. 2020;89(1):186–206. 28. Whitney NM, Pratt HL Jr, Pratt TC, Carrier JC. Identifying shark mating 6. Cotter CH. Early dead reckoning navigation. J Navig. 1978;31(1):20–8. behaviour using three‑ dimensional acceleration loggers. Endanger 7. Levi RW, Judd T. Dead reckoning navigational system using acceler‑ Species Res. 2010;10:71–82. ometer to measure foot impacts. In: Point Research Corporation, Santa 29. Lisovski S, Hewson CM, Klaassen RHG, Korner‑Nievergelt F, Kristensen Meijer, et al., Methods to assess physical activity with Ana, Calif., vol. MW, Hahn S. Geolocation by light: accuracy and precision affected by 5,583,776. 1996. p. 8. environmental factors. Methods Ecol Evol. 2012;3(3):603–12. 8. Beauregard S, Haas H. Pedestrian dead reckoning: a basis for personal 30. Miller PJO, Johnson MP, Madsen PT, Biassoni N, Quero M, Tyack PL. Using positioning. In: Proceedings of the 3rd workshop on positioning, navi‑ at‑sea experiments to study the effects of airguns on the foraging gation and communication; 2006. p. 27–35. behavior of sperm whales in the Gulf of Mexico. Deep Sea Res I Ocean‑ 9. Walker JS, Jones MW, Laramee RS, Holton MD, Shepard ELC, Williams ogr Res Pap. 2009;56(7):1168–81. HJ, Scantlebury DM, Marks NJ, Magowan EA, Maguire IE, Bidder OR, Di 31. Baumgartner MF, Freitag L, Partan J, Ball KR, Prada KE. Tracking large Virgilio A, Wilson RP. Prying into the intimate secrets of animal lives; soft‑ marine predators in three dimensions: the real‑time acoustic tracking ware beyond hardware for comprehensive annotation in ‘Daily Diary’ system. IEEE J Oceanic Eng. 2008;33(2):146–57. tags. Mov Ecol. 2015;3(1):29. 32. Williams LR, Fox DR, Bishop‑Hurley GJ, Swain DL. Use of radio frequency 10. Bidder OR, Walker JS, Jones MW, Holton MD, Urge P, Scantlebury DM, identification (RFID) technology to record grazing beef cattle water Marks NJ, Magowan EA, Maguire IE, Wilson RP. Step by step: reconstruc‑ point use. Comput Electron Agric. 2019;156:193–202. tion of terrestrial animal movement paths by dead‑reckoning. Mov Ecol. 33. Alexander JS, Zhang C, Shi K, Riordan P. A granular view of a snow 2015;3(1):23. leopard population using camera traps in Central China. Biol Conserv. 11. Wensveen PJ, Thomas L, Miller PJO. A path reconstruction method 2016;197:27–31. integrating dead‑reckoning and position fixes applied to humpback 34. English HM, Harvey L, Wilson RP, Gunner RM, Holton MD, Woodroffe R, whales. Mov Ecol. 2015;3(1):31. Börger L. Multi‑sensor biologgers and innovative training allow data 12. Andrzejaczek S, Gleiss AC, Lear KO, Pattiaratchi CB, Chapple TK, Meekan collection with high conservation and welfare value in zoos. https:// doi. MG. Biologging tags reveal links between fine ‑scale horizontal and org/ 10. 21203/ rs.3. rs‑ 562677/ v1. vertical movement behaviors in Tiger Sharks (Galeocerdo cuvier). Front 35. Fancy SG, Pank LF, Douglas DC, Curby CH, Garner GW. Satellite telem‑ Mar Sci. 2019;6:229. etry: a new tool for wildlife research and management. Washington, 13. Bowditch N. The American practical navigator, Bicentennial edition. D.C: Fish and Wildlife Service Washington DC; 1988. Bethesda: National Imagery and Mapping Agency; 2002. p. 879. 36. Soutullo A, Cadahía L, Urios V, Ferrer M, Negro JJ. Accuracy of light‑ 14. Wilson R, Wilson M‑P. Dead reckoning: a new technique for determining weight satellite telemetry: a case study in the Iberian Peninsula. J Wildl penguim movements at sea. Meeresforschung. 1988;32(2):155–8. Manag. 2007;71(3):1010–5. 15. Wilson RP, Wilson M‑PT, Link R, Mempel H, Adams NJ. Determination of 37. Catipovic JA. Performance limitations in underwater acoustic telemetry. movements of African penguins Spheniscus demersus using a compass IEEE J Oceanic Eng. 1990;15(3):205–16. system: dead reckoning may be an alternative to telemetry. J Exp Biol. 38. Hofman MPG, Hayward MW, Heim M, Marchand P, Rolandsen CM, 1991;157(1):557–64. Mattisson J, Urbano F, Heurich M, Mysterud A, Melzheimer J, Morellet 16. Wilson RP, Liebsch N, Davies IM, Quintana F, Weimerskirch H, Storch S, N, Voigt U, Allen BL, Gehr B, Rouco C, Ullmann W, Holand Ø, Jørgensen Lucke K, Siebert U, Zankl S, Müller G, Zimmer I, Scolaro A, Campagna C, NH, Steinheim G, Cagnacci F, Kroeschel M, Kaczensky P, Buuveibaatar Plötz J, Bornemann H, Teilmann J, McMahon CR. All at sea with animal B, Payne JC, Palmegiani I, Jerina K, Kjellander P, Johansson Ö, LaPoint tracks; methodological and analytical solutions for the resolution of S, Bayrakcismith R, Linnell JDC, Zaccaroni M, Jorge MLS, Oshima JEF, movement. Deep Sea Res II Top Stud Oceanogr. 2007;54(3):193–210. Songhurst A, Fischer C, Mc Bride RT Jr, Thompson JJ, Streif S, Sandfort R, 17. Mitani Y, Watanabe Y, Sato K, Cameron MF, Naito Y. 3D diving behavior Bonenfant C, Drouilly M, Klapproth M, Zinner D, Yarnell R, Stronza A, Wil‑ of Weddell seals with respect to prey accessibility and abundance. Mar mott L, Meisingset E, Thaker M, Vanak AT, Nicoloso S, Graeber R, Said S, Ecol Prog Ser. 2004;281:275–81. Boudreau MR, Devlin A, Hoogesteijn R, May‑ Junior JA, Nifong JC, Odden 18. Elkaim GH, Decker EB, Oliver G, Wright B. Marine mammal marker J, Quigley HB, Tortato F, Parker DM, Caso A, Perrine J, Tellaeche C, Zieba (MAMMARK) dead reckoning sensor for In‑Situ environmental monitor ‑ F, Zwijacz‑Kozica T, Appel CL, Axsom I, Bean WT, Cristescu B, Périquet ing. In: Proceedings of IEEE/ION PLANS 2006, San Diego, CA; 2006. p. S, Teichman KJ, Karpanty S, Licoppe A, Menges V, Black K, Scheppers 976–87. TL, Schai‑Braun SC, Azevedo FC, Lemos FG, Payne A, Swanepoel LH, 19. Wilson AD, Wikelski M, Wilson RP, Cooke SJ. Utility of biological sensor Weckworth BV, Berger A, Bertassoni A, McCulloch G, Šustr P, Athreya V, tags in animal conservation. Conserv Biol. 2015;29(4):1065–75. Gunner et al. Anim Biotelemetry (2021) 9:23 Page 34 of 37 Bockmuhl D, Casaer J, Ekori A, Melovski D, Richard‑Hansen C, van de 59. Goldbogen JA, Calambokidis J, Shadwick RE, Oleson EM, McDonald MA, Vyver D, Reyna‑Hurtado R, Robardet E, Selva N, Sergiel A, Farhadinia MS, Hildebrand JA. Kinematics of foraging dives and lunge‑feeding in fin Sunde P, Portas R, Ambarli H, Berzins R, Kappeler PM, Mann GK, Pyritz whales. J Exp Biol. 2006;209(7):1231–44. L, Bissett C, Grant T, Steinmetz R, Swedell L, Welch RJ, Armenteras D, 60. Zimmer WM, Tyack PL, Johnson MP, Madsen PT. Three‑ dimensional Bidder OR, González TM, Rosenblatt A, Kachel S, Balkenhol N. Right on beam pattern of regular sperm whale clicks confirms bent ‑horn track? Performance of satellite telemetry in terrestrial wildlife research. hypothesis. J Acoust Soc Am. 2005;117(3):1473–85. PLoS ONE. 2019;14(5):e0216223. 61. Wilson RP, Ropert‑ Coudert Y, Kato A. Rush and grab strategies in forag‑ 39. Newey S, Davidson P, Nazir S, Fairhurst G, Verdicchio F, Irvine RJ, van der ing marine endotherms: the case for haste in penguins. Anim Behav. Wal R. Limitations of recreational camera traps for wildlife manage‑ 2002;63(1):85–95. ment and conservation research: a practitioner’s perspective. Ambio. 62. Gabaldon J, Turner EL, Johnson‑Roberson M, Barton K, Johnson M, 2015;44(4):624–35. Anderson EJ, Shorter KA. Integration, calibration, and experimental 40. Leonardo M, Noss AJ, Erika C, Damián IR. Ocelot (Felis pardalis) popula‑ verification of a speed sensor for swimming animals. IEEE Sens J. tion densities, activity, and ranging behaviour in the dry forests of east‑ 2019;19(10):3616–25. ern Bolivia: data from camera trapping. J Trop Ecol. 2005;21(3):349–53. 63. Denny M. Air and water: the biology and physics of life’s media. Prince‑ 41. Lewis JS, Rachlow JL, Garton EO, Vierling LA. Eec ff ts of habitat on GPS ton University Press; 1993. collar performance: using data screening to reduce location error. J 64. Altynay K, Khan Mohammed A, Marengo M, Swanepoel L, Przybysz A, Appl Ecol. 2007;44(3):663–71. Muller C, Fahlman A, Buttner U, Geraldi NR, Wilson RP, Duarte CM, Kosel 42. Kaur M, Sandhu M, Mohan N, Sandhu PS. RFID technology principles, J. Wearable multifunctional printed graphene sensors. NPJ Flex Electron. advantages, limitations & its applications. Int J Comput Electr Eng. 2019;3(1):1–10. 2011;3(1):151. 65. Williams H, Shepard E, Duriez O, Lambertucci SA. Can accelerometry 43. Davis RW, Fuiman LA, Williams TM, Le Boeuf BJ. Three‑ dimensional be used to distinguish between flight types in soaring birds? Anim movements and swimming activity of a northern elephant seal. Comp Biotelem. 2015;3(1):1–11. Biochem Physiol A Mol Integr Physiol. 2001;129(4):759–70. 66. Williams HJ, Shepard E, Holton MD, Alarcón P, Wilson R, Lambertucci S. 44. Wilson RP. Movements in Adélie Penguins foraging for chicks at Ardley Physical limits of flight performance in the heaviest soaring bird. Proc Island, Antarctica: circles within spirals, wheels within wheels. Polar Natl Acad Sci. 2020;117(30):17884–90. Biosci. 2002;15:75–87. 67. Wilson RP, Börger L, Holton MD, Scantlebury DM, Gómez‑Laich A, 45. Johnson MP, Tyack PL. A digital acoustic recording tag for measuring Quintana F, Rosell F, Graf PM, Williams H, Gunner R, Hopkins L, Marks the response of wild marine mammals to sound. IEEE J Oceanic Eng. N, Geraldi NR, Duarte CM, Scott R, Strano MS, Robotka H, Eizaguirre 2003;28(1):3–12. C, Fahlman A, Shepard ELC. Estimates for energy expenditure in free‑ 46. Shiomi K, Sato K, Mitamura H, Arai N, Naito Y, Ponganis PJ. Eec ff t of living animals using acceleration proxies: a reappraisal. J Anim Ecol. ocean current on the dead‑reckoning estimation of 3‑D dive paths of 2020;89(1):161–72. emperor penguins. Aquat Biol. 2008;3(3):265–70. 68. Bidder OR, Qasem LA, Wilson RP. On higher ground: how well can 47. Laplanche C, Marques TA, Thomas L. Tracking marine mammals in 3D dynamic body acceleration determine speed in variable terrain? PLoS using electronic tag data. Methods Ecol Evol. 2015;6(9):987–96. ONE. 2012;7(11):e50556. 48. Adachi T, Costa DP, Robinson PW, Peterson SH, Yamamichi M, Naito 69. Bidder OR, Soresina M, Shepard ELC, Halsey LG, Quintana F, Gómez‑ Y, Takahashi A. Searching for prey in a three‑ dimensional environ‑ Laich A, Wilson RP. The need for speed: testing acceleration for estimat‑ ment: hierarchical movements enhance foraging success in northern ing animal travel rates in terrestrial dead‑reckoning systems. Zoology. elephant seals. Funct Ecol. 2017;31(2):361–9. 2012;115(1):58–64. 49. Bras YL, Jouma’a J, Guinet C. Three‑ dimensional space use during the 70. Markovska M, Svensson R. Evaluation of drift correction strategies for bottom phase of southern elephant seal dives. Mov Ecol. 2017;5(1):18. an inertial based dairy cow positioning system: a study on tracking the 50. Andrzejaczek S, Gleiss AC, Pattiaratchi CB, Meekan MG. First Insights position of dairy cows using a foot mounted IMU with drift correction Into the fine ‑scale movements of the Sandbar Shark, Carcharhinus from ZUPT or sparse RFID locations. Student thesis. 2019. plumbeus. Front Mar Sci. 2018;5:483. 71. di Virgilio A, Morales JM, Lambertucci SA, Shepard EL, Wilson RP. Multi‑ 51. Aoki K, Amano M, Mori K, Kourogi A, Kubodera T, Miyazaki N. Active dimensional precision livestock farming: a potential toolbox for sustain‑ hunting by deep‑ diving sperm whales: 3D dive profiles and maneuvers able rangeland management. PeerJ. 2018;6:e4867. during bursts of speed. Mar Ecol Prog Ser. 2012;444:289–301. 72. Shepard EL, Wilson RP, Halsey LG, Quintana F, Laich AG, Gleiss AC, Lieb‑ 52. Narazaki T, Sato K, Abernathy K, Marshall G, Miyazaki N. Sea turtles sch N, Myers AE, Norman B. Derivation of body motion via appropriate compensate deflection of heading at the sea surface during directional smoothing of acceleration data. Aquat Biol. 2008;4(3):235–41. travel. J Exp Biol. 2009;212(24):4019–26. 73. Noda T, Okuyama J, Koizumi T, Arai N, Kobayashi M. Monitoring attitude 53. Benoit‑Bird KJ, Battaile BC, Nordstrom CA, Trites AW. Foraging behavior and dynamic acceleration of free‑moving aquatic animals using a of northern fur seals closely matches the hierarchical patch scales of gyroscope. Aquat Biol. 2012;16(3):265–76. prey. Mar Ecol Prog Ser. 2013;479:283–302. 74. Wilson JW, Mills MG, Wilson RP, Peters G, Mills ME, Speakman JR, 54. Ware C, Friedlaender AS, Nowacek DP. Shallow and deep lunge feeding Durant SM, Bennett NC, Marks NJ, Scantlebury M. Cheetahs, Acinonyx of humpback whales in fjords of the West Antarctic Peninsula. Mar jubatus, balance turn capacity with pace when chasing prey. Biol Lett. Mamm Sci. 2011;27(3):587–605. 2013;9(5):20130620. 55. Wright BM, Ford JK, Ellis GM, Deecke VB, Shapiro AD, Battaile BC, Trites 75. McNarry MA, Wilson RP, Holton MD, Griffiths IW, Mackintosh KA. Inves‑ AW. Fine‑scale foraging movements by fish‑ eating killer whales (Orcinus tigating the relationship between energy expenditure, walking speed orca) relate to the vertical distributions and escape responses of salmo‑ and angle of turning in humans. PLoS ONE. 2017;12(8):e0182333. nid prey (Oncorhynchus spp.). Mov Ecol. 2017;5(1):1–18. 76. Kaniewski P, Kazubek J. Integrated system for heading determination. 56. Wensveen PJ, Isojunno S, Hansen RR, von Benda‑Beckmann AM, Klei‑ Acta Phys Pol Ser A Gen Phys. 2009;116(3):325. vane L, van IJsselmuide S, Lam FPA, Kvadsheim PH, DeRuiter SL, Curé C, 77. Noda T, Kawabata Y, Arai N, Mitamura H, Watanabe S. Animal‑mounted Narazaki T, Tyack PL, Miller PJO. Northern bottlenose whales in a pristine gyroscope/accelerometer/magnetometer: in situ measurement of the environment respond strongly to close and distant navy sonar signals. movement performance of fast‑start behaviour in fish. J Exp Mar Biol Proc R Soc B Biol Sci. 1899;2019(286):20182592. Ecol. 2014;451:55–68. 57. Ropert‑ Coudert Y, Wilson RP. Trends and perspectives in animal‑ 78. Patonis P, Patias P, Tziavos IN, Rossikopoulos D, Margaritis KG. A fusion attached remote sensing. Front Ecol Environ. 2005;3(8):437–44. method for combining low‑ cost IMU/magnetometer outputs for use in 58. Francisco FA, Nührenberg P, Jordan A. High‑resolution, non‑invasive applications on mobile devices. Sensors. 2018;18(8):2616. animal tracking and reconstruction of local environment in aquatic 79. Diebel J. Representing attitude: Euler angles, unit quaternions, and ecosystems. Mov Ecol. 2020;8(1):27. rotation vectors. Technical Report, vol. 58. Stanford: Stanford University; 2006. p. 1–35. Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 35 of 37 80. Han S, Wang J. A novel method to integrate IMU and magnetometers in DM, Duarte CM. Give the machine a hand: a Boolean time‑based attitude and heading reference systems. J Navig. 2011;64(4):727–38. decision‑tree template for rapidly finding animal behaviours in multi‑ 81. Säll J, Merkel J. Indoor navigation using accelerometer and magnetom‑ sensor data. Methods Ecol Evol. 2018;9(11):2206–15. eter. Student thesis. 2011. 102. Rong L, Zhiguo D, Jianzhong Z, Ming L. Identification of individual 82. Alaimo A, Artale V, Milazzo C, Ricciardello A. Comparison between walking patterns using gait acceleration. In: 2007 1st international con‑ Euler and quaternion parametrization in UAV dynamics. AIP Conf Proc. ference on bioinformatics and biomedical engineering, 6–8 July 2007; 2013;1558(1):1228–31. 2007. p. 543–6. 83. Valenti RG, Dryanovski I, Xiao J. Keeping a good attitude: a 103. Wilson RP, Hustler K, Ryan PG, Burger AE, Noldeke EC. Diving birds in quaternion‑based orientation filter for IMUs and MARGs. Sensors. cold water: do Archimedes and Boyle determine energetic costs? Am 2015;15(8):19302–30. Nat. 1992;140(2):179–200. 84. Feng K, Li J, Zhang X, Shen C, Bi Y, Zheng T, Liu J. A new quaternion‑ 104. Aoki K, Watanabe YY, Crocker DE, Robinson PW, Biuw M, Costa DP, based Kalman filter for real‑time attitude estimation using the two ‑step Miyazaki N, Fedak MA, Miller PJO. Northern elephant seals adjust gliding geometrically‑intuitive correction algorithm. Sensors. 2017;17(9):2146. and stroking patterns with changes in buoyancy: validation of at‑sea 85. Chiella ACB, Teixeira BOS, Pereira GAS. Quaternion‑based robust metrics of body density. J Exp Biol. 2011;214(17):2973–87. attitude estimation using an adaptive unscented Kalman filter. Sensors. 105. Williams HJ, Holton MD, Shepard ELC, Largey N, Norman B, Ryan PG, 2019;19(10):2372. Duriez O, Scantlebury M, Quintana F, Magowan EA, Marks NJ, Alagaili 86. Gunner RM, Wilson RP, Holton MD, Scott R, Hopkins P, Duarte CM. A new AN, Bennett NC, Wilson RP. Identification of animal movement patterns direction for differentiating animal activity based on measuring angular using tri‑axial magnetometry. Mov Ecol. 2017;5(1):6. velocity about the yaw axis. Ecol Evol. 2020;10(14):7872–86. 106. Whitford M, Klimley AP. An overview of behavioral, physiological, and 87. Wiltschko R. Magnetic orientation in animals, vol. 33. Berlin: Springer environmental sensors used in animal biotelemetry and biologging Science & Business Media; 2012. studies. Anim Biotelem. 2019;7(1):1–24. 88. Sequeira MM, Rickenbach M, Wietlisbach V, Tullen B, Schutz Y. 107. Shamoun‑Baranes J, Bom R, van Loon EE, Ens BJ, Oosterbeek K, Bouten Physical activity assessment using a pedometer and its comparison W. From sensor data to animal behaviour: an oystercatcher example. with a questionnaire in a large population survey. Am J Epidemiol. PLoS ONE. 2012;7(5):e37997. 1995;142(9):989–99. 108. Wilson R. The Jackass Penguin (Spheniscus demersus) as a pelagic preda‑ 89. Kerdok AE, Biewener AA, McMahon TA, Weyand PG, Herr HM. Energetics tor. Mar Ecol Prog Ser. 1985;25(3):219–27. and mechanics of human running on surfaces of different stiffnesses. J 109. Culik BM, Wilson RP, Dannfeld R, Adelung D, Spairani HJ, Coria NRC. Appl Physiol. 2002;92(2):469–78. Pygoscelid penguins in a swim canal. Polar Biol. 1991;11(4):277–82. 90. Halsey LG, Shepard ELC, Hulston CJ, Venables MC, White CR, Jeukend‑ 110. Bethge P, Nicol S, Culik BM, Wilson RP. Diving behaviour and rup AE, Wilson RP. Acceleration versus heart rate for estimating energy energetics in breeding little penguins (Eudyptula minor). J Zool. expenditure and speed during locomotion in animals: tests with an 1997;242(3):483–502. easy model species, Homo sapiens. Zoology. 2008;111(3):231–41. 111. Tobalske B, Dial K. Flight kinematics of black‑billed magpies and 91. Miwa M, Oishi K, Nakagawa Y, Maeno H, Anzai H, Kumagai H, Okano K, pigeons over a wide range of speeds. J Exp Biol. 1996;199(2):263–80. Tobioka H, Hirooka H. Application of overall dynamic body accelera‑ 112. Sato K, Sakamoto KQ, Watanuki Y, Takahashi A, Katsumata N, Bost C‑A, tion as a proxy for estimating the energy expenditure of grazing farm Weimerskirch H. Scaling of soaring seabirds and implications for flight animals: relationship with heart rate. PLoS ONE. 2015;10(6):e0128042. abilities of giant pterosaurs. PLoS ONE. 2009;4(4):e5400. 92. Chapman Jason W, Klaassen Raymond HG, Drake VA, Fossette S, Hays 113. Scharold J, Lai NC, Lowell W, Graham J. Metabolic rate, heart rate, and Graeme C, Metcalfe Julian D, Reynolds Andrew M, Reynolds Don R, tailbeat frequency during sustained swimming in the leopard shark Alerstam T. Animal orientation strategies for movement in flows. Curr Triakis semifasciata. Exp Biol. 1989;48(4):223–30. Biol. 2011;21(20):R861–70. 114. Kawabe R, Kawano T, Nakano N, Yamashita N, Hiraishi T, Naito Y. Simul‑ 93. R Development Core Team. R—a language and environment for statis‑ taneous measurement of swimming speed and tail beat activity of free‑ tical computing. Vienna: R Foundation for Statistical Computing; 2021. swimming rainbow trout Oncorhynchus mykiss using an acceleration 94. The R project for statistical computing. https:// www.r‑ proje ct. org/. data‑logger. Fish Sci. 2003;69(5):959–65. Accessed 07 Mar 2021. 115. Steinhausen MF, Steffensen JF, Andersen NG. Tail beat frequency as 95. Gundog.Tracks GitHub database. Available at https:// github. com/ Richa a predictor of swimming speed and oxygen consumption of saithe rd6195/ Dead‑ recko ning‑ animal‑ movem ents‑ in‑R. Accessed 09 June (Pollachius virens) and whiting (Merlangius merlangus) during forced 2021 swimming. Mar Biol. 2005;148(1):197–204. 96. Dowle M, Srinivasan A, Gorecki J, Chirico M, Stetsenko P, Short T, 116. Miller PJO, Johnson MP, Tyack PL, Terray EA. Swimming gaits, passive Lianoglou S, Antonyan E, Bonsch M, Parsonage H. Package ‘data.table’. drag and buoyancy of diving sperm whales Physeter macrocephalus. J Extension of ‘data frame; 2019. Exp Biol. 2004;207(11):1953–67. 97. Gunner RM, Wilson RP, Holton MD, Hopkins P, Bell SH, Marks NJ, Bennett 117. Laich AG, Wilson RP, Quintana F, Shepard EL. Identification of imperial NC, Ferreira S, Govender D, Viljoen P, Bruns A, Schalkwyk Lv, Bertelsen cormorant Phalacrocorax atriceps behaviour using accelerometers. MF, Duarte CM, Rooyen MCv, Tambling CJ, Goppert A, Diesel D, Scantle‑ Endanger Species Res. 2008;10:29–37. bury DM. Decision rules for determining terrestrial movement and the 118. Chakravarty P, Cozzi G, Ozgul A, Aminian K. A novel biomechanical consequences for filtering high‑resolution GPS tracks—a case study approach for animal behaviour recognition using accelerometers. using the African Lion (Panthera leo). Mov Ecol. (in review). Methods Ecol Evol. 2019;10(6):802–14. 98. Dunford CE, Marks NJ, Wilmers CC, Bryce CM, Nickel B, Wolfe LL, 119. Chakravarty P, Maalberg M, Cozzi G, Ozgul A, Aminian K. Behavioural Scantlebury DM, Williams TM. Surviving in steep terrain: a lab‑to ‑field compass: animal behaviour recognition using magnetometers. Mov assessment of locomotor costs for wild mountain lions (Puma concolor). Ecol. 2019;7(1):28. Mov Ecol. 2020;8(1):34. 120. Hochscheid S. Why we mind sea turtles’ underwater business: a review 99. Wilson RP, Rose KA, Gunner R, Holton M, Marks NJ, Bennett NC, Bell SH, on the study of diving behavior. J Exp Mar Biol Ecol. 2014;450:118–36. Twining JP, Hesketh J, Duarte CM, Bezodis N, Scantlebury DM. Forces 121. Constandache I, Bao X, Azizyan M, Choudhury RR. Did you see Bob? experienced by instrumented animals depend on lifestyle. bioRxiv. human localization using mobile phones. In: Proceedings of the 2020. https:// doi. org/ 10. 1101/ 2020. 08. 20. 258756. sixteenth annual international conference on Mobile computing and 100. Willener AST, Handrich Y, Halsey LG, Strike S. Eec ff t of walking speed networking; Chicago, Illinois, USA. Association for Computing Machin‑ on the gait of king penguins: an accelerometric approach. J Theor Biol. ery; 2010. p. 149–60. 2015;387:166–73. 122. Symington A, Trigoni N. Encounter based sensor tracking. In: Proceed‑ 101. Wilson RP, Holton MD, di Virgilio A, Williams H, Shepard ELC, Lamber‑ ings of the thirteenth ACM international symposium on mobile ad hoc tucci S, Quintana F, Sala JE, Balaji B, Lee ES, Srivastava M, Scantlebury Gunner et al. Anim Biotelemetry (2021) 9:23 Page 36 of 37 networking and computing, Hilton Head, South Carolina, USA. Associa‑ 146. Wilson R, Griffiths I, Legg P, Friswell M, Bidder O, Halsey L, Lambertucci tion for Computing Machinery; 2012. p. 15–24. SA, Shepard E. Turn costs change the value of animal search paths. Ecol 123. Harja YD, Sarno R. Determine the best option for nearest medical Lett. 2013;16(9):1145–50. services using Google maps API, Haversine and TOPSIS algorithm. In: 147. Potts J, Börger L, Scantlebury DM, Bennett NC, Alagaili A, Wilson RP. 2018 international conference on information and communications Finding turning‑points in ultra‑high‑resolution animal movement data. technology (ICOIACT ), 6–7 March 2018; 2018. p. 814–9. Methods Ecol Evol. 2018;9(10):2091–101. 124. Great circle distances and bearings between two locations. https:// 148. Narazaki T, Nakamura I, Aoki K, Iwata T, Shiomi K, Luschi P, Suganuma dtcen ter. org/ met/ users/ docs/ write_ ups/ gc_ simple. pdf. H, Meyer CG, Matsumoto R, Bost CA, Handrich Y, Amano M, Okamoto 125. Robusto CC. The Cosine‑Haversine formula. Am Math Mon. R, Mori K, Ciccione S, Bourjea J, Sato K. Similar circling movements 1957;64(1):38–40. observed across marine megafauna taxa. iScience. 2021;24:102221. 126. Zeileis A, Grothendieck G, Ryan JA, Andrews F, Zeileis MA. Package ‘zoo’. 149. Kunčar A, Sysel M, Urbánek T. Calibration of low‑ cost triaxial magnetom‑ 2020. eter. In: MATEC web of conferences 20th international conference on 127. Jockers ML, Thalken R. Introduction to dplyr. In: Text analysis with r for circuits, systems, communications and computers (CSCC 2016), 2016. students of literature. Cham: Springer; 2020. p. 121–32. EDP Sciences. 128. Calculating Azimuth, distance, and altitude from a pair of GPS locations 150. Sodhi R, Prunty J, Hsu G, Oh B. Automatic calibration of a three‑axis in JavaScript. https:// medium. com/ javas cript‑ in‑ plain‑ engli sh/ calcu lat ‑ magnetic compass. U.S. Patent 7,451,549 PNI Corporation (Santa Rosa, ing‑ azimu th‑ dista nce‑ and‑ altit ude‑ from‑a‑ pair‑ of‑ gps‑ locat ions‑ 36b43 CA, US); 2008. 25d8a b0. 151. Renaudin V, Afzal MH, Lachapelle G. Complete triaxis magnetometer 129. Handcock RN, Swain DL, Bishop‑Hurley GJ, Patison KP, Wark T, Valencia calibration in the magnetic domain. J Sens. 2010;2010:967245. P, Corke P, O’Neill CJ. Monitoring Animal behaviour and environmental 152. Tabatabaei SAH, Gluhak A, Tafazolli R. A fast calibration method for interactions using wireless sensor networks, GPS collars and satellite triaxial magnetometers. IEEE Trans Instrum Meas. 2013;62(11):2929–37. remote sensing. Sensors. 2009;9(5):3586–603. 153. Vitali A. Ellipsoid or sphere fitting for sensor calibration, Dt0059. ST 130. Gužvica G, Bošnjak I, Bielen A, Babić D, Radanović‑ Gužvica B, Šver L. Microelectronics, Design Tip; 2016. Comparative analysis of three different methods for monitoring the use 154. Simple and effective magnetometer calibration. https:// github. com/ of green bridges by wildlife. PLoS ONE. 2014;9(8):e106194.krisw iner/ MPU60 50/ wiki/ Simple‑ and‑ Eecff tive‑ Magne tomet er‑ Calib 131. Patel A, Stocks B, Fisher C, Nicolls F, Boje E. Tracking the cheetah ration. tail using animal‑borne cameras, GPS, and an IMU. IEEE Sens Lett. 155. Premerlani W, Bizard P. Direction cosine matrix imu: theory. Diy Drone: 2017;1(4):1–4. Usa; 2009:13–15. 132. Latham ADM, Latham MC, Anderson DP, Cruz J, Herries D, Hebblewhite 156. Grygorenko V. Sensing‑magnetic compass with tilt compensation. M. The GPS craze: six questions to address before deciding to deploy Cypress perform, semiconductor. Application Notes, AN2272. 2011. GPS technology on wildlife. N Z J Ecol. 2015;39(1):143–52. 157. Gleiss AC, Wilson RP, Shepard ELC. Making overall dynamic body 133. Hebblewhite M, Haydon DT. Distinguishing technology from biology: a acceleration work: on the theory of acceleration as a proxy for energy critical review of the use of GPS telemetry data in ecology. Philos Trans expenditure. Methods Ecol Evol. 2011;2(1):23–33. R Soc B Biol Sci. 2010;365(1550):2303–12. 158. Choi S, Youn I‑H, LeMay R, Burns S, Youn J‑H. Biometric gait recognition 134. Gibb R, Shoji A, Fayet AL, Perrins CM, Guilford T, Freeman R. Remotely based on wireless acceleration sensor using k‑nearest neighbor clas‑ sensed wind speed predicts soaring behaviour in a wide‑ranging sification. In: 2014 international conference on computing, networking pelagic seabird. J R Soc Interface. 2017;14(132):20170262. and communications (ICNC), 3–6 Feb 2014; 2014. p. 1091–5. 135. Swain DL, Wark T, Bishop‑Hurley GJ. Using high fix rate GPS data to 159. Sato K, Mitani Y, Cameron MF, Siniff DB, Naito Y. Factors affecting strok ‑ determine the relationships between fix rate, prediction errors and ing patterns and body angle in diving Weddell seals under natural patch selection. Ecol Model. 2008;212(3):273–9. conditions. J Exp Biol. 2003;206(9):1461–70. 136. Ryan PG, Petersen SL, Peters G, Grémillet D. GPS tracking a marine 160. Li W, Wang J. Eec ff tive adaptive Kalman filter for MEMS‑IMU/mag‑ predator: the effects of precision, resolution and sampling rate on netometers integrated attitude and heading reference systems. J Navig. foraging tracks of African Penguins. Mar Biol. 2004;145(2):215–23. 2013;66(1):99–113. 137. Bouvet D, Garcia G. GPS latency identification by Kalman filtering. 161. Pedley M. Tilt sensing using a three‑axis accelerometer. Freescale semi‑ Robotica. 2000;18(5):475–85. conductor. Application Note 2461. 2013;1(Rev. 6):1–22. 138. Frair JL, Fieberg J, Hebblewhite M, Cagnacci F, DeCesare NJ, Pedrotti L. 162. Shiomi K, Narazaki T, Sato K, Shimatani K, Arai N, Ponganis PJ, Resolving issues of imprecise and habitat‑biased locations in ecologi‑ Miyazaki N. Data‑processing artefacts in three ‑ dimensional dive path cal analyses using GPS telemetry data. Philos Trans R Soc B Biol Sci. reconstruction from geomagnetic and acceleration data. Aquat Biol. 2010;365(1550):2187–200. 2010;8(3):299–304. 139. Péron G, Calabrese JM, Duriez O, Fleming CH, García‑ Jiménez R, 163. Evans P. Rotations and rotation matrices. Acta Crystallogr Sect D. Johnston A, Lambertucci SA, Safi K, Shepard ELC. The challenges of 2001;57(10):1355–9. estimating the distribution of flight heights from telemetry or altimetry 164. Wilson RP, Holton MD, Walker JS, Shepard ELC, Scantlebury DM, Wilson data. Anim Biotelem. 2020;8(1):5. VL, Wilson GI, Tysse B, Gravenor M, Ciancio J, McNarry MA, Mackintosh 140. Poulin M‑P, Clermont J, Berteaux D. Extensive daily movement rates KA, Qasem L, Rosell F, Graf PM, Quintana F, Gomez‑Laich A, Sala J‑E, measured in territorial arctic foxes. Ecol Evol. 2021;00:1–12. Mulvenna CC, Marks NJ, Jones MW. A spherical‑plot solution to linking 141. Marcus Rowcliffe J, Carbone C, Kays R, Kranstauber B, Jansen PA. Bias acceleration metrics with animal performance, state, behaviour and in estimating animal travel distance: the effect of sampling frequency. lifestyle. Mov Ecol. 2016;4(1):22. Methods Ecol Evol. 2012;3(4):653–62. 165. Caruso MJ. Applications of magnetic sensors for low cost compass 142. Belant JL. Eec ff ts of antenna orientation and vegetation on global systems. In: IEEE 2000 position location and navigation symposium (Cat positioning system telemetry collar performance. Northeast Nat. No00CH37062): 13–16 March 2000; 2000. p. 177–84. 2009;16(4):577–84, 578. 166. Boelter KDH. Aircraft magnetic declinator system, vol. 2020/0182619. 143. Graf PM, Wilson RP, Qasem L, Hackländer K, Rosell F. The use of accel‑ Chicago: The Boeing Company; 2020. eration to code for animal behaviours; a case study in free‑ranging 167. World magnetic model 2020 calculator. http:// www. geomag. bgs. ac. Eurasian beavers castor fiber. PLoS ONE. 2015;10(8):e0136751.uk/ data_ servi ce/ models_ compa ss/ wmm_ calc. html. Accessed 07 Mar 144. Tonini MH, Palma ED. Tidal dynamics on the North Patagonian Argen‑ 2021. tinean Gulfs. Estuar Coast Shelf Sci. 2017;189:115–30. 168. Qasem L, Cardew A, Wilson A, Griffiths I, Halsey LG, Shepard ELC, Gleiss 145. Wilson RP, Kreye JM, Lucke K, Urquhart H. Antennae on transmitters AC, Wilson R. Tri‑axial dynamic acceleration as a proxy for animal energy on penguins: balancing energy budgets on the high wire. J Exp Biol. expenditure; should we be summing values or calculating the vector? 2004;207(15):2649–62. PLoS ONE. 2012;7(2):e31187. Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 37 of 37 169. Pewsey A, Neuhäuser M, Ruxton GD. Circular statistics in R. OUP Oxford; 172. Movable type scripts. https:// www. movab le‑ type. co. uk/ scrip ts/ latlo ng. 2013. html. Accessed 07 Mar 2021. 170. Farrahi V, Niemelä M, Kangas M, Korpelainen R, Jämsä T. Calibration 173. Grolemund G, Wickham H. Dates and times made easy with lubridate. J and validation of accelerometer‑based activity monitors: a systematic Stat Softw. 2011;40(i03):1–25. review of machine‑learning approaches. Gait Posture. 2019;68:285–99. 171. Ropert‑ Coudert Y, Kato A, Baudat J, Bost C‑A, Le Maho Y, Naito Y. Time/ Publisher’s Note depth usage of Adélie penguins: an approach based on dive angles. Springer Nature remains neutral with regard to jurisdictional claims in pub‑ Polar Biol. 2001;24(6):467–70. lished maps and institutional affiliations. Re Read ady y to to submit y submit your our re researc search h ? Choose BMC and benefit fr ? Choose BMC and benefit from om: : fast, convenient online submission thorough peer review by experienced researchers in your field rapid publication on acceptance support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Animal Biotelemetry Springer Journals

Loading next page...
 
/lp/springer-journals/dead-reckoning-animal-movements-in-r-a-reappraisal-using-gundog-tracks-2dioCzw0El
Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2021
eISSN
2050-3385
DOI
10.1186/s40317-021-00245-z
Publisher site
See Article on Publisher Site

Abstract

Background: Fine‑scale data on animal position are increasingly enabling us to understand the details of animal movement ecology and dead‑reckoning, a technique integrating motion sensor ‑ derived information on heading and speed, can be used to reconstruct fine ‑scale movement paths at sub ‑second resolution, irrespective of the environ‑ ment. On its own however, the dead‑reckoning process is prone to cumulative errors, so that position estimates quickly become uncoupled from true location. Periodic ground‑truthing with aligned location data (e.g., from global positioning technology) can correct for this drift between Verified Positions ( VPs). We present step ‑by‑step instruc‑ tions for implementing Verified Position Correction ( VPC) dead‑reckoning in R using the tilt ‑ compensated compass method, accompanied by the mathematical protocols underlying the code and improvements and extensions of this technique to reduce the trade‑ off between VPC rate and dead‑reckoning accuracy. These protocols are all built into a user‑friendly, fully annotated VPC dead‑reckoning R function; Gundog.Tracks, with multi‑functionality to reconstruct animal movement paths across terrestrial, aquatic, and aerial systems, provided within the Additional file 4 as well as online (GitHub). Results: The Gundog.Tracks function is demonstrated on three contrasting model species (the African lion Panthera leo, the Magellanic penguin Spheniscus magellanicus, and the Imperial cormorant Leucocarbo atriceps) moving on land, in water and in air. We show the effect of uncorrected errors in speed estimations, heading inaccuracies and infrequent VPC rate and demonstrate how these issues can be addressed. Conclusions: The function provided will allow anyone familiar with R to dead‑reckon animal tracks readily and accu‑ rately, as the key complex issues are dealt with by Gundog.Tracks. This will help the community to consider and imple‑ ment a valuable, but often overlooked method of reconstructing high‑resolution animal movement paths across diverse species and systems without requiring a bespoke application. Keywords: Animal behaviour, Animal movement, Global Positioning System, R (programming language), Track integration, Tri‑axial accelerometers, Tri‑axial magnetometers Background Reconstructing animal movement paths is an important *Correspondence: richard.m.g@hotmail.com 1 tool in ecology, providing insights into animal space-use, Swansea Lab for Animal Movement, Department of Biosciences, College of Science, Swansea University, Singleton Park, Swansea SA2 8PP, Wales, behaviour and habitat selection [1–3]. However, accurate UK estimation of paths at fine temporal scales has proved a Full list of author information is available at the end of the article © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Gunner et al. Anim Biotelemetry (2021) 9:23 Page 2 of 37 persistent challenge [4, 5]. Dead-reckoning is a method whilst the utility of transmission telemetry is restricted used to reconstruct animal movement paths, based on to periodic resurfacing events [58]. Moreover, speed can sequentially integrating the vector of travel from a pre- be more easily measured or approximated in water, with determined position using estimates of heading (also previous studies obtaining estimates via acoustic flow termed ‘bearing’ or ‘yaw’) and speed (and displacement noise [e.g., 59], passive sonar [e.g., 60], pitch and change about the vertical axis for 3-D movements), over an in depth [e.g., 11] and speed sensors [cf. 16, 61, 62]. The elapsed time interval [6–9]. In its most advanced form, it efficacy of such techniques diminishes within the aerial can provide positional data with sub-second resolution, environment, principally, due to the marked difference irrespective of the environment [e.g., 10, 11, 12] and it between water and air density [cf. 63] and the current therefore has huge potential for providing data that can speed and volatility of wind [cf. 64, 65]. Indeed, this may elucidate many fundamental behavioural and ecological explain why, in part (to our knowledge), only one study issues related to space-use. to date has dead-reckoned a volant species [66]. More The concept of dead-reckoning (also termed ‘track recently, dynamic body acceleration (DBA, see Wilson integration’) originated to aid nautical navigation [6, 13], et al. [67], for recent review) has been validated as a proxy though its utility to reconstruct uninterrupted fine-scale of speed for terrestrial animals [68, 69] although there are (in time and space) animal movement paths by inte- still very few studies that use the dead-reckoning method grating different sensors in animal-attached tags was in terrestrial animals [e.g., 9, 10, 26, 34, 70, 71]. suggested over three decades ago [14, 15]. Today, this We suggest that another reason that Verified Position typically involves the simultaneous deployment of tri- Correction (VPC) dead-reckoning has been little used axial accelerometers and magnetometers [e.g., 9, 10, 16, relates to the apparent difficulty and poor accessibility of 17–20], utilising the tilt-compensated compass method the analytical processes involved. With this in mind, the [21–24] (see “Glossary” for a definition of dead-reckon - primary aim of this paper is to provide potential users ing-related terminology used throughout). with a clear, concise roadmap for implementing dead- The utility of dead-reckoning depends on the accu - reckoning protocols. Specifically, we revisit the dead- racy of speed and heading estimates (see Table  1) and, reckoning methodology, from calibrating magnetometry due to the nature of vector integration, dead-reckoned data and deriving heading (tilt-compensated compass tracks accumulate errors (commonly termed ‘drift’) over method), to VPC dead-reckoning within both terrestrial time [15, 25, 26]. As a result, periodic ground-truthing and fluid media. We provide simplistic conceptual expla - by a secondary source is important for maintaining the nations and mathematical protocols and describe the accuracy of animal paths with all its underlying scales pitfalls within the procedure that can increase error. We and tortuosity of movement [9, 10, 27]. For this reason, also translate the relevant equations into complementary dead-reckoning data are normally enhanced by combin- R code [cf. 93, available at 94] throughout the text, with ing it with other methods for providing Verified Positions fully annotated scripts deposited in Additional files 2, 3, (VPs). These are primarily; direct observation [e.g., 28], 4, 5 and GitHub [available at 95]. light intensity-based geolocation [e.g., 29], VHF—[e.g., In addition to the above, we outline recent advances 30], acoustic—[e.g., 31] and GPS telemetry [e.g., 26]. to the VPC dead-reckoning technique. For use in ter- Other, less utilised, systems that may also have merit restrial environments, this includes implementing step at sites frequented by the tagged animals, include radio counts as a distance measure, by-passing dynamic body frequency identification (RFID) stations [cf. 32], camera acceleration (DBA) as a proxy of speed, and assessing the traps [cf. 33] and video footage, such as closed-circuit tel- value of ‘reverse dead-reckoning’ (useful when VPs are evision (CCTV) surveillance [e.g., 34]. Although all these concentrated to the latter end of an animal’s trajectory). systems are subject to a number of issues that can make For marine and aerial environments, we demonstrate their positional fixes temporally widely spaced [e.g., 4, 35, the value of integrating tidal-/air current data with dead- 36], inaccurate [e.g., 37, 38, 39] or impossible [e.g., 40, reckoned vectors (hereafter termed ‘current integration’) 41, 42], they can be critical in providing ground-truthed to reduce errors attributed to drift [cf. 46, 92]. Across all positions, even infrequently, with which to reset accumu- three media of travel (land, water and air), we show the lated drift [9, 26]. value of incorporating different speed coefficients accord - Of the above VP options, GPS-corrected dead-reck- ing to behaviour types. In addition, we provide examples oning is the most widely used and there is a marked bias of the various methods by which VP data can be under- towards marine studies [e.g., 10–12, 15–17, 19, 27, 43, sampled (relevant for high-res GPS datasets) prior to 44–56]. This is likely for logistical reasons, with many correcting dead-reckoned tracks and discuss the scales aquatic animals being larger (and thus can carry larger/ at which users should consider VP correction (which more devices) than their terrestrial counterparts [57], depend on the details of species-specific movement and Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 3 of 37 ‑ ‑ ff Table 1 Possible system errors that can affect the utility of animal deadreckoning within the ‘tilt compensated compass’ framework System error Reasons for error Underlying causes References Possible mitigation measures Derived heading Errors in deriving static acceleration (postural) During bouts of high centripetal (turning) accel [65, 72–75] Gyrointegrated data [cf. 76, 77] estimates eration Freefalling behaviour Using Euler angles (angle of rotation about each The orientation of the device with respect to [21, 79–82] Quaternion estimated heading [cf. 83, 84, 85] axis of a given coordinate system) the earth’s frame of reference (cf. Additional file 1: Text S2) can only be defined reasonably at angles less than perpendicular or less than a 180° inversion (dependent on pitch and roll equations used—cf. “VPC deadreckoning procedure in R” section) from their longitudinal and lateral axes of ‘normal’ posture (otherwise, unstable measures arise from the Gimbal lock singularity complex [cf. 78], whilst x, y and z values can become inversed and/or represent different ‘surge’, ‘sway ’ and ‘heave’ planes) Tag placement/dislodgment In line with the above—range for accurate angular [19, 86] Ensure tag orientation is noted during deployment (pitch and roll) measures are restricted in one or and retrieval operations (and subsequently used more dimensions in corrections) Heading will be biased according to the degree of displacement about the zaxis Variations in the strength and declination of Animals that undertake long journeys (regionally/ [87] Ensure at least one magnetic calibration procedure magnetic fields globally) is carried out (see Additional file 1: Text S3 for Environmental and manmade magnetic noise details) and apply magnetic declination offset to (iron distortions) heading values where required Derived speed Deviations of the DBA ~ speed relationship Load bearing [12, 65, 68, 69, 88–91] Iteratively modulate the gradient and/or intercept Moving over a deformable substrate/changeable within the DBA–speed linear regression according incline to environmental circumstance and mode of Changing gait movement [cf. 68] Moving within fluid media Bypass DBA (e.g., use speed/acoustic sensors, step/ Gliding/thermalling behaviour tail/wing beat frequency, vertical speed, etc. [e.g., 47]) Both External forces (e.g., current vectors in air and Decreases the signalto noise ratio of motion [11, 15, 46, 92] Smooth postural (and pre derivative data)/DBA water flow) sensor data estimates [cf. 67] Aects the relationship between an animal’s Incorporate current flow vectors within the dead (longitudinal axis) direction of travel from their reckoning procedure true vector of travel Some animals do not always move in the same direction as their anterior–posterior axis Gunner et al. Anim Biotelemetry (2021) 9:23 Page 4 of 37 User functionality length of data acquisition). We specifically demonstrate Gundog.Tracks is an all-encompassing dead-reckoning the above using our R-functions (Gundog.Tracks being function that can be used to dead-reckon animal paths the primary function for dead-reckoning), providing travelling terrestrially or through fluid media. Table  2 examples of its utility across various scenarios. Lastly, details all the function’s input requirements/options. we highlight the relevance of heading and distance cor- rection factors (derived from the VPC procedure), which Reverse dead‑reckoning can also be used to interrogate the animal–environ- Dead-reckoning backwards is useful when the start ment interaction and biases stemming from animal tag position is unknown, but the finishing coordinates are performance. known. For example, central-place foraging, diving ani- To illustrate our approach, we use three model species mals returning to land from the sea may not acquire a (the African lion Panthera leo, the Magellanic penguin satellite fix for an appreciable period of time following Spheniscus magellanicus and the Imperial cormorant submersion in water which can make determining the Leucocarbo atriceps) that cover almost two orders of size start position difficult. So, when VPs are skewed to the magnitude in body mass and that operate in markedly dif- latter part of the track, it may be beneficial to start the ferent environments and at different scales of movement. iterative dead-reckoning process from that end. This To make this review more broadly applicable to research- involves reversing the order of data to be dead-reckoned ers of varying dead-reckoning and R knowledge, we have and changing heading values by 180 degrees prior to departed from a traditional article format and instead, dead-reckoning. split this body of work into two distinct sections: firstly, we provide an overview of the critical Gundog.Tracks Integrating current vectors function and provide a brief review of the conceptual Wind or ocean currents can change the relationship workflow (“Implementation of Gundog.Tracks” section). between an animal’s (longitudinal axis) bearing and With respect to this, we discuss the relevant strengths speed of travel from their true vector of travel [46, 92]. and limitations of the current dead-reckoning framework This drift can be incorporated within movement paths by and the key considerations involved. Secondly, we detail advancing each iterated dead-reckoned vector according each ‘potential’ stage of the VPC dead-reckoning proce- to the direction and speed of the current at that point in dure with exemplar mathematical equations and R syntax space and time (cf. Fig. 1). (“VPC dead-reckoning procedure in R” section). DBA–speed derivation Implementation of Gundog.Tracks Given the approximate linear relationship between We expand on key concepts in Additional file  1 and DBA [sensu 67] and terrestrial animal speed provide complimentary R scripts (outlined below) in [speed = (DBA·m) + c], DBA estimates can be multi- Additional files 2, 3, 4 and 5. We also supply an exam- plied by a gradient, m (the multiplicative coefficient) ple data set of a Magellanic penguin walking out to sea and summed with an intercept, c (the constant) to in Additional file  6, which can be used to trial each of derive speed [10, 26]. These values are typically substi - the provided R scripts and perform the full dead-reck- tuted with results from DBA–speed linear regression oning process. Mathematical equations are referred to estimates, such as from treadmill tests or using GPS- as ‘Eqs.  1–33’ and R syntax as ‘R ’, where ‘ ’ is the refer- x x derived speed [26, 69, 97, 98]. The m -coefficient should ence number. To simplify concepts, we use base R syntax be selected such that (uncorrected) dead-reckoned (wherever possible) and typically use vectors to dem- tracks accord with the apparent straight-line distance onstrate points made, though ‘df$’ directly before the between VPs. Importantly, the DBA–speed relation- variable name indexes data retained within data frame ship may be a function of terrain-type (e.g., sand vs. columns (assuming data frame is called ‘df’). We note, concrete), animal state (e.g., weight variation) and however, that more efficient code implementations are mode of movement (e.g., running vs. climbing) [cf. 68]. possible (e.g., data.table [96] and lapply()) than pre- For instance, a condor gliding within a thermal would sented here, especially for large data, but here wanted have high speeds, despite having negligible DBA, while to make the code as readable as possible in this manu- an Ibex traversing across different substrate types and script, especially to persons not familiar with complex gradients would impart varying magnitudes of accel- coding. More efficient code will be implemented through eration that may scale non-linearly with a change in updated GitHub versions of the functions. See Additional stride gait. It may be of value, therefore, to iteratively file  1: Text S1 for our model species’ device setup and change the supplied m (and possibly c) values between capture protocol and the glossary for a definition of dead- VPs according to behaviour and environment. The user reckoning related terminology. Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 5 of 37 Table 2 Gundog.Tracks input fields and description of their role Funcon input Descripon Ref TS Timestamp - POSIXct object. No missing data (NA’s) permied - o o h Heading (0 to 360 ) - No missing data (NA’s) permied v - DBA (g or m/s ) or speed (m/s) - No missing data (NA’s) permied Elevaƒon / depth data (m) - No missing data (NA’s) permied 0 elv Pitch ( ) – Only supply if user wants radial distance modulated according to pitch (cf. "VPC dead-reckoning p NULL ). No Current speed (m/s) - Supplied as a single value or vector/column of changeable values. NA’s are replaced with cs NULL the most recent non-NA prior to it (observaƒons carried forward) o o Current heading (0 to 360 ) - Supplied as a single value or vector/column of changeable values. NA’s are replaced ch NULL with the most recent non-NA prior to it (observaƒons carried forward) Mulƒplicaƒve coefficient (gradient) - If speed (m/s) supplied for v, then m must be 1. m 1 Supplied as a single value or vector/column of changeable values Constant (y-intercept) – If speed supplied for v, then c must be 0. c 0 Supplied as a single value or vector/column of changeable values Marked Events – 0 denotes periods of sta�onary behaviour and 1 (or any integer number > 0) denotes periods of ME traversing movement. ME overrides ini�al speed input / DBA-derived speed (calculated within the func�on itself). 1 NA’s and character values are replaced with zero. lo Star�ng longitude coordinate to advance dead-reckon track from – Decimal format, e.g., 26.31989 0 Star�ng la�tude coordinate to advance dead-reckon track from – Decimal format e.g., -06.11995 0 la VP longitude coordinates – Decimal format. Missing reloca�on data expressed as either NA’s or 0’s. First (or last VP.lon NULL if reverse dead-reckoning) element/row allocated as lo within the func�on VP la�tude coordinates – Decimal format. Missing reloca�on data expressed as either NA’s or 0’s. First (or last if VP.lat NULL reverse dead-reckoning) element/row allocated as la within the func�on TRUE = Supplied VPs removed at �mes when ME = 0 (relevant for high-res VP datasets, when loca�on error is high VP.ME FALSE during rest). Note, this does not remove the element/row allocated as lo/la How the func�on under-samples VPs prior to correc�on (subsequent to the VP.ME subset, if set to TRUE) – “divide” = Fix kept every x (thresh) segments of supplied VPs, based on row number. The first and last fixes are always included “ me” (s) = F ix kept every x (thresh) accumulated seconds ( or t he n ext available fix a œer a period o f missing method NULL loca�onal data (≥ thresh). The first and last fixes are always included “distance” = Fix kept every x (thresh) p ropor�onal s egments of t he t otal a ccumulated distance ( m) b etween supplied VPs (using the stepping interval ‘dist.step’). The first and last fixes are always included “all” = Every supplied VP kept (irrespec�ve of thresh value) Threshold - Degree of V P under-sampling prior to d ead-reckon c orrec�on. The f requency of under sampling thresh 1 depends on the method selected The stepping interval used for calcula�ng distance b etween V Ps, both w ithin the VP s ummary d istance m etrics (see Table. 3) and within the ‘method = distance' VP under-sampling protocol prior to VPC. For example, dist.step dist.step 1 = 5 computes distance between every 5th VP (irrespec�ve of the �me difference between them) TRUE = VPC dead-reckoning is bounded by the first and last VP present bound FALSE = VPC dead-reckoning is unbounded by the last available VP. The last dead-reckoned track segment inherits TRUE the previous correc�on factors TRUE = ‘normal’ dead-reckoning procedure Outgoing FALSE = R everse dead-reckoning. N ote la a nd lo p osi�ons should n ow b e the finishing l ongitude a nd la�tude TRUE coordinates, respec�vely FALSE = No summary plots TRUE = R graphics window ini�alized: VPC = 4 summary plots / no VPC = 1 summary plot (cf. top leœ ) Top le) Uncorrected dead-reckoned track (blue) and VP track (red). If currents are supplied, the blue track has currents integrated and an addi�onal green track with no current integra�on is plo¤e d Plot Top right) VPC dead-reckoned track (blue) in rela�on to VP track (red) TRUE Bo�om le) Net error (m) between VPs and dead-reckoned posi�ons (un-corrected = red and corrected = black). If currents are supplied, then uncorrected w ith no c urrent i ntegra�on = green a nd uncorrected w ith current integra�on = red. Bo�om right) VPC corrected dead-reckoned track Ref. refers to the default value when no input is stated. Red shading represents required user inputs and green and orange shading reflect optional inputs (the latter change when using VPC). Note that if speed estimates (v) are directly inputted into the function then m and c (and possibly ME) defaults should not be changed. If either one of the VP.lon, VP.lat or method inputs is specified as NULL, then no VPC will occur Gunner et al. Anim Biotelemetry (2021) 9:23 Page 6 of 37 DBA is a weak proxy of speed for many marine animals because overall body tissue density changes with depth when air is associated so that speed may be invariant of the movement kinematics [cf. 103, 104]. DBA is also a weak proxy for flying animals that glide at constant veloc - ity, use thermals or bank [cf. 105]. One of the most com- mon methods for determining animal speed in water is via devices that estimate flow or resistance rate [16, 19, 64, 106]. These often have appreciable limitations, with currents, biofouling, blockage and turbulence affecting performance [64], and many of these issues are appli- cable to volant species, so that bird speed measures are typically restricted to GPS-derived estimates of ground Fig. 1 Schematic representation of a current flow vector (orange) speed [cf. 107]. In the absence of a reliable motion sen- (due to its speed and direction) being integrated to a given travel sor-derived speed proxy, previous reported approximated vector (blue). The x, y reflect the initial location of a dead‑reckoned speed estimates according to movement modes and/or track, x2 and y2 are the resultant location following the integration of topological whereabouts can be used [cf. 30]. For exam- a travel vector (prior to current integration) and xxx and yyy advance ple, for various diving animals such as penguins, a sim- these x 2 and y2 values a step further in the direction of the current flow vector. The dashed lines indicate the magnitude of the x and y ple depth threshold may prove effective to differentiate dimensions of travel (both pre‑ and post ‑ current flow integration) between various previous reported modal ‘surface-rest- and the green line reflects the actual travel vector ing’ and ‘underwater-commuting’ speeds [61, 108–110]. For volant species, whilst wingbeat frequency or ampli- tude does not scale reliably with air speed [cf. 111], the may also opt to supply a ‘marked events’ (ME) vector (a interplay between both can decipher various flight modes marked event is a term we use that refers to a number of (e.g., ‘cruising speed’ vs taking off/landing) [65, 112]. sequential (in time) data points within a dataset coded Furthermore, tail beat frequency has been shown to be by integer values) to ensure dead-reckoned tracks are a good predictor of swimming speed for various fish spe - not advanced with non-animal movement behaviours. cies [113–115]. For diving animals, a proxy for horizontal Within Gundog.Tracks, ME values of one or greater speed can be obtained based on animal pitch and rate of reflect progressive movement, and zero values code for change of depth [11, 116]. Specifically, the rate change of stationary behaviour—dead-reckoned tracks are not depth is divided by the tangent of the body pitch. advanced when ME = 0 (irrespective of the allocated In any case, when high-resolution VP data are available speed). For example, in its simplest form, ME could (e.g., 0.01–10  Hz GPS), for instance, during short-term be filled with binary 0  s and 1  s as governed by a DBA trial deployments, speed estimates can be compared threshold (labelling the ME vector 0 in sleep and resting alongside those derived between VPs and approximated behaviour). according to behaviour type (elucidated from, for exam- ple, accelerometer—[e.g., 117, 118], magnetometry— Pre‑determining speed [e.g., 105, 119], depth- [e.g., 120] or altitude—[e.g., 65] For terrestrial species (specifically bipeds and quadru - data), and uncorrected dead-reckoned tracks can be peds), the interplay between peak heave acceleration compared alongside VPs to determine where biases amplitude and periodicity may be a useful indicator for may occur visually. Furthermore, the correction factors the movement gait adopted [99], which may help decide obtained from the VPC process are viable comparators the m-coefficient in the DBA–speed relationship [69]. for detecting consistent under- or over-estimations of There may be times, however, when DBA is an unreliable speed and/or heading offsets (e.g., due to tag placement). proxy of terrestrial speed [cf. 68]. At this time, given that Essentially, when empirical speed evidence is unavail- the stride cycle can be easily detected by cyclic peaks in able, the user can ad hoc iteratively adjust allocated speed a given acceleration channel [e.g., 100, 101, 102], peak values or the underlying DBA–speed coefficients until periodicity (and amplitude) may be used as a proxy of uncorrected dead-reckoned track segments scale propor- distance moved by providing a distance per step estimate tionately to their aligned ground-truthed positions (pre- (assuming constant distance travelled between step gaits VPC). Within Gundog.Tracks, the user can modulate m, if only concerning step periodicity—cf. “VPC dead-reck- c and ME values to switch between predetermined speed oning procedure in R” section and Additional file  1: Text (m = 1, c = 0), DBA-derived speed (m > 0, c ≥ 0) and sta- S4). tionary behaviour (ME = 0). Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 7 of 37 VPC procedure with input modulated according to the animal in ques- Ground-truthing dead-reckoned tracks typically involves tion and data available (see Fig. 2). the linear drift correction method [cf. 26, 46], outlined The function outputs a data frame containing vari - in Constandache et al. [121] and Symington and Trigoni ous descriptive columns which, depending on the input, [122]. In essence, a shift vector aligns the starting dead- includes (but is not limited to): reckoned path segment with the VP at time point one, after which the difference between the VP and dead-reck - • The correction factors used oned path segment at time point two is calculated to pro- • Heading and radial distance estimates (both pre- and vide a correction vector that is applied linearly between post-current integration and/or VPC) time point one and time point two. Our method follows • Distance moved and speed estimates (both in 2-D the protocols outlined by Walker et  al. [9], whereby the and 3-D when elevation/depth data supplied) underlying correction coefficients (hereafter termed • Net error between dead-reckoned positions and VPs ‘factors’) for both heading and (radial) distance are cal- (both pre- and post-correction) culated—adjusting the length and heading at each dead- • Various VP summaries including notation of when reckoned path segment until the end points align to each VPs are present and which fixes were used in the cor - VP along the path. This process requires the trigonomet - rection process. ric ‘as the crow flies’ Haversine formulae [123–125] which allows one to translate a distance across the curvature of When specified, 2-D summary plots demonstrating the the Earth’s surface (detailed within “VPC dead-reckoning relationship between dead-reckoned positions and VPs procedure in R” section). The advantage of this method (both pre- and post-current integration and/or VPC) are is that, whilst correction factors are constant between provided (e.g., Fig.  3). Table  3 details all the function’s VPs, it does not assume that the dead-reckoned path available outputs (modulated according to input). Gun- deviates linearly over time from the true path because dog.Tracks uses the na.locf() function from the ‘zoo’ pack- (radial) distance is multiplied by the distance correction age [126] and the slice() function from the ‘dplyr’ package factor. This ensures that parts of track where the animal [127] (both are checked as dependencies and installed is determined to be stationary (e.g., ME = 0) are left unal- when required within this function). Output 2-D dis- tered. The function’s method of VPC, automatically han - tance/speed estimates are calculated with the Haversine dles NaN and Infinite (Inf ) values which can arise during formula. When depth/elevation data are supplied (and the derivation of the distance correction factors (when changes between sets of coordinates) 3-D distance/speed no dead-reckoned movement occurs between successive estimates are calculated with a variant of the Euclidean VPs—detailed within “VPC dead-reckoning procedure Formula—converting x, y, z from polar to Cartesian coor- in R” section). It is worth noting that even animals that dinates, and incorporating the Earth’s oblate spheroid [cf. travel in 3-D can be subject to the 2-D dead-reckoning World Geodetic System (WGS84)], via conversion from formulae and Haversine computation of distance correc- Geodetic- to Geocentric-latitude [cf. 128]. tion factors because we typically assume that both dead- The interplay between numerical precision in R, cor - reckoned- and VP positions are aligned in vertical space rection rate and net error can make more than one (assuming reliable pressure—[60]/depth [13] data) and round of adjustment necessary for dead-reckoning fixes attempt to control for the horizontal component of speed to accord exactly with ground-truthed locations (cf. [e.g., “VPC dead-reckoning procedure in R” section— Fig. 4a), particularly given that slight discrepancies accu- Eqs. (25, 27)] pre-correction. Although not covered here, mulate over time. Each iteration of the correction pro- we acknowledge that various state–space modelling tech- cess produces a tighter adherence between estimated and niques have also been developed to georeference dead- ground-truthed positions [cf. 9]. Typically, this does not reckoned tracks [e.g., 11, 47]. involve more than two rounds of VPC to achieve a maxi- mum net error of 0.01 m (the threshold used within Gun- Default inputs for calculations and outputs dog.Tracks) across a ca. (1 Hz) 2-week-long track. For an Gundog.Tracks default input takes the form: indication of processing times see Additional file  1: Text S6, Fig. S4; for example dead-reckoning a lion at 1  Hz Gundog.Tracks(TS, h, v, elv = 0, p = NULL, for 7 (continuous) days (with plot = TRUE, dist.step = 5, cs = NULL, ch = NULL, m = 1, c = 0, ME = 1, VP.ME = TRUE, method = “time” and thresh = 3600) lo = 0, la = 0, VP.lon = NULL, VP.lat = NULL, took 25 s to compute (on a MSI GP72 7RD Leopard lap- VP.ME = FALSE, method = NULL, thresh = 1, top with intel core i7 processor). Logically, the net error dist.step = 1, bound = TRUE, Outgo- between VPs and (corrected) dead-reckoned positions is ing = TRUE, plot = FALSE), positively correlated to the time between corrections (cf. Gunner et al. Anim Biotelemetry (2021) 9:23 Page 8 of 37 Fig. 2 Schematic of the conceptual workflow involved when dead‑reckoning using Gundog.Tracks—elaborated within “ VPC dead‑reckoning procedure in R” section. Note Gundog.Peaks (Additional file 3) is a peak finder function that locates peaks based on local signal maxima and Gundog.Compass (Additional file 2) is a function to correct iron distortions from tri‑axial magnetometry data and subsequently compute tilt‑ compensated heading. Both functions are elaborated within “VPC dead‑reckoning procedure in R ” section and in Additional file 1. The direction of workflow and key questions asked follow from green—(pre ‑processing and data alignment) to purple—(computing heading) sections, before splitting into blue—(air/water) and brown—(land) sections (computing speed) and culminating at the red section (final pre ‑ dead reckoning checks/data formats and post‑ dead‑reckoning checks/plots) in conjunction with the process of using Gundog.Tracks in R (yellow) Fig.  4b) [cf. 46], although the rate of net error ‘drop-off ’ according to the permissiveness of the environment, such is dependent on the accuracy of the initial (uncorrected) as high-density shrub or submersion in water [e.g., 38, dead-reckoned track (cf. Fig.  5), itself, modulated by the 129, 130]. GPS technology is arguably the most popular extent of system errors (Table  1) and initial user-defined and widely used method for determining estimates of track scaling. free-ranging animal movement [cf. 131, 132, 133]. This Within this process, people assume VPs to be perfect, is because inspection of data is less complex and time- however, across all VP determining methods, the rate consuming than some of the alternatives, whilst improve- and accuracy of data acquisition is highly moderated ments in design and battery longevity have enabled Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 9 of 37 Fig. 3 Dead‑reckoned (DR) movement path of lion as provided by Gundog.Tracks summary plots (within the initialised R graphics window). This is an approximate 2‑ week trajectory over an approximated total travel (DR) distance of > 142 km. (Pre‑filtered) GPS (red) was sampled at 1 Hz and derived heading and speed measurements were sub‑sampled to 1 Hz (initial acceleration/magnetometry data were recorded at 40 Hz). The VPC dead‑reckoned track (blue) was constructed using DBA–GPS‑ derived speed regression estimates and corrected approx. every 6 h. Note, for dead‑reckoning within fluid media, an additional green dead‑reckoned track with current integration and its associated distance estimates are also plotted (pre‑ correction) when wind/ocean currents are supplied (cf. Fig. 6). Accumulated 3‑D DR distance is shown when elevation/depth data are supplied GPS units to be attached to a plethora of animals (up to appreciably more depending upon the propagation of sig- almost four orders of magnitude in size and mass [cf. 11, nal quality and/or receiver reception capability [38, 138, 134]) and record at high frequencies (e.g., ≥ 1  Hz [131, 139]. As such, VP error becomes more relevant at smaller 135]). Consequently, GPS units are unparalleled for pro- scales of assessed movement and this is the reason why viding such detailed quantification of space-use out - VP distance-moved estimates can go from being typically side of the VPC dead-reckoning framework, and are the underestimated at low frequencies (due to linear interpo- most utilised VPC method within (including the case lation of tortuous movements) [26, 140, 141] to overesti- study datasets within this study). However, locational mated at high frequencies [97, 136] and result in highly accuracy (excepting precision error radius [cf. 136] and variable correction factors within the VPC dead-reckon- variable latency [cf. 137]) can vary by a few metres or be ing process [cf. 10]. Indeed, judicious selection of VPC Gunner et al. Anim Biotelemetry (2021) 9:23 Page 10 of 37 Table 3 Gundog.Tracks data frame output names and their parameters Funcon output Descripon Ref Row.number Row number Supplied mestamp - POSIXct object Timestamp DR.seconds Accumulated me (s) based on the supplied mestamp o o Heading Supplied heading (0 to 360 ) Marked.events Supplied Marked events (or replicated default) DBA.or.speed Supplied DBA (g or m/s ) or speed (m/s) Pitch Supplied pitch ( ) The calculated q-coefficient (prior to VPC) (see “VPC dead-reckoning procedure in R” secon—Eq. (26)) # Radial.distance Elevaon Supplied elevaon / depth (m) Elevaon.diff Rate change of supplied elevaon/depth (m/s) - (elevaon difference / me difference between rows) Supplied current speed (m/s) Current.strength o o Current.heading Supplied current heading (0 to 360 ) o o Heading.current.integrated # Updated heading (0 to 360 ) following addion of current vectors (prior to VPC) Radial.distance.current.integrated Updated q-coefficient following addion of current vectors (prior to VPC) # DR.longitude Dead-reckoned longitude coordinates – Decimal format (prior to VPC) DR.la tude Dead-reckoned latude coordinates – Decimal format (prior to VPC) Corrected dead-reckoned longitude coordinates – Decimal format (post VPC) DR.longitude.corr DR.la tude.corr Corrected dead-reckoned latude coordinates – Decimal format (post VPC) Dist.corr.factor Distance correcon factor (observaons carried forward) ↑ # o o Head.corr.factor ↑ # Heading correcon factor (0 to 360 ) (observaons carried forward) o o Heading.corr # Corrected heading (0 to 360 ) (post VPC) Corrected q-coefficient (post VPC) # Radial.distance.corr Distance (m) between uncorrected dead-reckoned posions and VPs (observaons carried forward), Distance.error.before.correc on ↑ subsequent to sub-sampling according to ME, if VP.ME = TRUE Distance (m) between corrected dead-reckoned posions and VPs (observaons carried forward), Distance.error.a er.correc on ↑ subsequent to sub-sampling according to ME, if VP.ME = TRUE DR.distance.2D Two-dimensional distance moved (m) between dead-reckoned fixes * DR.distance.3D Three-dimensional distance moved (m) between dead-reckoned fixes * DR.cumula ve.distance.2D Accumulated two-dimensional distance moved (m) between dead-reckoned fixes * DR.cumula ve.distance.3D Accumulated three-dimensional distance moved (m) between dead-reckoned fixes * Two-dimensional (straight-line) distance moved (m) from starng posion DR.distance.from.start.2D * Three-dimensional (straight-line) distance moved (m) from the starng posion * DR.distance.from.start.3D DR.speed.2D Horizontal speed (m/s) (DR.distance.2D / me difference between rows) * DR.speed.3D Total speed (m/s) (DR.distance.3D / me difference between rows) * Accumulated me (s) between supplied VPs (observaons carried forward) VP.seconds VP.longitude Supplied VP longitude values (observaons carried forward), sub-sampled according to ME, if VP.ME = TRUE ↑ VP.la tude Supplied VP latude values (observaons carried forward), sub-sampled according to ME, if VP.ME = TRUE ↑ Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 11 of 37 Table 3 (continued) Denotes when a fix was present (1) or absent (0), subsequent to sub-sampling according to ME, if VP.ME = VP.fix.present TRUE Denotes which VPs were used to correct (1) and which VPs were ignored (0) VP.used.to.correct Number.of.VPCs Increments by 1 each ƒme a VP was used to correct (observaƒons carried forward) ↑ Replicates the thresh value set (or default) or warns the user that addiƒonal VP under-sampling was VP.thresh required if ‘Inf’ values produced Two-dimensional distance moved (m) between VPs, subsequent to sub-sampling according to ME, if VP.ME VP.distance.2D = TRUE and using the stepping interval ‘dist.step’ Accumulated two-dimensional distance moved (m) between VPs, subsequent to sub-sampling according to VP.cumulave.distance.2D ME, if VP.ME = TRUE and using the stepping interval ‘dist.step’ The shading of Ref refers to when the outputs occur; red shading = always, purple shading = when pitch data are supplied, green shading = when elevation/depth data are supplied, blue = when current data are supplied and orange = when the user opts to undertake VPC. The symbol * demonstrates that the metrics will be derived from VPC tracks when correction is initialised. Note that, subsequent to reverse dead-reckoning, the data frame is reverted (back to original time order), though as a result, observations appear to be carried backwards in some instances, indicated by ↑. Due to the nature of reverse dead-reckoning (cf. Additional file 4), some input fields are shifted forward one row following the initial inversion of data. As such, fields indicated by #, are one row further forward in time (this is important when relating Head.corr.factor and Heading.corr output to the equivalent (uncorrected) Heading output. However, when currents are integrated, the Head. corr.factor and Heading.corr outputs refer to Heading.current.integrated and these are synchronised row-wise. All heading related data are rotated back 180 degrees following reverse dead-reckoning Fig. 4 Net error between (GPS‑ corrected) dead‑reckoned and GPS positions for a track from 5 African lions. a Maximum net error (m) between ground‑truthed GPS and time ‑matched dead‑reckoned positions after one iteration of correction, both as a function of GPS correction rate [one correction per 1—(red), 12—(green) and 24—(blue) hours] and underlying m‑ coefficient used to determine the DBA‑ derived speed. Data from 5 lions (individual denoted by symbol shape) over a period of 12 days. Note that the difference in error varies according to individual, initial speed estimate and the rate of correction. b Net error between dead‑reckoned positions and all available GPS fixes, according to correction rate (data from the same 5 lions), subsequent to the iterative procedure of GPS correction (maximum distance between GPS fix used in correction procedure and according dead‑reckoned position < 0.01 m). Boxes denote the median and 25–75% interquartile range with a blue ‘loess’ smooth line. Whiskers extend to 1.5 * interquartile range in both directions rate is critical in maximising dead-reckoned track accu- It was suggested by Bidder et  al. [10], that the next racy when relocation data are taken at fine spatial- and stage in this work is to derive a standardised set of rules temporal resolutions [26] (cf. Table  2—‘VP.ME’, ‘method’, to maximise the value of both GPS (though this applies to ‘thresh’ and ‘dist.step’ inputs to aid in modulating VPC any VP method) and dead-reckoned data in line with the rate). Likewise, the initial screening for location anoma- questions being asked. We argue that consistent trends lies, across all VP methods and sampling intervals, is in the magnitude and/or bias of correction factors can be important so as to prevent incorrect distortion of tracks. used as a diagnostic tool for elucidating: (i) VP inaccu- Put simply, the higher the quality of VP data input, the racy (e.g., possibly manifested by extremely high distance greater the robustness of the VPC dead-reckoning and heading correction factors), (ii) required alterations output. to the DBA–speed relationship [e.g., due to traversing Gunner et al. Anim Biotelemetry (2021) 9:23 Page 12 of 37 across different substrates (e.g., Fig.  5)] and (iii) drift due incorporated to some degree. Whilst not illustrated here, to current vectors [cf. 16, 46] (e.g., Fig. 6). advancing tracks without a VeDBA threshold would incur greater error still. Lastly, in this section, the dis- The case‑studies tance correction factor was consistently high (Fig. 5, inset An important question to address is how often to do VP b ) as the lion travelled along the Botswana fence bound- correction. This is obviously dependent upon the scales of ary, perhaps as a result of the animal walking on the com- movement elicited and the medium in/on which the ani- pact dirt road at this location (Fig. 5, inset a ), altering the mal in question navigates. Put simply, one should VP cor- VeDBA–speed relationship. Such patterns in correction rect as little as possible, but as much as is necessary and factors (whether consistent or highly variable) can high- we elaborate on this using our model species operating in light issues with the underlying track scaling. different media. Within Fig.  5, the 1 Hz GPS track (blue) Where animals move in water or air, obtaining accu- is plotted alongside two different dead-reckoned tracks; rate estimates of speed is more difficult without the use [(a) = uncorrected and (b) = corrected approx. every of speed sensors. Naturally, the resolution and accuracy 30  min (method = “time”)] from 12  days of data acquisi- of initial dead-reckoning track scaling (pre-VPC) reduces tion of one lion. There were two variations in the method when speed has to be approximated using constant values of scaling the dead-reckoned tracks; a track based on a according to behaviour type (a strategy used here). There Vectorial Dynamic Body Acceleration (VeDBA) threshold is a balance between initial dead-reckoning accuracy (red), and a track advanced based on periods of identified and required VPC. The lower the initial track accuracy, movement (purple). The m-coefficient and c-constant the more frequent it should be corrected, and additional values were determined from the VeDBA–GPS speed drift caused by external-force vectors compounds this relationship (Fig.  5, inset a ) and the Movement Verified issue. Within Fig.  6, we illustrate the value that current Filtering (MVF) protocol outlined by Gunner et  al. [97] correction, dependent on current information, brings to was used to depict movement and anomalous GPS fixes the VPC procedure if the derived track is to be superim- (green) and to compute reasonable GPS-derived speed posed on the environment. Here, one Magellanic pen- estimates. This case study demonstrates three impor - guin was dead-reckoned with and without tidal vector tant points. Firstly, on its own, dead-reckoning is subject integration (instantaneous tidal currents were deduced to substantial drift and so VPC is essential for resetting from a 3-D numerical model validated in the region this error. The more frequent a user corrects, the more [144], at hourly, 1  km grid nodes). Commuting speed accurate the dead-reckon track becomes (relative to VPs), was allocated 2.1 m/s [cf. 61, 145] and changed according though VP error can also be substantial, especially dur- to “VPC dead-reckoning procedure in R” section—R . ing rest behaviour (see Gunner et al. [97] for demonstra- Surface period ‘rest’ speed were allocated 0.416  m/s [cf. tion of this). For collared animals, heading measurements 108]. VP accuracy improved considerably both pre- and can become inaccurate at times of erratic collar roll (cf. post- VPC when currents were integrated which points Table 1) and conjointly, GPS performance is also reduced to the value of acquiring current data if possible, particu- when antenna position becomes compromised [e.g., 142]. larly if VPs are sparse. Notably the combination of dead- Secondly, and in conjunction to the above, irrespective reckoning and VP estimation of both movements relative of VPC rate, the initial allocation of speed is important. to the ground and fluid, may detail specific orientation Here, only dead-reckoning identified movement periods strategies used and thus can have value for assessing the resulted in greater accuracy than just advancing tracks ability of drift compensation in aquatic or volant animals based on a VeDBA threshold. This is because even sta - [46, 92] tionary behaviours can impart appreciable DBA [e.g., For all our case study animals, GPS units were set to 143] (beyond the threshold), and thus wrongly advance record at 1  Hz. With this temporal resolution (which tracks. The false patterns of tortuosity created from is not always possible anyway due to the high-power this, whilst scaled and possibly rotated with VPC (cf. requirements of the GPS), the value of dead-reckoning “VPC dead-reckoning procedure in R” section), remain would seem questionable. However, dead-reckoning can: (See figure on next page.) Fig. 5 Dead‑reckoned lion track in relation to GPS positions [a uncorrected and b corrected—approx. every 30 min (black circles)]. The start of the track (lo and la) is denoted with a black x. Three corresponding sections of each track are denoted with the same number and the finishing positions denoted with a circle (coloured according to its reference track). Note that the horizontal straight‑line sections (cf. yellow arrow) result from the lion following the Botswana boundary fence (which this individual eventually crossed). Mean net error between (corrected) dead‑reckoned positions and all available GPS fixes was higher for tracks resolved using a VeDBA threshold (0.11 g), than for tracks advanced only at times of depicted movement using the MVF method Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 13 of 37 Gunner et al. Anim Biotelemetry (2021) 9:23 Page 14 of 37 Fig. 6 One Magellanic penguin’s dead‑reckoned foraging trip at sea, lasting approximately 9 h (yellow arrow denotes the trajectory direction over time. Black track = GPS. Fifteen corrections (black circles) were made (method = “divide”). For comparison, the grey dotted track is the GPS‑ corrected dead‑reckoned track with current integration approx. every 1 min (where possible ‑method = “time”) (a). Note the difference of net error between dead‑reckoned positions and all available GPS fixes across the various tracks [insert = grey track] (b). Both uncorrected and corrected dead‑reckoned tracks had less error after current integration (black arrows vector every 5 min) and this was reflected in the direction and magnitude of heading correction factors required per unit time (c). Heading correction factors obtained from the track corrected approx. every 1 min; the colour of the scale bar indicates the extent of the heading correction factor required) Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 15 of 37 (i) work when GPS cannot—such as when an animal is error (relative to VPs). As such, and akin with VP under- underwater [e.g., 18] or in thick forest [cf. 2] and it can sampling, choice of under-sampling data to be dead-reck- (ii) by-pass the issues arising from GPS inaccuracies such oned may have implications to the resultant accuracy, as ‘jitter’ [cf. 97], allowing for more accurate and finer and this will be moderated according to the scales (and scale delineations of movement. This is illustrated in media) of movement elicited by the animal in question. Fig. 7, in which 12 outgoing (green) and incoming (blue) Beyond this, Fig.  7 also demonstrates the importance dead-reckoned trajectories from Magellanic penguins of initial track advancement, with three variants used, walking to and from their nest are plotted. Incoming including step counts instead of DBA. tracks were reverse-dead-reckoned (Outgoing = FALSE, Finally, obtaining accurate estimates of altitude or bound = FALSE), because the GPS did not always reg- depth allow users to plot and investigate scales of contin- ister fixes for minutes after birds left the water and uous movement in three dimensions and at times when because nest coordinates were known (Fig. 7, inset). This VP success rate fails completely (such as underwater). We explains why the blue tracks extend into the sea rather demonstrate this using the Imperial cormorant in Fig.  8. than encroach further inland when speed was overesti- After visual inspection of data, uncorrected tracks were mated. What is evident is that even ‘accurate’ GPS paths scaled according to the following speeds: periods of flying are coarsely resolved due to precision errors. Indeed, allocated 12 m/s, surface ‘rest’ periods allocated 0.1 m/s, even with little or no GPS error, this can greatly com- bottom phase of dives allocated 0.4 m/s and descent and promise movement estimates [cf. 136]. Conversely, the ascent speeds modulated according to “VPC dead-reck- precision of the dead-reckoned tracks is only limited by oning procedure in R” section—Eq.  (25). Note that ele- the amount of initial motion sensor data under-sampling vation was not resolved during flying periods (although (usually required in some capacity to make datasets more flying periods were dead-reckoned). Regardless of the manageable and less computationally expensive). Such current limitations, the VPC dead-reckoning proce- fine-scale estimates can therefore (with suitable VPC) dure represents a substantial advance for resolving, and allow users to define movement in space with unprec - thereby allowing investigation of, continuous, fine-scale, edented resolution. The benefit of this is that such reso - free-ranging 2- or 3-D space-use with all its underlying lution can resolve important metrics of movement, such scales of tortuosity and distances moved (e.g., Figs. 7 and as step duration [cf. 146] and the number and extent of 8). turns made [cf. 147]; useful parameters for investigat- ing navigation and foraging strategies according to envi- VPC dead‑reckoning procedure in R ronmental circumstance—though, such parameters are Preparing the three axes of rotation for derivation also useful without superimposing on the environment. of heading Moreover, even dead-reckoned tracks that are sparsely The tilt-compensated compass method is a well-known corrected or never corrected can detail important move- practice for deriving heading [e.g., 21, 22, 81]. Correct ment-specific behaviours [12], for example, circling coordinate system axis alignment and suitable calibra- behaviour [148]. tion of tri-axial magnetometry data [cf. 149] are crucial Ultimately, the higher the frequency at which dead- pre-processors, without which, heading estimates would reckoning is undertaken, the better the resolution and likely incorporate substantial error [cf. 21, 149]. The tilt- detail of reconstructed tracks. However, accuracy only compensated compass method described below (fol- improves up to a point because extrapolated travel vec- lowing the framework outlined by Pedley [21]), requires tors (heading and speed estimates) nearly always com- the aerospace (x-North, y-East, z-Down) (right-handed) prise some degree of error (no matter how small) and so, coordinate system, or ‘NED’ (cf. Additional file  1: Text with very high frequencies (> 1  Hz), more error is accu- S2, Fig. S1). We provide examples of axis alignment, out- mulated per unit time [cf. 16, 44]. In particular, when the line the importance of transforming between coordinate temporal resolution of dead-reckoning results in a spa- frames (relative to the Earths fixed frame) and recom - tial resolution dominated more by sensor noise than by mend a universal configuration calibration procedure to ‘actual’ movement of the animal in question, dead-reck- aid correct axis alignment within Additional file  1: Text oning accuracy will begin to decrease (at least pre-VPC). S2. The extent of this will depend on the size, speed and life - Multiple mathematically sophisticated algorithms style of the animal in question. For example, the benefits have been developed to correct distortions from each of dead-reckoning a lion at 40  Hz rather than 1  Hz are magnetometer channel’s output [e.g., 23, 149, 150, questionable (how often does a lion turn substantially 151, 152]. We provide an annotated R script—Gundog. within a second?), particularly given the additional com- Compass (Additional file 2) that corrects both soft and putation time (cf. Additional file  1: Text S6) and possible hard iron distortions from tri-axial magnetometry data Gunner et al. Anim Biotelemetry (2021) 9:23 Page 16 of 37 Fig. 7 Twelve outgoing (green) and incoming (blue) dead‑reckoned trajectories from Magellanic penguins walking to and from their nest. Three variants of track advancement were used: a A VeDBA threshold (0.1 g) and constant m‑ coefficient (1.4) (b), depicted movement periods using the LoCoD method to identify steps (cf. Wilson et al. 2018) and constant m‑ coefficient (1.4) and c depicted individual steps within depicted movement periods, from which a constant distance estimate (0.16 m) was multiplied by step frequency (x̄ no. steps/s) (full details within Additional file 1: Text S4) (c). Note that the accuracy with respect to the radial distance can be evaluated by examining the track stops in relation to the shore‑line. DBA‑ derived speed estimates were typically overestimated for incoming tracks, due to the birds being heavier (and thus impart greater DBA per stide cycle) after foraing. Tracks (from c) were GPS‑ corrected (d) (method = “ distance” , dist.step = 5, VP .ME = TRUE, thresh = between 8 and 15 (depending on track length) approx. every 50 m). A portion of the GPS‑ corrected dead‑reckoned tracks (bottom panel) are magnified (2 iterations) to show the difference in resolution of movement tortuosity, between GPS and dead‑reckoned tracks Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 17 of 37 Fig. 8 GPS‑ corrected dead‑reckoned tracks of Imperial cormorants foraging at sea: a 15 birds (blue = male, red = female). b Shows one of these tracks illustrated in 3‑D. Note gaps between dives are either associated with current drift, while the bird is resting at the sea surface, or periods of flight. c and d show the descent, bottom phase and ascent of a given dive in both 2‑D (c) and 3‑D, respectively and subsequently computes tilt-compensated heading Tilt‑compensated heading derivation (0° to 360°). Within this function, there are two main Device orientation is expressed in terms of a sequence of methods of correction to choose from, based on the Euler angle [roll (Φ), pitch (θ), yaw (Ψ)] rotations about mathematical protocols outlined by Vitali [153]—least- the x-, y- and z-axes, respectively, relative to the (iner- square error approximation (constructing an ellipsoid tial) Earths fixed frame of reference (e.g., Earth-Centre, rotation matrix) and Winer [154]—scale biases with Earth-Fixed (ECEF) system) [155]. Being a vector field simple orthogonal rescaling (avoiding matrices alto- sensor with two degrees of rotational freedom, acceler- gether). We expand on this user-defined function- ometers are insensitive to rotations about the gravity vec- ality, as well as outlining the causes of soft and hard tor and thus discerning heading requires the arctangent iron distortions and the initial calibration procedure of the ratio between the x- and y-orthogonal magnetom- required to correct such distortions within Additional eter measurements [156]. For the correct computation of file 1: Text S3. heading, these two channels need to be aligned parallel to Gunner et al. Anim Biotelemetry (2021) 9:23 Page 18 of 37 the earth’s surface. This is achieved by correcting any ori - animals, high frequency of body posture changes could entation (de-rotation) according to pitch and roll angles cause discrepancy between static acceleration data and (postural offsets) which can be deduced from accelera - magnetism data, which could consequently affect head - tion. These angles are typically approximated by deriving ing estimation [162]. Although this effect would not gravity-based (static) acceleration [see 72, 157] from each change general shapes of movement paths, we suggest channel by employing one of four approaches using: (i) that prior to the normalisation process (and magnetic a running mean [e.g., 72, 86], (ii) a Fast-Fourier transfor- calibration procedure), it may be of value to initially mation [e.g., 158], (iii) a high-pass filter [e.g., 159] or (iv) smooth out (see Eq. 1, R ) small deviations within mag- 1:4 a Kalman-filter [e.g., 160]. Here, we use a computation- netometry data, both to avoid this type of error and to ally simple running mean over 2 s [72] (Eq. 1): reduce the magnitude of anomalous spikes in magnetic inference. We used a smoothing window of 10 events for i+ the 40 Hz datasets used in this study. G = A x,y,z x,y,z (1)     j=i− NG G x x     NG � G = · y y G · G + G · G + G · G where w is an integer specifying the window size and x x y y z z NG G z z G and A represents the smoothed and raw compo- x,y,z x,y,z (2) nents of acceleration, respectively. In the absence of lin-     NM M x x ear (dynamic) acceleration [see 157, 161], values of G x,y,z 1     NM � M = · y y reflect the device orientation with respect to the earth’s M · M + M · M + M · M x x y y z z NM M z z reference frame (though see Table. 1), reading approx. (3) + 1  g when orientated directly towards the gravity vec- tor (down), − 1  g against the gravity vector (up) and 0  g (R5) NGx = Gx / sqrt(Gx^2 + Gy^2 + at perpendicular to it (horizontal). In R, the ‘zoo’ package Gz^2) [126] provides useful wrappers to apply arithmetic opera- (R6) NGy = Gy / sqrt(Gx^2 + Gy^2 + tions in a rolling fashion (R ). 1:4 Gz^2) (R7) NGz = Gz / sqrt(Gx^2 + Gy^2 + (R1) install.packages("zoo") ; Gz^2) library(zoo) (R8) NMx = Mx / sqrt(Mx^2 + My^2 + (R2) Gx = rollapply(Ax, width=w, Mz^2) FUN=mean, align="center", (R9) NMy = My / sqrt(Mx^2 + My^2 + fill="extend") Mz^2) (R3) Gy = rollapply(Ay, width=w, (R10) NMz = Mz / sqrt(Mx^2 + My^2 + FUN=mean, align="center", Mz^2) fill="extend") (R4) Gz = rollapply(Az, width=w, Depending on deployment position, the device-carried FUN=mean, align="center", NED coordinate frame (the x-, y-, z-axes) may not corre- fill="extend") spond with the animal’s body-carried NED frame. When this occurs, prior to deriving animal orientation, the nor- Here, w should be replaced with the window width malised gravity and magnetic vectors are required to be of choice (e.g., for 20  Hz data and a smoothing of 2  s corrected so that their measurements are expressed rela- required, replace w with 40). We use a centre-aligned tive to the body frame of the animal [45]. This requires index (compared to the rolling window of observations), three rotation sequences, using 3 by 3 rotation matrices with “extend” to indicate repetition of the leftmost or (Eqs. 4, 6) and involves two intermediate frames. The aer - rightmost non-NA value (though fill can equally be set as ospace sequence used here is as follows: NA, 0, etc.). Importantly, for correct trigonometric formulae out- 1. A right-handed rotation (C ), about the z-axis axis of put within the tilt-compensated compass method, the the device’s frame ( D ), through angle Ψ (Eq. 4), to get vectorial sum of static acceleration (G ) and calibrated x,y,z to the first intermediate frame (F ). magnetometry (M ) measurements across all three x,y,z 1 2. A right-handed rotation (C ) about the y-axis at F , spatial-dimensions must be normalised (to a unit vector) 1 through angle θ (Eq. 5), to get to the second interme- with a scaled magnitude (radius) of one (Eqs.  2, 3, R ). 5:10 diate frame ( F ). It was previously demonstrated that, for fast moving 2 Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 19 of 37 3. A right-handed rotation (C ) about the x-axis at F , The product of the conventionally used aerospace rota - through angle Φ (Eq.  6), to get to the animal’s body tion sequence outlined above (to get from the tag frame frame ( B). to the animal’s body frame) can be expressed as (Eq. 7).   (�,θ,�) (�) (θ) (�) cos (�) sin (�) 0 C = C · C · C (7) B/D B/F F /F F /D 2 2 1 1 (�)   − sin (�) cos (�) 0 C = (4) F /D 0 01 When matrix multiplied out, this yields (Eq.  8)—often referred to as a Direction Cosine Matrix (DCM). The   composition of this DCM varies according the (six pos- cos (θ) 0 − sin (θ) (θ) sible) orderings of the three rotation matrices (Eqs.  4,   01 0 C = (5) F /F 2 1 6) and the direction of intended rotation relative to the sin (θ) 0 cos (θ) direction of measured g within the NED system (see Additional file 1: Text S2).   cos (�) · cos (θ) sin (�) · cos (θ) − sin (θ) (�,θ ,�)   (8) C = cos (�) · sin (θ) · sin (�) − sin (�) · cos (�) sin (�) · sin (θ) · sin (�) + cos (�) · cos (�) cos (θ) · sin (�) B/D cos (�) · sin (θ) · cos (�) + sin (�) · sin (�) sin (�) · sin (θ) · cos (�) − cos (�) · sin (�) cos (θ) · cos (�) Note the left-handed rule of reading the vectorial nota-   (�,θ ,�) tion of ordered rotations, for example C means 10 0 B/D (�)   going from the device frame to the animal’s body frame, C = 0 cos (�) sin (�) (6) B/F by first rotating about the z-axis (though angle  ), fol- 0 − sin (�) cos (�) lowed by the y-axis (though angle θ) and then lastly the Note the right-handed rule of rotation; a positive x-axis (though angle  ). The device offset can be esti - reflects a clockwise rotation of the anterior–poste - mated from direct observation or deduced using photo- rior axis (relative to North), a positive θ reflects a nose- graphs or from the tag data itself. For example, assuming upward tilt of this axis and a positive  reflects a bank that ‘normal animal posture’ has no pitch and roll angle angle tilt to the right about this axis. Reversing the offset, then a tri-axial spherical plot of static acceleration direction of two axes causes a 180° inversion about the [164] would show a densely populated band of datapoints remaining axis and interchanging two axes (e.g., x with y) at the cross-sectional origin of 0  g about the x- and or reversing the direction of one or all three axes reverses y-axes, respectively, when the tag and body NED axes are the ‘handedness’ of rotation [right-handed—‘counter- in alignment. clockwise’ vs. left-handed—‘clockwise’ (when viewed NG and NM are pre-multiplied by the DCM to x,y,z x,y,z from the tip of the z-axis)]. Rotation matrices are orthog- compensate for offset. However, device offset is often onal (unitary), with every row and column being linearly parametrised by roll, pitch and/or yaw angles relative independent and normal to every other row and column. to the animal’s body frame and thus, the device actually The consequence of this is that the inverse of a rotation requires de-rotation (switching the ‘handedness’ of rota- matrix is its transpose [163] [which essentially reverses tion) according to these values. For example, a + 45° yaw the direction of rotation, and within (Eqs.  4, 6), this is offset requires an inverse rotation about the z-axis by achieved by negating the sign of the sines]. Importantly, − 45°, rather than a further + 45° rotation. This simply because rotation matrices are not symmetric, the order of involves taking the transpose of the DCM (Eq.  9), which matrix multiplication is important [45] (otherwise, Euler is the same as the transpose of each of the individual angles are without meaning for describing orientation). rotation matrices (Eqs. 10, 11).   cos � · cos θ cos � · sin θ · sin � − sin � · cos � cos � · sin θ · cos � + sin � · sin � ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (�,θ,�)   (9) sin � · cos θ sin � · sin θ · sin � + cos � · cos � sin � · sin θ · sin � − cos � · sin � C ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) B/D − sin (θ) cos (θ) · sin (�) cos (θ) · cos (�) Gunner et al. Anim Biotelemetry (2021) 9:23 Page 20 of 37     B D PitchSinAngle * RollSinAngle - YawSi- NGb NG x x T T T (�) (θ) (�) nAngle * RollCosAngle) + NMz *     NGb = C · C · C · NG (10) y y B/F F /F F /D 2 2 1 1 (YawCosAngle * PitchSinAngle * RollCo- NGb NG x x sAngle + YawSinAngle * RollSinAn- gle)     B D NMb NM (R21) NMby = NMx * YawSinAngle * Pitch- x x T T T (�) (θ) (�)     NMb = C · C · C · NM y y CosAngle + NMy * (YawSinAngle * B/F F /F F /D 2 2 1 1 NMb NM x x PitchSinAngle * RollSinAngle - YawSi- (11) nAngle * RollCosAngle) + NMz * (YawSinAngle * PitchSinAngle * RollSi- where T is the matrix transpose and resultant NGb x,y,z nAngle - YawCosAngle * RollSinAn- and NMb vectors are expressed in the animal’s body- x,y,z gle) carried NED frame. The input of these gravity and mag - (R22) NMbz = -NMx * PitchSinAngle + NMy netic vectors are supplied as 3 by 1 column matrices for * PitchCosAngle * RollSinAngle + true matrix multiplication, and when expanding out NMz * PitchCosAngle * RollCosAngle (Eq. 11), this results in (Eq. 12) [substituting NM with NG Here, Roll, Pitch and Yaw inputs denote the angular off - expands out (Eq. 10)].     B D NMb NM · cos (�) · cos (θ) + NM · (cos (�) · sin (θ) · sin (�) − sin (�) · cos (�)) + NM · (cos (�) · sin (θ) · cos (�) + sin (�) · sin (�)) x x y z      NMb   NM · sin (�) · cos (θ) + NM · (sin (�) · sin (θ) · sin (�) + cos (�) · cos (�)) + NM · (sin (�) · sin(θ) · sin(�) − cos(�) · sin(�))  y x y z NMb −NM · sin (θ) + NM · cos (θ) · sin (�) + NM · cos (θ) · cos(�) x x y z (12) In R then, the alignment of device to body axes for both set of the device, relative to the animal body frame. Note, gravity and magnetic vectors can be performed using the standard trigonometric functions operate in radians, not following procedure (R ). degrees. In base R, π = pi. Multiplying values by pi/180 11:22 coverts degrees into radians, whilst multiplying values (R11) RollSinAngle = sin(Roll * pi/180) by 180/pi does the reverse. This rotation correction pro - (R12) RollCosAngle = cos(Roll * pi/180) cedure is implemented within Gundog.Compass when (R13) PitchSinAngle = sin(Pitch * pitch, roll and/or yaw offsets are supplied (Additional pi/180) file 2). (R14) PitchCosAngle = cos(Pitch * Following the alignment of device and body axes, pitch pi/180) and roll of the animal are calculated from the DCM, and (R15) YawSinAngle = sin(Yaw * pi/180) because there are multiple variations in the order that (R16) YawCosAngle = cos(Yaw * pi/180) rotation sequences can be composed and applied, there (R17) NGbx = NGx * YawCosAngle * Pitch- are also different valid equations that output different CosAngle + NGy * (YawCosAngle * pitch and roll angle estimates, for equivalent static accel- PitchSinAngle * RollSinAngle - YawSi- eration input. The convention is to use formulae that nAngle * RollCosAngle) + NGz * have no dependence on yaw rotation and restrict either (YawCosAngle * PitchSinAngle * RollCo- the pitch or the roll angles within the range − 90° to + 90° sAngle + YawSinAngle * RollSinAn- (but not both), with the other axis of rotation able to lie gle) between − 180° and 180°, thereby eliminating duplicate (R18) NGby = NGx * YawSinAngle * Pitch- solutions at multiples of 360°. Multiplying (Eq.  8) by the CosAngle + NGy * (YawSinAngle * measured Earth’s gravitational field vector (+ 1  g when PitchSinAngle * RollSinAngle + YawCo- initially aligned downwards along the z-axis) simplifies sAngle * RollCosAngle) + NGz * down to (Eq.  13). The accelerometer output for this aer - (YawSinAngle * PitchSinAngle * RollSi- ospace rotation sequence is thus only dependent on the nAngle - YawCosAngle * RollSinAn- roll and pitch angles which can be solved (Eqs.  14, 15), gle) allowing roll angles the greater freedom [161]. This is (R19) NGbz = -NGx * PitchSinAngle + NGy relevant for studies using collar-mounted tags, whereby * PitchCosAngle * RollSinAngle + collar may roll > 90° in either direction from default NGz * PitchCosAngle * RollCosAngle orientation. (R20) NMbx = NMx * YawCosAngle * Pitch- CosAngle + NMy * (YawCosAngle * Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 21 of 37         0 cos (�) · cos (θ) sin (�) · cos (θ) − sin(θ) 0 − sin (θ)         NGb = · =  0   cos (�) · sin (θ) · sin (�) − sin (�) · cos (�) sin (�) · sin (θ) · sin (�) + cos (�) · cos (�) cos(θ) · sin(�)   0   cos (θ) · sin (�)  xyz 1 cos (�) · sin (θ) · cos (�) + sin (�) · sin (�) sin (�) · sin (θ) · cos (�) − cos (�) · sin (�) cos(θ) · cos(�) 1 cos (θ) · cos (�) (13)   � � � � −NGb 180   tan θ = � ⇒ θ = atan2 −NGb , NGb · NGb + NGb · NGb · (14) xyz x y y z z NGb + NGb 2 2 y z containing 1  s and − 1  s according to the direction of NGb 180 tan � = ⇒ � = atan2 NGb , NGb · xyz y z measured g from NGb (R ). NGb π 23 The magnetic vector of the device is then de-rotated (15) to the Earth frame (tilt-corrected) by pre-multiplying by The equation for roll (Eq.  15), however, has a region of the product of the inverse roll multiplied by inverse pitch instability at obtuse pitch angles (e.g., for NED systems, rotation matrix (Eq. 17), which when expanded out gives the x-axis points directly up or down, with respect to the (Eq. 18). Earth’s frame of reference). Whilst there is no ‘gold stand-     ard’ solution to this problem of singularity (using Euler NMbf cos (θ) sin (θ) · sin (�) sin (θ) · cos (�) angles), an attractive circumvention (detailed within     NMbf = 0 cos (�) − sin (�) [161]) is to modify (Eq. 15) and add a very small percent- NMbf − sin (θ) cos (θ) · sin (�) cos (θ) · cos (�)   age ( µ ) of the NGb reading into the denominator, pre- NMb venting it ever being zero and thus driving roll angles to   · NMb zero when pitch approaches −/+90° for stability (Eq. 16). NMb (17)     NMbf NMb · cos (θ) + NMb · sin (θ) · sin (�) + NMb · sin (θ) · cos (�) x x y z     NMbf = NM · cos (�) − NMb · sin (�)     y y z NMbf −NMb · sin (θ) + NMb · cos (θ) · sin (�) + NMb · cos (θ) · cos (�) z x y z (18) Here NMbf are the calibrated, normalised mag- x,y,z � = atan2 NGb , sign(NGb ) y z netometry data (expressed in the animal’s body-carried NED frame) after tilt-correction. Finally, yaw (ψ) (head- · (NGb · NGb + µ · NGb · NGb ) · z z x x ing—now defined by the compass convention, relative to (16) magnetic North) can be computed from the NMbf and where sign(NGb ) is allocated the value + 1 when NGb z z NMbf (Eq. 19) via; is non-negative and − 1, when NGb is negative (recovers directionality of NGb , subsequent to the square-root). ψ = atan2 −NMbf , NMbf · (19) y x Taken together then, in R, pitch and roll are computed according to, (R ) with outputs within the range of 24:25 We outline the R code for this procedure below (R ). 26:34 − 90° to + 90° for pitch and − 180° to + 180° for roll, and this is the formula we use in the tilt-compensated method (R26) RollSinAngle = sin(Roll * pi/180) outlined below (and within Additional file 2). (R27) RollCosAngle = cos(Roll * pi/180) (R28) PitchSinAngle = sin(Pitch * (R23) mu = 0.01 ; sign = ifelse(NGbz >= pi/180) 0, 1, -1) (R29) PitchCosAngle = cos(Pitch * (R24) Pitch = atan2(-NGbx, sqrt(NGby^2 pi/180) + NGbz^2)) * 180/pi (R30) NMbfx = NMbx * PitchCosAngle + (R25) Roll = atan2(NGby, sign * NMby * PitchSinAngle * RollSinAn- sqrt(NGbz^2 + mu * NGbx^2)) * gle + 180/pi NMbz * PitchSinAngle * RollCosAngle (R31) NMbfy = NMby * RollCosAngle – Here, prior to the derivation of pitch and roll, μ is allo- NMbz * RollSinAngle cated the value 0.01 and a vector termed ‘sign’ is created, Gunner et al. Anim Biotelemetry (2021) 9:23 Page 22 of 37 (R32) NMbfz = -NMbx * PitchSinAngle + acceleration (cf. Eq.  1, R ) from their raw equivalent 1:4 NMby * PitchCosAngle + RollSinAn- (R ). gle + NMbz * PitchCosAngle * RollCosAngle (R37) v = sqrt((Ax - Gx)^2 + (Ay - (R33) Yaw = atan2(-NMfby, NMfbx) * 180/ Gy)^2 + (Az - Gz)^2), pi (R34) Yaw = ifelse(Yaw < 0, Yaw + 360, where Ax, Ay, Az and Gx, Gy, Gz are the raw and static Yaw) (smoothed) values of each channel’s recorded acceleration. Note, yaw output from (R ) uses the scale − 180° to We recommend implementing a running mean (cf. + 180°. (R ) converts to the scale 0° to 360 (specifi - Eq.  1, R ) to raw VeDBA values to ensure that both 34 1:4 cally, 0° to 359°). This is also achieved by using a modulus acceleration and deceleration components of a stride (mod) operator (Eq.  20, R ), which in base R takes the cycle are incorporated together per unit time and to form %%. reduce the magnitude of small temporal spikes (likely not attributable to the scale of movement elicited [cf. 97]. ψ = mod(360 + ψ, 360) (20) Choice of smoothing window size is dependent on the scale of movement being investigated, though as a basic (R35) Yaw = (360 + Yaw) %% 360 rule, we suggest 1 to 2  s. For similar reasons, it is also worth post-smoothing raw pitch, roll and heading out- Magnetic declination is defined as the angle on the puts, although heading requires a circular mean (Eqs. 22, horizontal plane between magnetic north and true 23) [cf. 169]: north [165]. Prior to dead-reckoning, magnetic decli-   nation should be summed to heading values to convert n n � � � � � � 1 π 1 π from magnetic to true North [166]. There are many   θ = a tan 2 sin h · , cos h · p j j n 180 n 180 online sources to calculate the magnetic declination of j=i j=i an area [e.g., 167]. Notably, logical corrections may need (22) to be performed to ensure data does not exceed either circular direction after applying magnetic declination h = mod 360 + θ · , 360 p (23) (R ). where h and h are the unsmoothed and smoothed (R36) h = ifelse(h < 0, h + 360, h) ; h heading values, θ the arithmetic mean after converting = ifelse(h > 360, h - 360, h), degrees to cartesian coordinates and mod refers to the modulo operator. where h refers to the vector containing the heading In R, the above formula can be made into a function data. Should the user not correct for axis align- (R ), to be applied within the ‘rollapply’ wrapper (replac- ment between the device and animal body frame ing ‘FUN = mean’ with ‘FUN = Circ.Avg’) (cf. R ). 1:4 (cf. Eqs. 4–12, R ) then a reasonable post-cor- 11:22 rection for small discrepancies about the yaw axis (R38) Circ.Avg = function(x){ would be to subtract the difference to h values at this point. H.East = mean(sin(x * pi / 180)) Preparing speed estimates H.North = mean(cos(x * pi / 180)) The vectorial dynamic body acceleration (VeDBA) (Eq. 21) [cf. 67, 168] was our choice of DBA-based speed MH =(atan2(H.East, H.North)) * 180/pi proxy for terrestrial dead-reckoning purposes. This is given by: MH = (360 + MH) %% 360 2 2 2 v = D + D + D (21) x y z return(MH) where v represents VeDBA, D , D and D are the x y z } dynamic acceleration values from each axis, themselves Speed ( s ) can be estimated from VeDBA ( v ) via (Eq. 24). obtained by subtracting each axis’ static component of s = (v · m) + c (24) Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 23 of 37 where m is the multiplicative coefficient and c is a con- thus should only be calculated at times when the animal stant [10, 69]. Here, a user can define various bouts is travelling ‘ballistically’ (at considerable vertical speed). of movement from motion sensor data (e.g., via vari- (R41) s = ifelse(abs(p) >= 10, abs(RCD ous machine-learning approaches (for review see Far- / tan(p * pi/180)), s) rahi et al. [170]) or the Boolean-based LoCoD method [101]) and/or substrate condition (e.g., via GPS), to be In the above example (R ), nominal speed values are cross-referenced when allocating variants of the speed overwritten with the trigonometric formula output coefficients. As a simple example, in R, should walk- (Eq.  25) at times of ‘appreciable’ pitch (10°) [cf. 171], ing (coded for as 1) and running (coded for as 2) be where RCD is the rate change of depth and p is the pitch teased apart from all other (non-moving) data (coded (in radians). An upper limit should be imposed on speed for as 0) within a Marked Events vector (ME), then values derived in this way because values can become ME can be used to allocate various m (and if applica- highly inflated when the pitch angle is particularly acute. ble, c ) values using simple ‘ifelse’ statements (R ). 39:40 Converting speed to a distance coefficient (R39) m = ifelse(ME == 1, 1.5, Speed (s) estimates are multiplied by the time differ - ifelse(ME == 2, 3.5, 0)) ence between the values (TD ) to give a distance estimate (R40) c = ifelse(ME > 0, 0.1, 0) (units in metres) which, in turn, standardises coefficient comparisons across datasets sampled at different rates. These distance values are then divided by the approxi - Here, walking is given an arbitrary coefficient of 1.5 mate radius of the earth (R = 6,378,137 m) to give a radial and running, 3.5 with a value of 0.1 for their constants. distance coefficient ( q ) [se e 172] (Eq. 26)” All other ME values are given a 0 coefficient and 0 con - stant, which results in no speed at such times, regardless s · TD of DBA magnitude. q = (26) By‑passing DBA as a speed proxy Assuming that high-resolution depth data are not avail- Dividing the number of steps detected within a given able, but ‘absolute’ speed estimates have been obtained, rolling window length (cf. R ), by the window length then an alternative to Eq.  25, (in accordance with the 1:4 (s) gives an estimated step count per second. This can equal pitch assumption) is to derive horizontal distance be converted to speed by multiplying by a distance per estimates by multiplying the absolute distance by the step estimate (assuming constant distance travelled cosine of the pitch ( θ ) (converted from degrees to radi- between step gaits). We review this further in Addi- ans), which can equally be performed on the radial dis- tional file  1: Text S4, including a simple peak finder tance (Eq. 27): function—Gundog.Peaks (Additonal file 3) that locates peaks based on local signal maxima, using a given roll- q = q · cos θ · (27) ing window, with each candidate peak filtered accord- ing to whether it surpassed a threshold height (in In R, to determine accurate lengths of time between conjunction with other potential user-defined thresh- values, it is best to save date and time variables together olds). Note, this method can equally be applied to as POSIX class [173]. Creating timestamp (TS) objects non-terrestrial species, using flipper/tail beats instead, with POSIXct class enables greater control and manipula- where appropriate. tion of time data. This makes computing the rolling time For diving animals, a proxy for horizontal speed can difference ( TD ; units in seconds) between data points be obtained based on animal pitch and rate change in simple (R ): depth [47, 116]. Specifically, rate change of depth ( d ) (R42) TD = c(0, difftime(TS, lag(TS), (units in m/s) is divided by the tangent of pitch ( θ ) units = "secs")[-1]) (converted from degrees to radians) (Eq. 25): We detail how to create timestamp objects of POSIXct �d s = class within Additional file  1: Text S5, including format- (25) tan θ · ting with decimal seconds (important for infra-second datasets) and various codes useful for manipulating data Here, resultant speed values need to be made absolute to be dead-reckoned based on time. (positive). This calculation is only valid when the direction In R then, following the computation of TD , q is of movement is the same as the direction of the animal’s obtained via (R ). 43:44 longitudinal axis (equal pitch assumption) [cf. 47] and Gunner et al. Anim Biotelemetry (2021) 9:23 Page 24 of 37 (R43) s = (v * m) + c DR.lon[i] = DR.lon[i-1] + (R44) q = (s * TD) / 6378137 atan2(sin(h[i]) * sin(q[i]) * cos(DR.lat[i-1]), cos(q[i]) - sin(DR. Note, if a negative c intercept is used (e.g., to allow for lat[i-1]) * sin(DR.lat[i])) some body movement without translation), then any neg- ative speed values would need to be equated to zero as an additional step. Reverse dead‑reckoning As previously mentioned, the ME vector (progressive For this, firstly, the time difference is computed as usual movement coded by integer values greater than zero (R ) and the dimensions of each vector required in the (e.g., 1) and stationary behaviour coded by zero) can be dead-reckoning calculation are reversed. We bind all used to ensure q (essentially the distance moved) is zero relevant vectors into a data frame (df ) (R ), subsequent when ME reads zero, ensuring dead-reckoned tracks are to reversing data frame dimensions (R ); the last row not advanced at such times, regardless of the computed becomes the first row, second to last row becomes the speed (R ). second, etc. Note, this can equally be achieved by using the rev() function within base R, on each individual vec- (R45) q = ifelse(ME == 0, 0, q) tor. These reversed columns are now restored as vectors (R ) and shifted forward by one element (R ). This is 53 54 Derivation of coordinates required for correct alignment in time so that dead-reck- Once q and h are obtained, coordinates are advanced oning works in exactly the opposite manner to ‘forward’ using (Eqs. 28, 29); dead-reckoning. Lat = asin(sin Lat · cos q + cos Lat · sin q · cos h) i 0 0 (R50) TD = c(0, difftime(TS, lag(TS), (28) units = "secs")[-1]) Lon =Lon + atan2((sin h · sin q · cos Lat ), i 0 0 (R51) df = data.frame(TD, h, v, m, c, (cos q − sin Lat · sin Lat )) ME) 0 i (29) (R52) df = df[dim(df)[1]:1, ] where Lat , Lat and Lon , Lon are the previous and pre- 0 i 0 i (R53) TD = df[, ’TD’] ; h = df[, ’h’] ; sent latitude and longitude coordinates, respectively (in v = df[, ’v’] ; radians), h is the (present) heading (in radians) and q is m = df[, ’m’] ; c = df[, ’c’] ; ME = the (present) distance coefficient. df[, ’ME’] In R, the above can be performed iteratively within a (R54) TD = c(NA, TD[-length(TD)]) ; h = for-loop (iteration of code repeated per consecutive c(NA, h[-length(h)]) ; ith element of data; R ). Initialising the output latitude v = c(NA, v[-length(v)]) ; m = c(NA, (DR.lat) and longitude (DR.lon) variables to the required m[- length(m)]) ; length (e.g., as governed by the vector length of other c = c(NA, c[-length(c)]) ; ME = c(NA, input data (heading, speed, etc.) speeds up processing ME[-length(ME)]) time (R ). Within the trigonometric dead-reckoning for- The next step is to rotate heading 180° and correct for mulae, the starting latitude (la) and longitude (lo) coordi- its circular nature (R ). nates and heading ( h ) values must be supplied in radians (R ). The la and lo values are saved as the first elements (R55) h = h - 180 ; h = ifelse(h < 0, h of the DR.lat and DR.lon vectors to be advanced, respec- + 360, h) tively (R ). Lastly, is determined and DR.lon and DR.lat are advanced based on the dead-reckoning formula (cf. (R46) DR.lat = rep(NA, length(h)) ; R ), except in this instance, the first element of DR. 46:49 DR.lon = rep(NA, length(h)) lon and DR.lat needs to be supplied by the ‘known’ last lo (R47) la = la * pi/180 ; lo = lo * and la coordinates. pi/180 ; h = h * pi/180 (R48) DR.lat[1] = la DR.lon[1] = lo Integrating current vectors (R49) for(i in 2:length(DR.lat)) { In R, current vectors can be added according to (R ). 56:60 DR.lat[i] = asin(sin(DR.lat[i-1]) * Current speed (cs) is in m/s (ensure values are absolute) cos(q[i]) n+ cos(DR.lat[i-1]) * and current heading (ch) uses the scale 0° to 360°. Note sin(q[i]) * cos(h[i])) the use of ‘yy’ and ‘xx’ vectors, storing the previous DR. lat and DR.lon coordinates prior to implementing the Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 25 of 37 next ‘current drift’ vector per iteration. The current speed estimates from motion data—though they are essentially is also standardised according to the time period length the same). Lat − Lat Lon − Lon i 0 i 0 −1 2 2 (30) d = 2·R·sin sin + cos (Lat ) · cos (Lat ) · sin 0 i 2 2 where R is the Earth’s radius and d , the output in metres. sin (�Lon) · cos (Lat ), and Earth’s radius (analogous to the derivation of q ). i b = atan2 · cos (Lat ) · sin (�Lat) · cos (Lat ) · cos (�Lon) π 0 i When reverse dead-reckoning, it is important to ensure (31) that cs and ch are included in the steps outlined above (R ). where Lon represents Lon − Lon , Lat represents i 0 50:55 Lat − Lat and b output is in the scale − 180° to + 180°. To i 0 convert b to the conventional 0° to 360° scale, 360 should (R56) DR.lat = rep(NA, length(h)) ; be added to values < 0. DR.lon = rep(NA, length(h)) For each VP, the distance is divided by the dead-reck- (R57) xx <- rep(NA, length(cs)) ; yy <- oned distance providing a distance correction factor rep(NA, length(cs)) (ratio; Eq.  32). The heading correction factor is com - (R58) la = la * pi/180 ; lo = lo * puted by subtracting the dead-reckoned bearing from pi/180 ; the VP bearing (Eq.  33). To ensure that difference does h = h * pi/180 ; ch = ch * pi/180 not exceed 180° in either circular direction, 360 should (R59) DR.lat[1] = la DR.lon[1] = lo be added to values < − 180 and 360 subtracted from val- (R60) for(i in 2:length(DR.lat)) { ues > 180. A simple example of why this is relevant can be DR.lat[i] = asin(sin(DR.lat[i-1]) * illustrated by subtracting a dead-reckoned bearing value cos(q[i]) + cos(DR.lat[i-1]) * of 359° from a VP bearing value of 1°—post-correction, sin(q[i]) * cos(h[i])) the difference is + 2°. yy[i] = DR.lat[i] DR.lon[i] = DR.lon[i-1] + Distance VP atan2(sin(h[i]) * sin(q[i]) * Distance = corr.factor (32) Distance DR cos(DR.lat[i-1]), cos(q[i]) - sin(DR. lat[i-1]) * sin(DR.lat[i])) Heading = Bearing − Bearing (33) xx[i] = DR.lon[i] corr.factor VP DR DR.lat[i] = asin(sin(yy[i]) * All intermediate q values are multiplied by the distance cos((cs[i] * TD[i]) / 6378137) + correction factor and the heading correction factor is cos(yy[i]) * sin((cs[i] * TD[i]) / added to all intermediate h values (ensuring that h values 6378137) * cos(ch[i])) are in degrees). To ensure circular range is maintained DR.lon[i] = xx[i] + atan2(sin(ch[i]) * between 0° and 360°, 360 should be subtracted from val- sin((cs[i] * TD[i]) / ues > 360 and added to values < 0. 6378137) * cos(yy[i]), cos((cs[i] * Specifically, we follow the protocol illustrated within TD[i]) / 6378137) - sin(yy[i]) * Fig.  9 for intermediate values. Note the formulae to cal- sin(DR.lat[i])) culate both distance ( d ; Eq.  30) and bearing ( b ; Eq.  31) between two points, are also used to recalculate both the heading ( h ) and radial distance ( q ) between current-inte- VPC procedure grated dead-reckoned fixes (pre-VPC; cf. R ). Note that Specifically, this method entails calculating the difference the Haversine distance is required to be converted back of Haversine distance (net error) and bearing (from true to radial distance by dividing by R (R = 6,378,137). North) between consecutive VPs and the correspond- In R, the formulae to calculate the great-circle dis- ing time-matched dead-reckoned track positions. The tance and great circular bearing are saved within the trigonometric Haversine formulae (Eqs.  30, 31) are used disty (R ) and beary (R ) functions, respectively, where 61 62 to calculate the great-circle distance ( d ) and great circu- lon1, lat1, long2 and lat2 represent longitude and lati- lar bearing ( b ) between consecutive VPs and consecu- tude positions (decimal format) at t and t , ( t repre- i i+1 tive (time-matched) dead-reckoned positions (note we senting time). use the term ‘bearing’ to differentiate between heading Gunner et al. Anim Biotelemetry (2021) 9:23 Page 26 of 37 Fig. 9 Schematic diagram illustrating the order of fixes used when calculating the Distance and Heading (difference of both GPS corr.factor corr.factor and dead‑reckoned (DR) positions between arrow heads). Note the discrepancy with the order at which these correction factors are applied to intermediate DR positions (as denoted by colour shading). Known starting position denoted with asterisk (R61) disty = function(long1, lat1, d = 6378137 * c long2, lat2) { return(d) long1 = long1 * pi/180 ; long2 = long2 * pi/180 ; lat1 = lat1 * } (R62) beary = function(long1, lat1, pi/180 ; lat2 = lat2 * pi/180 long2, lat2) { long1 = long1 * pi/180 ; long2 = long2 * pi/180 ; lat1 = lat1 * pi/1 a = sin((lat2 - lat1) / 2) * sin((lat2 80 ; lat2 = lat2 * pi/180 - lat1) / 2) + cos(lat1) * a = sin(long2 - long1)*cos(lat2) b = cos(lat1) * sin(lat2) - sin(lat1) cos(lat2) * sin((long2 - long1) / 2) * * cos(lat2) * cos(long2 - long1) sin((long2 - long1) / 2) c = ((atan2(a, b) / pi)*180) return(c) c = 2 * atan2(sqrt(a), sqrt(1 - a)) } Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 27 of 37 Below, we outline an example of VPC in R and assume distances between VPs are stored within the column VP coordinates (decimal format) are aligned in the same termed VP.distance (R ). The VP.distance is divided by length vectors/columns as motion sensor-derived data, the DR.distance to provide the distance correction factor, e.g., heading, DBA/speed, etc., with the corresponding termed Dist.corr.factor (R ). Importantly here, an ifelse indexed (element-/row-wise) time. Typically, motion sen- statement is incorporated so that Dist.corr.factor defaults sor data are recorded at much higher frequency so that to zero at times when both VP.distance and DR.distance there are many dead-reckoned fixes between sequential are zero (otherwise dividing zero by zero in R produces VPs. As such, in the example below, we assume NAs are NaN’s). expressed in the VP longitude and latitude fields at times of missing locational data. This approach of synchronis - (R66) df.sub$DR.loni = c(df.sub[-1, ing VP—with motion sensor data also applies when inte- ’DR.longitude’], NA) grating current data; assuming ch and cs are element/ row-wise matched to the relevant VP grid node. df.sub$DR.lati = c(df.sub[-1, ’DR. Firstly, an indexing row number (Row.number) vector, latitude’], NA) the length of the data used in the dead-reckoning opera- tion (e.g., h ) is created (R ), which is relevant for merg- df.sub$VP.loni = c(df.sub[-1, ’VP. ing full-sized and under-sampled data frames together longitude’], NA) (seen later). Together, the row number, (uncorrected) dead-reckoned longitude and latitude coordinates, VP df.sub$VP.lati = c(df.sub[-1, ’VP. longitude and latitude coordinates, heading and the latitude’], NA) radial distance vectors are inputted column-wise into a (R67) df.sub$DR.distance= disty(df. ‘main’ data frame, termed ‘df ’ (R ; user-assigned column 64 sub$DR.longitude, names of each vector are within quotation marks). This df.sub$DR.latitude, df.sub$DR.loni, data frame is then filtered removing rows with missing df.sub$DR.lati) VP data and stored as df.sub (R ). This under-sampled 65 (R68) df.sub$VP.distance= disty(df. data frame thus, row-wise, contains the time-matched sub$VP.longitude, dead-reckoned and ground-truthed positions. The VPC df.sub$VP.latitude, df.sub$VP.loni, process is analogous for reverse dead-reckoned tracks— df.sub$VP.lati) although VP.lon and VP.lat must also be reversed [Row. (R69) df.sub$Dist.corr.factor = number remains in ascending order (not reversed)]. The ifelse(df.sub$VP.distance == 0 & first element of VP.lon and VP.lat must be the lo and la, df.sub$DR.distance == 0, 0, df.sub$VP. respectively (or for reverse dead-reckoning, the last ele- distance / df.sub$DR.distance) ment prior to reversing these vectors). Analogous to the distance correction, the bearings between consecutive dead-reckoned estimates are stored (R63) Row.number = rep(1:length(h)) within the column termed DR.head (R ) and the cor- (R64) df = data.frame(Row.number, ’DR. responding bearings between VPs are stored within the longitude’ = DR.lon, column termed VP.head (R ). Logical corrections are ’DR.latitude’ = DR.lat, ’VP.longitude’ performed to convert both to the 0° to 360° scale (R ), = VP.lon, DR.head is subtracted from VP.head providing the head- ’VP.latitude’ = VP.lat, h, q) ing correction factor, termed Head.corr.factor (R ) (R65) df.sub = df[!with(df, is.na(VP. and further logical corrections are performed to ensure longitude) | is.na(VP.latitude)) a minimum and maximum difference range between ,] − 180° to + 180 (R ). Both sets of dead-reckoned and VP coordinates are (R70) df.sub$DR.head = beary(df.sub$DR. shifted backwards one row within new columns termed; longitude, DR.loni, DR.lati, VP.loni, VP.lati (R ). Row-wise, these columns represent the consecutive fix at t with their i+1 df.sub$DR.latitude, df.sub$DR.loni, originals being t . This provides the correct format for df.sub$DR.lati) the inputs required within the disty (cf. R ) and beary 61 (R71) df.sub$VP.head = beary(df.sub$VP. (cf. R ) functions. The distances between consecu - 62 longitude, tive dead-reckoned estimates are stored within the col- df.sub$VP.latitude, df.sub$VP.loni, umn termed DR.distance (R ) and the corresponding 67 df.sub$VP.lati) Gunner et al. Anim Biotelemetry (2021) 9:23 Page 28 of 37 (R72) df.sub$DR.head = ifelse(df. coordinates, heading and radial distance each time) until sub$DR.head < 0, dead-reckoning fixes accord ‘exactly’ (Gundog.Tracks uses df.sub$DR.head + 360, df.sub$DR.head) a threshold of 0.01 m) with ground-truthed locations. An df.sub$VP.head = ifelse(df.sub$VP.head important pitfall of the correction process to consider is < 0, that dividing ‘any’ value (e.g., > 0) by 0 results in infinite df.sub$VP.head + 360, df.sub$VP.head) (Inf) values in R. This can arise during the correction (R73) df.sub$Head.corr.factor = process when there is a given distance between con- df.sub$VP.head - df.sub$DR.head secutive VPs, but no displacement between the accord- (R74) df.sub$Head.corr.factor = ing dead-reckoned positions. This can be a consequence ifelse(df.sub$Head.corr.factor < of ground-truthing too frequently (typically relevant to -180, high-res GPS studies), where positional noise is more (df.sub$Head.corr.factor + 360), apparent during rest periods [cf. 97] and/or wrongly df.sub$Head.corr.factor) assigned speed estimates/ME values. Gundog.Tracks df.sub$Head.corr.factor = ifelse(df. automatically resamples VPC rate where necessary to sub$Head.corr.factor > 180, avoid Inf values, essentially by changing the VPC rate to df.sub$Head.corr.factor - 360), avoid using successive VPs at times of no dead-reckoned df.sub$Head.corr.factor) track advancement. Lastly, Gundog.Tracks outputs mes- Only the relevant columns; Row.number, Dist.corr.fac- sages to the user’s console, detailing up to six stages of tor and Head.corr.factor are preserved (R ) and merged dead-reckoning progression, which includes reporting back into the main data frame (df ) based on the matching the maximum distance (units in metres) between dead- row numbers (R ). Both Dist.corr.factor and Head.corr. reckoned- and ground-truthed positions (used within factor express NA’s between VPs. These are replaced with the VPC procedure) at each iteration of correction and the most recent non-NA (observations carried forwards; whether automatic VPC resampling due to Inf values R ). Dist.corr.factor and Head.corr.factor values are occurred. shifted forward by one row (R ) for correct alignment purposes with respect to h and q values to be adjusted Conclusion (cf. Fig.  9; R ). A logical correction is performed to We have provided a comprehensive, fully integrated 79:80 ensure that a 0° to 360° circular scale is maintained after application of the dead-reckoning procedure within the the heading correction (R ). Note, the na.locf() function framework of the programming language, R, from pre- is required from the ‘zoo’ package, to replace NA values processing raw tri-axial accelerometery and magnetom- with the last non-NA value. etry data to VPC dead-reckoning. We have highlighted important considerations to increase the accuracy of the (R75) df.sub = df.sub[, c(’Row.number’, analytical procedure and to avoid misinterpretation of ’Dist.corr.factor’, ’Head.corr. error. We have also supplied extensive Additional files factor’)] 1, 2, 3, 4 and 5 and supporting functions to aid the pro- (R76) df = merge(df, df.sub, by = "Row. cess of deriving fine-scale movement paths, including number", all = TRUE) the protocols to correct magnetometry data and derive (R77) df$Dist.corr.factor = (tilt-compensated) heading. Importantly, we have dem- na.locf(df$Dist.corr.factor) onstrated the value of Gundog.Tracks; a multi-functional df$Head.corr.factor = na.locf(df$Head. and user-friendly tool to derive animal movement paths corr.factor) across all media of travel, with detailed input flexibil - (R78) df$Dist.corr.factor = c(NA, ity and output summaries. We suggest the next phase in df$Dist.corr.factor[-nrow(df)]) advancing the utility of animal dead-reckoning includes df$Head.corr.factor = c(NA, df$Head. looking for ‘track signatures’ that may signify a particu- corr.factor[-nrow(df)]) lar behaviour or reference a particular ‘ground-truthed’ (R79) q = (df$q * df$Dist.corr.factor) location. Lastly to advance the utility of Gundog.Tracks, (R80) h = (df$h + df$Head.corr.factor) we aim to optimise future iterations of the online code to (R81) h = ifelse(h > 360, h - 360, h) ; speed up computation time on larger datasets (e.g., sub- h = ifelse(h < 0, h + 360, h) second data collected over many months). Supplementary Information These updated coefficients are substituted into the The online version contains supplementary material available at https:// doi. dead-reckoning formula (cf. R ) and this process is 46:49 org/ 10. 1186/ s40317‑ 021‑ 00245‑z. repeated iteratively (using the updated dead-reckoned Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 29 of 37 Declarations Additional file 1. Methods expanded. Text S1. Device set up and cap ‑ ture protocol. Text S2. The importance of having the correct coordinate Ethics approval and consent to participate system and axis alignment. Text S3. Magnetometer calibration, rotation We thank the Conservation Agency from the Chubut Province, Argentina, correction and deriving yaw (heading)—Gundog.Compass() explained. for the permits to work at Punta León and Península Valdés protected areas Text S4. Step counts as a distance estimate—Gundog.Peaks() explained. (Disp No. 047/19‑SsCyAP). All penguin and cormorant handling procedures Text S5. Time Data in R (POSIXct). Text S6. VPC dead‑reckoning—Gundog. were reviewed and approved by the Dirección de Fauna y Flora Silvestre y el Tracks()—explained. Ministerio de Turismo y Áreas Protegidas de la Provincia de Chubut (permits Additional file 2. Gundog.Compass() (.R file). to work at San Lorenzo and Punta León, No. 060/19‑DFyFS‑MP and No. 047‑ SsCy/19). Ethical approval was also given by Animal Welfare Ethical Review Additional file 3. Gundog.Peaks() (.R file). Body (AWERB), approval number: SU‑Ethics‑Student ‑260919/1894, reference: Additional file 4. Gundog.Tracks() (.R file). IP‑1819‑30. Conditions and approvals for lion fieldwork were granted by the Animals Scientific Procedures Act (ASPA) at Queens University of Belfast (QUB‑ Additional file 5. Step by step guide of using Gundog.Tracks (.R file) (to BS‑AREC‑18‑006) and Pretoria University (NAS061‑19), permit authorisation use in conjunction with below). was given by South African National Parks (Permit Number SCAM 1550). Additional file 6. Raw sensor and GPS data frame (.txt) of a penguin walk ‑ ing out to sea from its nest. Consent for publication Not applicable. Acknowledgements Competing interests We thank South African National and the Department of Wildlife and National The authors declare no conflict of interest. Parks, Botswana, for allowing our research in the Kgalagadi Transfrontier Park. We are grateful to support and kind assistance of the staff and Rangers at the Kgalagadi National Park who were involved with this work, especially Steven Glossary Smith, Christa von Elling, Wayne Oppel and Corera Links. Pertaining to the field Acceleration The first der ivative (rate of change) work carried out in Argentina, we express our gratitude to Andrea Benvenuti, of an object’s velocity with respect to Fabian Gabelli, Monserrat Del Caño, La Chola, Miguel, Estancia El Pedral and Estancia San Lorenzo for assistance in various aspects of the research. We also time. Units are expressed as metres per thank the Instituto de Biología de Organismos Marinos (IBIOMAR‑ CONICET ) for second squared (m/s ) or in G-forces logistical support. HME is funded by an Irish Research Council Government of Ireland postgraduate scholarship. (g). A single G-force on Earth (though this does vary slightly with elevation) Authors’ contributions is 9.81 m/s . Tri-axial accelerometers The authors declare no conflict of interest and are in agreement to submit to Animal Biotelemetry. RMG conceived the study and RMG and RW wrote the measure acceleration in three orthogo- initial draft. PH constructed tag housings for all model species used. Data col‑ nal planes (surge—‘anterior–posterior’, lection for the lions was led by SF, DG, PV, LVS and AB and assisted by CJT, MFB, sway—‘medio-lateral’ and heave— DMS, SB, MVR, PH and RMG. Data collection for the penguins was led by FQ and data collection for the cormorants was led by AGL, with assistance from ‘dorsal–ventral’). Under non-moving KY, TY, and RMG. LB, MDH, RPW and RMG conceptualised the key considera‑ conditions, relative to gravity, the tions underlying the R code procedures and associated case‑studies, and RMG wrote the Gundog scripts and conducted the analysis of the case‑studies. device tilt (pitch and roll) can be cal- MHT supplied data for the Instantaneous tidal currents of the San Lorenzo culated directly from raw accelerom- region. All authors contributed to manuscript revision and LB, HJW, HME, AGL etery values since they are composed and LVS contributed to the testing and revision of the R syntax. All authors read and approved the final manuscript. entirely of the static force (gravity). Under linear acceleration, ‘moving’ Funding This research contributes to the CAASE project funded by King Abdullah forces applied to the device (e.g., due University of Science and Technology (KAUST ) under the KAUST Sensor to the animal moving) are superim- Initiative. Fieldwork in the Kgalagadi Transfrontier Park was supported in part posed to static readings and as such by a Department for Economy Global Challenges Research Fund grant to MS. Fieldwork within the Chubut Province was supported in part by the National measured animal acceleration is typi- Agency for Scientific and Technological Promotion of Argentina (PICT 2017‑ cally comprised of both a static and 1996 and PICT 2018‑1480), and the Grants‑in‑Aid for Scientific Research from the Japan Society for the Promotion of Science (16K18617). dynamic component. Barometric Pressure with the Earth’s atmosphere, Availability of data and materials pressure that is a measure of force per unit We provide a step‑by step example R script and an example data file of a Magellanic penguin walking out to sea (with time, raw acceleration and mag‑ area, often expressed as standard netometry data, marked events and aligned GPS positions) to demonstrate atmosphere (symbol: atm), defined as some of the key concepts outlined within “VPC dead‑reckoning procedure 101,325 Pa (1013.25 mbar; 1 Pa = 1 N/ in R” section when dead‑reckoning using Gundog.Tracks (including the initial calibration of magnetometry data with Gundog.Compass). All scripts, and the m ). The Earths mean sea-level atmos - example penguin data file have been uploaded to GitHub [95] and will be pheric pressure is approx. 1 atm. Baro- made available if the manuscript is accepted for publication. Online scripts will be continually updated and any queries, suggestions and/or reported bugs metric pressure decreases with eleva- should be emailed to the corresponding author. tion and increases with depth. Centripetal Inertial force caused by circular acceleration motion because an object is always Gunner et al. Anim Biotelemetry (2021) 9:23 Page 30 of 37 accelerating when either its direction (ECEF) system the Earth (this is often simplified to or magnitude (speed) changes, and in ‘Earth frame of reference’ or ‘Earths circular motion, the direction changes fixed frame’ in text). Its origin is fixed instantaneously. This can cause the at the Earth’s centre (the x-axis points animal to ‘pull g’, such as at times of towards the intersection of the Earth’s banking and cornering very fast. Greenwich Meridian and equatorial Coordinate In 3-D space, this is a set of three vectors plane, the y-axis pointing 90 degrees frame (x-, y-, z-axes) of unit length, perpen- East of the x-axis and the z pointing dicular (orthogonal) to each other. north, along the Earth’s rotation axis). Current flow (or ‘external’ current flow vectors). Note, this is different to the Earth- vectors The heading and speed of tidal-/ Centred Inertial (ECI) system, which air-currents. is non-rotating (and the x-axis instead Current Adding current flow vec tors to dead- always points towards the vernal integration reckoned travel vectors. equinox). De-rotation W ithin the tilt-compensated compass Equal pitch The animal moves in the same direction framework, this is the conversion of the assumption and angle as its anterior–posterior axis magnetic vector values through multi- (relative to North and the gravity vec- plying by the transpose (inverse) of the tor, respectively). pitch and then roll rotation matrices. Georeference W ithin the dead-reckoning frame- Distance The 2-D Haversine distance ratio work this is another term used for correction between successive Verified positions carrying out VPC-dead-reckoning, factor (VPs) (used in the VP-correction pro- or drift-correction or GPS-corrected cedure) and corresponding dead-reck- dead-reckoning. oned positions. This is multiplied to all Gimbal lock This is the loss of a degree of freedom intermediate (between VPs) radial dis- in 3-D, when two axes become paral- tance (q) values. lel to each other (locked in the same Drift The accumulation of spatial errors attitude, reflecting the same rotation). relative to a Verified Position, arising For example if the anterior–posterior from integrating incorrect dimensions axis (‘surge’ or ‘forward-back’—x-axis of travel. for NED coordinate frames) points in Dynamic The dynamic component of acceleration, the plane of the gravity vector (pitched body which is typically induced by the limb 90 degrees up or down), then the dor- acceleration and/or spine kinematics of the animal sal–ventral (‘heave’ or ‘up-down— (DBA) (and thus the attached accelerometer). z-axis for NED coordinate frames) and Generally, more mechanical work (via the medio-lateral axis (‘sway’ or ‘side- muscular contraction), corresponds to-side’—y-axis for NED coordinate to higher metabolic rate and greater frames) become parallel to each other, magnitudes in accelerometery read- and changes about the yaw can no ings (dependent on tag deployment longer be compensated for [changes site). Typically, dynamic values from in the roll (or ‘bank’) is equivalent to each multi-axial channel are integrated changing the heading]. into an overall metric, such as ‘Overall GPS-derived The Haversine distance calculated Dynamic Body Acceleration’ [ODBA = speed b etween successive GPS coordinates, │DBAx│ + │DBAy│ + │DBAz│] divided by the time taken between or ‘Vectorial Dynamic Body Accelera- locations. 2 2 tion’ [VeDBA = (DBAx + DBAy + D Ground- Empirical evidence (often information 2 0.5 BAz ) ]. Such derivatives have been truthing obtained by direct observation), as demonstrated as useful proxies for opposed to inference for validating movement-based power. something under investigation. Within Earth-Centre, This defines a non-inertial reference the dead-reckoning framework, VPs Earth-Fixed coordinate frame that rotates with such as GPS locations are used to Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 31 of 37 periodically ground-truth an animal’s system frame, the origin affixed as the devices position. centre of gravity and its axes oriented Haversine Computes the great-circular distance along the geodetic directions defined formula (units in metres) between two loca- by the Earth surface (the x- and y-axis tions (using their longitude and lati- pointing true north and East, respec- tude coordinates) on a sphere, apply- tively, parallel to the geoid surface ing trigonometry to map a triangle to and the z-axis pointing downwards the surface of a unit sphere. This for - towards the Earth’s surface). mula is only an approximation because Pitch The angle of device’s anterior–poste- the Earth is not a perfect sphere, with rior inclination or declination, relative numerical errors also arising at the to the horizontal plane of the Earth’s antipodal regions. surface. Pitch is often expressed as an Heading The difference of heading (or ‘bearing’) Euler angle, which describe the atti- correction factor f rom true North between consecutive tude and rotations of a device via a VPs (used in the VPC procedure) and given Euler angle sequence (yaw, pitch temporally aligned dead-reckoned posi- and roll) of rotations (using rotation tions. This is summed to all intermedi- matrices). Pitch can be derived from ate (between VPs) heading (h) values. the static component of accelera- Inf Results from numerical calculations tion. Assuming an NED system, pitch which are mathematically infinite (e.g., defines the degree of rotation about in R, dividing any value by zero results the y-axis. in Inf ). Radial Within the dead-reckoning framework, Linear drift At each path segment, the dead-reck- distance (q) this refers to a progression distance correction oned path is shifted to the position of accounting for the approximate cur- method the first VP encounter using a shift vature of the Earth (longlat projection vector. A correction vector then adds approximating the geoid to a sphere; the difference between the VP and radius (R) = 6,378,137 m). estimated dead-reckoned end points Right-handed The direction in which the linearly over this path segment period. coordinates/ fingers curl when pointing the right Multiplicative rotations thumb along the positive direction (m)-coefficient Within the dead-reckoning frame- (+ 1 g) of the z-axis (e.g., down for work, this refers to the gradient of a NED coordinates), reflect the direction linear regression, e.g., [speed =​ (Ve of rotation to be applied about each DBA· m) + c], where m is the multi- axis (for a given Euler angle sequence), plicative factor of VeDBA and c is the with the index finger representing the subsequently summed constant value x-axis and the middle finger repre- (reflecting the y-intercept). senting the y-axis, respectively, when NaN Non-n umeric (un-defined) values (e.g., splayed out at right angles to the in R, diving zero by zero results in thumb. NaN). Roll The angle of rotation about the device’s Net error Here, net error reflects the 2-D Haver - anterior–posterior axis. Roll is often sine distance (units in metres) between expressed as an Euler angle, which VPs and temporally aligned dead-reck- describe the attitude and rotations of a oned positions. device via a given Euler angle sequence Non-movement Behaviour performed while stationary, (yaw, pitch and roll) of rotations (using behaviours whereby the animal may be moving, rotation matrices). Roll is thus derived e.g., feeding on the spot, but there is after rotating by yaw and pitch and can no locomotion (not moving to a differ - be derived from the static component ent position in 3-D space). of acceleration. Assuming an NED sys- North–East– Often used in flight mechanics, tem, roll defines the degree of rotation Down (NED) this defines a non-inertial 3-D coordinate about the x-axis (also termed ‘bank Gunner et al. Anim Biotelemetry (2021) 9:23 Page 32 of 37 angle’). integration together. Assuming Cartesian coor- Static body The static component of dinates, vector addition is performed acceleration acceleration, due to gravity. Apart by adding the corresponding com- (SBA) from being used to calculate the angle ponents of the vectors together. E.g., of device tilt, increased inertial (cen- [A + B = (a + b , a + b , . . . , a + b )]. 1 1 2 2 n n tripetal) acceleration, e.g., when the Vertical speed Di stance travelled vertically up (at alti- animal ‘pulls g’, can be captured more tude) or down (at depth) divided by fully with static measures (rather than the time period between values. DBA estimates), and analogous to World Geodetic The typical model of the System 1984 VeDBA, the computation of the Vec- (WGS-84) E arth’s shape (standard for maps and torial Static Body Acceleration [VeSB satellite navigation), defining a coor - 2 2 2 0.5 A = (SBAx + SBAy + SBAz ) ] has dinate system that accounts for the been considered as a proxy of power. oblate spheroid. Tilt-compensated The compass heading compass Yaw The orientation of the device, gen- method (estimated using the arctangent ratio erally, with respect to true North between two orthogonal components (assuming any required magnetic dec- of the magnetic vector) is only accu- lination offset has been applied). Yaw, rate if the magnetometer outputs [typi- also termed ‘heading’ or ‘bearing’, is cally x, y channels—assuming the NED often expressed as an Euler angle, coordinate system is used (Additional which describe the attitude and rota- file  1: Text S2)] are taken when the tions of a device via a given Euler angle compass is level. Assuming the accel- sequence (yaw, pitch and roll) of rota- erometer-magnetometer approach, tions (using rotation matrices). Yaw static acceleration measures are used can be derived from the static com- to calculate the angles between the ponent of acceleration. Assuming an tag’s gravity (and thus magnetic) vec- NED system, yaw defines the degree of tor and the Earths frame of reference rotation about the z-axis. Yaw requires (e.g., Earth-Centered, Earth-Fixed the tilt-compensated compass method (ECEF) coordinate system). These to compute. angles are typically expressed as pitch Author details and roll Euler angles which are used to Swansea Lab for Animal Movement, Department of Biosciences, College of Science, Swansea University, Singleton Park, Swansea SA2 8PP, Wales, UK. compensate for variations in the mag- School of Biological Sciences, Queen’s University Belfast, Belfast, 19 Chlorine netometer output due to device tilt. Gardens, Belfast BT9 5DL, Northern Ireland, UK. Department of Agriculture, The tilt-compensated compass method Land Reform and Rural Development, Government of South Africa, Preto‑ ria 001, South Africa. Department of Migration, Max Planck Institute of Animal covers the procedures of adjusting the Behavior, 78315 Radolfzell, Germany. Department of Veterinary Tropical Dis‑ coordinate frame of the device to cor- eases, Faculty of Veterinary Science, University of Pretoria, Onderstepoort 0110, South Africa. School of Biology and Environmental Science, University College respond with a level inclination and Dublin, Belfield, Dublin, Ireland. Instituto de Biología de Organismos Marinos subsequently compute the compass (IBIOMAR), CONICET, Boulevard Brown 2915, U9120ACD Puerto Madryn, heading from the adjusted magnetom- Chubut, Argentina. Departamento de Ecología, Genética y Evolución & Insti‑ tuto de Ecología, Genética y Evolución de Buenos Aires (IEGEBA), CONICET, etry values. Pabellón II Ciudad Universitaria, C1428EGA Buenos Aires, Argentina. Centre Tortuosity The straight-line distance between for Biomathematics, College of Science, Swansea University, Swansea SA2 8PP, the start and end positions of a given UK. Graduate School of Environmental Studies, Nagoya University, Furo‑cho, Chikusa‑ku, Nagoya, Japan. Organization for the Strategic Coordination path segment, divided by the sum of of Research and Intellectual Properties, Meiji University, Nakano, Tokyo, Japan. the consecutive intermediate indi- Savanna and Grassland Research Unit, South African National Parks, Scientific Services Skukuza, Kruger National Park, Skukuza 1350, South Africa. Veteri‑ vidual distance steps that constituted nary Wildlife Services, South African National Parks, 97 Memorial Road, Old the total path segment’s length. Values Testing Grounds, Kimberley 8301, South Africa. Mammal Research Institute, closer to 0 (or conversely values closer Department of Zoology and Entomology, University of Pretoria, Pretoria 002, South Africa. Instituto Andino Patagónico de Tecnologías Biológicas y to 1 if subtracting the resultant ‘tortu- Geoambientales, Grupo GEA, IPATEC‑UNCO ‑ CONICET, San Carlos de Bariloche, osity’ value from 1) reflect more twists Río Negro, Argentina. Red Sea Research Centre, King Abdullah University of Science and Technology, Thuwal 23955, Saudi Arabia. Center for Zoo and turns in the movement path. and Wild Animal Health, Copenhagen Zoo, Roskildevej 38, 2000 Frederiksberg, Vector Adding vectors (of travel) Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 33 of 37 Denmark. Department of Zoology and Entomology, University of Fort Hare 20. Wilson RP, Shepard E, Liebsch N. Prying into the intimate details of Alice Campus, Ring Road, Alice 5700, South Africa. animal lives: use of a daily diary on animals. Endanger Species Res. 2008;4(1–2):123–37. Received: 8 March 2021 Accepted: 26 May 2021 21. Pedley M. eCompass‑build and calibrate a tilt ‑ compensating electronic compass. Circuit Cellar Mag Comput Appl. 2012;265:1–6. 22. Li Z, Li X, Wang Y. A calibration method for magnetic sensors and accelerometer in tilt‑ compensated digital compass. In: 2009 9th inter‑ national conference on electronic measurement & instruments, 16–19 Aug 2009; 2009. p. 2‑868‑862‑871. References 23. Ozyagcilar T. Implementing a tilt‑ compensated eCompass using 1. Browning E, Bolton M, Owen E, Shoji A, Guilford T, Freeman R. Predict‑ accelerometer and magnetometer sensors. Freescale semiconductor. ing animal behaviour using deep learning: GPS data alone accurately Application Note, AN4248; 2012. predict diving in seabirds. Methods Ecol Evol. 2018;9(3):681–92. 24. Gheorghe MV, Bodea MC. Calibration optimization study for tilt‑ 2. McClune DW. Joining the dots: reconstructing 3D environments compensated compasses. IEEE Trans Instrum Meas. 2018;67(6):1486–94. and movement paths using animal‑borne devices. Anim Biotelem. 25. Liu Y, Battaile BC, Trites AW, Zidek JV. Bias correction and uncertainty 2018;6(1):5. characterization of dead‑reckoned paths of marine mammals. Anim 3. Schlägel UE, Signer J, Herde A, Eden S, Jeltsch F, Eccard JA, Dammhahn Biotelem. 2015;3(1):51. M. Estimating interactions between individuals from concurrent animal 26. Dewhirst OP, Evans HK, Roskilly K, Harvey RJ, Hubel TY, Wilson AM. movements. Methods Ecol Evol. 2019;10(8):1234–45. Improving the accuracy of estimates of animal path and travel 4. Cagnacci F, Boitani L, Powell RA, Boyce MS. Animal ecology meets GPS‑ distance using GPS drift‑ corrected dead reckoning. Ecol Evol. based radiotelemetry: a perfect storm of opportunities and challenges. 2016;6(17):6210–22. Philos Trans R Soc B Biol Sci. 2010;365(1550):2157–62. 27. Mitani Y, Sato K, Ito S, Cameron MF, Siniff DB, Naito Y. A method for 5. Williams HJ, Taylor LA, Benhamou S, Bijleveld AI, Clay TA, de Grissac reconstructing three‑ dimensional dive profiles of marine mammals S, Demšar U, English HM, Franconi N, Gómez‑Laich A. Optimizing using geomagnetic intensity data: results from two lactating Weddell the use of biologgers for movement ecology research. J Anim Ecol. seals. Polar Biol. 2003;26(5):311–7. 2020;89(1):186–206. 28. Whitney NM, Pratt HL Jr, Pratt TC, Carrier JC. Identifying shark mating 6. Cotter CH. Early dead reckoning navigation. J Navig. 1978;31(1):20–8. behaviour using three‑ dimensional acceleration loggers. Endanger 7. Levi RW, Judd T. Dead reckoning navigational system using acceler‑ Species Res. 2010;10:71–82. ometer to measure foot impacts. In: Point Research Corporation, Santa 29. Lisovski S, Hewson CM, Klaassen RHG, Korner‑Nievergelt F, Kristensen Meijer, et al., Methods to assess physical activity with Ana, Calif., vol. MW, Hahn S. Geolocation by light: accuracy and precision affected by 5,583,776. 1996. p. 8. environmental factors. Methods Ecol Evol. 2012;3(3):603–12. 8. Beauregard S, Haas H. Pedestrian dead reckoning: a basis for personal 30. Miller PJO, Johnson MP, Madsen PT, Biassoni N, Quero M, Tyack PL. Using positioning. In: Proceedings of the 3rd workshop on positioning, navi‑ at‑sea experiments to study the effects of airguns on the foraging gation and communication; 2006. p. 27–35. behavior of sperm whales in the Gulf of Mexico. Deep Sea Res I Ocean‑ 9. Walker JS, Jones MW, Laramee RS, Holton MD, Shepard ELC, Williams ogr Res Pap. 2009;56(7):1168–81. HJ, Scantlebury DM, Marks NJ, Magowan EA, Maguire IE, Bidder OR, Di 31. Baumgartner MF, Freitag L, Partan J, Ball KR, Prada KE. Tracking large Virgilio A, Wilson RP. Prying into the intimate secrets of animal lives; soft‑ marine predators in three dimensions: the real‑time acoustic tracking ware beyond hardware for comprehensive annotation in ‘Daily Diary’ system. IEEE J Oceanic Eng. 2008;33(2):146–57. tags. Mov Ecol. 2015;3(1):29. 32. Williams LR, Fox DR, Bishop‑Hurley GJ, Swain DL. Use of radio frequency 10. Bidder OR, Walker JS, Jones MW, Holton MD, Urge P, Scantlebury DM, identification (RFID) technology to record grazing beef cattle water Marks NJ, Magowan EA, Maguire IE, Wilson RP. Step by step: reconstruc‑ point use. Comput Electron Agric. 2019;156:193–202. tion of terrestrial animal movement paths by dead‑reckoning. Mov Ecol. 33. Alexander JS, Zhang C, Shi K, Riordan P. A granular view of a snow 2015;3(1):23. leopard population using camera traps in Central China. Biol Conserv. 11. Wensveen PJ, Thomas L, Miller PJO. A path reconstruction method 2016;197:27–31. integrating dead‑reckoning and position fixes applied to humpback 34. English HM, Harvey L, Wilson RP, Gunner RM, Holton MD, Woodroffe R, whales. Mov Ecol. 2015;3(1):31. Börger L. Multi‑sensor biologgers and innovative training allow data 12. Andrzejaczek S, Gleiss AC, Lear KO, Pattiaratchi CB, Chapple TK, Meekan collection with high conservation and welfare value in zoos. https:// doi. MG. Biologging tags reveal links between fine ‑scale horizontal and org/ 10. 21203/ rs.3. rs‑ 562677/ v1. vertical movement behaviors in Tiger Sharks (Galeocerdo cuvier). Front 35. Fancy SG, Pank LF, Douglas DC, Curby CH, Garner GW. Satellite telem‑ Mar Sci. 2019;6:229. etry: a new tool for wildlife research and management. Washington, 13. Bowditch N. The American practical navigator, Bicentennial edition. D.C: Fish and Wildlife Service Washington DC; 1988. Bethesda: National Imagery and Mapping Agency; 2002. p. 879. 36. Soutullo A, Cadahía L, Urios V, Ferrer M, Negro JJ. Accuracy of light‑ 14. Wilson R, Wilson M‑P. Dead reckoning: a new technique for determining weight satellite telemetry: a case study in the Iberian Peninsula. J Wildl penguim movements at sea. Meeresforschung. 1988;32(2):155–8. Manag. 2007;71(3):1010–5. 15. Wilson RP, Wilson M‑PT, Link R, Mempel H, Adams NJ. Determination of 37. Catipovic JA. Performance limitations in underwater acoustic telemetry. movements of African penguins Spheniscus demersus using a compass IEEE J Oceanic Eng. 1990;15(3):205–16. system: dead reckoning may be an alternative to telemetry. J Exp Biol. 38. Hofman MPG, Hayward MW, Heim M, Marchand P, Rolandsen CM, 1991;157(1):557–64. Mattisson J, Urbano F, Heurich M, Mysterud A, Melzheimer J, Morellet 16. Wilson RP, Liebsch N, Davies IM, Quintana F, Weimerskirch H, Storch S, N, Voigt U, Allen BL, Gehr B, Rouco C, Ullmann W, Holand Ø, Jørgensen Lucke K, Siebert U, Zankl S, Müller G, Zimmer I, Scolaro A, Campagna C, NH, Steinheim G, Cagnacci F, Kroeschel M, Kaczensky P, Buuveibaatar Plötz J, Bornemann H, Teilmann J, McMahon CR. All at sea with animal B, Payne JC, Palmegiani I, Jerina K, Kjellander P, Johansson Ö, LaPoint tracks; methodological and analytical solutions for the resolution of S, Bayrakcismith R, Linnell JDC, Zaccaroni M, Jorge MLS, Oshima JEF, movement. Deep Sea Res II Top Stud Oceanogr. 2007;54(3):193–210. Songhurst A, Fischer C, Mc Bride RT Jr, Thompson JJ, Streif S, Sandfort R, 17. Mitani Y, Watanabe Y, Sato K, Cameron MF, Naito Y. 3D diving behavior Bonenfant C, Drouilly M, Klapproth M, Zinner D, Yarnell R, Stronza A, Wil‑ of Weddell seals with respect to prey accessibility and abundance. Mar mott L, Meisingset E, Thaker M, Vanak AT, Nicoloso S, Graeber R, Said S, Ecol Prog Ser. 2004;281:275–81. Boudreau MR, Devlin A, Hoogesteijn R, May‑ Junior JA, Nifong JC, Odden 18. Elkaim GH, Decker EB, Oliver G, Wright B. Marine mammal marker J, Quigley HB, Tortato F, Parker DM, Caso A, Perrine J, Tellaeche C, Zieba (MAMMARK) dead reckoning sensor for In‑Situ environmental monitor ‑ F, Zwijacz‑Kozica T, Appel CL, Axsom I, Bean WT, Cristescu B, Périquet ing. In: Proceedings of IEEE/ION PLANS 2006, San Diego, CA; 2006. p. S, Teichman KJ, Karpanty S, Licoppe A, Menges V, Black K, Scheppers 976–87. TL, Schai‑Braun SC, Azevedo FC, Lemos FG, Payne A, Swanepoel LH, 19. Wilson AD, Wikelski M, Wilson RP, Cooke SJ. Utility of biological sensor Weckworth BV, Berger A, Bertassoni A, McCulloch G, Šustr P, Athreya V, tags in animal conservation. Conserv Biol. 2015;29(4):1065–75. Gunner et al. Anim Biotelemetry (2021) 9:23 Page 34 of 37 Bockmuhl D, Casaer J, Ekori A, Melovski D, Richard‑Hansen C, van de 59. Goldbogen JA, Calambokidis J, Shadwick RE, Oleson EM, McDonald MA, Vyver D, Reyna‑Hurtado R, Robardet E, Selva N, Sergiel A, Farhadinia MS, Hildebrand JA. Kinematics of foraging dives and lunge‑feeding in fin Sunde P, Portas R, Ambarli H, Berzins R, Kappeler PM, Mann GK, Pyritz whales. J Exp Biol. 2006;209(7):1231–44. L, Bissett C, Grant T, Steinmetz R, Swedell L, Welch RJ, Armenteras D, 60. Zimmer WM, Tyack PL, Johnson MP, Madsen PT. Three‑ dimensional Bidder OR, González TM, Rosenblatt A, Kachel S, Balkenhol N. Right on beam pattern of regular sperm whale clicks confirms bent ‑horn track? Performance of satellite telemetry in terrestrial wildlife research. hypothesis. J Acoust Soc Am. 2005;117(3):1473–85. PLoS ONE. 2019;14(5):e0216223. 61. Wilson RP, Ropert‑ Coudert Y, Kato A. Rush and grab strategies in forag‑ 39. Newey S, Davidson P, Nazir S, Fairhurst G, Verdicchio F, Irvine RJ, van der ing marine endotherms: the case for haste in penguins. Anim Behav. Wal R. Limitations of recreational camera traps for wildlife manage‑ 2002;63(1):85–95. ment and conservation research: a practitioner’s perspective. Ambio. 62. Gabaldon J, Turner EL, Johnson‑Roberson M, Barton K, Johnson M, 2015;44(4):624–35. Anderson EJ, Shorter KA. Integration, calibration, and experimental 40. Leonardo M, Noss AJ, Erika C, Damián IR. Ocelot (Felis pardalis) popula‑ verification of a speed sensor for swimming animals. IEEE Sens J. tion densities, activity, and ranging behaviour in the dry forests of east‑ 2019;19(10):3616–25. ern Bolivia: data from camera trapping. J Trop Ecol. 2005;21(3):349–53. 63. Denny M. Air and water: the biology and physics of life’s media. Prince‑ 41. Lewis JS, Rachlow JL, Garton EO, Vierling LA. Eec ff ts of habitat on GPS ton University Press; 1993. collar performance: using data screening to reduce location error. J 64. Altynay K, Khan Mohammed A, Marengo M, Swanepoel L, Przybysz A, Appl Ecol. 2007;44(3):663–71. Muller C, Fahlman A, Buttner U, Geraldi NR, Wilson RP, Duarte CM, Kosel 42. Kaur M, Sandhu M, Mohan N, Sandhu PS. RFID technology principles, J. Wearable multifunctional printed graphene sensors. NPJ Flex Electron. advantages, limitations & its applications. Int J Comput Electr Eng. 2019;3(1):1–10. 2011;3(1):151. 65. Williams H, Shepard E, Duriez O, Lambertucci SA. Can accelerometry 43. Davis RW, Fuiman LA, Williams TM, Le Boeuf BJ. Three‑ dimensional be used to distinguish between flight types in soaring birds? Anim movements and swimming activity of a northern elephant seal. Comp Biotelem. 2015;3(1):1–11. Biochem Physiol A Mol Integr Physiol. 2001;129(4):759–70. 66. Williams HJ, Shepard E, Holton MD, Alarcón P, Wilson R, Lambertucci S. 44. Wilson RP. Movements in Adélie Penguins foraging for chicks at Ardley Physical limits of flight performance in the heaviest soaring bird. Proc Island, Antarctica: circles within spirals, wheels within wheels. Polar Natl Acad Sci. 2020;117(30):17884–90. Biosci. 2002;15:75–87. 67. Wilson RP, Börger L, Holton MD, Scantlebury DM, Gómez‑Laich A, 45. Johnson MP, Tyack PL. A digital acoustic recording tag for measuring Quintana F, Rosell F, Graf PM, Williams H, Gunner R, Hopkins L, Marks the response of wild marine mammals to sound. IEEE J Oceanic Eng. N, Geraldi NR, Duarte CM, Scott R, Strano MS, Robotka H, Eizaguirre 2003;28(1):3–12. C, Fahlman A, Shepard ELC. Estimates for energy expenditure in free‑ 46. Shiomi K, Sato K, Mitamura H, Arai N, Naito Y, Ponganis PJ. Eec ff t of living animals using acceleration proxies: a reappraisal. J Anim Ecol. ocean current on the dead‑reckoning estimation of 3‑D dive paths of 2020;89(1):161–72. emperor penguins. Aquat Biol. 2008;3(3):265–70. 68. Bidder OR, Qasem LA, Wilson RP. On higher ground: how well can 47. Laplanche C, Marques TA, Thomas L. Tracking marine mammals in 3D dynamic body acceleration determine speed in variable terrain? PLoS using electronic tag data. Methods Ecol Evol. 2015;6(9):987–96. ONE. 2012;7(11):e50556. 48. Adachi T, Costa DP, Robinson PW, Peterson SH, Yamamichi M, Naito 69. Bidder OR, Soresina M, Shepard ELC, Halsey LG, Quintana F, Gómez‑ Y, Takahashi A. Searching for prey in a three‑ dimensional environ‑ Laich A, Wilson RP. The need for speed: testing acceleration for estimat‑ ment: hierarchical movements enhance foraging success in northern ing animal travel rates in terrestrial dead‑reckoning systems. Zoology. elephant seals. Funct Ecol. 2017;31(2):361–9. 2012;115(1):58–64. 49. Bras YL, Jouma’a J, Guinet C. Three‑ dimensional space use during the 70. Markovska M, Svensson R. Evaluation of drift correction strategies for bottom phase of southern elephant seal dives. Mov Ecol. 2017;5(1):18. an inertial based dairy cow positioning system: a study on tracking the 50. Andrzejaczek S, Gleiss AC, Pattiaratchi CB, Meekan MG. First Insights position of dairy cows using a foot mounted IMU with drift correction Into the fine ‑scale movements of the Sandbar Shark, Carcharhinus from ZUPT or sparse RFID locations. Student thesis. 2019. plumbeus. Front Mar Sci. 2018;5:483. 71. di Virgilio A, Morales JM, Lambertucci SA, Shepard EL, Wilson RP. Multi‑ 51. Aoki K, Amano M, Mori K, Kourogi A, Kubodera T, Miyazaki N. Active dimensional precision livestock farming: a potential toolbox for sustain‑ hunting by deep‑ diving sperm whales: 3D dive profiles and maneuvers able rangeland management. PeerJ. 2018;6:e4867. during bursts of speed. Mar Ecol Prog Ser. 2012;444:289–301. 72. Shepard EL, Wilson RP, Halsey LG, Quintana F, Laich AG, Gleiss AC, Lieb‑ 52. Narazaki T, Sato K, Abernathy K, Marshall G, Miyazaki N. Sea turtles sch N, Myers AE, Norman B. Derivation of body motion via appropriate compensate deflection of heading at the sea surface during directional smoothing of acceleration data. Aquat Biol. 2008;4(3):235–41. travel. J Exp Biol. 2009;212(24):4019–26. 73. Noda T, Okuyama J, Koizumi T, Arai N, Kobayashi M. Monitoring attitude 53. Benoit‑Bird KJ, Battaile BC, Nordstrom CA, Trites AW. Foraging behavior and dynamic acceleration of free‑moving aquatic animals using a of northern fur seals closely matches the hierarchical patch scales of gyroscope. Aquat Biol. 2012;16(3):265–76. prey. Mar Ecol Prog Ser. 2013;479:283–302. 74. Wilson JW, Mills MG, Wilson RP, Peters G, Mills ME, Speakman JR, 54. Ware C, Friedlaender AS, Nowacek DP. Shallow and deep lunge feeding Durant SM, Bennett NC, Marks NJ, Scantlebury M. Cheetahs, Acinonyx of humpback whales in fjords of the West Antarctic Peninsula. Mar jubatus, balance turn capacity with pace when chasing prey. Biol Lett. Mamm Sci. 2011;27(3):587–605. 2013;9(5):20130620. 55. Wright BM, Ford JK, Ellis GM, Deecke VB, Shapiro AD, Battaile BC, Trites 75. McNarry MA, Wilson RP, Holton MD, Griffiths IW, Mackintosh KA. Inves‑ AW. Fine‑scale foraging movements by fish‑ eating killer whales (Orcinus tigating the relationship between energy expenditure, walking speed orca) relate to the vertical distributions and escape responses of salmo‑ and angle of turning in humans. PLoS ONE. 2017;12(8):e0182333. nid prey (Oncorhynchus spp.). Mov Ecol. 2017;5(1):1–18. 76. Kaniewski P, Kazubek J. Integrated system for heading determination. 56. Wensveen PJ, Isojunno S, Hansen RR, von Benda‑Beckmann AM, Klei‑ Acta Phys Pol Ser A Gen Phys. 2009;116(3):325. vane L, van IJsselmuide S, Lam FPA, Kvadsheim PH, DeRuiter SL, Curé C, 77. Noda T, Kawabata Y, Arai N, Mitamura H, Watanabe S. Animal‑mounted Narazaki T, Tyack PL, Miller PJO. Northern bottlenose whales in a pristine gyroscope/accelerometer/magnetometer: in situ measurement of the environment respond strongly to close and distant navy sonar signals. movement performance of fast‑start behaviour in fish. J Exp Mar Biol Proc R Soc B Biol Sci. 1899;2019(286):20182592. Ecol. 2014;451:55–68. 57. Ropert‑ Coudert Y, Wilson RP. Trends and perspectives in animal‑ 78. Patonis P, Patias P, Tziavos IN, Rossikopoulos D, Margaritis KG. A fusion attached remote sensing. Front Ecol Environ. 2005;3(8):437–44. method for combining low‑ cost IMU/magnetometer outputs for use in 58. Francisco FA, Nührenberg P, Jordan A. High‑resolution, non‑invasive applications on mobile devices. Sensors. 2018;18(8):2616. animal tracking and reconstruction of local environment in aquatic 79. Diebel J. Representing attitude: Euler angles, unit quaternions, and ecosystems. Mov Ecol. 2020;8(1):27. rotation vectors. Technical Report, vol. 58. Stanford: Stanford University; 2006. p. 1–35. Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 35 of 37 80. Han S, Wang J. A novel method to integrate IMU and magnetometers in DM, Duarte CM. Give the machine a hand: a Boolean time‑based attitude and heading reference systems. J Navig. 2011;64(4):727–38. decision‑tree template for rapidly finding animal behaviours in multi‑ 81. Säll J, Merkel J. Indoor navigation using accelerometer and magnetom‑ sensor data. Methods Ecol Evol. 2018;9(11):2206–15. eter. Student thesis. 2011. 102. Rong L, Zhiguo D, Jianzhong Z, Ming L. Identification of individual 82. Alaimo A, Artale V, Milazzo C, Ricciardello A. Comparison between walking patterns using gait acceleration. In: 2007 1st international con‑ Euler and quaternion parametrization in UAV dynamics. AIP Conf Proc. ference on bioinformatics and biomedical engineering, 6–8 July 2007; 2013;1558(1):1228–31. 2007. p. 543–6. 83. Valenti RG, Dryanovski I, Xiao J. Keeping a good attitude: a 103. Wilson RP, Hustler K, Ryan PG, Burger AE, Noldeke EC. Diving birds in quaternion‑based orientation filter for IMUs and MARGs. Sensors. cold water: do Archimedes and Boyle determine energetic costs? Am 2015;15(8):19302–30. Nat. 1992;140(2):179–200. 84. Feng K, Li J, Zhang X, Shen C, Bi Y, Zheng T, Liu J. A new quaternion‑ 104. Aoki K, Watanabe YY, Crocker DE, Robinson PW, Biuw M, Costa DP, based Kalman filter for real‑time attitude estimation using the two ‑step Miyazaki N, Fedak MA, Miller PJO. Northern elephant seals adjust gliding geometrically‑intuitive correction algorithm. Sensors. 2017;17(9):2146. and stroking patterns with changes in buoyancy: validation of at‑sea 85. Chiella ACB, Teixeira BOS, Pereira GAS. Quaternion‑based robust metrics of body density. J Exp Biol. 2011;214(17):2973–87. attitude estimation using an adaptive unscented Kalman filter. Sensors. 105. Williams HJ, Holton MD, Shepard ELC, Largey N, Norman B, Ryan PG, 2019;19(10):2372. Duriez O, Scantlebury M, Quintana F, Magowan EA, Marks NJ, Alagaili 86. Gunner RM, Wilson RP, Holton MD, Scott R, Hopkins P, Duarte CM. A new AN, Bennett NC, Wilson RP. Identification of animal movement patterns direction for differentiating animal activity based on measuring angular using tri‑axial magnetometry. Mov Ecol. 2017;5(1):6. velocity about the yaw axis. Ecol Evol. 2020;10(14):7872–86. 106. Whitford M, Klimley AP. An overview of behavioral, physiological, and 87. Wiltschko R. Magnetic orientation in animals, vol. 33. Berlin: Springer environmental sensors used in animal biotelemetry and biologging Science & Business Media; 2012. studies. Anim Biotelem. 2019;7(1):1–24. 88. Sequeira MM, Rickenbach M, Wietlisbach V, Tullen B, Schutz Y. 107. Shamoun‑Baranes J, Bom R, van Loon EE, Ens BJ, Oosterbeek K, Bouten Physical activity assessment using a pedometer and its comparison W. From sensor data to animal behaviour: an oystercatcher example. with a questionnaire in a large population survey. Am J Epidemiol. PLoS ONE. 2012;7(5):e37997. 1995;142(9):989–99. 108. Wilson R. The Jackass Penguin (Spheniscus demersus) as a pelagic preda‑ 89. Kerdok AE, Biewener AA, McMahon TA, Weyand PG, Herr HM. Energetics tor. Mar Ecol Prog Ser. 1985;25(3):219–27. and mechanics of human running on surfaces of different stiffnesses. J 109. Culik BM, Wilson RP, Dannfeld R, Adelung D, Spairani HJ, Coria NRC. Appl Physiol. 2002;92(2):469–78. Pygoscelid penguins in a swim canal. Polar Biol. 1991;11(4):277–82. 90. Halsey LG, Shepard ELC, Hulston CJ, Venables MC, White CR, Jeukend‑ 110. Bethge P, Nicol S, Culik BM, Wilson RP. Diving behaviour and rup AE, Wilson RP. Acceleration versus heart rate for estimating energy energetics in breeding little penguins (Eudyptula minor). J Zool. expenditure and speed during locomotion in animals: tests with an 1997;242(3):483–502. easy model species, Homo sapiens. Zoology. 2008;111(3):231–41. 111. Tobalske B, Dial K. Flight kinematics of black‑billed magpies and 91. Miwa M, Oishi K, Nakagawa Y, Maeno H, Anzai H, Kumagai H, Okano K, pigeons over a wide range of speeds. J Exp Biol. 1996;199(2):263–80. Tobioka H, Hirooka H. Application of overall dynamic body accelera‑ 112. Sato K, Sakamoto KQ, Watanuki Y, Takahashi A, Katsumata N, Bost C‑A, tion as a proxy for estimating the energy expenditure of grazing farm Weimerskirch H. Scaling of soaring seabirds and implications for flight animals: relationship with heart rate. PLoS ONE. 2015;10(6):e0128042. abilities of giant pterosaurs. PLoS ONE. 2009;4(4):e5400. 92. Chapman Jason W, Klaassen Raymond HG, Drake VA, Fossette S, Hays 113. Scharold J, Lai NC, Lowell W, Graham J. Metabolic rate, heart rate, and Graeme C, Metcalfe Julian D, Reynolds Andrew M, Reynolds Don R, tailbeat frequency during sustained swimming in the leopard shark Alerstam T. Animal orientation strategies for movement in flows. Curr Triakis semifasciata. Exp Biol. 1989;48(4):223–30. Biol. 2011;21(20):R861–70. 114. Kawabe R, Kawano T, Nakano N, Yamashita N, Hiraishi T, Naito Y. Simul‑ 93. R Development Core Team. R—a language and environment for statis‑ taneous measurement of swimming speed and tail beat activity of free‑ tical computing. Vienna: R Foundation for Statistical Computing; 2021. swimming rainbow trout Oncorhynchus mykiss using an acceleration 94. The R project for statistical computing. https:// www.r‑ proje ct. org/. data‑logger. Fish Sci. 2003;69(5):959–65. Accessed 07 Mar 2021. 115. Steinhausen MF, Steffensen JF, Andersen NG. Tail beat frequency as 95. Gundog.Tracks GitHub database. Available at https:// github. com/ Richa a predictor of swimming speed and oxygen consumption of saithe rd6195/ Dead‑ recko ning‑ animal‑ movem ents‑ in‑R. Accessed 09 June (Pollachius virens) and whiting (Merlangius merlangus) during forced 2021 swimming. Mar Biol. 2005;148(1):197–204. 96. Dowle M, Srinivasan A, Gorecki J, Chirico M, Stetsenko P, Short T, 116. Miller PJO, Johnson MP, Tyack PL, Terray EA. Swimming gaits, passive Lianoglou S, Antonyan E, Bonsch M, Parsonage H. Package ‘data.table’. drag and buoyancy of diving sperm whales Physeter macrocephalus. J Extension of ‘data frame; 2019. Exp Biol. 2004;207(11):1953–67. 97. Gunner RM, Wilson RP, Holton MD, Hopkins P, Bell SH, Marks NJ, Bennett 117. Laich AG, Wilson RP, Quintana F, Shepard EL. Identification of imperial NC, Ferreira S, Govender D, Viljoen P, Bruns A, Schalkwyk Lv, Bertelsen cormorant Phalacrocorax atriceps behaviour using accelerometers. MF, Duarte CM, Rooyen MCv, Tambling CJ, Goppert A, Diesel D, Scantle‑ Endanger Species Res. 2008;10:29–37. bury DM. Decision rules for determining terrestrial movement and the 118. Chakravarty P, Cozzi G, Ozgul A, Aminian K. A novel biomechanical consequences for filtering high‑resolution GPS tracks—a case study approach for animal behaviour recognition using accelerometers. using the African Lion (Panthera leo). Mov Ecol. (in review). Methods Ecol Evol. 2019;10(6):802–14. 98. Dunford CE, Marks NJ, Wilmers CC, Bryce CM, Nickel B, Wolfe LL, 119. Chakravarty P, Maalberg M, Cozzi G, Ozgul A, Aminian K. Behavioural Scantlebury DM, Williams TM. Surviving in steep terrain: a lab‑to ‑field compass: animal behaviour recognition using magnetometers. Mov assessment of locomotor costs for wild mountain lions (Puma concolor). Ecol. 2019;7(1):28. Mov Ecol. 2020;8(1):34. 120. Hochscheid S. Why we mind sea turtles’ underwater business: a review 99. Wilson RP, Rose KA, Gunner R, Holton M, Marks NJ, Bennett NC, Bell SH, on the study of diving behavior. J Exp Mar Biol Ecol. 2014;450:118–36. Twining JP, Hesketh J, Duarte CM, Bezodis N, Scantlebury DM. Forces 121. Constandache I, Bao X, Azizyan M, Choudhury RR. Did you see Bob? experienced by instrumented animals depend on lifestyle. bioRxiv. human localization using mobile phones. In: Proceedings of the 2020. https:// doi. org/ 10. 1101/ 2020. 08. 20. 258756. sixteenth annual international conference on Mobile computing and 100. Willener AST, Handrich Y, Halsey LG, Strike S. Eec ff t of walking speed networking; Chicago, Illinois, USA. Association for Computing Machin‑ on the gait of king penguins: an accelerometric approach. J Theor Biol. ery; 2010. p. 149–60. 2015;387:166–73. 122. Symington A, Trigoni N. Encounter based sensor tracking. In: Proceed‑ 101. Wilson RP, Holton MD, di Virgilio A, Williams H, Shepard ELC, Lamber‑ ings of the thirteenth ACM international symposium on mobile ad hoc tucci S, Quintana F, Sala JE, Balaji B, Lee ES, Srivastava M, Scantlebury Gunner et al. Anim Biotelemetry (2021) 9:23 Page 36 of 37 networking and computing, Hilton Head, South Carolina, USA. Associa‑ 146. Wilson R, Griffiths I, Legg P, Friswell M, Bidder O, Halsey L, Lambertucci tion for Computing Machinery; 2012. p. 15–24. SA, Shepard E. Turn costs change the value of animal search paths. Ecol 123. Harja YD, Sarno R. Determine the best option for nearest medical Lett. 2013;16(9):1145–50. services using Google maps API, Haversine and TOPSIS algorithm. In: 147. Potts J, Börger L, Scantlebury DM, Bennett NC, Alagaili A, Wilson RP. 2018 international conference on information and communications Finding turning‑points in ultra‑high‑resolution animal movement data. technology (ICOIACT ), 6–7 March 2018; 2018. p. 814–9. Methods Ecol Evol. 2018;9(10):2091–101. 124. Great circle distances and bearings between two locations. https:// 148. Narazaki T, Nakamura I, Aoki K, Iwata T, Shiomi K, Luschi P, Suganuma dtcen ter. org/ met/ users/ docs/ write_ ups/ gc_ simple. pdf. H, Meyer CG, Matsumoto R, Bost CA, Handrich Y, Amano M, Okamoto 125. Robusto CC. The Cosine‑Haversine formula. Am Math Mon. R, Mori K, Ciccione S, Bourjea J, Sato K. Similar circling movements 1957;64(1):38–40. observed across marine megafauna taxa. iScience. 2021;24:102221. 126. Zeileis A, Grothendieck G, Ryan JA, Andrews F, Zeileis MA. Package ‘zoo’. 149. Kunčar A, Sysel M, Urbánek T. Calibration of low‑ cost triaxial magnetom‑ 2020. eter. In: MATEC web of conferences 20th international conference on 127. Jockers ML, Thalken R. Introduction to dplyr. In: Text analysis with r for circuits, systems, communications and computers (CSCC 2016), 2016. students of literature. Cham: Springer; 2020. p. 121–32. EDP Sciences. 128. Calculating Azimuth, distance, and altitude from a pair of GPS locations 150. Sodhi R, Prunty J, Hsu G, Oh B. Automatic calibration of a three‑axis in JavaScript. https:// medium. com/ javas cript‑ in‑ plain‑ engli sh/ calcu lat ‑ magnetic compass. U.S. Patent 7,451,549 PNI Corporation (Santa Rosa, ing‑ azimu th‑ dista nce‑ and‑ altit ude‑ from‑a‑ pair‑ of‑ gps‑ locat ions‑ 36b43 CA, US); 2008. 25d8a b0. 151. Renaudin V, Afzal MH, Lachapelle G. Complete triaxis magnetometer 129. Handcock RN, Swain DL, Bishop‑Hurley GJ, Patison KP, Wark T, Valencia calibration in the magnetic domain. J Sens. 2010;2010:967245. P, Corke P, O’Neill CJ. Monitoring Animal behaviour and environmental 152. Tabatabaei SAH, Gluhak A, Tafazolli R. A fast calibration method for interactions using wireless sensor networks, GPS collars and satellite triaxial magnetometers. IEEE Trans Instrum Meas. 2013;62(11):2929–37. remote sensing. Sensors. 2009;9(5):3586–603. 153. Vitali A. Ellipsoid or sphere fitting for sensor calibration, Dt0059. ST 130. Gužvica G, Bošnjak I, Bielen A, Babić D, Radanović‑ Gužvica B, Šver L. Microelectronics, Design Tip; 2016. Comparative analysis of three different methods for monitoring the use 154. Simple and effective magnetometer calibration. https:// github. com/ of green bridges by wildlife. PLoS ONE. 2014;9(8):e106194.krisw iner/ MPU60 50/ wiki/ Simple‑ and‑ Eecff tive‑ Magne tomet er‑ Calib 131. Patel A, Stocks B, Fisher C, Nicolls F, Boje E. Tracking the cheetah ration. tail using animal‑borne cameras, GPS, and an IMU. IEEE Sens Lett. 155. Premerlani W, Bizard P. Direction cosine matrix imu: theory. Diy Drone: 2017;1(4):1–4. Usa; 2009:13–15. 132. Latham ADM, Latham MC, Anderson DP, Cruz J, Herries D, Hebblewhite 156. Grygorenko V. Sensing‑magnetic compass with tilt compensation. M. The GPS craze: six questions to address before deciding to deploy Cypress perform, semiconductor. Application Notes, AN2272. 2011. GPS technology on wildlife. N Z J Ecol. 2015;39(1):143–52. 157. Gleiss AC, Wilson RP, Shepard ELC. Making overall dynamic body 133. Hebblewhite M, Haydon DT. Distinguishing technology from biology: a acceleration work: on the theory of acceleration as a proxy for energy critical review of the use of GPS telemetry data in ecology. Philos Trans expenditure. Methods Ecol Evol. 2011;2(1):23–33. R Soc B Biol Sci. 2010;365(1550):2303–12. 158. Choi S, Youn I‑H, LeMay R, Burns S, Youn J‑H. Biometric gait recognition 134. Gibb R, Shoji A, Fayet AL, Perrins CM, Guilford T, Freeman R. Remotely based on wireless acceleration sensor using k‑nearest neighbor clas‑ sensed wind speed predicts soaring behaviour in a wide‑ranging sification. In: 2014 international conference on computing, networking pelagic seabird. J R Soc Interface. 2017;14(132):20170262. and communications (ICNC), 3–6 Feb 2014; 2014. p. 1091–5. 135. Swain DL, Wark T, Bishop‑Hurley GJ. Using high fix rate GPS data to 159. Sato K, Mitani Y, Cameron MF, Siniff DB, Naito Y. Factors affecting strok ‑ determine the relationships between fix rate, prediction errors and ing patterns and body angle in diving Weddell seals under natural patch selection. Ecol Model. 2008;212(3):273–9. conditions. J Exp Biol. 2003;206(9):1461–70. 136. Ryan PG, Petersen SL, Peters G, Grémillet D. GPS tracking a marine 160. Li W, Wang J. Eec ff tive adaptive Kalman filter for MEMS‑IMU/mag‑ predator: the effects of precision, resolution and sampling rate on netometers integrated attitude and heading reference systems. J Navig. foraging tracks of African Penguins. Mar Biol. 2004;145(2):215–23. 2013;66(1):99–113. 137. Bouvet D, Garcia G. GPS latency identification by Kalman filtering. 161. Pedley M. Tilt sensing using a three‑axis accelerometer. Freescale semi‑ Robotica. 2000;18(5):475–85. conductor. Application Note 2461. 2013;1(Rev. 6):1–22. 138. Frair JL, Fieberg J, Hebblewhite M, Cagnacci F, DeCesare NJ, Pedrotti L. 162. Shiomi K, Narazaki T, Sato K, Shimatani K, Arai N, Ponganis PJ, Resolving issues of imprecise and habitat‑biased locations in ecologi‑ Miyazaki N. Data‑processing artefacts in three ‑ dimensional dive path cal analyses using GPS telemetry data. Philos Trans R Soc B Biol Sci. reconstruction from geomagnetic and acceleration data. Aquat Biol. 2010;365(1550):2187–200. 2010;8(3):299–304. 139. Péron G, Calabrese JM, Duriez O, Fleming CH, García‑ Jiménez R, 163. Evans P. Rotations and rotation matrices. Acta Crystallogr Sect D. Johnston A, Lambertucci SA, Safi K, Shepard ELC. The challenges of 2001;57(10):1355–9. estimating the distribution of flight heights from telemetry or altimetry 164. Wilson RP, Holton MD, Walker JS, Shepard ELC, Scantlebury DM, Wilson data. Anim Biotelem. 2020;8(1):5. VL, Wilson GI, Tysse B, Gravenor M, Ciancio J, McNarry MA, Mackintosh 140. Poulin M‑P, Clermont J, Berteaux D. Extensive daily movement rates KA, Qasem L, Rosell F, Graf PM, Quintana F, Gomez‑Laich A, Sala J‑E, measured in territorial arctic foxes. Ecol Evol. 2021;00:1–12. Mulvenna CC, Marks NJ, Jones MW. A spherical‑plot solution to linking 141. Marcus Rowcliffe J, Carbone C, Kays R, Kranstauber B, Jansen PA. Bias acceleration metrics with animal performance, state, behaviour and in estimating animal travel distance: the effect of sampling frequency. lifestyle. Mov Ecol. 2016;4(1):22. Methods Ecol Evol. 2012;3(4):653–62. 165. Caruso MJ. Applications of magnetic sensors for low cost compass 142. Belant JL. Eec ff ts of antenna orientation and vegetation on global systems. In: IEEE 2000 position location and navigation symposium (Cat positioning system telemetry collar performance. Northeast Nat. No00CH37062): 13–16 March 2000; 2000. p. 177–84. 2009;16(4):577–84, 578. 166. Boelter KDH. Aircraft magnetic declinator system, vol. 2020/0182619. 143. Graf PM, Wilson RP, Qasem L, Hackländer K, Rosell F. The use of accel‑ Chicago: The Boeing Company; 2020. eration to code for animal behaviours; a case study in free‑ranging 167. World magnetic model 2020 calculator. http:// www. geomag. bgs. ac. Eurasian beavers castor fiber. PLoS ONE. 2015;10(8):e0136751.uk/ data_ servi ce/ models_ compa ss/ wmm_ calc. html. Accessed 07 Mar 144. Tonini MH, Palma ED. Tidal dynamics on the North Patagonian Argen‑ 2021. tinean Gulfs. Estuar Coast Shelf Sci. 2017;189:115–30. 168. Qasem L, Cardew A, Wilson A, Griffiths I, Halsey LG, Shepard ELC, Gleiss 145. Wilson RP, Kreye JM, Lucke K, Urquhart H. Antennae on transmitters AC, Wilson R. Tri‑axial dynamic acceleration as a proxy for animal energy on penguins: balancing energy budgets on the high wire. J Exp Biol. expenditure; should we be summing values or calculating the vector? 2004;207(15):2649–62. PLoS ONE. 2012;7(2):e31187. Gunner  et al. Anim Biotelemetry (2021) 9:23 Page 37 of 37 169. Pewsey A, Neuhäuser M, Ruxton GD. Circular statistics in R. OUP Oxford; 172. Movable type scripts. https:// www. movab le‑ type. co. uk/ scrip ts/ latlo ng. 2013. html. Accessed 07 Mar 2021. 170. Farrahi V, Niemelä M, Kangas M, Korpelainen R, Jämsä T. Calibration 173. Grolemund G, Wickham H. Dates and times made easy with lubridate. J and validation of accelerometer‑based activity monitors: a systematic Stat Softw. 2011;40(i03):1–25. review of machine‑learning approaches. Gait Posture. 2019;68:285–99. 171. Ropert‑ Coudert Y, Kato A, Baudat J, Bost C‑A, Le Maho Y, Naito Y. Time/ Publisher’s Note depth usage of Adélie penguins: an approach based on dive angles. Springer Nature remains neutral with regard to jurisdictional claims in pub‑ Polar Biol. 2001;24(6):467–70. lished maps and institutional affiliations. Re Read ady y to to submit y submit your our re researc search h ? Choose BMC and benefit fr ? Choose BMC and benefit from om: : fast, convenient online submission thorough peer review by experienced researchers in your field rapid publication on acceptance support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions

Journal

Animal BiotelemetrySpringer Journals

Published: Jul 1, 2021

Keywords: Animal behaviour; Animal movement; Global Positioning System; R (programming language); Track integration; Tri-axial accelerometers; Tri-axial magnetometers

References