dc: The depth-contour (DC) algorithm

In edwardlavender/flapper: Routines for the Analysis of Passive Acoustic Telemetry Data

The depth-contour (DC) algorithm

Description

This function implements the depth-contour (DC) algorithm. This algorithm relates one-dimensional depth time series to a two-dimensional bathymetry surface to determine the extent to which different parts of an area might have (or have not) been used, or effectively represent occupied depths, over time. Given a sequence of depth observations (archival) from a tagged animal and a measurement error parameter (calc_depth_error), at each time step the function determines the cells on a bathymetry raster (bathy) that match the observed depth*. Across all time steps, matches are summed to produce a single map representing the number of occasions when the depth in each cell matched the observed depth.

Usage

dc(
  archival,
  bathy,
  plot_ts = TRUE,
  calc_depth_error = function(...) matrix(c(-2.5, 2.5), nrow = 2),
  check_availability = TRUE,
  normalise = TRUE,
  save_record_spatial = 1L,
  write_record_spatial_for_pf = NULL,
  save_args = TRUE,
  verbose = TRUE,
  con = "",
  split = NULL,
  cl = NULL,
  varlist = NULL
)

Arguments

archival	A dataframe of depth time series (for a single individual). At a minimum, this should contain a column named ‘depth’ with depth observations. Depth should be recorded using absolute values in the same units as the bathymetry (bathy, see below).
bathy	A raster of the bathymetry in an area within which the animal is likely to have been located over the study. Bathymetry values should be recorded as absolute values and in the same units as for depths (see archival).
plot_ts	A logical input that defines whether or not to the depth time series before the algorithm is initiated.
calc_depth_error	A function that returns the depth errors around a vector of depths. The function should accept vector of depths (from archival$depth) and return a matrix, with one row for each (lower and upper) error and one one column for each depth (if the error varies with depth). For each depth, the two numbers are added to the observed depth to define the range of depths on the bathymetry raster (bathy) that the individual could plausibly have occupied at any time. Since the depth errors are added to the individual's depth, the first number should be negative (i.e., the individual could have been shallower that observed) and the second positive (i.e., the individual could have been deeper than observed). The appropriate form for calc_depth_error depends on the species (pelagic versus demersal/benthic species), the measurement error for the depth observations in archival and bathymetry (bathy) data, as well as the tidal range (m) across the area (over the duration of observations). For example, for a pelagic species, the constant function calc_depth_error = function(...) matrix(c(-2.5, Inf) implies that the individual could have occupied bathymetric cells that are deeper than the observed depth + (-2.5) m and shallower than Inf m (i.e. the individual could have been in any location in which the depth was deeper than the shallow depth limit for the individual). In contrast, for a benthic species, the constant function calc_depth_error = function(...) matrix(c(-2.5, 2.5), nrow = 2) implies that the individual could have occupied bathymetric cells whose depth lies within the interval defined by the observed depth + (-2.5) and + (+2.5) m.
check_availability	A logical input that defines whether or not to record explicitly, for each time step, whether or not there were any cells on bathy that matched the observed depth (within the bounds defined by calc_depth_error).
normalise	A logical variable that defines whether or not to normalise the map of possible locations at each time step so that their scores sum to one.
save_record_spatial	An integer vector that defines the time steps for which to return time step-specific and cumulative maps of the individual's possible locations. save_record_spatial = 0 suppresses the return of this information and save_record_spatial = NULL returns this information for all time steps.
write_record_spatial_for_pf	(optional) A named list, passed to writeRaster, to save the raster of the individual's possible positions at each time step to file. The ‘filename’ argument should be the directory in which to save files. Files are named by archival time steps as ‘arc_1’, ‘arc_2’ and so on.
save_args	A logical input that defines whether or not to save the list of function inputs in the returned object.
verbose	A logical variable that defines whether or not to print messages to the console or to file to relay function progress. If con = "", messages are printed to the console; otherwise, they are written to file (see below).
con	If verbose = TRUE, con is character string that defines how messages relaying function progress are returned. If con = "", messages are printed to the console (unless redirected by sink). Otherwise, con defines the full pathway to a .txt file (which can be created on-the-fly) into which messages are written to relay function progress.
split, cl, varlist	(optional) Parallelisation options. split is an integer which, if supplied, splits the archival dataframe every nth row into chunks†. The algorithm is applied sequentially within each chunk (if applicable) and chunk-wise maps can be summed afterwards to create a single map of space use. The advantage of this approach is that chunks can be analysed in parallel, via cl and varlist, while memory use is minimised. cl is (a) a cluster object from makeCluster or (b) an integer that defines the number of child processes. varlist is a character vector of variables for export (see cl_export). Exported variables must be located in the global environment. If a cluster is supplied, the connection to the cluster is closed within the function (see cl_stop). For further information, see cl_lapply and flapper-tips-parallel.

