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..bc6aecd 100755 --- a/R/Detectability.R +++ b/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)) 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") +})