Quickstart example

The methods provided in this package are applicable to various types of automated radio-telemetry data, including—but not limited to—the Motus Wildlife Tracking System.

library(movetrack)

## Warning: The 'fix' argument is deprecated and will be removed in a future
## release.

library(ggplot2)
theme_set(theme_bw(base_size = 15))

# Load example data
data(motusData)

Modelling

Parse your preprocessed raw data to track() and adjust additional arguments as needed. Input arguments are unique track IDs, timestamps, station locations, antenna bearings, and signal strength for each detection. By default, these input variable names align with standard Motus variable conventions. Use aRange to specify theoretical antenna detection ranges in kilometres. This can be a single integer value or a named list of values for different antenna types defined in the aType column. By default, approximations for different Yagi-antennas, as suggested in the Motus Docs, are used:

Antenna type	Theoretical range
3-element Yagi	~5 km
4-element Yagi	~6 km
5-element Yagi	~8 km
6-element Yagi	~10 km
9-element Yagi	~15 km

With dTime, specify the time interval in minutes for which positions are estimated. The default is 2 min, but this should be adjusted according to the tag’s burst interval and the expected flight speed—longer burst intervals or faster movements typically require longer intervals. Specify the number of states in the model using the states argument, where each state represents a distinct movement behaviour such as migratory flight, local movement, or stopover (default: 2). Optionally, the argument i_lambda can be used to fix the lag-correlation parameter across all tracks (see vignette("hmm")). Set model = FALSE to extract the internally calculated coarse raw positions and measurement errors for each time interval without modelling.

Here, we model the animals’ flight paths from the example data using track() with default arguments. The number of printed screen updates is reduced using the refresh argument passed to cmdstanr::sample(). See there for additional sampling options.

fit <- track(motusData, refresh = 1e3)

## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 2 finished in 184.6 seconds.
## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 4 finished in 188.5 seconds.
## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 3 finished in 220.0 seconds.
## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 1 finished in 259.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 213.1 seconds.
## Total execution time: 259.3 seconds.

Inspect results

track() returns a movetrack object containing complete posterior distributions for positions at each time step, along with corresponding distances in m and speeds between positions in m/s. We can peek at the results by simply running

fit

## # A tibble: 161 × 6
##       ID time                  lon   lat distance speed
##    <int> <dttm>              <dbl> <dbl>    <dbl> <dbl>
##  1 49237 2023-08-20 20:16:00  8.80  54.8      NA   NA  
##  2 49237 2023-08-20 20:18:00  8.80  54.8    1319.  11.0
##  3 49237 2023-08-20 20:20:00  8.79  54.8    1531.  12.8
##  4 49237 2023-08-20 20:22:00  8.79  54.8    1494.  12.5
##  5 49237 2023-08-20 20:24:00  8.78  54.8    1484.  12.4
##  6 49237 2023-08-20 20:26:00  8.77  54.8    1637.  13.6
##  7 49237 2023-08-20 20:28:00  8.76  54.8    1984.  16.5
##  8 49237 2023-08-20 20:30:00  8.75  54.8    2462.  20.5
##  9 49237 2023-08-20 20:32:00  8.73  54.7    2280.  19.0
## 10 49237 2023-08-20 20:34:00  8.73  54.7    1951.  16.3
## # ℹ 151 more rows

and transform it into a data.frame with

df <- as.data.frame(fit)

We can summarise one or multiple variables and inspect convergence measures using summary.movetrack():

summary(fit, var = "distance")

## # A tibble: 161 × 9
##       ID time                variable median    q5   q95  rhat ess_bulk ess_tail
##    <int> <dttm>              <chr>     <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 49237 2023-08-20 20:16:00 distance    NA    NA    NA  NA         NA       NA 
##  2 49237 2023-08-20 20:18:00 distance  1319.  356. 3501.  1.01    1153.    2098.
##  3 49237 2023-08-20 20:20:00 distance  1531.  377. 3623.  1.00    1838.    2985.
##  4 49237 2023-08-20 20:22:00 distance  1494.  394. 3546.  1.00    1796.    2373.
##  5 49237 2023-08-20 20:24:00 distance  1484.  389. 3398.  1.00    3191.    2877.
##  6 49237 2023-08-20 20:26:00 distance  1637.  424. 3564.  1.00    3681.    3744.
##  7 49237 2023-08-20 20:28:00 distance  1984.  575. 3954.  1.00    3294.    3452.
##  8 49237 2023-08-20 20:30:00 distance  2462.  836. 4634.  1.00    1935.    2336.
##  9 49237 2023-08-20 20:32:00 distance  2280.  687. 4485.  1.00    2751.    2163.
## 10 49237 2023-08-20 20:34:00 distance  1951.  530. 3964.  1.00    4457.    3088.
## # ℹ 151 more rows

To easily extract a defined number of posterior draws from the model for each time interval, run:

getDraws(fit, n = 10)

##      ID                time chain iter      lon      lat       tID
## 1 49237 2023-08-20 20:16:00     3  206 8.845960 54.81631 49237_206
## 2 49237 2023-08-20 20:18:00     3  206 8.816311 54.81043 49237_206
## 3 49237 2023-08-20 20:20:00     3  206 8.805651 54.80459 49237_206
## 4 49237 2023-08-20 20:22:00     3  206 8.795321 54.77523 49237_206
## 5 49237 2023-08-20 20:24:00     3  206 8.772312 54.77892 49237_206
## 6 49237 2023-08-20 20:26:00     3  206 8.752631 54.79020 49237_206
## 7 49237 2023-08-20 20:28:00     3  206 8.767458 54.79145 49237_206
##  [ reached 'max' / getOption("max.print") -- omitted 1603 rows ]

Plot results

We can plot the results per individual and output variable using plot():

plot(fit, vars = "lon", id = 49237)

Estimated median longitude together with 90% highest posterior density intervals over time.

Map flight paths

We can visualise the results on a static map using mapTrack():

mapTrack(fit) +
  geom_point(aes(recvDeployLon, recvDeployLat), data = motusData)

Modelled movement trajectories per individual. Posterior means are shown together with 50 posterior draws, circles on the map indicate receiver locations with detections of the animals.

We could also create an interactive Leaflet map using the following code:

library(sfheaders)
library(leaflet)

# Extract draws
draws <- getDraws(fit) |>
  sf_linestring("lon", "lat", linestring_id = "tID")

# Leaflet map
fit |>
  as.data.frame() |>
  sf_linestring("lon", "lat", linestring_id = "ID") |>
  leaflet() |>
  addTiles() |>
  addPolylines(data = draws, color = "grey", weight = 1, opacity = 0.2) |>
  addPolylines(color = ~ c("orange", "blue"))