Details

This algorithm can be applied to both pelagic species (which must be in an area where the depth is at least as deep as the observed depth) and benthic/demersal species (which must be in an area relatively close in depth to the observed depth). However, the results are likely to be much more precise for benthic species because depth observations for benthic animals are typically more informative about a tagged animal's possible location(s) than for pelagic species.

*Location probability could be weighted by the proximity between the observed depth and the depths of cells within the range defined by the observed depth and the measurement error (i.e., archival$depth[1] + calc_depth_error(archival$depth[1])[1], archival$depth[1] + calc_depth_error(archival$depth[1])[2]), but this is not currently implemented.

†Under the default options (split = NULL), the function starts with a blank map of the area and iterates over each time step, adding the ‘possible positions’ of the individual to the map at each step. By continuously updating a single map, this approach is slow but minimises memory requirements. An alternative approach is to split the time series into chunks, implement an iterative approach within each chunk, and then join the maps for each chunk. This is implemented by split (plus acdc_simplify).

Value

The function returns an acdc_archive-class object. If a connection to write files has also been specified, messages are also written to file relaying function progress.

Author(s)

Edward Lavender

See Also

acdc_simplify simplifies the outputs of the algorithm. acdc_plot_trace, acdc_plot_record and acdc_animate_record provide plotting routines. dcq implements a faster version of this algorithm termed the ‘quick depth-contour’ (DCQ) algorithm. Rather than considering the depth interval that the individual could have occupied at each time step, the DCQ algorithm considers a sequence of depth bins (e.g., 10 m bins), isolates these on the bathymetry raster (bathy) and counts the number of matches in each cell. The DCPF algorithm (see pf) extends the DC algorithm via particle filtering to reconstruct possible movement paths over bathy. The ACDC algorithm (see acdc) extends the depth-contour algorithm by integrating information from acoustic detections of individuals at each time step to restrict the locations in which depth contours are identified.

Examples

#### Define depth time series for examples
# We will use a sample depth time series for one individual
# We will select a small sample of observations for example speed
depth <- dat_archival[dat_archival$individual_id == 25, ][1:100, ]

#### Example (1): Implement algorithm with default options
## Implement algorithm
out_dc <- dc(archival = depth, bathy = dat_gebco)
## dc() returns a named list that is an 'acdc_archive' class object
summary(out_dc)
# ... archive elements contain the key information from the algorithm
# ... ts_by_chunk contains the time series for each chunk
# ... time contains a dataframe of the algorithm's progression
# ... args contains a list of user inputs
## Simplify outputs via acdc_simplify()
dc_summary <- acdc_simplify(out_dc, type = "dc")
summary(dc_summary)
## Examine time-specific maps via acdc_plot_record()
acdc_plot_record(dc_summary)
## Examine overall map via a raster* plotting function
# Each cell shows expected proportion of time steps spent in each cell
prettyGraphics::pretty_map(
  add_rasters = list(x = dc_summary$map),
  add_polys = list(x = dat_coast)
)
# Convert counts on map to percentages
prettyGraphics::pretty_map(
  add_rasters =
    list(
      x = dc_summary$map / nrow(depth) * 100,
      zlim = c(0, 100)
    ),
  add_polys = list(x = dat_coast)
)
# Check for occasions when the individual's depth was not consistent
# ... with the depth data for the area and the depth error e.g., possibly
# ... due to movement beyond this area:
any(do.call(rbind, lapply(dc_summary$record, function(elm) elm$dat))$availability == FALSE)

