From a640e40ab3a1b28f8c0cd89beb05504acec528e9 Mon Sep 17 00:00:00 2001 From: Martin Jung Date: Thu, 28 Mar 2024 22:29:53 +0100 Subject: [PATCH] Further fixes to derivates and scenario variable preparations #106 --- NEWS.md | 2 +- R/add_constraint.R | 27 ++++ R/add_predictors.R | 45 ++++-- R/class-distributionmodel.R | 22 ++- R/class-predictors.R | 21 ++- R/ibis.iSDM-package.R | 2 +- R/utils-predictors.R | 130 ++++++++++++------ R/utils-scenario.R | 20 ++- man/PredictorDataset-class.Rd | 22 ++- man/add_constraint_dispersal.Rd | 19 +++ man/add_predictors.Rd | 6 +- tests/testthat/test_Scenarios.R | 27 +++- .../articles/01_data_preparationhelpers.Rmd | 26 ++++ vignettes/articles/03_integrate_data.Rmd | 3 +- 14 files changed, 290 insertions(+), 82 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6ea3fe84..fcfc88b5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ * Small fix to parameter in `train()` #102 @jeffreyhanson * Small helper function for combining 2 different formula objects `combine_formulas()` * Small bug fixes dealing with `scenario()` projections and limits, plus unit tests #104 -* Fixed bug with adding `predictor_derivate()` to scenario predictors and added unit tests #106 +* Bug fixes with adding `predictor_derivate()` to scenario predictors and added unit tests #106 # ibis.iSDM 0.1.2 diff --git a/R/add_constraint.R b/R/add_constraint.R index 9fb47b41..382ef45c 100644 --- a/R/add_constraint.R +++ b/R/add_constraint.R @@ -149,6 +149,16 @@ methods::setMethod( #' instance \code{`kissmig`}. #' #' @details +#' **Dispersal**: +#' Parameters for \code{'method'}: +#' * \code{sdd_fixed} - Applies a fixed uniform dispersal distance per modelling timestep. +#' * \code{sdd_nexpkernel} - Applies a dispersal distance using a negative exponential kernel from its origin. +#' * \code{kissmig} - Applies the kissmig stochastic dispersal model. Requires \code{`kissmig`} package. Applied at each modelling time step. +#' * \code{migclim} - Applies the dispersal algorithm MigClim to the modelled objects. Requires \code{"MigClim"} package. +#' +#' A comprehensive overview of the benefits of including dispersal constrains in +#' species distribution models can be found in Bateman et al. (2013). +#' #' The following additional parameters can bet set: #' * \code{pext}: [`numeric`] indicator for \code{`kissmig`} of the probability a #' colonized cell becomes uncolonised, i.e., the species gets locally extinct @@ -156,6 +166,10 @@ methods::setMethod( #' * \code{pcor}: [`numeric`] probability that corner cells are considered in the #' 3x3 neighbourhood (Default: \code{0.2}). #' +#' @references +#' * Bateman, B. L., Murphy, H. T., Reside, A. E., Mokany, K., & VanDerWal, J. (2013). +#' Appropriateness of full‐, partial‐and no‐dispersal scenarios in climate change impact +#' modelling. Diversity and Distributions, 19(10), 1224-1234. #' @family constraint #' @keywords scenario #' @@ -424,6 +438,19 @@ methods::setMethod( #' @param resistance A [`SpatRaster`] object describing a resistance surface or #' barrier for use in connectivity constrains (Default: \code{NULL}). #' +#' @details +#' * \code{hardbarrier} - Defines a hard barrier to any dispersal events. By +#' definition this sets all values larger +#' than \code{0} in the barrier layer to \code{0} in the projection. Barrier has +#' to be provided through the \code{"resistance"} parameter. +#' * \code{resistance} - Allows the provision of a static or dynamic layer that is +#' multiplied with the projection at each time step. Can for example be used to +#' reduce the suitability of any given area (using pressures not included in the model). +#' The respective layer(s) have to be provided through the \code{"resistance"} parameter. +#' Provided layers are incorporated as \code{abs(resistance - 1)} and multiplied with +#' the prediction. +#' +#' #' @family constraint #' @keywords scenario #' diff --git a/R/add_predictors.R b/R/add_predictors.R index 0d0d9823..57a87859 100644 --- a/R/add_predictors.R +++ b/R/add_predictors.R @@ -51,10 +51,10 @@ NULL #' #' Available options for creating derivates are: #' * \code{'none'} - No additional predictor derivates are created. -#' * \code{'quad'} - Adds quadratic transformed predictors. +#' * \code{'quad'} - Adds quadratic derivate predictors. #' * \code{'interaction'} - Add interacting predictors. Interactions need to be specified (\code{"int_variables"})! -#' * \code{'thresh'} - Add threshold transformed predictors. -#' * \code{'hinge'} - Add hinge transformed predictors. +#' * \code{'thresh'} - Add threshold derivate predictors. +#' * \code{'hinge'} - Add hinge derivate predictors. #' * \code{'bin'} - Add predictors binned by their percentiles. #' #' @note @@ -161,6 +161,7 @@ methods::setMethod( assertthat::assert_that( all( priors$varnames() %in% names(env) ) ) y <- y$set_priors(priors) } + # Harmonize NA values if(harmonize_na){ if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','green','Harmonizing missing values...') @@ -248,6 +249,7 @@ methods::setMethod( # Finally set the data to the BiodiversityDistribution object pd <- PredictorDataset$new(id = new_id(), data = env, + transformed = ifelse('none' %notin% transform, TRUE, FALSE ), ...) y$set_predictors(pd) } @@ -727,7 +729,7 @@ methods::setMethod( derivate_knots = 4, int_variables = NULL, harmonize_na = FALSE, ... ) { # Try and match transform and derivatives arguments transform <- match.arg(transform, c('none','pca', 'scale', 'norm', 'windsor', 'percentile'), several.ok = TRUE) - derivates <- match.arg(derivates, c('none','thresh', 'hinge', 'quadratic', 'bin'), several.ok = TRUE) + derivates <- match.arg(derivates, c('none','thresh', 'hinge', 'quadratic', 'bin', 'interaction'), several.ok = TRUE) assertthat::validate_that(inherits(env,'stars'), msg = 'Projection rasters need to be stars stack!') assertthat::assert_that(inherits(x, "BiodiversityScenario"), @@ -762,26 +764,38 @@ methods::setMethod( if('none' %notin% transform){ if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','green','Transforming predictors...') for(tt in transform) env <- predictor_transform(env, option = tt) + } else { + # Check regardless + try({ test <- obj$model$predictors_object$is_transformed() },silent = TRUE) + if(!inherits(test, 'try-error')){ + if(test) if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','red','Option to transform predictors set to None, yet transformed variables were used in trained model?') + } } - # # Calculate derivates if set + # Calculate derivates if set if('none' %notin% derivates){ # Get variable names - varn <- obj$get_coefficients()[['Feature']] + varn <- obj$get_coefficients()[,1] # Are there any derivates present in the coefficients? - if(any( length( grep("hinge__|bin__|quad__|thresh__", varn ) ) > 0 )){ - if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','green','Creating predictor derivates...') - for(dd in derivates){ - if(any(grep(dd, varn))){ - env <- predictor_derivate(env, option = dd, nknots = derivate_knots, - deriv = varn, int_variables = int_variables) - } else { - if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','red', paste0(derivates,' derivates should be created, but not found among coefficients!')) - } + if(any( length( grep("hinge_|bin_|quadratic_|thresh_|interaction_", varn ) ) > 0 )){ + if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','green','Creating predictor derivates...') + for(dd in derivates){ + if(any(grep(dd, varn))){ + env <- predictor_derivate(env, option = dd, nknots = derivate_knots, + deriv = grep(dd,varn, value=TRUE), int_variables = int_variables) + } else { + if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','red', paste0(derivates,' derivates should be created, but not found among model coefficients!')) } + } } else { if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','red','No derivates found among coefficients. None created for projection!') } + } else { + # Check regardless as security check + try({ test <- obj$model$predictors_object$has_derivates() },silent = TRUE) + if(!inherits(test, 'try-error')){ + if(test) if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','red','Option to create derivates set to None, but likely derivates found among coefficients?') + } } # Get, guess and format Time period @@ -800,6 +814,7 @@ methods::setMethod( # Finally set the data to the BiodiversityScenario object pd <- PredictorDataset$new(id = new_id(), data = env, + transformed = ifelse('none' %notin% transform, TRUE, FALSE ), timeperiod = timeperiod ) # Make a clone copy of the object diff --git a/R/class-distributionmodel.R b/R/class-distributionmodel.R index b35ac869..2eb306cd 100644 --- a/R/class-distributionmodel.R +++ b/R/class-distributionmodel.R @@ -450,8 +450,16 @@ DistributionModel <- R6::R6Class( #' @param x A [`character`] stating what should be returned. #' @return A [`SpatRaster`] object with the prediction. get_data = function(x = "prediction") { - if (!x %in% names(self$fits)) + rr <- names(self$fits) + if(!x %in% names(self$fits)){ + # Check if x is present in rr, if so print a message + if(length(grep(x,rr))>0){ + if(getOption('ibis.setupmessages', default = TRUE)){ + myLog('[Estimation]','yellow','Output not found, but found: ', grep(x,rr,value = TRUE)[1]) + } + } return(new_waiver()) + } return(self$fits[[x]]) }, @@ -678,7 +686,7 @@ DistributionModel <- R6::R6Class( assertthat::assert_that( is.character(fname), type %in% c('gtif','gtiff','tif','nc','ncdf'), - 'fits' %in% self$ls(), + 'fits' %in% names(self), dt %in% c('LOG1S','INT1S','INT1U','INT2S','INT2U','INT4S','INT4U','FLT4S','FLT8S') ) type <- tolower(type) @@ -688,10 +696,18 @@ DistributionModel <- R6::R6Class( # Get raster file in fitted object cl <- sapply(self$fits, class) - ras <- self$fits[[grep('SpatRaster', cl,ignore.case = T)]] + if(length( grep('SpatRaster', cl,ignore.case = T) )==0){ + # Security check in case for some reason there are no predictions + if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Output]','red','No predictions found?') + return(NULL) + } + ras <- self$fits[grep('SpatRaster', cl,ignore.case = T)] assertthat::assert_that(length(ras)>0, msg = "No prediction to save found.") + # If is a list (multiple SpatRaster) -> Combine + if(is.list(ras)) ras <- Reduce('c', ras) + # Check that no-data value is not present in ras assertthat::assert_that(any(!terra::global(ras, "min", na.rm = TRUE)[,1] <= -9999), msg = 'No data value -9999 is potentially in prediction!') diff --git a/R/class-predictors.R b/R/class-predictors.R index a6258a6c..246443e7 100644 --- a/R/class-predictors.R +++ b/R/class-predictors.R @@ -21,25 +21,30 @@ PredictorDataset <- R6::R6Class( #' @field id The id for this collection as [`character`]. #' @field data A predictor dataset usually as [`SpatRaster`]. #' @field name A name for this object. + #' @field transformed Saves whether the predictors have been transformed somehow. #' @field timeperiod A timeperiod field id = character(0), data = new_waiver(), name = character(0), + transformed = logical(0), timeperiod = new_waiver(), #' @description #' Initializes the object and creates an empty list #' @param id The id for this collection as [`character`]. #' @param data A predictor dataset usually as [`SpatRaster`]. + #' @param transformed A [`logical`] flag if predictors have been transformed. Assume not. #' @param ... Any other parameters found. #' @return NULL - initialize = function(id, data, ...){ + initialize = function(id, data, transformed = FALSE, ...){ assertthat::assert_that( - is.Id(id) || is.character(id) + is.Id(id) || is.character(id), + is.logical(transformed) ) self$name <- 'Biodiversity data' self$id <- id self$data <- data + self$transformed <- transformed # Get Dots and save too dots <- list(...) for(el in names(dots)){ @@ -296,7 +301,17 @@ PredictorDataset <- R6::R6Class( #' @return A [`logical`] flag. has_derivates = function(){ return( - base::length( grep("hinge__|bin__|quad__|thresh__", self$get_names() ) ) > 0 + base::length( grep("hinge_|bin_|quadratic_|thresh_|interaction_", self$get_names() ) ) > 0 + ) + }, + + #' @description + #' Predictors have been transformed? + #' @seealso [predictor_transform()] + #' @return A [`logical`] flag. + is_transformed = function(){ + return( + self$transformed ) }, diff --git a/R/ibis.iSDM-package.R b/R/ibis.iSDM-package.R index f3744ac9..24bdf49b 100644 --- a/R/ibis.iSDM-package.R +++ b/R/ibis.iSDM-package.R @@ -30,7 +30,7 @@ globalVariables(c("background", "band", "bi_class", "bias", #MJ: Added self here hoping that does not crash all methods. "self", "id", "included", "i", - "km", "vt", + "km", "vt", "V2", "limit", "lower", "layer", "median", "model", "mad", "name", diff --git a/R/utils-predictors.R b/R/utils-predictors.R index 8c3fae58..95077532 100644 --- a/R/utils-predictors.R +++ b/R/utils-predictors.R @@ -196,10 +196,7 @@ predictor_transform <- function(env, option, windsor_props = c(.05,.95), pca.var # Sample covariance from stack and fit PCA covMat <- terra::layerCor(env, fun = "cov", na.rm = TRUE) - pca <- stats::princomp(covmat = covMat[[1]], cor = FALSE) - # Add means and grid cells - pca$center <- covMat$mean - pca$n.obs <- terra::ncell(env) + pca <- terra::princomp(env) # Use Terra directly # Check how many components are requested: if(pca.var<1){ @@ -208,7 +205,7 @@ predictor_transform <- function(env, option, windsor_props = c(.05,.95), pca.var nComp <- length( which(props <= pca.var) ) } # Predict principle components - out <- terra::predict(env, pca,na.rm = TRUE, index = 1:nComp) + out <- terra::predict(env, pca, index = 1:nComp) names(out) <- paste0("PC", 1:nComp) return(out) @@ -328,19 +325,19 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab # If deriv is not set, assume all variables should be derivated and set if(is.null(deriv)){ - deriv <- paste0(option, "__", names(env)) + deriv <- paste0(option, "_", names(env)) } # If stars see if we can convert it to a stack if(inherits(env, 'stars')){ assertthat::assert_that(!is.null(deriv),msg = "Derivate names could not be found!") # Decompose derivate variable names if set - deriv <- grep(paste0(option, "__"), deriv, value = TRUE) + deriv <- grep(paste0(option, "_"), deriv, value = TRUE) if(length(deriv)==0){ if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','red','Predictors with derivates not found!') return(NULL) } - cutoffs <- do.call(rbind,strsplit(deriv, "__")) |> as.data.frame() + cutoffs <- do.call(rbind,strsplit(deriv, "_")) |> as.data.frame() cutoffs$deriv <- deriv lyrs <- names(env) # Names of predictors @@ -367,12 +364,12 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab } else { new_env <- terra::app(env, function(x) I(x^2)) } - names(new_env) <- paste0('quad__', names(env)) + names(new_env) <- paste0('quadratic_', names(env)) } else { # Convert to quadratic transform (without terra::app as this uses the time dimension) new_env <- lapply(env_list, function(z){ o <- (z^2) - names(o) <- paste0('quad__', names(o)) + names(o) <- paste0('quadratic_', names(o)) return(o) } ) } @@ -423,7 +420,7 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab o[] <- hingeval(env_list[[val]][[k]][], first, last) terra::time(o) <- terra::time( env_list[[val]][[k]] ) - names(o) <- paste0(val, "__",first,"_",last) + names(o) <- paste0(val, "_",first,"_",last) suppressWarnings( nn <- append(nn, o) ) rm(o) } @@ -432,7 +429,7 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab new_env[[val]] <- nn rm(nn) } - } + } invisible(gc()) } @@ -457,14 +454,30 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab if(any(substr(cu,1, 1)==".")){ cu[which(substr(cu,1, 1)==".")] <- gsub("^.","",cu[which(substr(cu,1, 1)==".")]) } - cu <- as.numeric(cu) - assertthat::assert_that(!anyNA(cu), is.numeric(cu)) + suppressWarnings(cu <- as.numeric(cu) ) + if(is.na(cu[1]) || is.na(cu[2])){ + # Use the default knots for cutting and get knots + o <- makeThresh(env_list[[val]][[1]], n = val, nknots = nknots, cutoffs = NULL) + suppressWarnings( cu <- strsplit(names(o), "_") |> + unlist() |> + as.numeric()) + cu <- subset(cu, stats::complete.cases(cu)) + cu <- matrix(cu, byrow = FALSE) # Convert to pmin and pmax + } + nn <- terra::rast() for(k in 1:terra::nlyr(env_list[[val]])){ - o <- emptyraster(env_list[[val]][[k]]) - o[] <- thresholdval(env_list[[val]][[k]][], cu) - env_list[[val]][[k]] <- o - rm(o) + # For each cutoff + for(i in seq(1, length(cu))){ + o <- emptyraster(env_list[[val]][[k]]) + o[] <- thresholdval(env_list[[val]][[k]][], cu[i]) + terra::time(o) <- terra::time( env_list[[val]][[k]] ) + names(o) <- paste0(val, "_",cu[i]) + suppressWarnings( nn <- append(nn, o) ) + rm(o) + } } + new_env[[val]] <- nn + rm(nn) } invisible(gc()) } @@ -490,16 +503,32 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab if(any(substr(cu,1, 1)==".")){ cu[which(substr(cu,1, 1)==".")] <- gsub("^.","",cu[which(substr(cu,1, 1)==".")]) } - cu <- as.numeric(cu) - assertthat::assert_that(!anyNA(cu), is.numeric(cu)) + suppressWarnings( cu <- as.numeric(cu) ) + nn <- terra::rast() + # If NA, set + if(is.na(cu[1]) || is.na(cu[2])){ + # Use the default knots for cutting and get knots + o <- makeBin(env_list[[val]][[1]], n = val, nknots = nknots, cutoffs = NULL) + suppressWarnings( cu <- strsplit(names(o), "_") |> + unlist() |> + as.numeric()) + cu <- subset(cu, stats::complete.cases(cu)) + cu <- matrix(cu,ncol=2,byrow = TRUE) # Convert to pmin and pmax + cu <- cbind(cu, 1:nrow(cu)) + } for(k in 1:terra::nlyr(env_list[[val]])){ - o <- emptyraster(env_list[[val]][[k]]) - o <- terra::classify(env_list[[val]][[k]], cu) - o[is.na(o)] <- 0 - o <- terra::mask(o, env_list[[val]][[k]] ) - env_list[[val]][[k]] <- o - rm(o) + newcu <- cu + # Set smallest and largest value to global minimum/maximum to account for rounding issues + newcu[1] <- terra::global(env_list[[val]][[k]], "min",na.rm=T)[,1] + newcu[nrow(newcu)*2] <- terra::global(env_list[[val]][[k]], "max",na.rm=T)[,1] + + o <- terra::classify(env_list[[val]][[k]], newcu, include.lowest=TRUE) + terra::time(o) <- terra::time( env_list[[val]][[k]] ) + names(o) <- paste0(val, "_",nrow(newcu)) + suppressWarnings( nn <- append(nn, o) ) + rm(o, newcu) } + new_env[[val]] <- nn } invisible(gc()) } @@ -511,34 +540,57 @@ predictor_derivate <- function(env, option, nknots = 4, deriv = NULL, int_variab if(is.null(int_variables)){ int_variables <- attr(env, "int_variables") } - assertthat::assert_that(is.vector(int_variables)) + # If interaction check that vector is not null and exactly 2 + assertthat::assert_that(!is.null(int_variables), + is.vector(int_variables), + msg = "Provide a vector of 2 variables for an interaction!") - if(is.Raster(env)){ - # Make unique combinations - ind <- utils::combn(int_variables, 2) + # Make unique combinations + ind <- utils::combn(int_variables, 2) + if(is.Raster(env)){ # Now for each combination build new variable new_env <- terra::rast() for(i in 1:ncol(ind)){ # Multiply first with second entry o <- env[[ind[1,i]]] * env[[ind[2,i]]] - names(o) <- paste0('inter__', ind[1, i],".", ind[2, i]) + names(o) <- paste0('interaction_', ind[1, i],".", ind[2, i]) suppressWarnings( new_env <- c(new_env, o) ) rm(o) } } else { # Stars processing - stop("Not yet implemented!") + cu <- subset(cutoffs, V2 %in% ind) # Get only relevant variables + + for(i in 1:ncol(ind)){ + # For each layer + first <- cu$deriv[cu[,2]==ind[1,i]] + last <- cu$deriv[cu[,2]==ind[2,i]] + # Multiply first with second entry + o <- env_list[[first]] * env_list[[last]] + val <- paste0('interaction_', ind[1, i],".", ind[2, i]) # Interaction term + names(o) <- rep(paste0(val),terra::nlyr(o)) + new_env[[val]] <- o + rm(o) + } } } # If stars convert results back to stars object if(inherits(env, 'stars')){ # Add the original layers back - if(option == "quadratic"){ - o <- stars::st_as_stars(new_env[[cutoffs$deriv]], crs = sf::st_crs(env)) - names(o) <- cutoffs$deriv + if(option %in% c("quadratic")){ + o <- lapply(new_env, function(z) { + stars::st_as_stars(z, crs = sf::st_crs(env)) + }) + o <- Reduce("c", o) + } else if(option == "interaction"){ + o <- list() + for(val in names(new_env)){ + o[[val]] <- stars::st_as_stars(new_env[[val]]) + } + o <- Reduce("c", o) } else { o <- list() for(val in unique(names(new_env[[cutoffs$deriv]])) ){ @@ -690,8 +742,8 @@ makeHinge <- function(v, n, nknots = 4, cutoffs = NULL){ lh <- outer(v[] |> as.vector(), utils::head(k, -1), function(w, h) hingeval(w,h, v.max)) # Hinge starting from min rh <- outer(v[] |> as.vector(), k[-1], function(w, h) hingeval(w, v.min, h)) - colnames(lh) <- paste0("hinge__",n,'__', round( utils::head(k, -1), 2),'_', round(v.max, 2)) - colnames(rh) <- paste0("hinge__",n,'__', round( v.min, 2),'_', round(k[-1], 2)) + colnames(lh) <- paste0("hinge_",n,'_', round( utils::head(k, -1), 2),'_', round(v.max, 2)) + colnames(rh) <- paste0("hinge_",n,'_', round( v.min, 2),'_', round(k[-1], 2)) o <- as.data.frame( cbind(lh, rh) ) @@ -738,7 +790,7 @@ makeThresh <- function(v, n, nknots = 4, cutoffs = NULL){ } if(length(k)<=1) return(NULL) f <- outer(v[] |> as.vector(), k, function(w, t) ifelse(w >= t, 1, 0)) - colnames(f) <- paste0("thresh__", n, "__", round(k, 2)) + colnames(f) <- paste0("thresh_", n, "_", round(k, 2)) f <- as.data.frame(f) attr(f, "deriv.thresh") <- k return(f) @@ -793,7 +845,7 @@ makeBin <- function(v, n, nknots, cutoffs = NULL){ cu.brk <- gsub(",","_",cu.brk) cu.brk <- gsub("\\(|\\]", "", cu.brk) # names(out) <- paste0("bin__",n, "__", gsub(x = names(cu)[-1], pattern = "\\D", replacement = ""),"__", cu.brk ) - names(out) <- paste0("bin__",n, "__", cu.brk ) + names(out) <- paste0("bin_",n, "_", cu.brk ) for(i in 1:terra::nlyr(out)){ attr(out[[i]], "deriv.bin") <- cu[i:(i+1)] } diff --git a/R/utils-scenario.R b/R/utils-scenario.R index cc5ae245..d3d55142 100644 --- a/R/utils-scenario.R +++ b/R/utils-scenario.R @@ -27,7 +27,9 @@ interpolate_gaps <- function(env, date_interpolation = "annual"){ check_package("dplyr") date_interpolation <- match.arg(date_interpolation, - c("none", "yearly", "annual", "monthly", "daily"), + c("none", "yearly", "year", "annual", + "month","monthly", + "day","daily"), several.ok = FALSE) if(date_interpolation=="none") return(env) @@ -39,16 +41,19 @@ interpolate_gaps <- function(env, date_interpolation = "annual"){ assertthat::assert_that(tzone != "", length(times)>=2) # Interpolate time steps inc <- switch (date_interpolation, - "yearly" = "year", + "yearly" = "year", "year" = "year", "annual" = "year", - "monthly" = "month", - "daily" = "day" + "monthly" = "month", "month" = "month", + "daily" = "day", "day" = "day" ) # Create new time layer new_times <- seq.Date(from = as.Date(times[1],tz = tzone), to = as.Date(times[length(times)],tz = tzone), by = inc) - new_times <- to_POSIXct(new_times) + # Check time zone and convert if POSIXct + if(inherits( terra::time(stars_to_raster(env[1])[[1]]), "POSIXct")){ + new_times <- to_POSIXct(new_times) + } ori.dims <- names(stars::st_dimensions(env)) # Now for each variable, interpolate @@ -61,7 +66,7 @@ interpolate_gaps <- function(env, date_interpolation = "annual"){ o <- Reduce(c, stars_to_raster(env[v]) ) # Create empty copies per times - nt <- new_times[new_times %notin% terra::time(o)] + nt <- new_times[new_times %notin% terra::time(o)] # Ignore existing values oo <- rep(emptyraster(o, vals = NA), length(nt)) terra::time(oo) <- nt oo <- c(o,oo) @@ -82,7 +87,8 @@ interpolate_gaps <- function(env, date_interpolation = "annual"){ # Checks assertthat::assert_that( is.stars(out), - all( names(out) %in% names(env) ) + all( names(out) %in% names(env) ), + length( stars::st_get_dimension_values(out, 3) ) == length(new_times) ) return(out) } diff --git a/man/PredictorDataset-class.Rd b/man/PredictorDataset-class.Rd index 6023648c..fbe0bcee 100644 --- a/man/PredictorDataset-class.Rd +++ b/man/PredictorDataset-class.Rd @@ -9,6 +9,8 @@ This class describes the PredictorDataset and is used to store covariates within } \seealso{ \code{\link[=predictor_derivate]{predictor_derivate()}} + +\code{\link[=predictor_transform]{predictor_transform()}} } \keyword{classes} \section{Public fields}{ @@ -20,6 +22,8 @@ This class describes the PredictorDataset and is used to store covariates within \item{\code{name}}{A name for this object.} +\item{\code{transformed}}{Saves whether the predictors have been transformed somehow.} + \item{\code{timeperiod}}{A timeperiod field} } \if{html}{\out{}} @@ -44,6 +48,7 @@ This class describes the PredictorDataset and is used to store covariates within \item \href{#method-PredictorDataset-show}{\code{PredictorDataset$show()}} \item \href{#method-PredictorDataset-summary}{\code{PredictorDataset$summary()}} \item \href{#method-PredictorDataset-has_derivates}{\code{PredictorDataset$has_derivates()}} +\item \href{#method-PredictorDataset-is_transformed}{\code{PredictorDataset$is_transformed()}} \item \href{#method-PredictorDataset-length}{\code{PredictorDataset$length()}} \item \href{#method-PredictorDataset-ncell}{\code{PredictorDataset$ncell()}} \item \href{#method-PredictorDataset-plot}{\code{PredictorDataset$plot()}} @@ -56,7 +61,7 @@ This class describes the PredictorDataset and is used to store covariates within \subsection{Method \code{new()}}{ Initializes the object and creates an empty list \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{PredictorDataset$new(id, data, ...)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{PredictorDataset$new(id, data, transformed = FALSE, ...)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -66,6 +71,8 @@ Initializes the object and creates an empty list \item{\code{data}}{A predictor dataset usually as \code{\link{SpatRaster}}.} +\item{\code{transformed}}{A \code{\link{logical}} flag if predictors have been transformed. Assume not.} + \item{\code{...}}{Any other parameters found.} } \if{html}{\out{}} @@ -344,6 +351,19 @@ Indication if there are any predictors that are derivates of outers \if{html}{\out{
}}\preformatted{PredictorDataset$has_derivates()}\if{html}{\out{
}} } +\subsection{Returns}{ +A \code{\link{logical}} flag. +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PredictorDataset-is_transformed}{}}} +\subsection{Method \code{is_transformed()}}{ +Predictors have been transformed? +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{PredictorDataset$is_transformed()}\if{html}{\out{
}} +} + \subsection{Returns}{ A \code{\link{logical}} flag. } diff --git a/man/add_constraint_dispersal.Rd b/man/add_constraint_dispersal.Rd index eec8257d..24d77aec 100644 --- a/man/add_constraint_dispersal.Rd +++ b/man/add_constraint_dispersal.Rd @@ -32,6 +32,18 @@ constraints.} Add dispersal constraint to an existing \code{scenario} } \details{ +\strong{Dispersal}: +Parameters for \code{'method'}: +\itemize{ +\item \code{sdd_fixed} - Applies a fixed uniform dispersal distance per modelling timestep. +\item \code{sdd_nexpkernel} - Applies a dispersal distance using a negative exponential kernel from its origin. +\item \code{kissmig} - Applies the kissmig stochastic dispersal model. Requires \code{`kissmig`} package. Applied at each modelling time step. +\item \code{migclim} - Applies the dispersal algorithm MigClim to the modelled objects. Requires \code{"MigClim"} package. +} + +A comprehensive overview of the benefits of including dispersal constrains in +species distribution models can be found in Bateman et al. (2013). + The following additional parameters can bet set: \itemize{ \item \code{pext}: \code{\link{numeric}} indicator for \code{`kissmig`} of the probability a @@ -41,6 +53,13 @@ colonized cell becomes uncolonised, i.e., the species gets locally extinct 3x3 neighbourhood (Default: \code{0.2}). } } +\references{ +\itemize{ +\item Bateman, B. L., Murphy, H. T., Reside, A. E., Mokany, K., & VanDerWal, J. (2013). +Appropriateness of full‐, partial‐and no‐dispersal scenarios in climate change impact +modelling. Diversity and Distributions, 19(10), 1224-1234. +} +} \seealso{ Other constraint: \code{\link{add_constraint}()}, diff --git a/man/add_predictors.Rd b/man/add_predictors.Rd index 46186b67..6008379f 100644 --- a/man/add_predictors.Rd +++ b/man/add_predictors.Rd @@ -163,10 +163,10 @@ remove extreme outliers. Available options for creating derivates are: \itemize{ \item \code{'none'} - No additional predictor derivates are created. -\item \code{'quad'} - Adds quadratic transformed predictors. +\item \code{'quad'} - Adds quadratic derivate predictors. \item \code{'interaction'} - Add interacting predictors. Interactions need to be specified (\code{"int_variables"})! -\item \code{'thresh'} - Add threshold transformed predictors. -\item \code{'hinge'} - Add hinge transformed predictors. +\item \code{'thresh'} - Add threshold derivate predictors. +\item \code{'hinge'} - Add hinge derivate predictors. \item \code{'bin'} - Add predictors binned by their percentiles. } } diff --git a/tests/testthat/test_Scenarios.R b/tests/testthat/test_Scenarios.R index d5d502f7..e652c377 100644 --- a/tests/testthat/test_Scenarios.R +++ b/tests/testthat/test_Scenarios.R @@ -61,21 +61,34 @@ test_that('Testing data prep functions for spatial-temporal data in stars', { new <- interpolate_gaps(pred_future, date_interpolation = "annual") ) expect_length(new, 9) - expect_length(stars::st_get_dimension_values(new, 3), 81) + expect_length(stars::st_get_dimension_values(new, "time"), 81) # --- # # Create derivates of stars data for testing test <- pred_future[1] expect_length(test, 1) + # Nothing expect_no_error(new <- predictor_derivate(test, option = "none")) - expect_length(new, 1) + expect_length(new, 1); expect_true(utils::hasName(new,"bio01")) + # Quadratic transform expect_no_error(new <- predictor_derivate(test, option = "quad")) - expect_length(new, 2) + expect_length(new, 2); expect_true(utils::hasName(new,"bio01")) + # Hinge transform expect_no_error(new <- predictor_derivate(test, option = "hinge", nknots = 4)) - expect_length(new, 5) - # expect_no_error(new <- predictor_derivate(test, option = "thresh",nknots = 4)) - # expect_length(new, 2) - # expect_no_error(new <- predictor_derivate(test, option = "bin",nknots = 4)) + expect_length(new, 5); expect_true(utils::hasName(new,"bio01")) + # Threshold transform + expect_no_error(new <- predictor_derivate(test, option = "thresh",nknots = 4)) + expect_length(new, 4);expect_true(utils::hasName(new,"bio01")) + # Binning + expect_no_error(new <- predictor_derivate(test, option = "bin",nknots = 4)) + expect_length(new, 2);expect_true(utils::hasName(new,"bio01")) + # Interaction + test <- pred_future[1:2] + expect_length(test, 2) + expect_error(new <- predictor_derivate(test, option = "interaction")) + expect_no_error(new <- predictor_derivate(test, option = "interaction", int_variables = names(test))) + expect_length(new, 3) + expect_true(utils::hasName(new,"bio01"));expect_true(utils::hasName(new,"bio12")) # --- # # Create threshold diff --git a/vignettes/articles/01_data_preparationhelpers.Rmd b/vignettes/articles/01_data_preparationhelpers.Rmd index 03b4294c..186b4692 100644 --- a/vignettes/articles/01_data_preparationhelpers.Rmd +++ b/vignettes/articles/01_data_preparationhelpers.Rmd @@ -389,3 +389,29 @@ stars::st_get_dimension_values(new, 3) ``` +### Derivates of scenario predictors + +As with `SpatRaster` covariates, it is also possible to create transformed or +derivate versions of scenario predictors through the same helper function. + +```{r} +# Load some stars rasters +ll <- list.files(system.file('extdata/predictors_presfuture/', + package = 'ibis.iSDM', + mustWork = TRUE),full.names = TRUE) + +# Load the same files future ones +suppressWarnings( + pred_future <- stars::read_stars(ll[1:3]) |> dplyr::slice('Time', seq(1, 86, by = 10)) +) +sf::st_crs(pred_future) <- sf::st_crs(4326) + +# Scale +new <- predictor_transform(pred_future, option = "scale") + +# Add hinge transformed variables to the scenario object +new2 <- predictor_derivate(pred_future, option = "hinge") +# New variable names +names(new2) + +``` diff --git a/vignettes/articles/03_integrate_data.Rmd b/vignettes/articles/03_integrate_data.Rmd index a84c54ea..778146a5 100644 --- a/vignettes/articles/03_integrate_data.Rmd +++ b/vignettes/articles/03_integrate_data.Rmd @@ -340,8 +340,7 @@ mod1 <- distribution(background) |> # INLA requires the specification of a mesh which in this example is generated from the data. engine_inlabru() |> # Train - train(runname = "Combined prediction", only_linear = TRUE, - method_integration = "predictor") + train(runname = "Combined prediction", only_linear = TRUE) # The resulting object contains the combined prediction with shared coefficients among datasets. plot(mod1)