Version bump and small 🐛 fixes and edits

DESCRIPTION

@@ -1,21 +1,22 @@
 Package: insights
 Title: An R implementation of the InSiGHTS framework
-Version: 0.3
-Year: 2023
+Version: 0.4
+Year: 2024
 Authors@R: 
     person("Martin", "Jung", , "jung@iiasa.ac.at", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-7569-1390"))
 Description: The package provides an implementation of the InSiGHTS modelling framework for creating climate and land-use indicators. This package acts a simple wrapper to do an area-of-habitat refinements on top of climatic envelope models or species distribution models (SDMs). By default it is assumed that such models are obtained through the ibis.iSDM package and this package contains a simple wrapper to link the two approaches.
 License: CC BY 4.0
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 Imports: 
     sf,
     terra,
     ibis.iSDM,
     stars,
     assertthat,
+    lubridate,
     tibble,
     dplyr
 Remotes: 

NEWS.md

@@ -1,3 +1,7 @@
+# InSiGHTS 0.4
+
+* Bug fixes and more support for future clipping when `stars` is provided.
+
 # InSiGHTS 0.3
 
 * Full ibis.iSDM scenario support

R/data.R

@@ -7,12 +7,30 @@
 #'
 #' @note
 #' This layer is an example Species distribution modelling estimate of *Spermophilus citellus*
+#' Alternatively a layer of the frog *Bombina bombina* is loaded.
 #'
+#' @param timeperiod A [`character`] on which example layer is to be loaded. Options
+#' include \code{"current"} and \code{"future"}.
 #' @format A \code{SpatRaster} object in binary format providing a range.
 #' @author Martin Jung
 #' @keywords internal
 #' @noRd
-load_exampledata <- function(){
-  range <- terra::rast(system.file('extdata/example_range.tif', package='insights',mustWork = TRUE))
+load_exampledata <- function(timeperiod = "current"){
+  assertthat::assert_that(
+    is.character(timeperiod)
+  )
+  # Match
+  timeperiod <- match.arg(timeperiod, c("current", "future"),several.ok = FALSE)
+
+  if(timeperiod == "current"){
+    range <- terra::rast(system.file('extdata/example_range.tif', package='insights',mustWork = TRUE))
+  } else {
+    range <- stars::read_mdim(system.file('extdata/Bombina_bombina__ssp126.nc',
+                                          package = 'insights',mustWork = TRUE))
+    range <- range |> split()
+    # Focus on threshold only
+    range <- range |> dplyr::select(ensemble_threshold)
+  }
+  assertthat::assert_that(inherits(range, 'stars') || inherits(range, 'SpatRaster'))
   return(range)
 }

R/insights_fraction.R

@@ -9,6 +9,9 @@
 #' Optionally also a elevation (\code{elev}) layer and habitat condition (\code{condition}) can be provided to support refinements
 #' by elevational range or habitat condition.
 #'
+#' @note
+#' This function does not do the refinement of land-use fraction to relevant habitats.
+#' This needs to be done by the analyst a-priori.
 #' @param range A [`SpatRaster`] or temporal [`stars`] object describing the estimated distribution of a
 #' biodiversity feature (e.g. species). **Has to be in binary format!**
 #' Alternatively a \code{DistributionModel} fitted with \code{ibis.iSDM} package can be supplied.
@@ -21,6 +24,11 @@
 #' @author Martin Jung
 #' @importClassesFrom terra SpatRaster
 #' @importFrom ibis.iSDM is.Raster
+#' @examples
+#' \dontrun{
+#'  out <- insights_fraction(range, landuse)
+#' }
+#'
 #' @references
 #' * Rondinini, Carlo, and Piero Visconti. "Scenarios of large mammal loss in Europe for the 21st century." Conservation Biology 29, no. 4 (2015): 1028-1036.
 #' * Visconti, Piero, Michel Bakkenes, Daniele Baisero, Thomas Brooks, Stuart HM Butchart, Lucas Joppa, Rob Alkemade et al. "Projecting global biodiversity indicators under future development scenarios." Conservation Letters 9, no. 1 (2016): 5-13.
@@ -55,8 +63,8 @@ methods::setMethod(
     rr <- terra::global(range,"range",na.rm=TRUE)
     assertthat::assert_that(all(rr[["min"]]>=0 ),
                             all(rr[["max"]]<=1 ),
-                            all(apply(terra::unique(range), 2, function(z) length(which(!is.nan(unique(z))))) <= 2),
-                            msg = "Input range has to be in binary format!"
+                            # all(apply(terra::unique(range), 2, function(z) length(which(!is.nan(unique(z))))) <= 2),
+                            msg = "Input range has to be in binary or fractional format!"
                             )
     rm(rr)
     assertthat::assert_that(
@@ -76,11 +84,11 @@ methods::setMethod(
     if(ibis.iSDM::is.Raster(lu) && ibis.iSDM::is.Raster(range)){
       if(!(terra::same.crs(range, lu) && terra::compareGeom(range, lu, stopOnError = FALSE))){
         if(!terra::same.crs(range, lu)){
-          ibis.iSDM::myLog("Preparation", "yellow", "Reprojecting land-use layer to range crs.")
+          ibis.iSDM:::myLog("Preparation", "yellow", "Reprojecting land-use layer to range crs.")
           lu <- terra::project(lu, terra::crs(range))
         }
         if(!terra::compareGeom(range, lu, stopOnError = FALSE)){
-          ibis.iSDM::myLog("Preparation", "yellow", "Cropping and resampling land-use layer(s) to range.")
+          ibis.iSDM:::myLog("Preparation", "yellow", "Cropping and resampling land-use layer(s) to range.")
           lu <- terra::crop(lu, range)
           lu <- terra::resample(lu, range, method = "average", threads = TRUE)
         }
@@ -101,11 +109,11 @@ methods::setMethod(
 
       if(!(terra::same.crs(range, other) && terra::compareGeom(range, other, stopOnError = FALSE))){
         if(!terra::same.crs(range, other)){
-          ibis.iSDM::myLog("Preparation", "yellow", "Reprojecting other layers to range crs.")
+          ibis.iSDM:::myLog("Preparation", "yellow", "Reprojecting other layers to range crs.")
           other <- terra::project(other, terra::crs(range))
         }
         if(!terra::compareGeom(range, other, stopOnError = FALSE)){
-          ibis.iSDM::myLog("Preparation", "yellow", "Cropping and resampling other layer(s) to range.")
+          ibis.iSDM:::myLog("Preparation", "yellow", "Cropping and resampling other layer(s) to range.")
           other <- terra::crop(other, range)
           other <- terra::resample(other, range, method = "average", threads = TRUE)
         }
@@ -227,7 +235,7 @@ methods::setMethod(
     proj <- stars::st_as_stars(proj,
                                crs = sf::st_crs(range)
     )
-    names(proj)
+
     # Reset time dimension for consistency
     dims <- stars::st_dimensions(proj)
     names(dims)[3] <- "time"
@@ -271,7 +279,7 @@ methods::setMethod(
       # Reproejct and rewarp
       other <- other |> sf::st_transform(crs = sf::st_crs(range))
 
-      ibis.iSDM::myLog("[Reprojection]", "yellow", "Aligning other layers to range.")
+      ibis.iSDM:::myLog("[Reprojection]", "yellow", "Aligning other layers to range.")
       grid <- other |> sf::st_bbox() |> stars::st_as_stars()
       other <- other |>
         stars::st_warp(grid, cellsize = stars::st_res(range),use_gdal = FALSE)
@@ -337,6 +345,100 @@ methods::setMethod(
   }
 )
 
+#' @name insights_fraction
+#' @rdname insights_fraction
+#' @usage \S4method{insights_fraction}{stars,stars,ANY,character}(range,lu,other,outfile)
+methods::setMethod(
+  "insights_fraction",
+  methods::signature(range = "stars", lu = "SpatRaster"),
+  function(range, lu, other, outfile = NULL) {
+    assertthat::assert_that(
+      inherits(range, "stars"),
+      ibis.iSDM::is.Raster(lu),
+      missing(other) || (inherits(other, "stars") || ibis.iSDM::is.Raster(other)),
+      is.null(outfile) || is.character(outfile)
+    )
+
+    # Some check
+    assertthat::assert_that(
+      length(range)==1,
+      msg = "More than one layer in range found...?"
+    )
+
+    # Convert if needed
+    if(!missing(other)){
+      # In this case we recreate / warp the raster to dem
+      other <- stars::st_as_stars(other)
+      names(other) <- "other"
+      # Reproject and rewarp
+      other <- other |> sf::st_transform(crs = sf::st_crs(range))
+
+      ibis.iSDM:::myLog("[Reprojection]", "yellow", "Aligning other layers to range.")
+      grid <- other |> sf::st_bbox() |> stars::st_as_stars()
+      other <- other |>
+        stars::st_warp(grid, cellsize = stars::st_res(range),use_gdal = FALSE)
+    }
+
+    # Correct output file extension if necessary
+    if(!is.null(outfile)){
+      assertthat::assert_that(dir.exists(dirname(outfile)),
+                              msg = "Output file directory does not exist!")
+      # Correct output file name depending on type
+      if((inherits(lu, "stars") || ibis.iSDM:::is.Raster(lu))
+         && tolower(tools::file_ext(outfile))!="nc"){
+        outfile <- paste0(outfile, ".nc")
+      }
+    }
+
+    # SpatRaster assumed, sum if too many
+    if(terra::nlyr(lu)>1){
+      lu <- terra::app(lu, 'sum', na.rm = TRUE)
+    }
+
+    # Get dimensions from range
+    times <- stars::st_get_dimension_values(range, 3)
+    assertthat::assert_that(
+      is.numeric(times) || lubridate::is.Date(times) || lubridate::is.POSIXct(times)
+    )
+
+    # --- #
+    # Then convert each time step to a SpatRaster and pass to insights_fraction
+    proj <- terra::rast()
+    for(tt in 1:length(unique(times))){
+      # Make a slice
+      s <- range |> stars:::slice.stars('time', tt)
+      # Convert to raster
+      s <- terra::rast(s)
+      assertthat::assert_that(terra::global(s, "max", na.rm = TRUE)[,1] <=1,
+                              msg = "Values in range larger than 1 found?")
+      o <- insights_fraction(range = s,
+                             lu = lu,
+                             # other = other,
+                             outfile = NULL)
+      suppressWarnings(
+        proj <- c(proj, o)
+      )
+    }
+    # Finally convert to stars and rename
+    proj <- stars::st_as_stars(proj,
+                               crs = sf::st_crs(range)
+    )
+
+    # Reset time dimension for consistency
+    dims <- stars::st_dimensions(proj)
+    names(dims)[3] <- "time"
+    dims$time <- stars::st_dimensions(range)[[3]]
+    stars::st_dimensions(proj) <- dims
+
+    # Return result or write respectively
+    if(is.null(outfile)){
+      return(proj)
+    } else {
+      assertthat::assert_that(inherits(proj, "stars"))
+      stars::write_stars(proj, outfile)
+    }
+  }
+)
 
 #### Implementation for ibis.iSDM predictions and projections ####
 #' @name insights_fraction

README.Rmd

@@ -143,7 +143,7 @@ plot(o$suitability~o$band, type = "b",
 
 ## Citations
 
-Jung M (2023). _insights: An R implementation of the InSiGHTS framework_. R package version 0.2.
+See [CITATION.cff] file
 
 P. Visconti, M. Bakkenes, D. Baisero, T. Brooks, S.H.M. Butchart, L. Joppa, R. Alkemade, M. Di Marco, L. Santini, M. Hoffmann, C. Rondinini
 _Projecting global biodiversity indicators under future development scenarios_ Conserv. Lett., 9 (2016), pp. 5-13 [DOI](https://doi.org/10.1111/conl.12159)

tests/testthat/test_ibis.R

@@ -1,18 +1,13 @@
 # Test that package works and data can be loaded
 test_that('Train a ibis.iSDM model and apply inSights on it', {
 
-  skip_if_not_installed("glmnet")
   skip_if_not_installed("ibis.iSDM")
 
   suppressWarnings( requireNamespace("terra", quietly = TRUE) )
   suppressWarnings( requireNamespace("ibis.iSDM", quietly = TRUE) )
-  suppressWarnings( requireNamespace("glmnet", quietly = TRUE) )
   suppressPackageStartupMessages(
     require("ibis.iSDM")
   )
-  suppressPackageStartupMessages(
-    require("glmnet")
-  )
 
   options("ibis.setupmessages" = FALSE)
 
@@ -32,15 +27,16 @@ test_that('Train a ibis.iSDM model and apply inSights on it', {
 
   # Now train a small little model
   fit <- ibis.iSDM::distribution(background) |>
-    ibis.iSDM::add_biodiversity_poipo(virtual_points) |>
+    ibis.iSDM::add_biodiversity_poipo(virtual_points,field_occurrence = "Observed") |>
     ibis.iSDM::add_predictors(predictors) |>
-    ibis.iSDM::engine_glmnet() |>
+    ibis.iSDM::engine_glm() |>
     ibis.iSDM::train() |>
     ibis.iSDM::threshold(method = "perc", value = .33)
 
   expect_s3_class(fit, "DistributionModel")
   tr <- fit$get_data("threshold_percentile")
   expect_s4_class(tr, "SpatRaster")
+  if(terra::nlyr(tr)>1) tr <- tr[[1]]
 
   # --- #
   # Now apply insights
@@ -61,22 +57,19 @@ test_that('Train a ibis.iSDM model and apply inSights on it', {
   expect_equal(round(terra::global(out, "max", na.rm = TRUE)[,1],3),
                round(terra::global(out2, "max", na.rm = TRUE)[,1],3))
 
-  })
+})
 
 # Test that package works and data can be loaded
 test_that('Make a ibis.iSDM scenario projection and apply InSiGHTS on it', {
 
-  skip_if_not_installed("glmnet")
   skip_if_not_installed("ibis.iSDM")
 
   suppressWarnings( requireNamespace("terra", quietly = TRUE) )
   suppressWarnings( requireNamespace("ibis.iSDM", quietly = TRUE) )
-  suppressWarnings( requireNamespace("glmnet", quietly = TRUE) )
   suppressWarnings( requireNamespace("stars", quietly = TRUE) )
-  suppressPackageStartupMessages( require("glmnet"))
-  suppressPackageStartupMessages( require("terra"))
-  suppressPackageStartupMessages( require("stars"))
-  suppressPackageStartupMessages( require("ibis.iSDM"))
+  suppressWarnings(
+    suppressPackageStartupMessages( require("ibis.iSDM"))
+  )
 
   options("ibis.setupmessages" = FALSE)
 
@@ -104,7 +97,7 @@ test_that('Make a ibis.iSDM scenario projection and apply InSiGHTS on it', {
     ibis.iSDM::add_biodiversity_poipa(virtual_points,
                                       field_occurrence = 'Observed', name = 'Virtual points') |>
     ibis.iSDM::add_predictors(pred_current, transform = 'scale',derivates = "none") |>
-    ibis.iSDM::engine_glmnet()
+    ibis.iSDM::engine_glm()
 
   modf <- ibis.iSDM::train(x, runname = 'Null', verbose = FALSE) |>
     ibis.iSDM::threshold(method = 'percent',value = .3)
@@ -132,5 +125,4 @@ test_that('Make a ibis.iSDM scenario projection and apply InSiGHTS on it', {
   # Summarize
   o <- insights_summary(out)
   expect_s3_class(o, "data.frame")
-
 })