In paleolimbot/headings: Tools for Working with Compass Headings

library(tidyverse)
set.seed(39374)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  dpi = 300
)

headings

The goal of headings is to provide well-tested functions for working with orientation data (headings), including aggregation (e.g., mean, sd) and correction (e.g., declination).

Installation

You can install the development version from GitHub with:

# install.packages("devtools")
devtools::install_github("paleolimbot/headings")

Example

The headings package can calculate a proper mean and standard deviation of compass headings:

library(headings)

# totally bogus since 0 and 360 are identical!
mean(c(350, 355, 360, 0, 5, 10))
sd(c(350, 355, 360, 0, 5, 10))

# everybody wins!
hdg_mean(c(350, 355, 360, 0, 5, 10))
hdg_sd(c(350, 355, 360, 0, 5, 10))

For more statistical analyses of directions, you can use the circular package. A convenience constructor hdg_circular() is provided to construct objects in the way that the circular package expects for compass headings.

Functions to correct for magnetic declination are also provided for the IGRF13, WMM2020, and EMM2017 models.

hdg_decl(-64, 45, year = 2021)
hdg_true_from_magnetic(13.40, -64, 45, year = 2021)
hdg_magnetic_from_true(356.51, -64, 45, year = 2021)

As an example, heading aggregation can be used to summarize climate and/or current directions over time.

library(tidyverse)
data("kamloops2016")

kamloops2016 %>% 
  as_tibble() %>%
  filter(is.finite(wind_dir), is.finite(wind_spd)) %>% 
  group_by(month) %>% 
  summarise(
    mean_wtd = hdg_mean(wind_dir, weights = wind_spd),
    sd_wtd = hdg_sd(wind_dir, weights = wind_spd)
  )

The headings package also comes with a circle-aware kernel density estimator for visualizing vectors of headings. For circular kernel densities in a non-visualization context, use the circular package (circular::density.circular()).

plot(
  hdg_density(
    kamloops2016$wind_dir, 
    weights = kamloops2016$wind_spd, 
    na.rm = TRUE
  )
)

To use in ggplot2, you will need to extract the values from the density output:

kamloops2016 %>% 
  group_by(month) %>% 
  summarise(
    broom::tidy(hdg_density(wind_dir, weights = wind_spd))
  ) %>% 
  ggplot(aes(x, y)) +
  geom_line() +
  coord_polar() +
  scale_x_continuous(breaks = seq(0, 360, by = 90)) +
  facet_wrap(vars(month)) +
  theme_bw()