Merge pull request #135 from iiasa/dev

Merging of current `dev`

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: ibis.iSDM
 Type: Package
 Title: Modelling framework for integrated biodiversity distribution scenarios
-Version: 0.1.4
+Version: 0.1.5
 Authors@R: 
     c(person(given = "Martin",
            family = "Jung",
@@ -32,16 +32,17 @@ Imports:
     R6 (>= 2.5.0),
     assertthat (>= 0.2.0),
     doFuture (>= 0.12.2),
-    dplyr,
-    foreach,
     future (>= 1.23.0),
+    parallelly (>= 1.30.0),
+    parallel,
+    foreach,
+    dplyr,
     geodist,
     ggplot2,
     graphics,
     methods,
     Matrix,
     ncdf4,
-    parallel,
     posterior,
     sf (>= 1.0),
     stars (>= 0.5),
@@ -60,7 +61,6 @@ Suggests:
     cubelyr,
     dbarts (>= 0.9-22),
     deldir,
-    doParallel,
     ellipsis,
     glmnet (>= 4.1),
     glmnetUtils,

NAMESPACE

@@ -105,8 +105,11 @@ export(get_ngbvalue)
 export(get_priors)
 export(get_rastervalue)
 export(ibis_dependencies)
+export(ibis_enable_parallel)
 export(ibis_future)
 export(ibis_options)
+export(ibis_set_strategy)
+export(ibis_set_threads)
 export(interpolate_gaps)
 export(is.Id)
 export(is.Raster)
@@ -119,8 +122,11 @@ export(mask.BiodiversityDatasetCollection)
 export(mask.BiodiversityScenario)
 export(mask.DistributionModel)
 export(mask.PredictorDataset)
+export(modal)
 export(new_id)
 export(new_waiver)
+export(nicheplot)
+export(objects_size)
 export(partial)
 export(partial.DistributionModel)
 export(partial_density)
@@ -140,6 +146,7 @@ export(rm_limits)
 export(rm_offset)
 export(rm_predictors)
 export(rm_priors)
+export(run_parallel)
 export(run_stan)
 export(sanitize_names)
 export(scenario)
@@ -166,3 +173,4 @@ import(terra)
 importFrom(foreach,"%do%")
 importFrom(foreach,"%dopar%")
 importFrom(stats,effects)
+importFrom(utils,object.size)

NEWS.md

@@ -1,4 +1,18 @@
-# ibis.iSDM 0.1.4 (current dev branch)
+# ibis.iSDM 0.1.5 (current dev branch)
+
+#### New features
+* New visualization function `nicheplot()` to visualize suitability across 2 axes #87.
+* Support for 'modal' value calculations in `ensemble()`.
+* Support for 'superlearner' in `ensemble()`.
+* Support for 'kmeans' derived threshold calculation in `threshold()` and `predictor_derivate()`.
+* Support for future processing streamlined. See FAQ section for instructions #18.
+
+#### Minor improvements and bug fixes
+* Now overwriting temporary data by default in `predictor_transform()` and similar functions.
+* Minor :bug: fix related to misaligned thresholds and negative exponential kernels.
+* :fire: :bug: fix for scenario projections that use different grain sizes than for inference.
+
+# ibis.iSDM 0.1.4 
 
 #### New features
 * Support for carnying over latent spatial effects (`add_latent_spatial()`) to `scenario()` projections.

R/add_constraint.R

@@ -382,8 +382,11 @@ methods::setMethod(
   # Divide alpha values by 2
   alpha <- value/2
 
+  # Scale value for different projections
+  value_scale <- ifelse(terra::is.lonlat(baseline_threshold), terra::res(baseline_threshold)[1] * 10000, 1)
+
   # Grow baseline raster by using an exponentially weighted kernel
-  ras_dis <- terra::gridDist(baseline_threshold, target = 1)
+  ras_dis <- terra::gridDist(baseline_threshold, target = 1, scale = value_scale)
   # Normalized (with a constant) negative exponential kernel
   ras_dis <- terra::app(ras_dis, fun = function(x) (1 / (2 * pi * value ^ 2)) * exp(-x / value) )
   # Equivalent to alpha = 1/value and
@@ -847,12 +850,15 @@ methods::setMethod(
 
     # Rasterize the layer
     # First try and dig out a layer from a predictor dataset if found
-    if(inherits( mod$get_predictors(), "PredictorDataSet")){
-      ras <- mod$get_predictors()$get_data() |> stars_to_raster()
-      ras <- ras[[1]]
+    if(inherits( mod$get_predictors(), "PredictorDataset")){
+      ras <- mod$get_predictors()$get_data()
+      if(inherits(ras, 'stars')){
+        ras <- stars_to_raster(ras)[[1]]
+      }
     } else {
       # Try and get the underlying model and its predictors
       ras <- mod$get_model()$get_data()
+      if(is.null(ras)) ras <- emptyraster(mod$get_model()$model$predictors_object$get_data())
     }
     assertthat::assert_that(is.Raster(ras))
     bb <- try({ terra::rasterize(layer, ras, 1)}, silent = TRUE)

R/add_control_bias.R

@@ -77,7 +77,7 @@
 #' estimating spatial sampling effort and habitat suitability for multiple species
 #' from opportunistic presence‐only data. Methods in Ecology and Evolution, 12(5), 933-945.
 #'
-#' @seealso [add_control_extrapolation()]
+#' @seealso [add_limits_extrapolation()]
 #' @keywords bias offset control
 #' @concept The spatial bias weighting was inspired by code in the \code{enmSdmX} package.
 #'

R/add_predictors.R

@@ -60,6 +60,7 @@ NULL
 #' * \code{'interaction'} - Add interacting predictors. Interactions need to be specified (\code{"int_variables"})!
 #' * \code{'thresh'} - Add threshold derivate predictors.
 #' * \code{'hinge'} - Add hinge derivate predictors.
+#' * \code{'kmeans'} - Add k-means derived factors.
 #' * \code{'bin'} - Add predictors binned by their percentiles.
 #'
 #' @note
@@ -130,7 +131,7 @@ methods::setMethod(
     assertthat::assert_that(inherits(x, "BiodiversityDistribution"),
                             is.Raster(env),
                             all(transform == 'none') || all( transform %in% c('pca', 'scale', 'norm', 'windsor') ),
-                            all(derivates == 'none') || all( derivates %in% c('thresh', 'hinge', 'quadratic', 'bin', 'interaction') ),
+                            all(derivates == 'none') || all( derivates %in% c('thresh', 'hinge', 'quadratic', 'bin', 'kmeans', 'interaction') ),
                             is.vector(derivate_knots) || is.numeric(derivate_knots),
                             is.null(names) || assertthat::is.scalar(names) || is.vector(names),
                             is.logical(explode_factors),
@@ -247,7 +248,7 @@ methods::setMethod(
 
     # Mask predictors with existing background layer
     if(bgmask){
-      env <- terra::mask(env, mask = x$background)
+      env <- terra::mask(env, mask = x$background, overwrite = TRUE)
       # Reratify, work somehow only on stacks
       if(has_factors && any(is.factor(env)) ){
         new_env <- env
@@ -349,7 +350,7 @@ methods::setMethod(
       # If it is a raster
       if(is.Raster(x$background)){
         # Check that background and range align, otherwise raise error
-        if(is_comparable_raster(layer, x$background)){
+        if(!is_comparable_raster(layer, x$background)){
           warning('Supplied range does not align with background! Aligning them now...')
           layer <- alignRasters(layer, x$background, method = 'bilinear', func = mean, cl = FALSE)
         }
@@ -371,12 +372,12 @@ methods::setMethod(
     if(terra::global(ras1, "min", na.rm = TRUE) == terra::global(ras1, "max", na.rm = TRUE)){
       o <- ras2
       # Ensure that all layers have a minimum and a maximum
-      o[is.na(o)] <- 0; o <- terra::mask(o, x$background)
+      o[is.na(o)] <- 0; o <- terra::mask(o, x$background, overwrite = TRUE)
       names(o) <- c('elev_high')
     } else {
       o <- c(ras1, ras2)
       # Ensure that all layers have a minimum and a maximum
-      o[is.na(o)] <- 0; o <- terra::mask(o, x$background)
+      o[is.na(o)] <- 0; o <- terra::mask(o, x$background, overwrite = TRUE)
       names(o) <- c('elev_low', 'elev_high')
     }
     rm(ras1,ras2)
@@ -552,7 +553,8 @@ methods::setMethod(
     #     ras_range <- raster::rasterize(layer, temp, field = 1, background = NA)
     #   }
     # } else {
-    ras_range <- terra::rasterize(layer, temp, field = 1, background = 0)
+    ras_range <- terra::rasterize(layer, temp, field = 1,
+                                  background = 0, overwrite = TRUE)
     # }
 
     # -------------- #
@@ -564,8 +566,8 @@ methods::setMethod(
       names(dis) <- 'binary_range'
     } else if(method == 'distance'){
       # Calculate the linear distance from the range
-      dis <- terra::gridDist(ras_range, target = 1)
-      dis <- terra::mask(dis, x$background)
+      dis <- terra::gridDist(ras_range, target = 1, overwrite = TRUE)
+      dis <- terra::mask(dis, x$background, overwrite = TRUE)
       # If max distance is specified
       if(!is.null(distance_max) && !is.infinite(distance_max)){
         dis[dis > distance_max] <- NA # Set values above threshold to NA
@@ -580,7 +582,7 @@ methods::setMethod(
 
       # Set NA to 0 and mask again
       dis[is.na(dis)] <- 0
-      dis <- terra::mask(dis, x$background)
+      dis <- terra::mask(dis, x$background, overwrite = TRUE)
       names(dis) <- 'distance_range'
     }
 
@@ -753,7 +755,7 @@ methods::setMethod(
     # names = names = NULL; transform = 'none'; derivates = 'none'; derivate_knots = 4; int_variables = NULL;harmonize_na = FALSE; state = NULL
     # Try and match transform and derivatives arguments
     transform <- match.arg(transform, c('none','pca', 'scale', 'norm', 'windsor', 'percentile'), several.ok = FALSE) # Several ok set to FALSE as states are not working otherwise
-    derivates <- match.arg(derivates, c('none','thresh', 'hinge', 'quadratic', 'bin', 'interaction'), several.ok = TRUE)
+    derivates <- match.arg(derivates, c('none','thresh', 'hinge', 'quadratic', 'bin', 'kmeans', 'interaction'), several.ok = TRUE)
 
     assertthat::validate_that(inherits(env,'stars'), msg = 'Projection rasters need to be stars stack!')
     assertthat::assert_that(inherits(x, "BiodiversityScenario"),
@@ -767,8 +769,8 @@ methods::setMethod(
     assertthat::validate_that(length(env) >= 1)
 
     # Get model object
-    obj <- x$get_model()
-    assertthat::assert_that(!(is.null(obj) || is.Waiver(obj)),
+    obj <- x$get_model(copy = TRUE)
+    assertthat::assert_that(!(isFALSE(obj) || is.Waiver(obj)),
                             msg = "No model object found in scenario?")
     model <- obj$model
 
@@ -856,7 +858,7 @@ methods::setMethod(
       # Get variable names
       varn <- obj$get_coefficients()[,1]
       # Are there any derivates present in the coefficients?
-      if(any( length( grep("hinge_|bin_|quadratic_|thresh_|interaction_", varn ) ) > 0 )){
+      if(any( length( grep("hinge_|bin_|kmeans_|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))){

R/class-biodiversitydistribution.R

@@ -142,7 +142,7 @@ BiodiversityDistribution <- R6::R6Class(
     #' @description
     #' Specify new limits to the background
     #' @param x A [`list`] object with method and limit type.
-    #' @seealso [add_control_extrapolation()]
+    #' @seealso [add_limits_extrapolation()]
     #' @return This object.
     set_limits = function(x){
       # Specify list

R/class-biodiversityscenario.R

@@ -103,6 +103,12 @@ BiodiversityScenario <- R6::R6Class(
       assertthat::assert_that(is.Waiver(self$get_predictors()) || inherits(self$get_predictors(), "PredictorDataset"))
       assertthat::assert_that(is.Waiver(self$get_data()) || (inherits(self$get_data(), "stars") || is.Raster(self$get_data())) )
       assertthat::assert_that(is.Waiver(self$get_constraints()) || is.list(self$get_constraints()))
+      # Check predictor mismatch
+      if(!is.Waiver(self$get_predictors())){
+        ori <- x$get_projection()
+        test <- self$get_projection()
+        assertthat::validate_that(sf::st_crs(test) == sf::st_crs(ori),msg = "Predictor and fitted predictor projections mismatch!")
+      }
       invisible(self)
     },
 
@@ -351,6 +357,15 @@ BiodiversityScenario <- R6::R6Class(
       return(self[[what]])
     },
 
+    #' @description
+    #' Remove scenario predictions
+    #' @param what A [`character`] vector with names of what
+    #' @return Invisible
+    rm_data = function(){
+      self$scenarios <- new_waiver()
+      invisible()
+    },
+
     #' @description
     #' Set new data in object.
     #' @param x A new data object measuing scenarios.
@@ -382,7 +397,6 @@ BiodiversityScenario <- R6::R6Class(
     #' Get latent factors if found in object.
     #' @return A [`list`] with the latent settings
     get_latent = function(){
-      if(is.Waiver(self$latentfactors)) return('None')
       self$latentfactors
     },