diff --git a/NAMESPACE b/NAMESPACE
index 50f26890..fe75b39e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,7 @@ export(centroid_dyad)
 export(centroid_fusion)
 export(centroid_group)
 export(direction_group)
+export(direction_polarization)
 export(direction_step)
 export(direction_to_centroid)
 export(distance_to_centroid)
diff --git a/R/build_lines.R b/R/build_lines.R
index 7f7dfbc4..f3a3843c 100644
--- a/R/build_lines.R
+++ b/R/build_lines.R
@@ -2,7 +2,7 @@
 #'
 #'
 #' `build_lines` generates a simple feature collection with LINESTRINGs from a
-#' `data.table`. The function accepts a `data.table` with relocation data,
+#' `data.table`. The function expects a `data.table` with relocation data,
 #' individual identifiers, a sorting column and a `projection`. The relocation
 #' data is transformed into LINESTRINGs for each individual and, optionally,
 #' combination of columns listed in `splitBy`. Relocation data should be in two
diff --git a/R/build_polys.R b/R/build_polys.R
index d68fd458..3a0de2cd 100644
--- a/R/build_polys.R
+++ b/R/build_polys.R
@@ -1,7 +1,7 @@
 #' Build Polygons
 #'
 #' `build_polys` generates a simple feature collection with POLYGONs from a
-#' `data.table`. The function accepts a `data.table` with
+#' `data.table`. The function expects a `data.table` with
 #' relocation data, individual identifiers, a projection,
 #' home range type and parameters. The relocation
 #' data is transformed into POLYGONs using either [adehabitatHR::mcp] or
diff --git a/R/direction_polarization.R b/R/direction_polarization.R
new file mode 100644
index 00000000..1b03210c
--- /dev/null
+++ b/R/direction_polarization.R
@@ -0,0 +1,122 @@
+#' Polarization
+#'
+#' \code{direction_polarization} calculates the polarization of individual
+#' directions in each spatiotemporal group identified by \code{group_pts}. The
+#' function expects a \code{data.table} with relocation data appended with a
+#' \code{direction} column from \code{direction_step} and a \code{group} column
+#' from \code{group_pts}.
+#'
+#' The \code{DT} must be a \code{data.table}. If your data is a
+#' \code{data.frame}, you can convert it by reference using
+#' \code{\link[data.table:setDT]{data.table::setDT}} or by reassigning using
+#' \code{\link[data.table:data.table]{data.table::data.table}}.
+#'
+#' The \code{direction} and \code{group} arguments expect the names of columns
+#' in \code{DT} which correspond to the direction and group columns. The
+#' direction column is expected in units of radians and the polarization is
+#' calculated with [CircStats::r.test()].
+#'
+#' @inheritParams direction_group
+#'
+#' @return \code{direction_polarization} returns the input \code{DT} appended
+#'   with a \code{polarization} column representing the direction polarization
+#'   of all individuals in each spatiotemporal group.
+#'
+#'   The direction polarization is calculated using [CircStats::r.test()]
+#'   which expects units of radians.
+#'
+#'   A message is returned when the \code{polarization} columns already
+#'   exists in the input \code{DT}, because it will be overwritten.
+#'
+#' @export
+#' @seealso \code{\link{direction_step}}, \code{\link{group_pts}},
+#'   [CircStats::r.test()]
+#' @family Direction functions
+#'
+#' @references
+#' See example of using polarization:
+#'  * <https://doi.org/10.1016/j.cub.2017.08.004>
+#'  * <10.1371/journal.pcbi.1009437>
+#'  * <https://doi.org/10.7554/eLife.19505>
+#'
+#' @examples
+#' # Load data.table
+#' library(data.table)
+#' \dontshow{data.table::setDTthreads(1)}
+#'
+#' # Read example data
+#' DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+#'
+#' # Cast the character column to POSIXct
+#' DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+#'
+#' # Temporal grouping
+#' group_times(DT, datetime = 'datetime', threshold = '20 minutes')
+#'
+#' # Spatial grouping with timegroup
+#' group_pts(DT, threshold = 50, id = 'ID',
+#'           coords = c('X', 'Y'), timegroup = 'timegroup')
+#'
+#' # Calculate direction at each step
+#' direction_step(
+#'   DT = DT,
+#'   id = 'ID',
+#'   coords = c('X', 'Y'),
+#'   projection = 32736
+#' )
+#'
+#' # Calculate polarization
+#' direction_polarization(DT)
+direction_polarization <- function(
+    DT,
+    direction = 'direction',
+    group = 'group') {
+
+  if (is.null(DT)) {
+    stop('input DT required')
+  }
+
+  if (is.null(direction)) {
+    stop('direction column name required')
+  }
+
+  if (is.null(group)) {
+    stop('group column name required')
+  }
+
+  if (any(!(
+    c(direction, group) %in% colnames(DT)
+  ))) {
+    stop(paste0(
+      as.character(paste(setdiff(
+        c(direction, group),
+        colnames(DT)
+      ), collapse = ', ')),
+      ' field(s) provided are not present in input DT'
+    ))
+  }
+
+  if (any(!(DT[, vapply(.SD, is.numeric, TRUE), .SDcols = c(direction)]))) {
+    stop('direction must be numeric')
+  }
+
+  out <- 'polarization'
+
+  if (out %in% colnames(DT)) {
+    message(paste(out, 'column will be overwritten by this function'))
+    data.table::set(DT, j = out, value = NULL)
+  }
+
+  if (DT[, !inherits(.SD[[1]], 'units'), .SDcols = c(direction)] ||
+      DT[, units(.SD[[1]])$numerator != 'rad', .SDcols = c(direction)]) {
+    stop('units(DT$direction) is not radians, did you use direction_step?')
+  }
+
+  DT[, c(out) := CircStats::r.test(units::drop_units(.SD[[1]]), degree = FALSE)$r.bar,
+    by = c(group),
+    .SDcols = c(direction)]
+
+  return(DT[])
+}
+
+
diff --git a/R/edge_dist.R b/R/edge_dist.R
index 7201ab12..127041f2 100644
--- a/R/edge_dist.R
+++ b/R/edge_dist.R
@@ -2,7 +2,7 @@
 #'
 #'
 #' \code{edge_dist} returns edge lists defined by a spatial distance within the
