-
Notifications
You must be signed in to change notification settings - Fork 0
/
fine-scale-tracking.Rmd
444 lines (324 loc) · 19.9 KB
/
fine-scale-tracking.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
---
title: "Fine-scale tracking"
author: "tagtools project team"
date: "`r Sys.Date()`"
output:
html_document:
code_folding: hide
vignette: >
%\VignetteIndexEntry{fine-scale-tracking}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
require(tagtools)
```
Welcome to the `fine-scale-tracking` vignette! Thanks for taking some time to get to know our package; we hope it sparks joy.
This practical covers how to compute a dead-reckoned track from accelerometer and magnetometer data and compares this track estimate to GPS positions that were taken by the tag. This illustrates both the high resolution of dead-reckoned tracks and the **very** large errors that can accumulate. You will then explore how to merge the GPS and dead-reckoned tracks to produce a high-resolution track with lower error and try visualizing the track coloured by another parameter.
*Estimated time for this vignette: 40 minutes.*
These practicals all assume that you have R/Rstudio installed on your machine, some basic experience working with them, and can execute provided code, making some user-specific changes along the way (e.g. to help R find a file you downloaded). We will provide you with quite a few lines. To boost your own learning, you would do well to try and write them before opening what we give, using this just to check your work.
Additionally, be careful when copy-pasting special characters such as \_underscores\_ and 'quotes'. If you get an error, one thing to check is that you have just single, simple underscores, and `'straight quotes'`, whether `'single'` or `"double"` (rather than “smart quotes”).
# Load and Explore Data
If you want to run this example, download the "testset7.nc" file from https://github.com/animaltags/tagtools_data and change the file path to match where you've saved the files
```{r, echo = FALSE, eval = TRUE}
d7_file_path <- "nc_files/testset7.nc"
d7 <- load_nc(d7_file_path)
```
Load the data set `testset7.nc`. This dataset is built into the `tagtools` package, so you can access it using `system.file`.
Then, check the contents of the data structures to answer the following questions. (See the code below if and only if you need a hint.)
- What species did the data come from, and where did the data come from?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_quickie1"> _Show/Hide Answers_ </button>
<div id="ans_quickie1" class="collapse">
*A humpback whale, *Megaptera novaeangliae*.*
</div>
- What is the sampling rate of the accelerometer data?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_quickie2"> _Show/Hide Answers_ </button>
<div id="ans_quickie2" class="collapse">
*1 Hz.*
</div>
- What processing steps have already been applied to the magnetometer data?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_quickie3"> _Show/Hide Answers_ </button>
<div id="ans_quickie3" class="collapse">
*sens_struct, do_cal, tag2animal, crop_to, decdc(5)... So, it has been converted to a sensor structure, calibrated, put into the animal frame, cropped to a certain segment of data, and decimated by a factor of five (as a result, the sampling rate is now the same as the acceleration data, 1 Hz).*
</div>
- What is in the 3 columns in the POS (GPS position) data?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_quickie4"> _Show/Hide Answers_ </button>
<div id="ans_quickie4" class="collapse">
*The three columns are the time, the latitude, and the longitude.*
</div>
- Which frame is the accelerometer data in? Which frame is the magnetometer data in?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_quickie5"> _Show/Hide Answers_ </button>
<div id="ans_quickie5" class="collapse">
*The accelerometer data is oriented in the animal frame, and so is the magnetometer data.*
</div>
```{r, eval = TRUE, echo = TRUE}
str(d7$info, max.level = 1)
str(d7$Aa, max.level = 1)
str(d7$POS, max.level = 1)
```
Plot the pressure and accelerometer data with `plott()` to get a sense for what the animal might be doing in this data segment. Note that the code example below assumes you have called the data set `d7`:
```{r, eval = FALSE, echo = TRUE}
library(tagtools)
d7_file_path <- ""
d7 <- load_nc(d7_file_path)
sampling_rate <- d7$P$sampling_rate
plott(X = list(Depth = d7$P, Acceleration = d7$Aa), fsx = sampling_rate, interactive = TRUE)
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#plott_d7"> Show/Hide Results </button>
<div id="plott_d7" class="collapse">
```{r, eval = TRUE, echo = FALSE}
sampling_rate <- d7$P$sampling_rate
plott(X = list(Depth = d7$P, Acceleration = d7$Aa), fsx = sampling_rate, interactive = FALSE)
```
</div>
Then, plot the GPS positions:
```{r, echo = TRUE, eval = FALSE}
plot(d7$POS$data[,3], d7$POS$data[,2], type = 'b', pch = 20)
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#plott_d7_POS"> Show/Hide Results </button>
<div id="plott_d7_POS" class="collapse">
```{r, echo = FALSE, eval = TRUE}
plot(d7$POS$data[,3], d7$POS$data[,2], type = 'b', pch = 20)
```
</div>
# Dead-reckoning
## What is Dead-reckoning?
Hearken back to your earliest trigonometry course. Chances are, you used angles and distances to calculate other distances (on a right triangle, especially). Then, in calculus, you've probably used angles and rates of change to calculate distance or other rates of change. Dead-reckoning is essentially this---using just a heading and a forward speed in order to plot a forward path. The technique is fairly commonplace in sailing and aviation, and can be made accurate with the help of known position data. If you want more background, consult [this video](https://www.youtube.com/watch?v=doZ1k1Oeeh8) for an example with some simple trigonometry.
However, per [Wikipedia](https://en.wikipedia.org/wiki/Dead_reckoning#Errors) (which simply says it very well), dead-reckoning "is subject to significant errors of approximation. For precise positional information, both speed and direction must be accurately known at all times during travel. Most notably, dead reckoning does not account for directional drift during travel through a fluid medium. These errors tend to compound themselves over greater distances". We will see how this is true over the course of this vignette---water is seldom entirely still, and there is quite a current present in the data we will soon investigate.
## Why use Dead-reckoning?
The plot shows a mix of intensive and extensive movements, but the constraint of only getting positions when the animal is at the surface means we cannot infer much about the movement behaviour within dives. Dead-reckoning from accelerometer and magnetometer data is the only way to estimate movement within dives without requiring external tracking infrastructure.
## Estimating Animal Speed
Dead-reckoning uses the accelerometer and magnetometer to calculate the direction of travel. In what frame do these data need to be?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_frame_to_calc"> _Show/Hide Answers_ </button>
<div id="ans_frame_to_calc" class="collapse">
*They both need to be in the animal frame. Thankfully, they both already are!*
</div>
An estimate of the forward speed is also required. We don't have a speed sensor---doing this would actually be quite cumbersome, as compared to an acceleration sensor. Nevertheless, we can compute the vertical speed (i.e., the differential of the depth) during descents and ascents, which might be a good starting guess:
```{r, echo = TRUE, eval = FALSE}
v <- depth_rate(d7$P)
plott(X = list(Depth = d7$P$data, `Vertical Speed` = v),
fsx = sampling_rate)
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#plott_v"> Show/Hide Results </button>
<div id="plott_v" class="collapse">
```{r, echo = FALSE, eval = TRUE}
v <- depth_rate(d7$P)
plott(X = list(Depth = d7$P$data, `Vertical Speed` = v),
fsx = sampling_rate)
```
</div>
Set `interactive = TRUE` and zoom in to individual dives on this plot to get a rough idea of the descent and ascent speed of the whale.
Or, if the interactive figure drives you a bit crazy, just set the axis limits:
```{r, echo = TRUE, eval = FALSE}
plott(X = list(Depth = d7$P$data, `Vertical Speed` = v),
fsx = sampling_rate, xl = c(0,0.1))
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#plott_v2"> Show/Hide Results </button>
<div id="plott_v2" class="collapse">
```{r, echo = FALSE, eval = TRUE}
plott(X = list(Depth = d7$P$data, `Vertical Speed` = v),
fsx = sampling_rate, xl = c(0,0.1))
```
</div>
## Computing a Dead-Reckoned Track
Call your speed estimate speed `spd` and use it in the following line to compute the dead-reckoned track:
```{r, echo = TRUE, eval = FALSE}
spd <- #???? (fill in your estimate here)
DR <- ptrack(A = d7$Aa, M = d7$Ma, s = spd)
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track')
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#plott_deadreckon1"> Show/Hide Results </button>
<div id="plott_deadreckon1" class="collapse">
```{r, echo = FALSE, eval = TRUE}
spd <- 1.7
DR <- ptrack(A = d7$Aa, M = d7$Ma, s = spd)
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track')
```
</div>
A dead-reckoned track is a series of distances north (or south if negative) and east (or west if negative) of the starting point which is position (0, 0) in the plot. The first two columns of DR are these 'northing' and 'easting' values. The DR track is defined in a Local Level Frame (LLF) as opposed to latitude and longitude. An LLF is like a map---a region that is small enough that we can assume the earth is flat---centered on the starting point.
**How does the spatial resolution of the dead-reckoned track look compared to the surface GPS positions?**
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_spatial_rez"> _Show/Hide Answers_ </button>
<div id="ans_spatial_rez" class="collapse">
*At first glance, the shape is quite similar, so this might not concern us yet. Also, the spatial resolution looks much, much better! So, making this track seems like a good idea.*
</div>
## Adding GPS positions
To plot the GPS positions on the same plot, we need to first convert them into the same LLF. The first GPS position is only 0.8s into the dataset (this is `\$d7\$POS\$data[1,1]`) so we can say that the dead-reckoned track starts from this point. To convert latitude and longitude into the LLF use:
```{r, echo = TRUE, eval = FALSE}
POSLLF <- lalo2llf(trk = d7$POS$data[,c(2:3)])
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track',
yl = c(-1000, 5000))
lines(POSLLF$easting, POSLLF$northing, type = 'b', col = 'red', pch = 20)
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#plott_posllf"> Show/Hide Results </button>
<div id="plott_posllf" class="collapse">
```{r, echo = FALSE, eval = TRUE}
POSLLF <- lalo2llf(trk = d7$POS$data[,c(2:3)])
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track',
yl = c(-1000, 5000))
lines(POSLLF$easting, POSLLF$northing, type = 'b', col = 'red', pch = 20)
```
</div>
How well does the dead-reckoned track match up to the GPS positions?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_NOTMATCHINGATALL"> _Show/Hide Answers_ </button>
<div id="ans_NOTMATCHINGATALL" class="collapse">
*They do not match up much at all. They touch at the beginning... and that's about it. After that, the shapes are somewhat preserved, but the predicted track is way farther north, and quite a ways farther east, than the actual known positions.*
</div>
The dead-reckoned track is always computed with respect to the water, not the ground, but we are plotting it here with respect to the ground. A more accurate track would take into account the water current. Can you imagine what current direction would be needed to make the dead-reckoned track more closely match the GPS positions?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_SOUTHCURRENT"> _Show/Hide Answers_ </button>
<div id="ans_SOUTHCURRENT" class="collapse">
*The current would need to be moving quite a bit south, and also some distance west, in order to compensate for the large errors we are seeing.*
</div>
There are a number of ways to combine the GPS positions and the dead-reckoned track into a single track which has both the absolute accuracy of GPS and the high temporal resolution of dead-reckoning. A simple method is to force the dead-reckoned track to meet the GPS positions at each surfacing by adding a constant 'current' to each track point in the preceding dive. This is done by `fit_tracks`:
```{r, eval = FALSE, echo = TRUE}
FT <- fit_tracks(P = POSLLF, times = d7$POS$data[,1],
D = DR[,c(1:2)],
sampling_rate = d7$Aa$sampling_rate)
# add to plot
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track',
yl = c(-1000, 5000))
lines(POSLLF$easting, POSLLF$northing, type = 'b', col = 'red', pch = 20)
lines(FT$easting, FT$northing, col = 'darkgreen')
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#fit_tracks"> Show/Hide Results </button>
<div id="fit_tracks" class="collapse">
```{r, eval = TRUE, echo = FALSE}
FT <- fit_tracks(P = POSLLF, times = d7$POS$data[,1],
D = DR[,c(1:2)],
sampling_rate = d7$Aa$sampling_rate)
# add to plot
plot(DR$easting, DR$northing, type = 'l',
xlab = 'Easting, m', ylab = 'Northing, m',
main = 'Dead-Reckoned Track',
yl = c(-1000, 5000))
lines(POSLLF$easting, POSLLF$northing, type = 'b', col = 'red', pch = 20)
lines(FT$easting, FT$northing, col = 'darkgreen')
```
</div>
## A digression about graphics
**Skip to the next section if you're not interested in improving the plots.**
If you are interested in a nicer, zoomable plot and either have, or are willing to install, packages `ggformula` and `plotly`, give the following a try!
```{r, echo = TRUE, eval = FALSE}
# if you need to install:
# install.packages(pkgs = c('ggformula', 'plotly'))
library(ggformula)
library(plotly)
theme_set(theme_bw(base_size = 12))
track_fig <- gf_path(northing ~ easting, data = DR,
xlab = 'Easting, m', ylab = 'Northing, m') %>%
gf_path(northing ~ easting, data = POSLLF, color = 'darkred') %>%
gf_point(northing ~ easting, data = POSLLF, color = 'darkred') %>%
gf_path(northing ~ easting, data = FT, color = 'darkgreen')
track_fig
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#track_fig"> Show/Hide Results </button>
<div id="track_fig" class="collapse">
```{r, eval = TRUE, echo = FALSE}
library(ggformula)
library(plotly)
theme_set(theme_bw(base_size = 12))
track_fig <- gf_path(northing ~ easting, data = DR,
xlab = 'Easting, m', ylab = 'Northing, m') %>%
gf_path(northing ~ easting, data = POSLLF, color = 'darkred') %>%
gf_point(northing ~ easting, data = POSLLF, color = 'darkred') %>%
gf_path(northing ~ easting, data = FT, color = 'darkgreen')
track_fig
```
</div>
OK, that's the track figure. What about zooming and interaction?
```{r, echo = TRUE, eval = FALSE}
track_fig %>% ggplotly()
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#track_fig_ggplotly"> Show/Hide Results </button>
<div id="track_fig_ggplotly" class="collapse">
```{r, echo = FALSE, eval = TRUE}
track_fig %>% ggplotly()
```
</div>
Ahhh, so much nicer! But do note that `plotly()` figures render only in interactive R sessions or, if using Rmarkdown, in html output (not PDF).
This figure is also nicer in that we can update it by chaining. For example, to change the axis limits:
```{r, echo = TRUE, eval = FALSE}
track_fig_zoom <- track_fig %>%
gf_lims(x = c(-500, 1500), y = c(-1100, 100))
track_fig_zoom
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#track_fig_zoom"> Show/Hide Results </button>
<div id="track_fig_zoom" class="collapse">
```{r, echo = FALSE, eval = TRUE}
track_fig_zoom <- track_fig %>%
gf_lims(x = c(-500, 1500), y = c(-1100, 100))
track_fig_zoom
```
</div>
## Interpreting the Tracks
**Now, back to the main tutorial...**
Examine the plot to see how the green fitted track interpolates the red GPS positions. If the green track is to be believed, how effectively do the surface positions capture the movement of the animal?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_stillNotWell"> _Show/Hide Answers_ </button>
<div id="ans_stillNotWell" class="collapse">
*They don't capture the movement of the animal very well. So much happens between the positions that we're able to collect at the surface, so we are missing out if surface positions is all we collect/show.*
</div>
# Tracks Colored by a Variable
We often want to plot a track colored proportional to another variable of interest. For example, it can be useful to see where the animal is diving along the track. To colour the track by depth, use `col_line()` with `P` as the color information:
```{r, echo = TRUE, eval = FALSE}
CF <- col_line(y = northing, x = easting, data = FT, color = d7$P$data)
CF
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#col_line_static"> Show/Hide Results </button>
<div id="col_line_static" class="collapse">
```{r, echo = TRUE, eval = FALSE}
CF <- col_line(y = northing, x = easting, data = FT, color = d7$P$data)
CF
```
</div>
-->
What about an interactive version?
```{r, echo = TRUE, eval = FALSE}
col_line(x = FT$easting, y = FT$northing, color = d7$P$data)
```
Zoom in to see what the scale is of the track [tortuosity](https://en.wikipedia.org/wiki/Tortuosity): Is there tortuosity within individual dives, or is the tortuosity occurring across dives?
<button class="btn btn-primary" data-toggle="collapse" data-target="#ans_tortuous_errywher"> _Show/Hide Answers_ </button>
<div id="ans_tortuous_errywher" class="collapse">
*While most of the tortuosity occurs across dives, there is some tortuosity within dives as well. The following picture comes from the next chunk of code.*
```{r, echo=FALSE, eval = TRUE}
knitr::include_graphics('tortuosity-within-dive.png')
```
</div>
If you are done and want a challenge, try colouring the track by the absolute roll angle instead of the depth. You could use col_line3 to plot the 3-d positions (i.e., `FT$easting`, `FT$northing`, and `d7\$P\$data`), and then colour the plot by absolute roll angle remembering to convert from radians to degrees.
```{r, echo = TRUE, eval = FALSE}
pitch_roll <- a2pr(d7$Aa)
roll_deg <- pitch_roll$r/pi*180
col_line3(x = FT$easting, y = FT$northing,
z = d7$P$data, c = roll_deg)
```
<button class="btn btn-primary" data-toggle="collapse" data-target="#col_line3"> Show/Hide Results </button>
<div id="col_line3" class="collapse">
```{r, echo = FALSE, eval = TRUE}
pitch_roll <- a2pr(d7$Aa)
roll_deg <- pitch_roll$r/pi*180
col_line3(x = FT$easting, y = FT$northing,
z = d7$P$data, c = roll_deg)
```
</div>
# Review
What have you learned? What you can get out of---despite what problems can arise from---dead-reckoning.
And that's it! Fantastic work with this vignette.
*If you'd like to continue working through these vignettes, you could take a stab at `find-dives` or `dive-stats`. These deal with the functions find\_dives and dive\_stats, respectively.*
```{r, echo = TRUE, eval = FALSE}
vignette('find-dives')
vignette('dive-stats')
```
***
Animaltags home pages: http://animaltags.org/ (old), https://animaltags.netlify.app/ (new), https://github.com/stacyderuiter/TagTools (for latest beta source code), https://stacyderuiter.github.io/TagTools/articles/TagTools (vignettes overview)