TestVignette
John Fieberg
Purpose
- test whether you have everything set up for the workshop
- retrieve and explore animal movement data from Movebank
- illustrate how to work with functions in the animal movement tools package (amt) ## Preamble
Load libraries
Fist make sure all packages that are needed are available:
packages_needed <- c("knitr", "lubridate", "maptools", "raster", "move",
"amt", "tibble", "leaflet", "dplyr", "readr", "ggplot2",
"glmmTMB", "lme4", "tidyr", "purrr", "glue", "sf",
"here", "moveVis", "GGally", "devtools", "TwoStepCLogit",
"broom", "tictoc", "ezknitr", "moveVis", "maps", "rgeos",
"maptools")
new_packages <- packages_needed[!(packages_needed %in%
installed.packages()[,"Package"])]
if(length(new_packages))
install.packages(new_packages, repos = "https://cloud.r-project.org")Now load all packages.
library(knitr)
library(lubridate)
library(maptools)
library(raster)
library(move)
library(amt)
library(tibble)
library(leaflet)
library(dplyr)
library(readr)
library(ggplot2)
library(glmmTMB)
library(sf)
library(here)
options(width=165,digits.secs = 3)
opts_chunk$set(fig.width=12,fig.height=4.5, error=TRUE,cache = FALSE)Record time for running all code
ptm<-proc.time()Read in data from movebank
Create a login object for a user account at movebank.org
loginStored <- movebankLogin(username="MovebankWorkshop", password="genericuserforthisexercise")Get overview information about a Movebank study. Be sure to check the citation and license terms if not using your own data.
getMovebankStudy(study="Martes pennanti LaPoint New York", login=loginStored) # see study-level info## acknowledgements
## 1 This work was supported in part by NSF (grant 0756920 to RK), the New York State Museum, \nthe Max Planck Institute for Ornithology; the North Carolina Museum of Natural Sciences, and a National Geographic Society Waitt Grant (SDL). We thank the Albany Pine Bush Preserve, Dina Dechmann, Neil Gifford, Wolfgang Heidrich, Bart Kranstauber, Franz Kümmeth, Roger Powell, Kamran Safi, Marco Smolla, and Brad Stratton for their logistical support and valuable input.
## bounding_box
## 1 NA
## citation
## 1 LaPoint, S, Gallery P, Wikelski M, Kays R (2013) Animal behavior, cost-based corridor models, and real corridors. Landscape Ecology, v 28 i 8, p 1615–1630. doi:10.1007/s10980-013-9910-0
## comments enable_for_animal_tracker grants_used has_quota i_am_owner
## 1 NA NA National Science Foundation (#0756920 to RWK) and the National Geographic Society's Waitt Grant (to SDL) true false
## id
## 1 6925808
## license_terms
## 1 These data have been published by the Movebank Data Repository with DOI 10.5441/001/1.2tp2j43g. See www.datarepository.movebank.org/handle/10255/move.328. Please contact the PI prior to use of these data for any purposes.
## location_description main_location_lat main_location_long name number_of_deployments number_of_events number_of_individuals
## 1 NA 42.75205 -73.86246 Martes pennanti LaPoint New York NA NA NA
## number_of_tags principal_investigator_address principal_investigator_email principal_investigator_name
## 1 NA NA NA NA
## study_objective
## 1 To investigate the use of corridors by fisher within a suburban landscape and to validate cost-based corridor model predictions with both animal tracking data and camera traps.
## study_type suspend_license_terms timestamp_end timestamp_start i_can_see_data there_are_data_which_i_cannot_see
## 1 research false NA NA true falseYou may receive a message that you do not have access to data in a study. If you are not a data manager for a study, the owner may require that you read and accept the license terms once before you may download it. Currently license terms cannot be accepted directly from the move package. To view and accept them, go to movebank.org, search for the study and begin a download, and view and accept the license terms. After that you will be able to load from R. Load data from a study in Movebank and create a MoveStack object. For more details and options see https://cran.r-project.org/web/packages/move/index.html.
fisher.move <- getMovebankData(study = "Martes pennanti LaPoint New York",
login = loginStored)
head(fisher.move)## sensor_type_id behavioural_classification eobs_battery_voltage eobs_fix_battery_voltage eobs_horizontal_accuracy_estimate eobs_key_bin_checksum
## 5 653 NA NA NA NA
## 6 653 NA NA NA NA
## 7 653 NA NA NA NA
## 9 653 NA NA NA NA
## 15 653 NA NA NA NA
## 16 653 NA NA NA NA
## eobs_speed_accuracy_estimate eobs_start_timestamp eobs_status eobs_temperature eobs_type_of_fix eobs_used_time_to_get_fix ground_speed heading
## 5 NA NA NA NA NA 80 NA NA
## 6 NA NA NA NA NA 7 NA NA
## 7 NA NA NA NA NA 58 NA NA
## 9 NA NA NA NA NA 90 NA NA
## 15 NA NA NA NA NA 71 NA NA
## 16 NA NA NA NA NA 57 NA NA
## height_above_ellipsoid location_lat location_long manually_marked_outlier timestamp update_ts deployment_id event_id tag_id
## 5 NA 42.79528 -73.86015 NA 2011-02-11 18:06:13.999 1970-01-01 00:00:00.000 8623923 170635812 8623917
## 6 NA 42.79534 -73.86001 NA 2011-02-11 18:10:08.999 1970-01-01 00:00:00.000 8623923 170635813 8623917
## 7 NA 42.79528 -73.86013 NA 2011-02-11 18:20:58.997 1970-01-01 00:00:00.000 8623923 170635814 8623917
## 9 NA 42.79483 -73.86012 NA 2011-02-11 18:41:30.999 1970-01-01 00:00:00.000 8623923 170635816 8623917
## 15 NA 42.79509 -73.86037 NA 2011-02-11 19:41:10.997 1970-01-01 00:00:00.000 8623923 170635822 8623917
## 16 NA 42.79515 -73.86023 NA 2011-02-11 19:50:58.000 1970-01-01 00:00:00.000 8623923 170635823 8623917
## sensor_type
## 5 GPS
## 6 GPS
## 7 GPS
## 9 GPS
## 15 GPS
## 16 GPSNote: If there are duplicate animal-timestamp records in Movebank, you will get a warning. You can exclude duplicate records on import using removeDuplicatedTimestamps = TRUE. If you are a data manager for a study in Movebank you can also filter them out directly in the study so they are excluded by default in downloads (see https://www.movebank.org/node/27252). Create a data frame from the MoveStack object
fisher.dat <- as(fisher.move, "data.frame")Data cleaning
Delete observations where missing lat or long or a timestamp. There are no missing observations in this data set, but it is still good practice to check.
ind <- complete.cases(fisher.dat[, c("location_lat", "location_long", "timestamp")])The number of relocations with missing coordinates or timestamp (if any).
table(ind)## ind
## TRUE
## 32904fisher.dat <- fisher.dat %>% filter(ind)Check for duplicated observations (ones with same lat, long, timestamp, and individual identifier). There are no duplicate observations in this data set, but it is still good practice to check.
ind2 <- fisher.dat %>%
select(timestamp, location_long, location_lat, local_identifier) %>%
duplicated
sum(ind2) # no duplicates## [1] 0fisher.dat <- fisher.dat %>% filter(!ind2)Make timestamp a date/time variable
fisher.dat <- fisher.dat %>% mutate(timestamp = ymd_hms(timestamp, tz = "UTC"))Move package
Look at functions in the move package.
plot(fisher.move)
show(fisher.move)## class : MoveStack
## features : 32904
## extent : -73.94032, -73.38947, 42.70898, 42.851 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0
## variables : 24
## names : sensor_type_id, behavioural_classification, eobs_battery_voltage, eobs_fix_battery_voltage, eobs_horizontal_accuracy_estimate, eobs_key_bin_checksum, eobs_speed_accuracy_estimate, eobs_start_timestamp, eobs_status, eobs_temperature, eobs_type_of_fix, eobs_used_time_to_get_fix, ground_speed, heading, height_above_ellipsoid, ...
## min values : 653, , NA, NA, NA, NA, NA, NA, NA, NA, NA, 3, NA, NA, NA, ...
## max values : 653, traveling, NA, NA, NA, NA, NA, NA, NA, NA, NA, 120, NA, NA, NA, ...
## timestamps : 2009-02-11 12:16:45.000 ... 2011-05-28 09:10:24.999 Time difference of 836 days (start ... end, duration)
## sensors : GPS
## indiv. data : comments, death_comments, earliest_date_born, exact_date_of_birth, individual_id, latest_date_born, local_identifier, nick_name, ring_id, sex, taxon_canonical_name
## min ID Data : "Leroy", adult male, 5.6kg. Fixed temporal schedule GPS tag = 15 minutes., , NA, NA, 8623916, NA, F1, NA, NA, f, NA
## max ID Data : "Zissou", subadult female, 2.35kg. Accelerometer-informed GPS tag, temporal schedule = 2 to 60 minute interval., 7 March 2011 road-killed at 42.752050, -73.820170, NA, NA, 8630581, NA, M5, NA, NA, m, NA
## individuals : F1, F2, F3, M1, M2, M3, M4, M5
## unused rec. : 14443
## license : These data have been published by the Movebank Data Repository with DOI 10.5441/001/1.2tp2j43g. See www.datarepository.movebank.org/handle/10255/mo ...
## citation : LaPoint, S, Gallery P, Wikelski M, Kays R (2013) Animal behavior, cost-based corridor models, and real corridors. Landscape Ecology, v 28 i 8, p 1615–1630. doi:10.1007/s10980-013-9910-0
## study name : Martes pennanti LaPoint New York
## date created: 2019-03-08 15:06:28.398summary(fisher.move)## Object of class MoveStack
## Coordinates:
## min max
## location_long -73.94032 -73.38947
## location_lat 42.70898 42.85100
## Is projected: FALSE
## proj4string :
## [+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0]
## Number of points: 32904
## Data attributes:
## sensor_type_id behavioural_classification eobs_battery_voltage eobs_fix_battery_voltage eobs_horizontal_accuracy_estimate eobs_key_bin_checksum
## Min. :653 :30468 Mode:logical Mode:logical Mode:logical Mode:logical
## 1st Qu.:653 active : 551 NA's:32904 NA's:32904 NA's:32904 NA's:32904
## Median :653 hunting : 1521
## Mean :653 resting : 63
## 3rd Qu.:653 traveling: 301
## Max. :653
## eobs_speed_accuracy_estimate eobs_start_timestamp eobs_status eobs_temperature eobs_type_of_fix eobs_used_time_to_get_fix ground_speed heading
## Mode:logical Mode:logical Mode:logical Mode:logical Mode:logical Min. : 3.00 Mode:logical Mode:logical
## NA's:32904 NA's:32904 NA's:32904 NA's:32904 NA's:32904 1st Qu.: 7.00 NA's:32904 NA's:32904
## Median : 13.00
## Mean : 24.74
## 3rd Qu.: 37.00
## Max. :120.00
## height_above_ellipsoid location_lat location_long manually_marked_outlier timestamp update_ts deployment_id
## Mode:logical Min. :42.71 Min. :-73.94 Mode:logical Min. :2009-02-11 12:16:45.0 1970-01-01 00:00:00.000:32904 Min. :8623923
## NA's:32904 1st Qu.:42.76 1st Qu.:-73.90 NA's:32904 1st Qu.:2010-03-23 03:57:38.0 1st Qu.:8624627
## Median :42.83 Median :-73.82 Median :2011-02-08 09:33:13.5 Median :8625155
## Mean :42.81 Mean :-73.69 Mean :2010-11-15 07:50:30.5 Mean :8627068
## 3rd Qu.:42.84 3rd Qu.:-73.42 3rd Qu.:2011-04-17 08:12:59.5 3rd Qu.:8630587
## Max. :42.85 Max. :-73.39 Max. :2011-05-28 09:10:25.0 Max. :8630587
## event_id tag_id sensor_type
## Min. :170635812 Min. :8623917 GPS:32904
## 1st Qu.:170652431 1st Qu.:8624568
## Median :170664696 Median :8624887
## Mean :170671625 Mean :8626969
## 3rd Qu.:170696240 3rd Qu.:8630582
## Max. :170705307 Max. :8630582Plots of the data
Note: there are lots of ways to plot data using R. I have included code that illustrates how to create a map using the leaflet package with data for a single individual. This approach can become cumbersome and slow with too many observations.
Lets look at the data from F2. Note: When reading data from Movebank always be sure to refer to the animal-id and individual-local-identifer and not the tag-id or tag-local-identifier in order to exclude pre- and post-deployment records and correctly separate out individual animals when the study includes redeployments.
fisherF2 <- fisher.dat %>% filter(local_identifier == "F2")leaflet
leaflet lets you interact / move the map around which is nice.
leaflet(fisherF2) %>% addTiles()%>%
addCircles(fisherF2$location_long, fisherF2$location_lat)Using ggplot2 to plot the data
Use separate axes for each individual (add scales=“free” to facet_wrap)
ggplot(fisher.dat, aes(x = location_long,
y = location_lat)) + geom_point() +
facet_wrap(~local_identifier, scales = "free")
Now, all on 1 plot
ggplot(fisher.dat,
aes(location_long, location_lat, color = local_identifier,
group = local_identifier))+
geom_point() + coord_equal() +
theme(legend.position = "bottom")
Creating a track in amt
Before we can use the amt package to calculate step lengths, turn angles, and bearings for fisher data, we need to add a class (track) to the data. Then, we can summarize the data by individual, month, etc.
If we have a data set with projected coordinates, we can use: trk.temp <- make_track(fisher.dat, .x=utm.easting, .y=utm.northing, .t=timestamp, id = individual_local.identifier)
Note: we did not need to explicitly specify x, y and t (but probably good to do so). This would also have worked trk <- make_track(fisher.dat, utm.easting, utm.northing, timestamp, id = local_identifier).
We can also use lat, long, which will allow us to determine time of day.
trk <- make_track(fisher.dat, .x = location_long, .y = location_lat,
.t = timestamp, id = local_identifier, crs = CRS("+init=epsg:4326"))## .t found, creating `track_xyt`.Now it is easy to calculate day/night from the movement track
trk <- trk %>% time_of_day()Now, we can transform the geographic coordinates to projected coordinates (we use NAD83).
trk <- transform_coords(trk, CRS("+init=epsg:5070"))Movement Characteristics
Functions in the amt package:
- dir_abs will calculate absolute angles for steps
- dir_rel will calculate turning angles (relative angles)
- step_lengths will calculate distances between points
- nsd = Net Squared Displacement (distance from first point)
Arguments direction_abs:
- full_circle will calculate between 0 and 2 \(\pi\) (rather than \(-\pi\) and \(\pi\))
- zero gives the direction = 0 for absolute angle
Note: we have to calculate these characteristics separately for each individual (to avoid calculating a distance between say the last observation of the first individual and the first observation of the second individual).
To do this, we could loop through individuals, calculate these characteristics for each individual, then rbind the data back together. Or, use nested data frames and the map function in the purrr package to do this with very little code.
To see how nesting works, we can create a nested object by individual
nesttrk <- trk %>% group_by(id) %>% # here we say which variable (column) contains
# information about groups. Note, this could be more than one.
nest() # And now we nest the data frame, creating a list column.
nesttrk## # A tibble: 8 x 2
## id data
## <fct> <list>
## 1 M1 <tibble [919 x 4]>
## 2 M4 <tibble [8,958 x 4]>
## 3 F2 <tibble [3,004 x 4]>
## 4 M3 <tibble [2,436 x 4]>
## 5 F3 <tibble [1,501 x 4]>
## 6 M2 <tibble [1,638 x 4]>
## 7 F1 <tibble [1,349 x 4]>
## 8 M5 <tibble [13,099 x 4]>Each row contains data from an individual. For example, we can access data from the first individual using:
nesttrk$data[[1]]## # A tibble: 919 x 4
## x_ y_ t_ tod_
## * <dbl> <dbl> <dttm> <chr>
## 1 1782673. 2402297. 2009-02-11 12:16:45.000 day
## 2 1782680. 2402297. 2009-02-11 12:31:38.000 day
## 3 1782683. 2402292. 2009-02-11 12:45:48.996 day
## 4 1782686. 2402305. 2009-02-11 13:00:16.001 day
## 5 1782681. 2402297. 2009-02-11 13:15:19.000 day
## 6 1782671. 2402290. 2009-02-11 13:30:13.000 day
## 7 1782683. 2402298. 2009-02-11 13:45:37.000 day
## 8 1782686. 2402281. 2009-02-11 14:00:35.000 day
## 9 1782682. 2402290. 2009-02-11 14:15:49.000 day
## 10 1782681. 2402291. 2009-02-11 14:30:49.000 day
## # ... with 909 more rowsWe could calculate movement characteristics by individual using:
temp <- direction_rel(nesttrk$data[[1]])
head(temp)## [1] NA -0.9787252 2.2942174 2.7885018 -0.4359982 3.1269565or:
temp <- trk %>% filter(id=="M1") %>% direction_rel
head(temp)## [1] NA -0.9787252 2.2942174 2.7885018 -0.4359982 3.1269565Or, we can add a columns to each nested column of data using purrr::map
trk1 <- trk %>% group_by(id) %>% nest() %>%
mutate(dir_abs = map(data, direction_abs, full_circle = TRUE, zero = "N", clockwise = TRUE),
dir_rel = map(data, direction_rel),
sl = map(data, step_lengths),
nsd_=map(data, nsd)) %>% unnest()
trk1## # A tibble: 32,904 x 9
## id dir_abs dir_rel sl nsd_ x_ y_ t_ tod_
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dttm> <chr>
## 1 M1 1.54 NA 6.35 0 1782673. 2402297. 2009-02-11 12:16:45.000 day
## 2 M1 2.52 -0.979 5.63 40.3 1782680. 2402297. 2009-02-11 12:31:38.000 day
## 3 M1 0.225 2.29 12.7 112. 1782683. 2402292. 2009-02-11 12:45:48.996 day
## 4 M1 3.72 2.79 9.65 220. 1782686. 2402305. 2009-02-11 13:00:16.001 day
## 5 M1 4.16 -0.436 11.7 51.7 1782681. 2402297. 2009-02-11 13:15:19.000 day
## 6 M1 1.03 3.13 13.9 46.3 1782671. 2402290. 2009-02-11 13:30:13.000 day
## 7 M1 2.93 -1.90 17.4 84.6 1782683. 2402298. 2009-02-11 13:45:37.000 day
## 8 M1 5.84 -2.91 10.5 423. 1782686. 2402281. 2009-02-11 14:00:35.000 day
## 9 M1 5.46 0.376 1.57 112. 1782682. 2402290. 2009-02-11 14:15:49.000 day
## 10 M1 0.526 -1.35 7.51 80.9 1782681. 2402291. 2009-02-11 14:30:49.000 day
## # ... with 32,894 more rowsNow, calculate month, year, hour, week of each observation and append these to the dataset. Unlike the movement charactersitics, these calculations can be done all at once, since they do not utilize successive observations (like step lengths and turn angles do).
trk1 <- trk1 %>%
mutate(
week = week(t_),
month = month(t_, label=TRUE),
year = year(t_),
hour = hour(t_)
)Some plots of movement characteristics
Absolute angles (for each movement) relative to North
We could use a rose diagram (below) to depict the distribution of angles.
ggplot(trk1, aes(x = dir_abs, y = ..density..)) +
geom_histogram(breaks = seq(0, 2 * pi, len = 30))+
coord_polar(start = 0) + theme_minimal() +
scale_fill_brewer() + labs(y = "Density", title = "Angles Direct") +
scale_x_continuous(limits = c(0, 2 * pi),
breaks = c(0, pi/2, pi, 3 * pi/2),
labels = c("0", "pi/2", "pi", "3pi/2")) +
facet_wrap( ~ id)## Warning: Removed 8 rows containing non-finite values (stat_bin).
Turning angles
Note: a 0 indicates the animal continued to move in a straight line, a \(\pi\) indicates the animal turned around (but note, resting + measurement error often can make it look like the animal turned around).
ggplot(trk1, aes(x = dir_rel, y = ..density..)) +
geom_histogram(breaks = seq(-pi, pi, length = 20))+
coord_polar(start = 0) + theme_minimal() +
scale_fill_brewer() + ylab("Density") + ggtitle("Angles Direct") +
scale_x_continuous(limits = c(-pi, pi), breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c("-pi", "-pi/2", "0", "pi/2", "pi")) +
facet_wrap(~id)## Warning: Removed 16 rows containing non-finite values (stat_bin).
Net-squared displacement over time for each individual
ggplot(trk1, aes(x = t_, y = nsd_)) + geom_path()+
facet_wrap(~id, scales = "free")
Explore movement characteristics by (day/night, hour, month)
step length distribution by day/night
ggplot(trk1, aes(x = tod_, y = log(sl))) +
geom_boxplot() + geom_smooth() + facet_wrap(~id)
Space use (MCP or KDE) by week, month, and year
Note: this code will only work for your critters if you have multiple observations for each combination of (month, year). If you don’t have many observations, you could try: group_by(id, year).
mcps.week <- trk %>%
mutate(year = year(t_), month = month(t_), week = week(t_)) %>%
group_by(id, year, month, week) %>% nest() %>%
mutate(mcparea = map(data, ~ hr_mcp(., levels = c(0.95)) %>% hr_area)) %>%
select(id, year, month, week, mcparea) %>% unnest()ggplot(mcps.week, aes(x = week, y = area, colour = factor(year))) +
geom_point() +
geom_smooth() + facet_wrap(~id, scales="free")
Same for KDE
kde.week <- trk %>%
mutate(year = year(t_), month = month(t_), week = week(t_)) %>%
group_by(id, year, month, week) %>% nest() %>%
mutate(kdearea = map(data, ~ hr_kde(., levels=c(0.95)) %>% hr_area)) %>%
select(id, year, month, week, kdearea) %>% unnest()ggplot(kde.week, aes(x = week, y = kdearea, colour = factor(year))) +
geom_point() +
geom_smooth() + facet_wrap(~id, scales = "free")
Write out files for later analyses
write_rds(trk, path = here("data","trk.Rdata"))
write_rds(fisher.dat, path = here("data","fisher.Rdata"))Investigate step length and turn angle distributions further
We will use only one animal to illustrate these approaches, but it could be easily extended to multiple animals using the same strategies as above.
Different habitats
dat <- trk %>% filter(id == "F1")We first have to make sure that the sampling rate is constant. Otherwise step length and turning angles can not be compared. We first can check what the current sampling rate is:
summarize_sampling_rate(dat)## # A tibble: 1 x 9
## min q1 median mean q3 max sd n unit
## <S3: table> <S3: table> <S3: table> <S3: table> <S3: table> <S3: table> <dbl> <int> <chr>
## 1 3.916667 9.799996 10.01662 20.69303 10.29995 1650.317 95.2 1348 min10 minutes seems to be appropriate. We can resample our track to a 10 min sampling rate with a tolerance of 1 minute.
dat <- track_resample(dat, rate = minutes(10), tolerance = minutes(1)) %>%
filter_min_n_burst(min_n = 3)dat gained a new column: burst_. This column indicates for each point to which burst it belongs. A burst is a sequence of relocations with equidstant (in time) relocations.
steps <- dat %>% steps_by_burst()
steps %>% ggplot(aes(sl_)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Let us add to each step in which habitat it started. It might be good to group some of the categories. 1 = 21, 22, 23: Developed 4 = 90: Wetland
data("amt_fisher_lu")
lu <- amt_fisher_lu %in% c(41:43, 90)
plot(lu)
steps %>% extract_covariates(lu, where = "start") %>%
ggplot(aes(sl_)) + geom_histogram() + facet_wrap(~ layer)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
steps %>% extract_covariates(lu, where = "both") %>%
mutate(where = case_when(
layer_start == 1 & layer_end == 1 ~ "in forest",
layer_start == 0 & layer_end == 0 ~ "in other",
layer_start != layer_end ~ "transistion")) %>%
ggplot(aes(sl_)) + geom_density() +
facet_wrap(~ where)
Look at step turning angles for long and short steps
We first need to split steps into long and short steps and therefore determine a threshold.
steps %>% ggplot(aes(sl_)) + geom_histogram() +
geom_vline(xintercept = 120, col = "red")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
120 m seems to be reasonabLe.
steps %>%
mutate(step_length = case_when(sl_ > 120 ~ "long", TRUE ~ "short")) %>%
ggplot(aes(x = ta_, y = ..density..)) +
geom_histogram(breaks = seq(-pi, pi, length = 20))+
coord_polar(start = 0) + theme_minimal() +
scale_fill_brewer() + ylab("Density") + ggtitle("Angles Direct") +
scale_x_continuous(limits = c(-pi, pi), breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c("-pi", "-pi/2", "0", "pi/2", "pi")) +
facet_wrap(~ step_length)## Warning: Removed 66 rows containing non-finite values (stat_bin).
We can see that longer steps tend to be more directed (with turn angles near 0), and that short steps often have turn angles near -\(\pi\) or \(\pi\). Some of these angles associated with small steps are likely due to measurement error.
Here is another plot of turn angles for short and long steps.
steps %>%
mutate(step_length = case_when(sl_ > 120 ~ "long", TRUE ~ "short")) %>%
ggplot(aes(x = ta_, y = ..density..)) +
geom_histogram(breaks = seq(-pi, pi, length = 20))+
scale_x_continuous(limits = c(-pi, pi), breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c("-pi", "-pi/2", "0", "pi/2", "pi")) +
facet_wrap(~ step_length)## Warning: Removed 66 rows containing non-finite values (stat_bin).