From f496c0a1d04b804060a2f44a8f6fdfcf57818f13 Mon Sep 17 00:00:00 2001 From: hambrecht <137965368+hambrecht@users.noreply.github.com> Date: Tue, 5 May 2026 15:06:29 -0700 Subject: [PATCH 1/3] feat: add spatially-explicit raster covariate support to detectability - add cov.surface support to Detectability and make.detectability - apply raster covariates in calculate.scale.param - validate cov.surface in simulation checks - add raster covariate test coverage - update related simulation/covariate paths used by tests --- DESCRIPTION | 1 + R/ClassConstructors.R | 20 +- R/Detectability.R | 41 +++- R/calculate.scale.param.R | 26 +- R/check.simulation.R | 10 +- .../test-check_IndividualCovariates.R | 2 +- tests/testthat/test-check_RasterCovariate.R | 227 ++++++++++++++++++ 7 files changed, 319 insertions(+), 8 deletions(-) create mode 100644 tests/testthat/test-check_RasterCovariate.R diff --git a/DESCRIPTION b/DESCRIPTION index 67c9158..8d9dc72 100644 --- 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 index c669fbc..e99d8c4 100755 --- 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 index 3e0fe08..6f206b9 100755 --- 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) } ) diff --git a/R/calculate.scale.param.R b/R/calculate.scale.param.R index d3fea0e..42d79b4 100755 --- a/R/calculate.scale.param.R +++ b/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){ diff --git a/R/check.simulation.R b/R/check.simulation.R index 855b9ac..d8a8f07 100644 --- a/R/check.simulation.R +++ b/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")){ diff --git a/tests/testthat/test-check_IndividualCovariates.R b/tests/testthat/test-check_IndividualCovariates.R index 45c329e..f670024 100644 --- a/tests/testthat/test-check_IndividualCovariates.R +++ b/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)) diff --git a/tests/testthat/test-check_RasterCovariate.R b/tests/testthat/test-check_RasterCovariate.R new file mode 100644 index 0000000..3559d8a --- /dev/null +++ b/tests/testthat/test-check_RasterCovariate.R @@ -0,0 +1,227 @@ +context("Raster-based spatial covariate detectability (cov.surface)") + +# --------------------------------------------------------------------------- +# Helper: build a minimal single-stratum region covering [0,1] x [0,1] +# --------------------------------------------------------------------------- +make_test_region <- function(){ + pol <- sf::st_polygon(list(matrix(c(0,0, 0,1, 1,1, 1,0, 0,0), ncol = 2, byrow = TRUE))) + sfc <- sf::st_sfc(pol) + sf.pol <- sf::st_sf(strata = "main", geom = sfc) + make.region(region.name = "test", strata.name = "main", shape = sf.pol) +} + +# --------------------------------------------------------------------------- +# 1. Canopy raster reduces scale.param at higher cover values +# --------------------------------------------------------------------------- +test_that("Canopy raster covariate: higher cover yields smaller scale.param", { + skip_if_not_installed("terra") + + region <- make_test_region() + + # Two-zone raster: left half (x < 0.5) canopy = 0.2, right half = 0.8 + r <- terra::rast(nrows = 10, ncols = 10, xmin = 0, xmax = 1, ymin = 0, ymax = 1) + m <- matrix(c(rep(0.2, 50), rep(0.8, 50)), nrow = 10, ncol = 10, byrow = FALSE) + terra::values(r) <- as.vector(t(m)) + + rast_path <- tempfile(fileext = ".tif") + on.exit(unlink(rast_path), add = TRUE) + terra::writeRaster(r, rast_path, overwrite = TRUE) + + # Slope = -1 on log scale: higher canopy -> smaller sigma -> lower detectability + detect <- make.detectability( + key.function = "hn", + scale.param = 25, + cov.param = list(canopy = -1), + cov.surface = list(canopy = rast_path), + truncation = 50 + ) + + pop.data <- data.frame( + individual = 1:4, + x = c(0.1, 0.3, 0.6, 0.8), + y = c(0.5, 0.5, 0.5, 0.5), + Region.Label = "main" + ) + + result <- dsims:::calculate.scale.param(pop.data, detect, region) + + # Expected: sigma = exp(log(25) + (-1) * canopy_value) + low_scale <- exp(log(25) + (-1) * 0.2) # left-zone animals + high_scale <- exp(log(25) + (-1) * 0.8) # right-zone animals + + # Animals 1 & 2 are in low-canopy zone; 3 & 4 are in high-canopy zone + expect_lt(result$scale.param[3], result$scale.param[1], + label = "high-canopy scale.param < low-canopy scale.param") + expect_equal(result$scale.param[1], low_scale, tolerance = 0.05) + expect_equal(result$scale.param[2], low_scale, tolerance = 0.05) + expect_equal(result$scale.param[3], high_scale, tolerance = 0.05) + expect_equal(result$scale.param[4], high_scale, tolerance = 0.05) +}) + +# --------------------------------------------------------------------------- +# 2. File path and in-memory SpatRaster give identical results +# --------------------------------------------------------------------------- +test_that("File path and in-memory SpatRaster produce the same scale.param", { + skip_if_not_installed("terra") + + region <- make_test_region() + + r <- terra::rast(nrows = 5, ncols = 5, xmin = 0, xmax = 1, ymin = 0, ymax = 1) + terra::values(r) <- seq(0, 1, length.out = 25) + + rast_path <- tempfile(fileext = ".tif") + on.exit(unlink(rast_path), add = TRUE) + terra::writeRaster(r, rast_path, overwrite = TRUE) + + pop.data <- data.frame( + individual = 1:3, + x = c(0.1, 0.5, 0.9), + y = c(0.5, 0.5, 0.5), + Region.Label = "main" + ) + + detect_path <- make.detectability(cov.param = list(canopy = -0.5), + cov.surface = list(canopy = rast_path)) + detect_mem <- make.detectability(cov.param = list(canopy = -0.5), + cov.surface = list(canopy = r)) + + res_path <- dsims:::calculate.scale.param(pop.data, detect_path, region) + res_mem <- dsims:::calculate.scale.param(pop.data, detect_mem, region) + + expect_equal(res_path$scale.param, res_mem$scale.param, tolerance = 1e-6) +}) + +# --------------------------------------------------------------------------- +# 3. Validation: cov.surface name not in cov.param triggers an error +# --------------------------------------------------------------------------- +test_that("make.detectability errors when cov.surface name absent from cov.param", { + skip_if_not_installed("terra") + + rast_path <- tempfile(fileext = ".tif") + r <- terra::rast(nrows = 2, ncols = 2, xmin = 0, xmax = 1, ymin = 0, ymax = 1) + terra::values(r) <- 1:4 + on.exit(unlink(rast_path), add = TRUE) + terra::writeRaster(r, rast_path, overwrite = TRUE) + + expect_error( + make.detectability( + cov.param = list(size = log(1.1)), # 'size', not 'canopy' + cov.surface = list(canopy = rast_path) # 'canopy' has no cov.param slope + ), + regexp = "cov.surface entry 'canopy' has no matching" + ) +}) + +# --------------------------------------------------------------------------- +# 4. Validation: non-character / non-SpatRaster entry triggers an error +# --------------------------------------------------------------------------- +test_that("make.detectability errors when cov.surface element is wrong type", { + expect_error( + make.detectability( + cov.param = list(canopy = -0.5), + cov.surface = list(canopy = 42) # numeric scalar is invalid + ), + regexp = "must be a file path" + ) +}) + +# --------------------------------------------------------------------------- +# 5. Validation: non-existent file path triggers an error +# --------------------------------------------------------------------------- +test_that("make.detectability errors when cov.surface file path does not exist", { + expect_error( + make.detectability( + cov.param = list(canopy = -0.5), + cov.surface = list(canopy = "/does/not/exist/canopy.tif") + ), + regexp = "does not exist" + ) +}) + +# --------------------------------------------------------------------------- +# 6. Existing behaviour unchanged when cov.surface is empty (regression guard) +# --------------------------------------------------------------------------- +test_that("calculate.scale.param unchanged when no cov.surface is supplied", { + region <- make_test_region() + + detect <- make.detectability( + key.function = "hn", + scale.param = 25, + cov.param = list(size = log(1.02)), + truncation = 50 + ) + + pop.data <- data.frame( + individual = 1:3, + x = c(0.2, 0.5, 0.8), + y = c(0.5, 0.5, 0.5), + Region.Label = "main", + size = c(10L, 20L, 30L) + ) + + result <- dsims:::calculate.scale.param(pop.data, detect, region) + + expected <- exp(log(25) + log(1.02) * c(10, 20, 30)) + expect_equal(result$scale.param, expected, tolerance = 1e-8) +}) + +# --------------------------------------------------------------------------- +# 7. Integration: run.simulation uses raster covariate each replicate +# --------------------------------------------------------------------------- +test_that("run.simulation completes with cov.surface raster (2 reps)", { + skip_if_not_installed("terra") + skip_if_not_installed("Distance") + + region <- make_test_region() + + # Flat-gradient raster (constant value = 0.5, so slope effect is fixed) + r <- terra::rast(nrows = 10, ncols = 10, xmin = 0, xmax = 1, ymin = 0, ymax = 1) + terra::values(r) <- 0.5 + rast_path <- tempfile(fileext = ".tif") + on.exit(unlink(rast_path), add = TRUE) + terra::writeRaster(r, rast_path, overwrite = TRUE) + + density <- make.density(region = region, x.space = 0.1, constant = 200) + + popdesc <- make.population.description( + region = region, + density = density, + fixed.N = TRUE, + N = 200 + ) + + detect <- make.detectability( + key.function = "hn", + scale.param = 0.05, + cov.param = list(canopy = -0.5), + cov.surface = list(canopy = rast_path), + truncation = 0.1 + ) + + design <- make.design( + region = region, + transect.type = "line", + design = "systematic", + samplers = 5, + truncation = 0.1 + ) + + analysis <- make.ds.analysis( + dfmodel = list(~1), + key = "hn", + truncation = 0.1, + criteria = "AIC" + ) + + sim <- make.simulation( + reps = 2, + design = design, + population.description = popdesc, + detectability = detect, + ds.analysis = analysis + ) + + # Should run without error and return a Simulation object with results + result <- run.simulation(sim, run.parallel = FALSE, counter = FALSE) + expect_s4_class(result, "Simulation") +}) From 45d9abc78bfd2b14827ab0248d5ab64690cf8be3 Mon Sep 17 00:00:00 2001 From: hambrecht <137965368+hambrecht@users.noreply.github.com> Date: Tue, 5 May 2026 15:04:49 -0700 Subject: [PATCH 2/3] fix: allow plotting detectability with raster covariates Update Detectability plot validation so covariates supplied via cov.surface are treated as valid even when absent from population covariates. Add raster-specific plotting path that reads surface values and uses raster quantiles to generate representative detection curves, enabling mixed covariate plotting (for example size + canopy) without false missing-covariate errors. --- R/Detectability.R | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/R/Detectability.R b/R/Detectability.R index 6f206b9..d40c69b 100755 --- a/R/Detectability.R +++ b/R/Detectability.R @@ -194,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 @@ -222,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")){ @@ -236,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 = "") } @@ -312,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){ @@ -350,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)) From 971b41aa86cefcd58fc427d255f1cf91fc3ceb3e Mon Sep 17 00:00:00 2001 From: Leonard Hambrecht <137965368+hambrecht@users.noreply.github.com> Date: Tue, 5 May 2026 16:01:35 -0700 Subject: [PATCH 3/3] Potential fix for pull request finding Updated the @slot cov.param docs to reflect the actual allowed shapes (numeric vectors for continuous/discrete + data.frames for categorical levels). Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- R/Detectability.R | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/R/Detectability.R b/R/Detectability.R index d40c69b..bc6aecd 100755 --- a/R/Detectability.R +++ b/R/Detectability.R @@ -15,12 +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{"list"}; Named list of slope -#' coefficients (log scale) for individual-level covariates. +#' @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 entries in \code{cov.param}. +#' 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