udpate fn and add tests for MakeMaxEntAbundance()

Author: MargaretSiple-NOAA

R/MakeMaxEntAbundance.R

@@ -7,78 +7,88 @@
 #' @param land raster; a mask to be applied to the output usually landmasses
 #' @param mask raster; an additional mask to apply to the output, if needed
 #' @param type character; use "maxnet" for prob suitable habitat or "cloglog" for approximate abundance
-#' @param clamp Logical; just leave this as F; shoudl the covariates be restricted to the range observe in fitting the model?
-#' @param filename a filename to save results, "" writes to active memory
+#' @param clamp Logical; just leave this as F; should the covariates be restricted to the range observe in fitting the model?
+#' @param filename a filename to save results; "" writes to active memory
 #'
 #' @return a raster map with the desired prediction
 #' @export
 #' @importFrom terra values
 #' @importFrom terra setValues
 #'
 #' @examples
-MakeMaxEntAbundance<-function(model,
-                              maxent.stack,
-                              scale.fac=1,
-                              land=NULL,
-                              mask=NULL,
-                              type="cloglog",
-                              clamp=F,
-                              filename=""){
-
-  #correct a common mistake
-  if(is.null(filename)||is.na(filename)){filename<-""}
+MakeMaxEntAbundance <- function(model,
+                                maxent.stack,
+                                scale.fac = 1,
+                                land = NULL,
+                                mask = NULL,
+                                type = "cloglog",
+                                clamp = F,
+                                filename = "") {
+  # correct a common mistake
+  if (is.null(filename) || is.na(filename)) {
+    filename <- ""
+  }
 
   # Identify the type and make the main prediction
-  type=tolower(type)
+  type <- tolower(type)
   # Somewhat counterintuitive, but this is the type if using the cloglog link to make an abundance estimate
-  if(type=="cloglog"){
+  if (type == "cloglog") {
     # since ENMeval 2.0, they got rid of the useful function and I need to make the predictions the long way
-    dat<-terra::values(maxent.stack)
+    dat <- terra::values(maxent.stack)
 
     # using predict with maxnet will quietly remove the NAs, so need to track them manually
-    na.spots<-which(apply(X = dat,MARGIN = 1,FUN = function(x){return(any(is.na(x)))}))
-    dat.spots<-which(seq(1:nrow(dat))%in%na.spots==F)
+    na.spots <- which(apply(X = dat, MARGIN = 1, FUN = function(x) {
+      return(any(is.na(x)))
+    }))
+    dat.spots <- which(seq(1:nrow(dat)) %in% na.spots == F)
 
-    preds<-stats::predict(model,dat[dat.spots,],type="link")
-    preds2<-exp(preds+model$ent)*scale.fac
-    new.vals<-vector(length=nrow(dat))
-    new.vals[na.spots]<-NA
-    new.vals[dat.spots]<-preds2
-    habitat.prediction<-terra::setValues(x = terra::rast(maxent.stack[[1]]),values = new.vals) #terra needs to call just one raster in the stack
+    preds <- stats::predict(model, dat[dat.spots, ], type = "link")
+    preds2 <- exp(preds + model$ent) * scale.fac
+    new.vals <- vector(length = nrow(dat))
+    new.vals[na.spots] <- NA
+    new.vals[dat.spots] <- preds2
+    habitat.prediction <- terra::setValues(x = terra::rast(maxent.stack[[1]]), values = new.vals) # terra needs to call just one raster in the stack bc it's just setting dimensions.
   }
   # this makes a habitat suitability map from a maxnet model
-  if(type=="maxnet"){
-    dat<-terra::values(maxent.stack)
+  if (type == "maxnet") {
+    dat <- terra::values(maxent.stack)
 
     # using predict with maxnet will quietly remove the NAs, so need to track them manually
-    na.spots<-which(apply(X = dat,MARGIN = 1,FUN = function(x){return(any(is.na(x)))}))
-    dat.spots<-which(seq(1:nrow(dat))%in%na.spots==F)
+    na.spots <- which(apply(X = dat, MARGIN = 1, FUN = function(x) {
+      return(any(is.na(x)))
+    }))
+    dat.spots <- which(seq(1:nrow(dat)) %in% na.spots == F)
 
-    preds<-stats::predict(model,dat[dat.spots,],type="cloglog")
-    new.vals<-vector(length=nrow(dat))
-    new.vals[na.spots]<-NA
-    new.vals[dat.spots]<-preds
-    habitat.prediction<-terra::setValues(x = terra::rast(maxent.stack),values = new.vals)
+    preds <- stats::predict(model, dat[dat.spots, ], type = "cloglog")
+    new.vals <- vector(length = nrow(dat))
+    new.vals[na.spots] <- NA
+    new.vals[dat.spots] <- preds
+    habitat.prediction <- terra::setValues(x = terra::rast(maxent.stack[[1]]), values = new.vals)
   }
   # need to add a check to see about the strange problems with the EBS
-  if(is.null(land)==F){
-    habitat.prediction<-raster::extend(x=habitat.prediction,y=land)
+  if (is.null(land) == F) {
+    habitat.prediction <- raster::extend(x = habitat.prediction, y = land)
   }
 
   # For some reason, the crs info isn't always carrying over
   terra::crs(habitat.prediction) <- terra::crs(maxent.stack)
-  if(filename!=""){terra::writeRaster(x = habitat.prediction,filename = filename, overwrite = TRUE)}
+  if (filename != "") {
+    terra::writeRaster(x = habitat.prediction, filename = filename, overwrite = TRUE)
+  }
 
   # Apply additional masks if necessary
-  if(is.null(land)==F){
-    habitat.prediction<-terra::mask(habitat.prediction, land, inverse = T, overwrite = TRUE,
-                                    filename = filename)
+  if (is.null(land) == F) {
+    habitat.prediction <- terra::mask(habitat.prediction, land,
+      inverse = T, overwrite = TRUE,
+      filename = filename
+    )
   }
   # Apply additional masks if necessary
-  if(is.null(mask)==F){
-    habitat.prediction<-terra::mask(habitat.prediction, mask, overwrite = TRUE,
-                                    filename = filename)
+  if (is.null(mask) == F) {
+    habitat.prediction <- terra::mask(habitat.prediction, mask,
+      overwrite = TRUE,
+      filename = filename
+    )
   }
   return(habitat.prediction)
 }
-

man/MakeMaxEntAbundance.Rd

@@ -0,0 +1,40 @@
+library(terra)
+library(maxnet)
+
+# testing data -------------------------------------------------------------
+# Create some presence/absence data to fit a model
+species.data <- readRDS(system.file("test_files", "goa_data_logarea_folds.rds", package = "EFHSDM"))
+data(raster_stack) # not sure if this will work in context. This is GOA data only.
+raster.stack <- terra::rast(raster_stack)
+names(raster.stack) <- c("lon", "lat", "bdepth", "btemp", "slope", "sponge")
+
+maxnet.covars <- c("bdepth", "btemp", "slope")
+cofactors <- c("sponge")
+r.mult <- 0.5
+maxnet.model0 <- FitMaxnet(
+  data = species.data,
+  species = "dogfish",
+  vars = maxnet.covars,
+  facs = cofactors,
+  regmult = r.mult,
+  reduce = TRUE
+)
+
+# tests -------------------------------------------------------------------
+test_that("cloglog prediction works and returns a SpatRaster", {
+  result <- MakeMaxEntAbundance(model = maxnet.model0, maxent.stack = raster.stack,
+                                scale.fac = 1.619, type = "cloglog")
+  expect_s4_class(result, "SpatRaster")
+  expect_equal(ncell(result), ncell(raster.stack[[1]]))
+})
+
+test_that("maxnet prediction type works", {
+  result <- MakeMaxEntAbundance(model = maxnet.model0, maxent.stack = raster.stack,
+                                scale.fac = 1.619, type = "maxnet")
+  expect_s4_class(result, "SpatRaster")
+})
+
+test_that("does not error with NULL or NA filename", {
+  expect_no_error(MakeMaxEntAbundance(maxnet.model0, raster.stack, filename = NULL))
+  expect_no_error(MakeMaxEntAbundance(maxnet.model0, raster.stack, filename = NA))
+})