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 as of CmdStanR 1.0.0 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 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 1 finished in 187.5 seconds.
## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 2 finished in 190.1 seconds.
## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 3 finished in 230.6 seconds.
## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 4 finished in 231.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 209.8 seconds.
## Total execution time: 231.2 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    1363.  11.4
##  3 49237 2023-08-20 20:20:00  8.79  54.8    1548.  12.9
##  4 49237 2023-08-20 20:22:00  8.79  54.8    1543.  12.9
##  5 49237 2023-08-20 20:24:00  8.78  54.8    1515.  12.6
##  6 49237 2023-08-20 20:26:00  8.77  54.8    1647.  13.7
##  7 49237 2023-08-20 20:28:00  8.76  54.8    1943.  16.2
##  8 49237 2023-08-20 20:30:00  8.75  54.8    2446.  20.4
##  9 49237 2023-08-20 20:32:00  8.73  54.7    2301.  19.2
## 10 49237 2023-08-20 20:34:00  8.73  54.7    1946.  16.2
## # ℹ 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  1363.  348. 3560.  1.01     923.    1332.
##  3 49237 2023-08-20 20:20:00 distance  1548.  437. 3739.  1.00    1492.    2192.
##  4 49237 2023-08-20 20:22:00 distance  1543.  401. 3590.  1.00    1653.    1994.
##  5 49237 2023-08-20 20:24:00 distance  1515.  387. 3416.  1.00    3221.    3785.
##  6 49237 2023-08-20 20:26:00 distance  1647.  430. 3578.  1.00    3113.    2521.
##  7 49237 2023-08-20 20:28:00 distance  1943.  556. 4038.  1.00    3802.    3753.
##  8 49237 2023-08-20 20:30:00 distance  2446.  779. 4616.  1.00    2225.    2121.
##  9 49237 2023-08-20 20:32:00 distance  2301.  679. 4531.  1.00    2688.    2700.
## 10 49237 2023-08-20 20:34:00 distance  1946.  499. 4070.  1.00    4480.    2486.
## # ℹ 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     1 1000 8.808933 54.81710 49237_1000
## 2 49237 2023-08-20 20:18:00     1 1000 8.802298 54.81925 49237_1000
## 3 49237 2023-08-20 20:20:00     1 1000 8.789184 54.79849 49237_1000
## 4 49237 2023-08-20 20:22:00     1 1000 8.787332 54.78480 49237_1000
## 5 49237 2023-08-20 20:24:00     1 1000 8.792912 54.78192 49237_1000
## 6 49237 2023-08-20 20:26:00     1 1000 8.791950 54.77302 49237_1000
## 7 49237 2023-08-20 20:28:00     1 1000 8.787439 54.75661 49237_1000
##  [ 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.

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.

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"))