-#' user defined threshold. The function accepts a \code{data.table} with
+#' user defined threshold. The function expects a \code{data.table} with
 #' relocation data, individual identifiers and a threshold argument. The
 #' threshold argument is used to specify the criteria for distance between
 #' points which defines a group. Relocation data should be in two columns
diff --git a/R/edge_nn.R b/R/edge_nn.R
index a96f2818..c89bcd7e 100644
--- a/R/edge_nn.R
+++ b/R/edge_nn.R
@@ -2,7 +2,7 @@
 #'
 #'
 #' \code{edge_nn} returns edge lists defined by the nearest neighbour. The
-#' function accepts a \code{data.table} with relocation data, individual
+#' function expects a \code{data.table} with relocation data, individual
 #' identifiers and a threshold argument. The threshold argument is used to
 #' specify the criteria for distance between points which defines a group.
 #' Relocation data should be in two columns representing the X and Y
diff --git a/R/get_gbi.R b/R/get_gbi.R
index 792d35ad..05d02dbe 100644
--- a/R/get_gbi.R
+++ b/R/get_gbi.R
@@ -1,7 +1,7 @@
 #' Generate group by individual matrix
 #'
 #'
-#' \code{get_gbi} generates a group by individual matrix. The function accepts a
+#' \code{get_gbi} generates a group by individual matrix. The function expects a
 #' \code{data.table} with individual identifiers and a group column. The group
 #' by individual matrix can then be used to build a network using
 #' \code{\link[asnipe:get_network]{asnipe::get_network}}.
diff --git a/R/group_lines.R b/R/group_lines.R
index ece668d9..9336e43c 100644
--- a/R/group_lines.R
+++ b/R/group_lines.R
@@ -1,7 +1,7 @@
 #' Group Lines
 #'
 #' `group_lines` groups rows into spatial groups by generating LINESTRINGs and
-#' grouping based on spatial intersection. The function accepts a `data.table`
+#' grouping based on spatial intersection. The function expects a `data.table`
 #' with relocation data, individual identifiers and a distance threshold. The
 #' relocation data is transformed into sf LINESTRINGs using [build_lines] and
 #' intersecting LINESTRINGs are grouped. The threshold argument is used to
