DistanceDevelopment · hambrecht · May 5, 2026 · May 5, 2026 · May 5, 2026 · May 5, 2026
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -5,6 +5,7 @@ Imports:
     mrds,
     Distance,
     sf,
+    terra,
     ggplot2,
     purrr,
     dplyr,

diff --git a/R/ClassConstructors.R b/R/ClassConstructors.R
@@ -270,6 +270,22 @@ make.population.description <- make.pop.description <- function(region = make.re
 #' @param scale.param numeric vector with either a single value to be applied globally or a value for each strata. These should be supplied on the natural scale.
 #' @param shape.param numeric vector with either a single value to be applied globally or a value for each strata. These should be supplied on the natural scale.
 #' @param cov.param Named list with one named entry per individual level covariate. Covariate parameter values should be defined on the log scale (rather than the natural scale), this is the same scale as provided in the ddf output in mrds and also in the MCDS output in Distance. Cluster sizes parameter values can be defined here. Each list entry will either be a data.frame containing 2 or 3 columns: level, param and where desired strata. If the region has multiple strata but this column is omitted then the values will be assumed to apply globally. The cluster size entry in the list must be named 'size'. Alternatively the list element may a numeric vector with either a single value to be applied globally or a value for each strata.
+#' @param cov.surface Optional named list of raster surfaces providing spatially-explicit
+#'   covariate values based on each animal's simulated location. Each element name must
+#'   match a numeric entry in \code{cov.param} (which gives the log-scale slope applied to
+#'   the extracted raster value). Each element should be either a character file path to a
+#'   single-layer raster readable by \code{terra::rast()}, or an in-memory
+#'   \code{terra} SpatRaster object. At every simulation replicate the raster is sampled at
+#'   each animal's (x, y) location; the value is used as the covariate in the standard
+#'   log-linear scaling formula:
+#'   \eqn{\log(\sigma_i) = \log(\sigma_0) + \beta \cdot z_i}.
+#'   Use a negative \code{cov.param} slope to model lower detectability where raster values
+#'   are higher (e.g. denser canopy cover reduces detection range).
+#'   Animal locations outside the raster extent or with NA cell values are replaced with the
+#'   raster mean and a warning is issued.
+#'   \strong{Parallel note:} file paths are strongly preferred over in-memory SpatRaster
+#'   objects when using \code{run.parallel = TRUE}, as large objects may not serialise
+#'   efficiently to worker processes.
 #' @param truncation the maximum perpendicular (or radial) distance at which
 #'   objects may be detected from a line (or point) transect.
 #' @return \code{\link{Detectability-class}} object
@@ -323,9 +339,9 @@ make.population.description <- make.pop.description <- function(region = make.re
 #'
 #' plot(detect, popdesc)
 #'
-make.detectability <- function(key.function = "hn", scale.param = 25, shape.param = numeric(0), cov.param = list(), truncation = 50){
+make.detectability <- function(key.function = "hn", scale.param = 25, shape.param = numeric(0), cov.param = list(), cov.surface = list(), truncation = 50){
   # Passes all arguments to function to make a new instance of the class
-  detectability <- new(Class = "Detectability", key.function = key.function, scale.param = scale.param, shape.param = shape.param, cov.param = cov.param, truncation = truncation)
+  detectability <- new(Class = "Detectability", key.function = key.function, scale.param = scale.param, shape.param = shape.param, cov.param = cov.param, cov.surface = cov.surface, truncation = truncation)
   return(detectability)
 }
 

diff --git a/R/Detectability.R b/R/Detectability.R
@@ -15,8 +15,12 @@
 #' parameter for the detection function.
 #' @slot shape.param Object of class \code{"numeric"}; The shape
 #' parameter for the detection function.
-#' @slot cov.param Object of class \code{"numeric"}; The parameter
-#' values associated with the covariates. Not yet implemented
+#' @slot cov.param Object of class \code{"list"}; Named list of slope
+#' coefficients (log scale) for individual-level covariates.
+#' @slot cov.surface Object of class \code{"list"}; Optional named list of
+#' raster surfaces (character file paths or terra SpatRaster objects) providing
+#' spatially-explicit covariate values sampled at each animal's location. Names
+#' must match numeric entries in \code{cov.param}.
 #' @slot truncation Object of class \code{"numeric"}; The maximum
 #' distance at which objects may be detected.
 #' @keywords classes
@@ -26,12 +30,13 @@ setClass("Detectability", representation(key.function    = "character",
                                          scale.param     = "numeric",
                                          shape.param     = "numeric",
                                          cov.param       = "list",
+                                         cov.surface     = "list",
                                          truncation      = "numeric"))
 #' @importFrom methods validObject is
 setMethod(
   f="initialize",
   signature="Detectability",
-  definition=function(.Object, key.function = "hn", scale.param = 25, shape.param = numeric(0), covariates = character(0), cov.param = list(), truncation = 50){
+  definition=function(.Object, key.function = "hn", scale.param = 25, shape.param = numeric(0), covariates = character(0), cov.param = list(), cov.surface = list(), truncation = 50){
     #Input pre-processing
     # - this needs to be done here as cannot alter object inside validation method
     cov.names <- names(cov.param)
@@ -55,6 +60,7 @@ setMethod(
     .Object@scale.param  <- scale.param
     .Object@shape.param  <- shape.param
     .Object@cov.param    <- cov.param
+    .Object@cov.surface  <- cov.surface
     .Object@truncation   <- truncation
     #Check object is valid
     valid <- validObject(.Object, test = TRUE)
@@ -112,6 +118,35 @@ setValidity("Detectability",
                   }
                 }
               }
+              # Validate cov.surface
+              if(length(object@cov.surface) > 0){
+                surf.names <- names(object@cov.surface)
+                if(is.null(surf.names) || any(surf.names == "")){
+                  return("Not all elements of cov.surface are named. Please provide names matching entries in cov.param.")
+                }
+                cov.param.names <- names(object@cov.param)
+                for(sn in surf.names){
+                  if(!sn %in% cov.param.names){
+                    return(paste0("cov.surface entry '", sn, "' has no matching numeric entry in cov.param. ",
+                                  "Each raster surface needs a corresponding slope in cov.param."))
+                  }
+                  if(is.data.frame(object@cov.param[[sn]])){
+                    return(paste0("cov.param entry '", sn, "' paired with cov.surface must be a numeric slope, not a factor data.frame."))
+                  }
+                  surf <- object@cov.surface[[sn]]
+                  if(!is.character(surf) && !inherits(surf, "SpatRaster")){
+                    return(paste0("cov.surface entry '", sn, "' must be a file path (character) or a terra SpatRaster object."))
+                  }
+                  if(is.character(surf)){
+                    if(length(surf) != 1){
+                      return(paste0("cov.surface entry '", sn, "' must be a single file path."))
+                    }
+                    if(!file.exists(surf)){
+                      return(paste0("cov.surface raster file does not exist: '", surf, "'."))
+                    }
+                  }
+                }
+              }
               return(TRUE)
             }
 )
@@ -159,7 +194,9 @@ setMethod(
     cov.names <- names(object@cov.param)
     # Check if there are covariates in detectability that are not in pop.desc
     pop.covs <- names(pop.desc@covariates)
-    if(any(!cov.names %in% pop.covs)){
+    surface.covs <- names(object@cov.surface)
+    missing.covs <- setdiff(cov.names, c(pop.covs, surface.covs))
+    if(length(missing.covs) > 0){
       stop("You have defined detectability for covariates that are not included in the population description.", call. = FALSE)
     }
     # set mfrow storing old settings
@@ -187,7 +224,21 @@ setMethod(
       for(cov in seq(along = object@cov.param)){
         cov.params <- object@cov.param[[cov]]
         cov.dist <- pop.desc@covariates[[cov.names[cov]]]
-        if(is(object@cov.param[[cov]], "data.frame")){
+        if(is.null(cov.dist) && cov.names[cov] %in% surface.covs){
+          param.type <- "surface"
+          no.cov.strata <- max(1, length(cov.params))
+          cov.surface <- object@cov.surface[[cov.names[cov]]]
+          if(is.character(cov.surface)){
+            cov.surface <- terra::rast(cov.surface)
+          }
+          surf.values <- terra::values(cov.surface)[,1]
+          surf.values <- surf.values[!is.na(surf.values)]
+          if(length(surf.values) == 0){
+            quantiles <- c(0, 0, 0)
+          }else{
+            quantiles <- as.numeric(quantile(surf.values, c(0.025, 0.5, 0.975)))
+          }
+        }else if(is(object@cov.param[[cov]], "data.frame")){
           param.type = "categorical"
           no.cov.strata <- ifelse(is.null(cov.params$strata), 1, length(unique(cov.params$strata)))
         }else if(is(cov.dist[[1]], "data.frame")){
@@ -201,6 +252,8 @@ setMethod(
           plot.title <- paste("Covariate: ", cov.names[cov], " (factor)", sep = "")
         }else if(param.type == "discrete"){
           plot.title <- paste("Covariate: ", cov.names[cov], " (discrete)", sep = "")
+        }else if(param.type == "surface"){
+          plot.title <- paste("Covariate: ", cov.names[cov], " (raster)", sep = "")
         }else{
           plot.title <- paste("Covariate: ", cov.names[cov], " (continuous)", sep = "")
         }
@@ -277,6 +330,9 @@ setMethod(
                                     "poisson" = qpois(int, dist.param$lambda),
                                     "lognormal" = qlnorm(int, dist.param$meanlog, dist.param$sdlog))
               }
+            }else if(param.type == "surface"){
+              # quantiles were already computed from the raster values
+              quantiles <- quantiles
             }
             # get adjustment paramters
             if(length(cov.params) == no.strata){
@@ -315,6 +371,8 @@ setMethod(
           desc.ints <- "min,max"
         }else if(param.type == "continuous"){
           desc.ints <- "95%ints"
+        }else if(param.type == "surface"){
+          desc.ints <- "raster 95%ints"
         }
         if(param.type == "categorical"){
           no.levels <- length(unique(cov.params$level))

diff --git a/R/accumulate.PP.results.R b/R/accumulate.PP.results.R
@@ -1,15 +1,61 @@
+extract.rep.result <- function(results, i){
+  rep.result <- list(
+    individuals = list(
+      summary = results$individuals$summary[,,i, drop = FALSE],
+      N = results$individuals$N[,,i, drop = FALSE],
+      D = results$individuals$D[,,i, drop = FALSE]
+    ),
+    Detection = results$Detection[,,i, drop = FALSE]
+  )
+
+  if(!is.null(results$clusters)){
+    rep.result$clusters <- list(
+      summary = results$clusters$summary[,,i, drop = FALSE],
+      N = results$clusters$N[,,i, drop = FALSE],
+      D = results$clusters$D[,,i, drop = FALSE]
+    )
+    rep.result$expected.size <- results$expected.size[,,i, drop = FALSE]
+  }
+
+  rep.result
+}
+
+apply.rep.result <- function(results, rep.result, i){
+  results$individuals$summary[,,i] <- rep.result$individuals$summary[,,1]
+  results$individuals$N[,,i] <- rep.result$individuals$N[,,1]
+  results$individuals$D[,,i] <- rep.result$individuals$D[,,1]
+  results$Detection[,,i] <- rep.result$Detection[,,1]
+  if(!is.null(results$clusters) && !is.null(rep.result$clusters)){
+    results$clusters$summary[,,i] <- rep.result$clusters$summary[,,1]
+    results$clusters$N[,,i] <- rep.result$clusters$N[,,1]
+    results$clusters$D[,,i] <- rep.result$clusters$D[,,1]
+    results$expected.size[,,i] <- rep.result$expected.size[,,1]
+  }
+  results
+}
+
 accumulate.PP.results <- function(simulation, results){
   simulation@results$filename <- rep(NA, length(results))
   for(i in seq(along = results)){
-    simulation@results$individuals$summary[,,i] <- results[[i]]$individuals$summary[,,i]
-    simulation@results$individuals$N[,,i] <- results[[i]]$individuals$N[,,i]
-    simulation@results$individuals$D[,,i] <- results[[i]]$individuals$D[,,i]
-    simulation@results$Detection[,,i] <- results[[i]]$Detection[,,i]
-    if(!is.null(simulation@results$clusters)){
-      simulation@results$clusters$summary[,,i] <- results[[i]]$clusters$summary[,,i]
-      simulation@results$clusters$N[,,i] <- results[[i]]$clusters$N[,,i]
-      simulation@results$clusters$D[,,i] <- results[[i]]$clusters$D[,,i]
-      simulation@results$expected.size[,,i] <- results[[i]]$expected.size[,,i]
+    if(!is.null(results[[i]]$rep.result)){
+      simulation@results <- apply.rep.result(simulation@results,
+                                             results[[i]]$rep.result,
+                                             i)
+      if(!is.null(results[[i]]$filename) && length(results[[i]]$filename) > 0){
+        simulation@results$filename[i] <- results[[i]]$filename
+      }
+    }else{
+      # Backward-compatible path for legacy full-object result returns
+      simulation@results$individuals$summary[,,i] <- results[[i]]$individuals$summary[,,i]
+      simulation@results$individuals$N[,,i] <- results[[i]]$individuals$N[,,i]
+      simulation@results$individuals$D[,,i] <- results[[i]]$individuals$D[,,i]
+      simulation@results$Detection[,,i] <- results[[i]]$Detection[,,i]
+      if(!is.null(simulation@results$clusters)){
+        simulation@results$clusters$summary[,,i] <- results[[i]]$clusters$summary[,,i]
+        simulation@results$clusters$N[,,i] <- results[[i]]$clusters$N[,,i]
+        simulation@results$clusters$D[,,i] <- results[[i]]$clusters$D[,,i]
+        simulation@results$expected.size[,,i] <- results[[i]]$expected.size[,,i]
+      }
     }
   }
   return(simulation)

diff --git a/R/calc.perp.dists.R b/R/calc.perp.dists.R
@@ -1,7 +1,6 @@
 #' @importFrom graphics points
-#' @importFrom sf st_intersection st_drop_geometry st_crs
+#' @importFrom sf st_drop_geometry st_crs st_join st_intersects st_distance st_geometry st_coordinates
 #' @importFrom methods is
-#' @importFrom purrr reduce
 calc.perp.dists <- function(population, transects, plot = FALSE){
   # Calculates the possible detection distances to the transects
   # Arguments:
@@ -10,82 +9,34 @@ calc.perp.dists <- function(population, transects, plot = FALSE){
   # Returns:
   #   A data frame of possible detection distances
 
-  subset.buff.func <- function(i, sf.pop, samplers, cov.areas){
-    #returns the locations of the population within the truncation distance of transect i.
-    # Extract relevant sampler
-    sf.column.t <- attr(samplers, "sf_column")
-    samp <- samplers[[sf.column.t]][[i]]
-    #Extract associated covered area
-    cov.area <- cov.areas[cov.areas$transect == i,]
-    # Find the population in the covered area of transect i
-    pop.in.cov <- suppressWarnings(
-      st_intersection(sf.pop, cov.area))
-    #Turn into a data.frame
-    sub.pop.coords <- as.data.frame(sf::st_coordinates(pop.in.cov))
-    names(sub.pop.coords) <- c("x","y")
-    # Add other info back in
-    sub.pop.coords <- cbind(sub.pop.coords, st_drop_geometry(pop.in.cov))
-    #Find start and end point [note may be a multilinestring]
-    if(is(samp, "LINESTRING")){
-      start.X <- samp[1]
-      start.Y <- samp[3]
-      end.X <- samp[2]
-      end.Y <- samp[4]
-    }else if(is(samp, "MULTILINESTRING")){
-      index.final <- length(samp)
-      start.X <- samp[[1]][1,1]
-      start.Y <- samp[[1]][1,2]
-      end.X <- samp[[index.final]][2,1]
-      end.Y <- samp[[index.final]][2,2]
-    }else{
-      stop("sampler is not of type linestring or multilinestring", call. = TRUE)
-    }
-    #now calculate distances to transect
-    #find the angle between the transect and the vector from the animal to the start of the transect
-    transect.angle <- atan2(end.Y-start.Y, end.X-start.X)
-    animal.angle <- atan2(sub.pop.coords$y-start.Y, sub.pop.coords$x-start.X)
-    delta.angle <- abs(animal.angle-transect.angle)
-    delta.angle <- (ifelse(delta.angle > pi, 2*pi - delta.angle, delta.angle))
-    #calculate the distance from the transect start to the animal (the hypotenuse)
-    hyp <- sqrt((sub.pop.coords$y-start.Y)^2+(sub.pop.coords$x-start.X)^2)
-    #calculate the perpendicular distance (the opposite side of the RA triangle)
-    perp.dists  <- hyp*sin(delta.angle)
-    #Add perp distances
-    if(nrow(sub.pop.coords) > 0){
-      #Make new dataset
-      new.data <- cbind(sub.pop.coords,
-                        Sample.Label = rep(samplers$transect[i], nrow(sub.pop.coords)),
-                        distance = perp.dists)
-    }else{
-      new.data <- NULL
-    }
-    return(new.data)
-  }
-  #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  #Main function begins
   samplers <- transects@samplers
   covered.areas <- transects@cov.area.polys
   pop <- population@population
-  sf.pop <- st_as_sf(pop, coords = c("x", "y")) 
+  sf.pop <- sf::st_as_sf(pop, coords = c("x", "y"))
   sf::st_crs(sf.pop) <- sf::st_crs(covered.areas)
-  #get all possible detection distances
-  all.poss.detects <- lapply(1:nrow(samplers),
-                             FUN = subset.buff.func,
-                             sf.pop = sf.pop,
-                             samplers = samplers,
-                             cov.areas = covered.areas)
-
-  #Build up into a single data.frame
-  new.dataframe <- reduce(all.poss.detects, rbind)
-  if(!is.null(new.dataframe)){
-    # Order the data by individual
-    index <- order(new.dataframe$individual)
-    new.dataframe <- new.dataframe[index,]
-  }else{
+  # One indexed spatial join is significantly faster than N per-transect intersections
+  covered.lookup <- covered.areas[, c("transect")]
+  joined <- suppressWarnings(
+    sf::st_join(sf.pop, covered.lookup, join = sf::st_intersects, left = FALSE)
+  )
+
+  if(nrow(joined) == 0){
     new.dataframe <- data.frame()
+    return(new.dataframe)
   }
-  # remove duplicate / redundant cols
-  #index <- which(names(tmp4) %in% c("transect", "strata"))
-  #ordered.data <- ordered.data[,-index]
+
+  transect.ids <- as.character(joined$transect)
+  sampler.ids <- as.character(samplers$transect)
+  sampler.index <- match(transect.ids, sampler.ids)
+  sampler.geom <- sf::st_geometry(samplers)[sampler.index]
+  distances <- as.numeric(sf::st_distance(sf::st_geometry(joined), sampler.geom, by_element = TRUE))
+
+  coords <- as.data.frame(sf::st_coordinates(joined))
+  names(coords) <- c("x", "y")
+  new.dataframe <- cbind(coords, sf::st_drop_geometry(joined))
+  new.dataframe$Sample.Label <- joined$transect
+  new.dataframe$distance <- distances
+  index <- order(new.dataframe$individual)
+  new.dataframe <- new.dataframe[index,]
   return(new.dataframe)
 }
-Original file line number
+Diff line change
@@ Expand Up / @@ -5,6 +5,7 @@ Imports: @@
         mrds,
         Distance,
         sf,
+        terra,
         ggplot2,
         purrr,
         dplyr,
@@ Expand Down @@