#### Example (2): Implement depth error functions that depend on depth
## Here, we will define a calc_depth_error function that depends on:
# ... tag error
# ... tidal range
# ... a depth-dependent bathymetry error
cde <- function(depth) {
  e <- 4.77 + 2.5 + sqrt(0.5^2 + (0.013 * depth)^2)
  e <- matrix(c(-e, e), nrow = 2)
  return(e)
}
# Vectorise function over depths (essential)
cde <- Vectorize(cde)
## Implement algorithm with  depth-dependent error
out_dc <- dc(archival = depth, bathy = dat_gebco, calc_depth_error = cde)

#### Example (3): Write timestep-specific maps of allowed positions to file
out_dc <- dc(
  archival = depth, bathy = dat_gebco,
  write_record_spatial_for_pf = list(filename = tempdir())
)
list.files(tempdir())

#### Example (4): Implement the algorithm in parallel
## Trial different options for 'split' and compare speed
# Approach (1)
at1 <- Sys.time()
cl <- parallel::makeCluster(2L)
parallel::clusterEvalQ(cl = cl, library(raster))
out_dc <- dc(
  archival = depth,
  bathy = dat_gebco,
  plot_ts = FALSE,
  split = 1,
  cl = cl
)
at2 <- Sys.time()
# Approach (2)
bt1 <- Sys.time()
cl <- parallel::makeCluster(2L)
parallel::clusterEvalQ(cl = cl, library(raster))
out_dc <- dc(
  archival = depth,
  bathy = dat_gebco,
  plot_ts = FALSE,
  split = 5,
  cl = cl
)
bt2 <- Sys.time()
# Compare timings
difftime(at2, at1)
difftime(bt2, bt1)

#### Example (5): Write messages to file via con
out_dc <- dc(
  archival = depth, bathy = dat_gebco,
  con = paste0(tempdir(), "/dc_log.txt")
)
readLines(paste0(tempdir(), "/dc_log.txt"))

#### Example (6) Compare an automated vs. manual implementation of DC

## (A) Compare step-wise results for randomly select time steps
# Implement algorithm (using default calc_depth_error)
out_dc <- dc(
  archival = depth,
  bathy = dat_gebco,
  save_record_spatial = NULL
)
# Compare results for randomly selected time steps to manual implementation
# ... with the same simple calc_depth_error model
pp <- graphics::par(mfrow = c(5, 2))
for (i in sample(1:nrow(depth), 5)) {
  # Extract map created via dc()
  raster::plot(out_dc$archive[[1]]$record[[i]]$spatial[[1]]$map_timestep)
  # Compare to manual creation of map
  raster::plot(dat_gebco >= depth$depth[i] - 2.5 & dat_gebco <= depth$depth[i] + 2.5)
}
graphics::par(pp)

## (B) Compare the chunk-wise results and manual implementation
# Implement algorithm chunk-wise
out_dc <- dc(
  archival = depth,
  bathy = dat_gebco,
  save_record_spatial = NULL,
  split = 10L
)
# Get overall map via acdc_simplify()
map_1 <- acdc_simplify(out_dc, type = "dc", mask = dat_gebco)$map
# Get overall map via a simple custom approach
map_2 <- raster::setValues(dat_gebco, 0)
for (i in 1:nrow(depth)) {
  map_2 <- map_2 + (dat_gebco >= depth$depth[i] - 2.5 & dat_gebco <= depth$depth[i] + 2.5)
}
map_2 <- raster::mask(map_2, dat_gebco)
# Show that the two overall maps are identical
pp <- par(mfrow = c(1, 3))
raster::plot(map_1, main = "Automated")
raster::plot(map_2, main = "Custom")
raster::plot(map_1 - map_2, main = "Difference")
par(pp)

edwardlavender/flapper documentation built on Jan. 22, 2025, 2:44 p.m.