diff --git a/R/group_polys.R b/R/group_polys.R
index d82cbd6f..17db4582 100644
--- a/R/group_polys.R
+++ b/R/group_polys.R
@@ -1,7 +1,7 @@
 #' Group Polygons
 #'
 #' `group_polys` groups rows into spatial groups by overlapping polygons (home
-#' ranges). The function accepts a `data.table` with relocation data, individual
+#' ranges). The function expects a `data.table` with relocation data, individual
 #' identifiers and an `area` argument.  The relocation data is transformed into
 #' home range POLYGONs using [build_polys()] with [adehabitatHR::mcp] or
 #' [adehabitatHR::kernelUD]. If the `area` argument is `FALSE`, `group_polys`
diff --git a/R/group_pts.R b/R/group_pts.R
index 3b4065ca..90b25545 100644
--- a/R/group_pts.R
+++ b/R/group_pts.R
@@ -1,6 +1,6 @@
 #' Group Points
 #'
-#' \code{group_pts} groups rows into spatial groups. The function accepts a
+#' \code{group_pts} groups rows into spatial groups. The function expects a
 #' \code{data.table} with relocation data, individual identifiers and a
 #' threshold argument. The threshold argument is used to specify the criteria
 #' for distance between points which defines a group. Relocation data should be
diff --git a/R/group_times.R b/R/group_times.R
index 4fb19e96..7db4fe5b 100644
--- a/R/group_times.R
+++ b/R/group_times.R
@@ -1,6 +1,6 @@
 #' Group Times
 #'
-#' \code{group_times} groups rows into time groups. The function accepts date
+#' \code{group_times} groups rows into time groups. The function expects date
 #' time formatted data and a threshold argument. The threshold argument is used
 #' to specify a time window within which rows are grouped.
 #'
diff --git a/R/randomizations.R b/R/randomizations.R
index e1033402..353a477d 100644
--- a/R/randomizations.R
+++ b/R/randomizations.R
@@ -1,7 +1,7 @@
 #' Data-stream randomizations
 #'
 #' \code{randomizations} performs data-stream social network randomization. The
-#' function accepts a \code{data.table} with relocation data, individual
+#' function expects a \code{data.table} with relocation data, individual
 #' identifiers and a randomization \code{type}. The \code{data.table} is
 #' randomized either using \code{step} or \code{daily} between-individual
 #' methods, or within-individual daily \code{trajectory} method described by
diff --git a/man/build_lines.Rd b/man/build_lines.Rd
index 1cbb7455..dcaf5593 100644
--- a/man/build_lines.Rd
+++ b/man/build_lines.Rd
@@ -41,7 +41,7 @@ build a line.
 }
 \description{
 \code{build_lines} generates a simple feature collection with LINESTRINGs from a
-\code{data.table}. The function accepts a \code{data.table} with relocation data,
+\code{data.table}. The function expects a \code{data.table} with relocation data,
 individual identifiers, a sorting column and a \code{projection}. The relocation
 data is transformed into LINESTRINGs for each individual and, optionally,
 combination of columns listed in \code{splitBy}. Relocation data should be in two
diff --git a/man/build_polys.Rd b/man/build_polys.Rd
index 53a76ff1..1a1651f5 100644
--- a/man/build_polys.Rd
+++ b/man/build_polys.Rd
@@ -47,7 +47,7 @@ of the respective \code{hrType} \code{adehabitatHR} function.
 }
 \description{
 \code{build_polys} generates a simple feature collection with POLYGONs from a
-\code{data.table}. The function accepts a \code{data.table} with
+\code{data.table}. The function expects a \code{data.table} with
 relocation data, individual identifiers, a projection,
 home range type and parameters. The relocation
 data is transformed into POLYGONs using either \link[adehabitatHR:mcp]{adehabitatHR::mcp} or
diff --git a/man/direction_group.Rd b/man/direction_group.Rd
index 1797176f..839fc7b5 100644
--- a/man/direction_group.Rd
+++ b/man/direction_group.Rd
@@ -86,6 +86,7 @@ See example of using mean group direction:
 \code{\link[CircStats:circ.mean]{CircStats::circ.mean()}}
 
 Other Direction functions: 
+\code{\link{direction_polarization}()},
 \code{\link{direction_step}()}
 }
 \concept{Direction functions}
