Feat/use calc direction

NEWS.md

@@ -1,5 +1,12 @@
 # Development version
 
+Major changes:
+
+* use new internal `calc_direction` function in and provide sf interface for
+`direction_step`, `direction_to_centroid`, `direction_to_leader` and 
+`edge_direction` ([PR 125](https://github.com/ropensci/spatsoc/pull/125))
+
+
 New experimental functions:
 
 * `get_geometry` helper function for setting up an input DT with a sfc geometry

R/direction_step.R

@@ -1,39 +1,51 @@
 #' Direction step
 #'
-#' `direction_step` calculates the direction of movement steps in radians.
-#' The function expects a `data.table` with relocation data and individual
+#' `direction_step` calculates the direction of movement steps in radians. The
+#' function expects a `data.table` with relocation data and individual
 #' identifiers. Relocation data should be in two columns representing the X and
-#' Y coordinates. Note the order of rows is not modified by this function and
+#' Y coordinates, or in a geometry column prepared by the helper function
+#' [get_geometry()]. Note the order of rows is not modified by this function and
 #' therefore users must be cautious to set it explicitly. See example for one
 #' approach to setting order of rows using a datetime field.
 #'
-#' The `DT` must be a `data.table`. If your data is a
-#' `data.frame`, you can convert it by reference using
-#' [data.table::setDT()] or by reassigning using
+#' The `DT` must be a `data.table`. If your data is a `data.frame`, you can
+#' convert it by reference using [data.table::setDT()] or by reassigning using
 #' [data.table::data.table()].
 #'
-#' The `id`, `coords`, and optional `splitBy` arguments expect
-#' the names of a column in `DT` which correspond to the individual
-#' identifier, X and Y coordinates, and additional grouping columns.
+#' The `id`, and optional `splitBy` arguments expect the names of a column in
+#' `DT` which correspond to the individual identifier and additional grouping
+#' columns.
 #'
-#' The `crs` argument expects a character string or numeric defining
-#' the coordinate reference system to be passed to [sf::st_crs]. For example,
-#' for UTM zone 36S (EPSG 32736), the crs argument is
-#' `crs = "EPSG:32736"` or `crs = 32736`. See
-#' <https://spatialreference.org> for a list of EPSG codes.
+#' See below under Interface for details on providing coordinates.
 #'
 #' The `splitBy` argument offers further control over grouping. If within
 #' your `DT`, you have distinct sampling periods for each individual, you
 #' can provide the column name(s) which identify them to `splitBy`. The
 #' direction calculation by `direction_step` will only consider rows within
 #' each `id` and `splitBy` subgroup.
 #'
+#' @section Interface:
+#'  Two interfaces are available for providing coordinates:
+#'
+#'  1. Provide `coords` and `crs`. The `coords` argument expects the names of
+#'  the X and Y coordinate columns. The `crs` argument expects a character
+#'  string or numeric defining the coordinate reference system to be passed to
+#'  [sf::st_crs]. For example, for UTM zone 36S (EPSG 32736), the crs argument
+#'  is `crs = "EPSG:32736"` or `crs = 32736`. See <https://spatialreference.org>
+#'  for a list of EPSG codes.
+#'  2. (New!) Provide `geometry`. The `geometry` argument allows the user to
+#'  supply a `geometry` column that represents the coordinates as a simple
+#'  feature geometry list column. This interface expects the user to prepare
+#'  their input DT with [get_geometry()]. To use this interface, leave the
+#'  `coords` and `crs` arguments `NULL`, and the default argument for `geometry`
+#'  ('geometry') will be used directly.
+#'
 #' @return `direction_step` returns the input `DT` appended with
-#'  a direction column with units set to radians using the `units`
+#'  a `direction` column with units set to radians using the `units`
 #'  package.
 #'
 #'   This column represents the azimuth between the sequence of points for
-#'   each individual computed using `lwgeom::st_geod_azimuth`. Note, the
+#'   each individual computed using [lwgeom::st_geod_azimuth()]. Note, the
 #'   order of points is not modified by this function and therefore it is
 #'   crucial the user sets the order of rows to their specific question
 #'   before using `direction_step`. In addition, the direction column
@@ -43,14 +55,23 @@
 #'   A message is returned when a direction column are already exists in
 #'   the input `DT`, because it will be overwritten.
 #'
+#'   An error is returned if there are any missing values in coordinates for
+#'   the focal individual or the group leader, as the underlying direction
+#'   function ([lwgeom::st_geod_azimuth()]) does not accept missing values. An
+#'   error is also returned if the `crs` argument is NULL when `coords` are
+#'   provided.
+#'
 #'   See details for appending outputs using modify-by-reference in the
 #'   [FAQ](https://docs.ropensci.org/spatsoc/articles/faq.html).
 #'
 #' @inheritParams group_pts
 #' @inheritParams build_polys
+#' @param geometry simple feature geometry list column name, generated by
+#'   [get_geometry()]. Default 'geometry', see details under Interface
 #'
 #' @family Direction functions
-#' @seealso [amt::direction_abs()], [geosphere::bearing()]
+#' @seealso [lwgeom::st_geod_azimuth()], [amt::direction_abs()],
+#' [geosphere::bearing()]
 #' @export
 #'
 #' @examples
@@ -75,6 +96,10 @@
 #'   crs = 32736
 #' )
 #'
+#' # Or: sfc interface
+#' get_geometry(DT, coords = c('X', 'Y'), crs = 32736)
+#' direction_step(DT, id = 'ID')
+#'
 #' # Example result for East, North, West, South steps
 #' example <- data.table(
 #'   X = c(0, 5, 5, 0, 0),
@@ -91,10 +116,11 @@ direction_step <- function(
     coords = NULL,
     crs = NULL,
     splitBy = NULL,
+    geometry = 'geometry',
     projection = NULL) {
 
   # due to NSE notes in R CMD check
-  direction <- NULL
+  geo <- x <- y <- NULL
 
   if (!is.null(projection)) {
     warning('projection argument is deprecated, setting crs = projection')
@@ -104,35 +130,49 @@ direction_step <- function(
   assert_not_null(DT)
   assert_is_data_table(DT)
   assert_not_null(id)
+  assert_are_colnames(DT, c(id, splitBy))
 
-  check_cols <- c(id, splitBy)
-  assert_are_colnames(DT, check_cols)
+  out_colname <- 'direction'
 
-  assert_not_null(crs)
+  if (is.null(coords)) {
+    if (!is.null(crs)) {
+      message('crs argument is ignored when coords are null, using geometry')
+    }
 
-  assert_are_colnames(DT, coords)
-  assert_length(coords, 2)
-  assert_col_inherits(DT, coords, 'numeric')
+    assert_are_colnames(DT, geometry, ', did you run get_geometry()?')
+    assert_col_inherits(DT, geometry, 'sfc_POINT')
 
-  if ('direction' %in% colnames(DT)) {
-    message('direction column will be overwritten by this function')
-    data.table::set(DT, j = 'direction', value = NULL)
-  }
+    if ('direction' %in% colnames(DT)) {
+      message('direction column will be overwritten by this function')
+      data.table::set(DT, j = 'direction', value = NULL)
+    }
 
-  if (sf::st_is_longlat(crs)) {
-    DT[, direction := c(
-      lwgeom::st_geod_azimuth(
-        sf::st_as_sf(.SD, coords = coords, crs = crs)),
-      units::set_units(NA, 'rad')),
-      by = c(id, splitBy)]
-  } else if (!sf::st_is_longlat(crs)) {
-    DT[, direction := c(
-      lwgeom::st_geod_azimuth(
-        sf::st_transform(
-          sf::st_as_sf(.SD, coords = coords, crs = crs),
-          crs = 4326)
-        ),
-      units::set_units(NA, 'rad')),
-      by = c(id, splitBy)]
+    DT[, c(out_colname) := c(calc_direction(
+      geometry_a = geo
+    ), units::set_units(NA, 'rad')),
+    by = c(id, splitBy),
+    env = list(geo = geometry)
+    ]
+
+  } else {
+    assert_are_colnames(DT, coords)
+    assert_length(coords, 2)
+    assert_col_inherits(DT, coords, 'numeric')
+    assert_not_null(crs)
+
+    if (out_colname %in% colnames(DT)) {
+      message(out_colname, ' column will be overwritten by this function')
+      data.table::set(DT, j = 'direction', value = NULL)
+    }
+
+    DT[, c(out_colname) := c(calc_direction(
+      x_a = x,
+      y_a = y,
+      crs = crs
+    ), units::set_units(NA, 'rad')),
+    by = c(id, splitBy),
+    env = list(x = data.table::first(coords), y = data.table::last(coords))
+    ]
   }
+
 }

R/direction_to_centroid.R

@@ -1,41 +1,47 @@
 #' Direction to group centroid
 #'
-#' `direction_to_centroid` calculates the direction of each relocation to
-#' the centroid of the spatiotemporal group identified by `group_pts`. The
-#' function expects a `data.table` with relocation data appended with a
-#' `group` column from `group_pts` and centroid columns from
-#' `centroid_group`. Relocation data should be in planar coordinates
-#' provided in two columns representing the X and Y coordinates.
-#'
-#' The `DT` must be a `data.table`. If your data is a
-#' `data.frame`, you can convert it by reference using
-#' [data.table::setDT()] or by reassigning using
+#' `direction_to_centroid` calculates the direction of each relocation to the
+#' centroid of the spatiotemporal group identified by `group_pts`. The function
+#' expects a `data.table` with relocation data appended with a `group` column
+#' from `group_pts` and centroid columns from `centroid_group`. Relocation data
+#' should be in two columns representing the X and Y coordinates, or in a
+#' geometry column prepared by the helper function [get_geometry()].
+#'
+#' The `DT` must be a `data.table`. If your data is a `data.frame`, you can
+#' convert it by reference using [data.table::setDT()] or by reassigning using
 #' [data.table::data.table()].
 #'
-#' This function expects a `group` column present generated with the
-#' `group_pts` function and centroid coordinate columns generated with the
-#' `centroid_group` function. The `coords` and `group` arguments
-#' expect the names of columns in `DT` which correspond to the X and Y
-#' coordinates and group columns.
+#' This function expects a `group` column present generated with the `group_pts`
+#' function and centroid coordinates generated with the `centroid_group`
+#' function. The `group` argument expects the name of the column in `DT` which
+#' correspond to the group column.
+#'
+#' See below under Interface for details on providing coordinates.
 #'
+#' @inheritSection direction_step Interface
 #' @inheritParams distance_to_centroid
+#' @inheritParams direction_step
 #' @inheritParams group_pts
 #'
-#' @return `direction_to_centroid` returns the input `DT` appended
-#'   with a `direction_centroid` column indicating the direction to group
-#'   centroid in radians. The direction is measured in radians in the range
-#'   of 0 to 2 * pi from the positive x-axis.
+#' @return `direction_to_centroid` returns the input `DT` appended with a
+#'   `direction_centroid` column indicating the direction to the group centroid
+#'   in radians. A value of NaN is returned when the coordinates of the focal
+#'   individual equal the coordinates of the centroid.
 #'
 #'   A message is returned when `direction_centroid` column already exist
 #'   in the input `DT`, because they will be overwritten.
 #'
+#'   An error is returned if there are any missing values in coordinates for
+#'   the focal individual or the group leader, as the underlying direction
+#'   function ([lwgeom::st_geod_azimuth()]) does not accept missing values.
+#'
 #'   See details for appending outputs using modify-by-reference in the
 #'   [FAQ](https://docs.ropensci.org/spatsoc/articles/faq.html).
 #'
 #' @export
 #' @family Direction functions
 #' @family Centroid functions
-#' @seealso [centroid_group], [group_pts]
+#' @seealso [centroid_group], [group_pts], [lwgeom::st_geod_azimuth()]
 #' @references
 #' See example of using direction to group centroid:
 #'  * \doi{doi:10.1016/j.cub.2017.08.004}
@@ -62,48 +68,72 @@
 #' centroid_group(DT, coords = c('X', 'Y'), group = 'group', na.rm = TRUE)
 #'
 #' # Calculate direction to group centroid
-#' direction_to_centroid(DT, coords = c('X', 'Y'))
+#' direction_to_centroid(DT, coords = c('X', 'Y'), crs = 32736)
 direction_to_centroid <- function(
     DT = NULL,
-    coords = NULL) {
+    coords = NULL,
+    crs = NULL,
+    geometry = 'geometry') {
 
   # Due to NSE notes in R CMD check
-  direction_centroid <- NULL
+  geo <- cent <- x <- y <- x_centroid <- y_centroid <- NULL
 
   assert_not_null(DT)
   assert_is_data_table(DT)
 
-  assert_are_colnames(DT, coords)
-  assert_length(coords, 2)
-  assert_col_inherits(DT, coords, 'numeric')
+  out_colname <- 'direction_centroid'
+
+  if (is.null(coords)) {
+    if (!is.null(crs)) {
+      message('crs argument is ignored when coords are null, using geometry')
+    }
+
+    assert_are_colnames(DT, geometry, ', did you run get_geometry()?')
+    assert_col_inherits(DT, geometry, 'sfc_POINT')
+    centroid_col <- 'centroid'
+    assert_are_colnames(DT, centroid_col, ', did you run centroid_group?')
+    assert_col_inherits(DT, centroid_col, 'sfc_POINT')
+
+    if (out_colname %in% colnames(DT)) {
+      message(out_colname, ' column will be overwritten by this function')
+      data.table::set(DT, j = out_colname, value = NULL)
+    }
+
+    DT[, c(out_colname) := calc_direction(
+      geometry_a = geo,
+      geometry_b = cent
+    ),
+    env = list(geo = geometry, cent = centroid_col)]
+
+  } else {
+    assert_are_colnames(DT, coords)
+    assert_length(coords, 2)
+    assert_col_inherits(DT, coords, 'numeric')
+    assert_not_null(crs)
+
+    xcol <- data.table::first(coords)
+    ycol <- data.table::last(coords)
+    pre <- 'centroid_'
+    xcol_centroid <- paste0(pre, xcol)
+    ycol_centroid <- paste0(pre, ycol)
+    coords_centroid  <- c(xcol_centroid, ycol_centroid)
 
-  xcol <- data.table::first(coords)
-  ycol <- data.table::last(coords)
-  pre <- 'centroid_'
-  centroid_xcol <- paste0(pre, xcol)
-  centroid_ycol <- paste0(pre, ycol)
-  centroid_coords  <- c(centroid_xcol, centroid_ycol)
+    assert_are_colnames(DT, coords_centroid, ', did you run centroid_group?')
+    assert_col_inherits(DT, coords_centroid, 'numeric')
 
-  assert_are_colnames(DT, centroid_coords, ', did you run centroid_group?')
-  assert_col_inherits(DT, centroid_coords, 'numeric')
+    if (out_colname %in% colnames(DT)) {
+      message(out_colname, ' column will be overwritten by this function')
+      data.table::set(DT, j = out_colname, value = NULL)
+    }
 
-  if ('direction_centroid' %in% colnames(DT)) {
-    message('direction_centroid column will be overwritten by this function')
-    data.table::set(DT, j = 'direction_centroid', value = NULL)
+    DT[, c(out_colname) := calc_direction(
+      x_a = x, y_a = y,
+      x_b = x_centroid, y_b = y_centroid,
+      crs = crs
+    ),
+    env = list(x = xcol, y = ycol,
+               x_centroid = xcol_centroid, y_centroid = ycol_centroid)]
   }
-
-  DT[, direction_centroid := fifelse(
-    .SD[[xcol]] == .SD[[centroid_xcol]] &
-      .SD[[ycol]] == .SD[[centroid_ycol]],
-    units::as_units(NaN, 'rad'),
-    units::as_units(
-      atan2(.SD[[centroid_ycol]] - .SD[[ycol]],
-            (.SD[[centroid_xcol]] - .SD[[xcol]])),
-      'rad'
-    )
-  )]
-  DT[direction_centroid < units::as_units(0, 'rad'),
-     direction_centroid := direction_centroid + units::as_units(2 * pi, 'rad')]
 
   return(DT[])
 }