Further work latent effects and zoning support in projections #121

Author: Martin-Jung

DESCRIPTION

@@ -84,7 +84,7 @@ Suggests:
     xgboost
 URL: https://iiasa.github.io/ibis.iSDM/
 BugReports: https://github.com/iiasa/ibis.iSDM/issues
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Config/testthat/edition: 3
 Roxygen: list(markdown = TRUE)
 Biarch: true

NEWS.md

@@ -1,6 +1,7 @@
 # ibis.iSDM 0.1.4 (current dev branch)
 
 #### New features
+* Support for carnying over latent spatial effects (`add_latent_spatial()`) to `scenario()` projections.
 * Convenience functions to remove limits and controls `rm_limits()`/`rm_control()` #121
 * :fire: Enable stars and multi-temporal SpatRaster zones for `scenario()` and `distribution()` #121 
 

R/add_latent.R

@@ -110,6 +110,49 @@ methods::setMethod(
   }
 )
 
+#### Latent for scenarios ####
+#' Add latent spatial effect to the model equation
+#' @param layer A [`SpatRaster`] layer describing alternative latent effects to be used instead
+#' if \code{"reuse_latent"} is set to \code{FALSE}.
+#' @param reuse_latent A [`logical`] flag on whether any latent effects found
+#' in the fitted model should be reused (Default \code{TRUE}).
+#' @rdname add_latent_spatial
+methods::setMethod(
+  "add_latent_spatial",
+  methods::signature(x = "BiodiversityScenario"),
+    function(x, layer = NULL, reuse_latent = TRUE, ...) {
+      assertthat::assert_that(
+        inherits(x, "BiodiversityScenario"),
+        is.Raster(layer) || is.null(layer),
+        is.logical(reuse_latent)
+      )
+
+      # Check existing model if latent factors are set there
+      if(reuse_latent){
+        # Try and see if model has been defined?
+        fit <- x$get_model()
+        if(!is.Waiver(fit)){
+          # Latent-factors specified?
+          check <- fit$has_latent()
+          assertthat::assert_that(
+            check,
+            msg = "Latent factors are to be reused, but could not be found in fitted model?"
+          )
+        }
+      }
+
+      # Build latent object
+      lat <- list("layer" = layer, "reuse_latent" = reuse_latent)
+
+      # Make a clone copy of the object
+      y <- x$clone(deep = TRUE)
+
+      # Add to data to the BiodiversityDistribution object
+      y$set_latent(lat)
+    }
+)
+
+
 #' Function to remove a latent effect
 #'
 #' @description This is just a wrapper function for removing specified offsets
@@ -156,3 +199,23 @@ methods::setMethod(
     y$rm_latent()
   }
 )
+
+#' @rdname rm_latent
+methods::setMethod(
+  "rm_latent",
+  methods::signature(x = "BiodiversityScenario"),
+  function(x) {
+    assertthat::assert_that(inherits(x, "BiodiversityScenario") )
+    # If no offset can be found, just return proto object
+    if(is.Waiver(x$latentfactors)){ return(x) }
+
+    # Messenger
+    if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Setup]','yellow','Removing latent effects.')
+
+    # Make a clone copy of the object
+    y <- x$clone(deep = TRUE)
+
+    # Now remove the offset
+    y$rm_latent()
+  }
+)

R/add_predictors.R

@@ -395,7 +395,7 @@ methods::setMethod(
         r <- o[[n]]
         # If predictor transformation is specified, apply
         if(transform != "none") r <- predictor_transform(r, option = transform)
-        y$predictors <- y$predictors$set_data(n, r)
+        y$predictors <- y$predictors$set_data(r)
         rm(r)
       }
     }