diff --git a/man/direction_polarization.Rd b/man/direction_polarization.Rd
new file mode 100644
index 00000000..b4832f6a
--- /dev/null
+++ b/man/direction_polarization.Rd
@@ -0,0 +1,92 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/direction_polarization.R
+\name{direction_polarization}
+\alias{direction_polarization}
+\title{Polarization}
+\usage{
+direction_polarization(DT, direction = "direction", group = "group")
+}
+\arguments{
+\item{DT}{input data.table with distance column generated by
+\code{distance_step} and group column generated with \code{group_pts}}
+
+\item{direction}{character string of direction column name, default
+"direction"}
+
+\item{group}{character string of group column name, default "group"}
+}
+\value{
+\code{direction_polarization} returns the input \code{DT} appended
+with a \code{polarization} column representing the direction polarization
+of all individuals in each spatiotemporal group.
+
+The direction polarization is calculated using \code{\link[CircStats:r.test]{CircStats::r.test()}}
+which expects units of radians.
+
+A message is returned when the \code{polarization} columns already
+exists in the input \code{DT}, because it will be overwritten.
+}
+\description{
+\code{direction_polarization} calculates the polarization of individual
+directions in each spatiotemporal group identified by \code{group_pts}. The
+function expects a \code{data.table} with relocation data appended with a
+\code{direction} column from \code{direction_step} and a \code{group} column
+from \code{group_pts}.
+}
+\details{
+The \code{DT} must be a \code{data.table}. If your data is a
+\code{data.frame}, you can convert it by reference using
+\code{\link[data.table:setDT]{data.table::setDT}} or by reassigning using
+\code{\link[data.table:data.table]{data.table::data.table}}.
+
+The \code{direction} and \code{group} arguments expect the names of columns
+in \code{DT} which correspond to the direction and group columns. The
+direction column is expected in units of radians and the polarization is
+calculated with \code{\link[CircStats:r.test]{CircStats::r.test()}}.
+}
+\examples{
+# Load data.table
+library(data.table)
+\dontshow{data.table::setDTthreads(1)}
+
+# Read example data
+DT <- fread(system.file("extdata", "DT.csv", package = "spatsoc"))
+
+# Cast the character column to POSIXct
+DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+
+# Temporal grouping
+group_times(DT, datetime = 'datetime', threshold = '20 minutes')
+
+# Spatial grouping with timegroup
+group_pts(DT, threshold = 50, id = 'ID',
+          coords = c('X', 'Y'), timegroup = 'timegroup')
+
+# Calculate direction at each step
+direction_step(
+  DT = DT,
+  id = 'ID',
+  coords = c('X', 'Y'),
+  projection = 32736
+)
+
+# Calculate polarization
+direction_polarization(DT)
+}
+\references{
+See example of using polarization:
+\itemize{
+\item \url{https://doi.org/10.1016/j.cub.2017.08.004}
+\item <10.1371/journal.pcbi.1009437>
+\item \url{https://doi.org/10.7554/eLife.19505}
+}
+}
+\seealso{
+\code{\link{direction_step}}, \code{\link{group_pts}},
+\code{\link[CircStats:r.test]{CircStats::r.test()}}
+
+Other Direction functions: 
+\code{\link{direction_group}()},
+\code{\link{direction_step}()}
+}
+\concept{Direction functions}
diff --git a/man/direction_step.Rd b/man/direction_step.Rd
index 0b6df87c..07a56fba 100644
--- a/man/direction_step.Rd
+++ b/man/direction_step.Rd
@@ -99,6 +99,7 @@ direction_step(
 \code{\link[amt:steps]{amt::direction_abs()}}, \code{\link[geosphere:bearing]{geosphere::bearing()}}
 
 Other Direction functions: 
-\code{\link{direction_group}()}
+\code{\link{direction_group}()},
+\code{\link{direction_polarization}()}
 }
 \concept{Direction functions}
diff --git a/man/edge_dist.Rd b/man/edge_dist.Rd
index df362f3f..c9e3ba20 100644
--- a/man/edge_dist.Rd
+++ b/man/edge_dist.Rd
@@ -53,7 +53,7 @@ temporal with \code{group_times}) thresholds.
 }
 \description{
 \code{edge_dist} returns edge lists defined by a spatial distance within the
-user defined threshold. The function accepts a \code{data.table} with
+user defined threshold. The function expects a \code{data.table} with
 relocation data, individual identifiers and a threshold argument. The
 threshold argument is used to specify the criteria for distance between
 points which defines a group. Relocation data should be in two columns
diff --git a/man/edge_nn.Rd b/man/edge_nn.Rd
index 71e4a676..a6794a0b 100644
--- a/man/edge_nn.Rd
+++ b/man/edge_nn.Rd
@@ -50,7 +50,7 @@ neighbour.
 }
 \description{
 \code{edge_nn} returns edge lists defined by the nearest neighbour. The
-function accepts a \code{data.table} with relocation data, individual
+function expects a \code{data.table} with relocation data, individual
 identifiers and a threshold argument. The threshold argument is used to
 specify the criteria for distance between points which defines a group.
 Relocation data should be in two columns representing the X and Y
diff --git a/man/get_gbi.Rd b/man/get_gbi.Rd
index d9bbb6ce..cef044f5 100644
--- a/man/get_gbi.Rd
+++ b/man/get_gbi.Rd
@@ -25,7 +25,7 @@ but is more efficient thanks to
 \code{\link[data.table:dcast.data.table]{data.table::dcast}}.
 }
 \description{
-\code{get_gbi} generates a group by individual matrix. The function accepts a
+\code{get_gbi} generates a group by individual matrix. The function expects a
 \code{data.table} with individual identifiers and a group column. The group
 by individual matrix can then be used to build a network using
 \code{\link[asnipe:get_network]{asnipe::get_network}}.
diff --git a/man/group_lines.Rd b/man/group_lines.Rd
index 0803fbc2..531df3f3 100644
--- a/man/group_lines.Rd
+++ b/man/group_lines.Rd
@@ -62,7 +62,7 @@ input \code{DT}, because it will be overwritten.
 }
 \description{
 \code{group_lines} groups rows into spatial groups by generating LINESTRINGs and
-grouping based on spatial intersection. The function accepts a \code{data.table}
+grouping based on spatial intersection. The function expects a \code{data.table}
 with relocation data, individual identifiers and a distance threshold. The
 relocation data is transformed into sf LINESTRINGs using \link{build_lines} and
 intersecting LINESTRINGs are grouped. The threshold argument is used to
diff --git a/man/group_polys.Rd b/man/group_polys.Rd
index 6643b7c8..029871b8 100644
--- a/man/group_polys.Rd
+++ b/man/group_polys.Rd
@@ -61,7 +61,7 @@ through the \code{units} package.
 }
 \description{
 \code{group_polys} groups rows into spatial groups by overlapping polygons (home
-ranges). The function accepts a \code{data.table} with relocation data, individual
+ranges). The function expects a \code{data.table} with relocation data, individual
 identifiers and an \code{area} argument.  The relocation data is transformed into
 home range POLYGONs using \code{\link[=build_polys]{build_polys()}} with \link[adehabitatHR:mcp]{adehabitatHR::mcp} or
 \link[adehabitatHR:kernelUD]{adehabitatHR::kernelUD}. If the \code{area} argument is \code{FALSE}, \code{group_polys}
