Feature: spatially explicit detectability covariates through a new cov.surface interface

DESCRIPTION

@@ -5,6 +5,7 @@ Imports:
     mrds,
     Distance,
     sf,
+    terra,
     ggplot2,
     purrr,
     dplyr,

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)
 }
 

R/Detectability.R

@@ -15,8 +15,15 @@
 #' 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 covariate
+#' effect parameters (log scale) for individual-level covariates. Entries may
+#' be numeric vectors for continuous/discrete covariates, or data.frames for
+#' categorical covariates with columns such as \code{strata}, \code{level},
+#' and \code{param}.
+#' @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-vector 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 +33,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 +63,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 +121,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 +197,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 +227,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 +255,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 +333,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 +374,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))

R/calculate.scale.param.R

@@ -1,4 +1,5 @@
 #' @importFrom methods is
+#' @importFrom terra rast cellFromXY values
 calculate.scale.param <- function(pop.data, detectability, region){
 # This function calculates the scale parameters including any covariate effects
 # and adds these values to the population dataframe which is then returned. Also
@@ -32,6 +33,29 @@ calculate.scale.param <- function(pop.data, detectability, region){
       detectability@shape.param <- rep(shape[1], strata.no)
     }
   }
+  # Extract raster surface values at each animal's location and append as columns.
+  # This must happen before the covariate-name check so the columns are present.
+  if(length(detectability@cov.surface) > 0){
+    surf.names <- names(detectability@cov.surface)
+    pts <- cbind(pop.data$x, pop.data$y)
+    for(sn in surf.names){
+      surf <- detectability@cov.surface[[sn]]
+      if(is.character(surf)){
+        surf <- terra::rast(surf)
+      }
+      cell_ids <- terra::cellFromXY(surf, pts)
+      rast_vals <- terra::values(surf)[, 1]
+      vals <- rast_vals[cell_ids]
+      if(any(is.na(vals))){
+        warning(paste0("NA raster values for surface covariate '", sn,
+                       "' at some animal locations (outside raster extent). ",
+                       "Replacing with global raster mean."),
+                call. = FALSE, immediate. = TRUE)
+        vals[is.na(vals)] <- mean(vals, na.rm = TRUE)
+      }
+      pop.data[[sn]] <- vals
+    }
+  }
   # Check the covariate names in detectability are present in the population
   pop.cov.names <- names(pop.data)
   detect.cov.names <- names(detectability@cov.param)
@@ -44,7 +68,7 @@ calculate.scale.param <- function(pop.data, detectability, region){
   for(cov in seq(along = detectability@cov.param)){
     current.cov <- detectability@cov.param[[cov]]
     if(!is(current.cov, "data.frame")){
-      if(length(current.cov == 1)){
+      if(length(current.cov) == 1){
         # repeat it for the number of strata
         detectability@cov.param[[cov]] <- rep(current.cov, strata.no)
       }else if(length(current.cov) != strata.no){

R/check.simulation.R

@@ -28,7 +28,15 @@ check.simulation <- function(object){
       }
     }
   }
-
+  # Check that cov.surface names all have corresponding cov.param entries
+  surf.names <- names(detect@cov.surface)
+  if(length(surf.names) > 0){
+    for(sn in surf.names){
+      if(!sn %in% cov.names){
+        return(paste0("cov.surface entry '", sn, "' does not have a matching numeric entry in cov.param."))
+      }
+    }
+  }
   # DESIGN CHECKS
   design <- object@design
   if(any(design@edge.protocol == "plus")){

tests/testthat/test-check_IndividualCovariates.R

@@ -269,7 +269,7 @@ test_that("Factor level covariates", {
 
   pop.data <- data.frame(individual, x, y, Region.Label, size, sex)
 
-  check.scale.params <- calculate.scale.param(pop.data, detect, region)
+  check.scale.params <- dsims:::calculate.scale.param(pop.data, detect, region)
 
   ind1.scale <- exp(log(0.08) + log(1.02)*size[1] + log(1.5))
   ind2.scale <- exp(log(0.08) + log(1.02)*size[2] + log(1.5))