@@ -504,7 +504,8 @@ methods::setMethod(
     if(is.Waiver(x$predictors)){
       y <- add_predictors(y, env = layer, transform = 'none',derivates = 'none', priors)
     } else {
-      y$predictors <- y$predictors$set_data('range_distance', layer)
+      names(layer) <- "range_distance"
+      y$predictors <- y$predictors$set_data(layer)
       if(!is.null(priors)) {
         # FIXME: Ideally attempt to match varnames against supplied predictors vis match.arg or similar
         assertthat::assert_that( all( priors$varnames() %in% names(layer) ) )
@@ -608,7 +609,8 @@ methods::setMethod(
     if(is.Waiver(x$predictors)){
       y <- add_predictors(y, env = dis, transform = 'none',derivates = 'none')
     } else {
-      y$predictors <- y$predictors$set_data('range_distance', dis)
+      names(dis) <- "range_distance"
+      y$predictors <- y$predictors$set_data(dis)
     }
     return(y)
   }

R/class-biodiversityscenario.R

@@ -21,12 +21,14 @@ BiodiversityScenario <- R6::R6Class(
     #' @field limits A [`sf`] object used to constraint the prediction.
     #' @field predictors A predictor object for projection.
     #' @field constraints Any constraints set for projection.
+    #' @field latentfactors A [`list`] on whether latentfactors are used.
     #' @field scenarios The resulting [`stars`] objects.
     modelobject = new_waiver(),
     modelid = new_waiver(),
     limits = new_waiver(),
     predictors = new_waiver(),
     constraints = new_waiver(),
+    latentfactors = new_waiver(),
     scenarios = new_waiver(),
 
     #' @description
@@ -60,6 +62,10 @@ BiodiversityScenario <- R6::R6Class(
       # Thresholds
       tr <- self$get_threshold()
 
+      # Latent factors if set
+      lat <- self$get_latent()
+      if(!is.Waiver(lat))
+
       # Any other simulation outputs modules
       simmods <- self$get_simulation()
 
@@ -71,6 +77,7 @@ BiodiversityScenario <- R6::R6Class(
         "\n  Time period:    ", tp,
         ifelse(!is.Waiver(cs)||!is.Waiver(tr), "\n --------- ", ""),
         ifelse(is.Waiver(cs),"", paste0("\n  Constraints:    ", text_green(paste(paste0(names(cs),' (',cs,')'),collapse = ', ')) ) ),
+        ifelse(is.Waiver(lat),"", paste0("\n  Latent factors:    ", text_green("<Found>")) ),
         ifelse(is.Waiver(tr),"", paste0("\n  Threshold:      ", round(tr[1], 3),' (',names(tr[1]),')') ),
         ifelse(is.Waiver(simmods),"", paste0("\n  Simulations:    ", text_green(paste(paste0(names(simmods),' (',simmods[[1]][[1]],')'),collapse = ', ')) ) ),
         "\n --------- ",
@@ -350,6 +357,36 @@ BiodiversityScenario <- R6::R6Class(
       return( self )
     },
 
+    #' @description
+    #' Adding latent factors to the object.
+    #' @note
+    #' The latent factor is usually obtained from the fitted model object, unless
+    #' re-specified and added here to the list.
+    #' @param latent A [`list`] containing the data object.
+    #' @seealso [add_latent_spatial()]
+    #' @return This object.
+    set_latent = function(latent){
+      assertthat::assert_that(is.list(latent))
+      self$latentfactors <- latent
+      return(self)
+    },
+
+    #' @description
+    #' 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
+    },
+
+    #' @description
+    #' Remove latent factors if found in object.
+    #' @return This object.
+    rm_latent = function(){
+      self$latentfactors <- new_waiver()
+      return(self)
+    },
+
     #' @description
     #' Plot the predictions made here.
     #' @param what A [`character`] describing the layers to be plotted.

R/class-distributionmodel.R

@@ -655,6 +655,19 @@ DistributionModel <- R6::R6Class(
       }
     },
 
+    #' @description
+    #' Logical indication if the prediction has added latent factors.
+    #' @return A [`logical`] flag.
+    has_latent = function(){
+      # Check for settings here
+      settings <- self$settings
+      if(!is.Waiver(settings)){
+        return(
+          settings$get('has_latent')
+        )
+      }
+    },
+
     #' @description
     #' Has a offset been used?
     #' @return A [`logical`] flag.

R/class-predictors.R

@@ -279,14 +279,20 @@ PredictorDataset <- R6::R6Class(
 
     #' @description
     #' Add a new Predictor dataset to this collection
-    #' @param x [`character`] of the new name to be stored.
-    #' @param value A new [`SpatRaster`] object.
+    #' @param value A new [`SpatRaster`] or [`stars`] object.
     #' @return This object
-    set_data = function(x, value){
-      assertthat::assert_that(assertthat::is.string(x),
-                              is.Raster(value),
-                              is_comparable_raster(self$get_data(), value))
-      self$data <- suppressWarnings( c(self$get_data(), value) )
+    set_data = function(value){
+      assertthat::assert_that(is.Raster(value) || inherits(value, "stars"))
+      if(is.Raster(self$get_data()) && is.Raster(value)){
+        assertthat::assert_that(is_comparable_raster(self$get_data(), value))
+        self$data <- suppressWarnings( c(self$get_data(), value) )
+      } else if(inherits(self$get_data(), "stars") && is.Raster(value)) {
+        # Stars assumed
+        self$data <- st_add_raster(self$get_data(), value)
+      } else {
+        # Simply try combining them
+        self$data <- suppressWarnings( c(self$get_data(), value) )
+      }
       return(self)
     },
 

R/project.R

@@ -139,7 +139,34 @@ methods::setMethod(
       new <- interpolate_gaps(new_preds$get_data(),
                               date_interpolation = date_interpolation)
       # Set new data
-      new_preds$set_data(new)
+      new_preds <- new_preds$set_data(new)
+    }
+
+    # Check latent factors and add them to future if found
+    if(!is.Waiver(mod$latentfactors) || fit$has_latent()){
+      # Get predictor names
+      pn <- fit$model$predictors_names
+      ind <- grep("nearestpoint_|spatialtrend_|kde__coordinates", pn,value = TRUE)
+      if(length(mod$latentfactors)==2){
+        # Get the layer if found
+        if(!is.null(mod$latentfactors$layer) && length(ind)>0){
+          layer <- mod$latentfactors$layer
+          new_preds <- new_preds$set_data(layer)
+        } else if(length(ind)>0){
+          # Get the predictors
+          sps <- fit$model$predictors_object$get_data() |> terra::subset(ind)
+          new_preds <- new_preds$set_data(sps)
+        }
+        if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Scenario]','green', "Found latent factors and added them to future predictors stack.")
+      } else {
+        if(length(ind)>0){
+          # Get the predictors
+          sps <- fit$model$predictors_object$get_data() |> terra::subset(ind)
+          new_preds <- new_preds$set_data(sps)
+        }
+        if(getOption('ibis.setupmessages', default = TRUE)) myLog('[Scenario]','green', "Found latent factors and added them to future predictors stack.")
+      }
+      try({rm(pn,ind,sps)},silent = TRUE)
     }
 
     # Get limits if flagged as present in scenario object
@@ -149,7 +176,8 @@ methods::setMethod(
 
       # Clip new predictors
       # This also takes account of nearest time zones
-      new_preds$crop_data(limits, apply_time = TRUE)
+      new_preds$crop_data(limits,
+                          apply_time = ifelse(any(limits$time == "No set time"), FALSE, TRUE))
 
       # Filter to first time slot if present for remaining clipping
       if(utils::hasName(limits,"time")){

R/scenario.R

@@ -61,19 +61,23 @@ methods::setMethod(
     # Convert limits if provided
     if(!is.null(limits)){
       # Convert to polygon if raster
-      if(inherits(limits,'SpatRaster')){
+      if(is.Raster(limits)){
         # Remove 0 from ratified raster assuming this is no-data
         limits[limits == 0] <- NA
         limits <- terra_to_sf(limits)
       }
       # Ensure that limits has the same projection as background
       if(sf::st_crs(limits) != sf::st_crs(fit$model$background)) limits <- sf::st_transform(limits, fit$model$background)
       # Ensure that limits is intersecting the background
-      if(suppressMessages(length( sf::st_intersects(limits, fit$model$background)))==0) { limits <- NULL; warning('Provided limits do not intersect the background!') }
+      if(suppressWarnings( suppressMessages(length( sf::st_intersects(limits, fit$model$background)))==0)) {
+        limits <- NULL; warning('Provided limits do not intersect the background!')
+      }
 
       # Rename geometry just to be sure
-      limits <- rename_geometry(limits, "geometry")
-      names(limits)[1] <- "limit" # The first entry is the layer name
+      if(inherits(limits, 'sf')){
+        limits <- rename_geometry(limits, "geometry")
+        names(limits)[1] <- "limit" # The first entry is the layer name
+      }
     }
 
     # Also check if limits are to be reused if found

R/train.R

@@ -314,7 +314,10 @@ methods::setMethod(
           pred <- model$predictors_object$get_data(df = FALSE)
           new <-  fill_rasters(coords_poly, emptyraster(pred))
           for(val in names(new)){
-            model$predictors_object <- model$predictors_object$set_data(val, new[[val]] )
+            v <- new[[val]]
+            names(v) <- val
+            model$predictors_object <- model$predictors_object$set_data(v)
+            rm(v)
           }
           rm(pred, new)
         } else if(m == "kde") {
@@ -376,7 +379,13 @@ methods::setMethod(
 
         }
       }
-    } else { model[["latent"]] <- new_waiver() }# End of latent factor loop
+
+      # Make a flag to the settings
+      settings$set('has_latent', TRUE)
+    } else {
+      model[["latent"]] <- new_waiver()
+      settings$set('has_latent', FALSE)
+    }# End of latent factor loop
 
     # Set offset if existing
     if(!is.Waiver(x$offset)){
@@ -415,7 +424,7 @@ methods::setMethod(
           settings$set("bias_value", bias$bias_value )
           # Check that variable is already in the predictors object
           if(!(names(bias$layer) %in% model$predictors_names)){
-            model$predictors_object <- model$predictors_object$set_data(names(bias$layer), bias$layer)
+            model$predictors_object <- model$predictors_object$set_data(bias$layer)
             # Also set predictor names
             model[['predictors_names']] <- model$predictors_object$get_names()
             model[['predictors']] <- model$predictors_object$get_data(df = TRUE, na.rm = FALSE)
@@ -727,6 +736,8 @@ methods::setMethod(
                                          layer)
           )
         )
+        # Save the zones categories for later too!
+        settings$set("limits_zones_categories", unique(zones$limit))
         # Limit zones (using the full layer here!)
         zones <- subset(layer, limit %in% unique(zones$limit) )
       } else if(x$get_limits()$limits_method %in% c("nt2", "mess")){
@@ -818,13 +829,6 @@ methods::setMethod(
             point_in_polygon(poly = model$background, points = model$predictors[,c('x','y')] )[['limit']]
           )),model$predictors_names] <- NA # Fill with NA
         }
-        # The same with offset if specified, Note this operation below is computationally quite costly
-        # MJ: 18/10/22 Removed below as (re)-extraction further in the pipeline makes this step irrelevant
-        # if(!is.Waiver(x$offset)){
-        #   model$offset[which( is.na(
-        #     point_in_polygon(poly = zones, points = model$offset[,c('x','y')] )[['limit']]
-        #   )), "spatial_offset" ] <- NA # Fill with NA
-        # }
       }
       # Reset the zones, but save the created layer
       # Note that we are using the original saved layer here again!
@@ -837,6 +841,8 @@ methods::setMethod(
                 "mcp_buffer" = x$limits$mcp_buffer,
                 "limits_clip" = x$limits$limits_clip)
       settings$set("limits", l)
+      # Save the zones categories for later too!
+      settings$set("limits_zones_categories", unique(zones$limit))
       x <- x$set_limits(x = l)
       rm(zones)
     }