diff --git a/man/group_pts.Rd b/man/group_pts.Rd
index 0015239f..0ab9cc1d 100644
--- a/man/group_pts.Rd
+++ b/man/group_pts.Rd
@@ -44,7 +44,7 @@ A message is returned when a column named \code{group} already exists in
 the input \code{DT}, because it will be overwritten.
 }
 \description{
-\code{group_pts} groups rows into spatial groups. The function accepts a
+\code{group_pts} groups rows into spatial groups. The function expects a
 \code{data.table} with relocation data, individual identifiers and a
 threshold argument. The threshold argument is used to specify the criteria
 for distance between points which defines a group. Relocation data should be
diff --git a/man/group_times.Rd b/man/group_times.Rd
index 9cfcd4aa..c84e5dff 100644
--- a/man/group_times.Rd
+++ b/man/group_times.Rd
@@ -39,7 +39,7 @@ A message is returned when any of these columns already exist in the input
 \code{DT}, because they will be overwritten.
 }
 \description{
-\code{group_times} groups rows into time groups. The function accepts date
+\code{group_times} groups rows into time groups. The function expects date
 time formatted data and a threshold argument. The threshold argument is used
 to specify a time window within which rows are grouped.
 }
diff --git a/man/randomizations.Rd b/man/randomizations.Rd
index ebd21334..164df6b9 100644
--- a/man/randomizations.Rd
+++ b/man/randomizations.Rd
@@ -56,7 +56,7 @@ day relocations are swapped to. }
 }
 \description{
 \code{randomizations} performs data-stream social network randomization. The
-function accepts a \code{data.table} with relocation data, individual
+function expects a \code{data.table} with relocation data, individual
 identifiers and a randomization \code{type}. The \code{data.table} is
 randomized either using \code{step} or \code{daily} between-individual
 methods, or within-individual daily \code{trajectory} method described by
diff --git a/tests/testthat/test-direction-polarization.R b/tests/testthat/test-direction-polarization.R
new file mode 100644
index 00000000..fd74de23
--- /dev/null
+++ b/tests/testthat/test-direction-polarization.R
@@ -0,0 +1,82 @@
+# Test direction_polarization
+context('test direction_polarization')
+
+library(spatsoc)
+
+DT <- fread('../testdata/DT.csv')
+id <- 'ID'
+datetime <- 'datetime'
+timethreshold <- '20 minutes'
+threshold <- 50
+coords <- c('X', 'Y')
+timegroup <- 'timegroup'
+group <- 'group'
+projection <- 32736
+
+DT[, datetime := as.POSIXct(datetime, tz = 'UTC')]
+group_times(DT, datetime = datetime, threshold = timethreshold)
+group_pts(DT, threshold = threshold, id = id, coords = coords, timegroup = timegroup)
+direction_step(DT, id = id, coords = coords, projection = projection)
+
+clean_DT <- copy(DT)
+
+test_that('DT is required', {
+  expect_error(direction_polarization(DT = NULL))
+})
+
+test_that('arguments required, otherwise error detected', {
+  expect_error(direction_polarization(DT, group = NULL),
+               'group column name required')
+  expect_error(direction_polarization(DT, direction = NULL),
+               'direction column name required')
+})
+
+test_that('column names must exist in DT', {
+  expect_error(direction_polarization(DT, direction = 'potato'),
+               'potato field')
+  expect_error(direction_polarization(DT, group = 'potato'),
+               'potato field')
+})
+
+test_that('radians expected else error', {
+  expect_error(direction_polarization(DT, direction = 'X'),
+               'direction_step')
+})
+
+test_that('direction expected numeric', {
+  expect_error(direction_polarization(DT, direction = 'ID'),
+               'direction must be numeric')
+})
+
+
+test_that('polarization column succesfully detected', {
+  copyDT <- copy(clean_DT)[, polarization := 1]
+  expect_message(
+    direction_polarization(copyDT),
+    'polarization column will be overwritten'
+  )
+})
+
+test_that('no rows are added to the result DT', {
+  copyDT <- copy(clean_DT)
+
+  expect_equal(nrow(copyDT),
+               nrow(direction_polarization(copyDT)))
+})
+
+test_that('one column added to the result DT', {
+  copyDT <- copy(clean_DT)
+
+  expect_equal(ncol(copyDT) + 1,
+               ncol(direction_polarization(DT)))
+})
+
+test_that('column added to the result DT is numeric, between 0-1', {
+  expect_type(direction_polarization(DT)$polarization, 'double')
+  expect_gte(min(direction_polarization(DT)$polarization, na.rm = TRUE), 0)
+  expect_lte(max(direction_polarization(DT)$polarization, na.rm = TRUE), 1)
+})
+
+test_that('returns a data.table', {
+  expect_s3_class(direction_polarization(DT), 'data.table')
+})