fix #48 (#49)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: aoh
 Type: Package
-Version: 0.0.2.5
+Version: 0.0.2.6
 Title: Create Area of Habitat Data
 Description: Create Area of Habitat data to characterize species distributions.
     Data are produced following procedures outlined by Brooks et al. (2019)
@@ -44,6 +44,7 @@ Imports:
     raster (>= 3.5-15),
     rvest (>= 1.0.3)
 Suggests:
+    sp (>= 2.0.0),
     testthat (>= 2.0.1),
     knitr (>= 1.2.0),
     roxygen2 (>= 6.1.1),
@@ -52,12 +53,12 @@ Suggests:
     ggmap (>= 2.6.1),
     fields (>= 14.0),
     smoothr (>= 0.2.2),
-    rnaturalearth (>= 0.1.0),
+    rnaturalearth (>= 0.3.3),
     rgdal (>= 1.5.27),
     gdalUtilities (>= 1.2.1),
     archive (>= 1.1.2),
     link2GI (>= 0.5-0),
-    rgrass (>= 0.3-7),
+    rgrass (>= 0.3-8),
     prepr (>= 0.1.9000),
     vcr (>= 0.6.0)
 Depends:

NEWS.md

@@ -1,3 +1,19 @@
+# aoh 0.0.2.6
+
+- Update `create_spp_aoh_data()` and `create_spp_frc_data()` to have a
+  new `rasterize_touches` parameter (#48). This parameter can be toggled so
+  that when rasterizing species' range data, raster cells that partially
+  overlap with any part of the species' range are treated as covered by
+  the species' range. This functionality may be especially useful for
+  species with very small geographic ranges.
+- Fix bug in `create_spp_aoh_data()` and `create_spp_frc_data()` that causes
+  Python errors when using the GDAL engine and a `habitat_data` or
+  `elevation_data` raster that is stored only in memory (and not associated with
+  any file on disk).
+- Fix bug in `create_spp_aoh_data()` and `create_spp_frc_data()` that causes
+  the GRASS engine to throws errors.
+- Update package dependency versions.
+
 # aoh 0.0.2.5
 
 - Fix compatibility with updates to `terra::compareGeom()`.

R/create_spp_aoh_data.R

@@ -91,6 +91,22 @@ NULL
 #'   (see below for details).
 #'   Defaults to `"terra"`.
 #'
+#' @param rasterize_touches `logical` How should `x` (the species' range data)
+#'   be rasterized when overlapped with the elevation and habitat raster data?
+#'   If `rasterize_touches = FALSE`, the species' range data are treated as
+#'   overlapping with a raster cell, if the range data overlap with the
+#'   centroid of the raster cell.
+#'   If `rasterize_touches = TRUE`, then the species' range data are treated as
+#'   overlapping with a raster cell, if the range data overlap with any
+#'   part of the raster cell.
+#'   Since some species' ranges might be too small to overlap
+#'   with the centroid of any raster cells (meaning that the output
+#'   Area of Habitat map does not contain any suitable habitat for the species),
+#'   it may be preferable to use `rasterize_touches = TRUE`.
+#'   Note that `rasterize_touches = TRUE` is not compatible with the GRASS
+#'   engine.
+#'   Defaults to `FALSE` (following Lumbierres *et al.* 2022).
+#'
 #' @param verbose `logical` Should progress be displayed while processing data?
 #'  Defaults to `TRUE`.
 #'
@@ -206,6 +222,10 @@ NULL
 #' habitat of terrestrial vertebrates. *Conservation Biology*, 36, e13851.
 #' \doi{10.1111/cobi.13851}
 #'
+#' Lumbierres M, Dahal PR, Soria CD, Di Marco M, Butchart SHM, Donald PF, and
+#' Rondinini C (2022) Area of Habitat maps for the world’s terrestrial birds
+#' and mammals. *Scientific Data*, 9, 749.
+#'
 #' Robinson N, Regetz J, and Guralnick RP (2014) EarthEnv-DEM90: A nearly-
 #' global, void-free, multi-scale smoothed, 90m digital elevation model from
 #' fused ASTER and SRTM data.
@@ -278,6 +298,7 @@ create_spp_aoh_data <- function(x,
                                 n_threads = 1,
                                 cache_limit = 1000,
                                 engine = "terra",
+                                rasterize_touches = FALSE,
                                 verbose = TRUE) {
   create_spp_data(
     x = x,
@@ -293,6 +314,7 @@ create_spp_aoh_data <- function(x,
     n_threads = n_threads,
     cache_limit = cache_limit,
     engine = engine,
+    rasterize_touches = rasterize_touches,
     verbose = verbose
   )
 }

R/create_spp_data.R

@@ -25,6 +25,7 @@ create_spp_data <- function(x,
                             n_threads = 1,
                             cache_limit = 1000,
                             engine = "terra",
+                            rasterize_touches = FALSE,
                             verbose = TRUE) {
 
   # initialization
@@ -77,9 +78,20 @@ create_spp_data <- function(x,
     assertthat::is.string(engine),
     assertthat::noNA(engine),
     engine %in% c("terra", "gdal", "grass"),
+    assertthat::is.flag(rasterize_touches),
+    assertthat::noNA(rasterize_touches),
     assertthat::is.flag(verbose),
     assertthat::noNA(verbose)
   )
+  if (isTRUE(rasterize_touches)) {
+    assertthat::assert_that(
+      !identical(engine, "grass"),
+      msg = paste0(
+        "argument to \"rasterize_touches\" is TRUE, and so",
+        "argument to \"engine\" cannot be \"grass\""
+      )
+    )
+  }
   if (!is.null(res)) {
     assertthat::assert_that(
       assertthat::is.count(res),
@@ -112,6 +124,15 @@ create_spp_data <- function(x,
       msg = "can't use GRASS for processing because it's not available."
     )
   }
+  if (identical(engine, "grass")) {
+    assertthat::assert_that(
+      requireNamespace("sp", quietly = TRUE),
+      msg = paste(
+        "the \"sp\" package needs to be installed, use",
+        "install.packages(\"sp\")"
+      )
+    )
+  }
   if (identical(engine, "grass")) {
     assertthat::assert_that(
       is_gdal_calc_available(),
@@ -347,6 +368,7 @@ create_spp_data <- function(x,
     frc_template_data = frc_template,
     n_threads = n_threads,
     cache_limit = cache_limit,
+    rasterize_touches = rasterize_touches,
     verbose = verbose
   )
 

R/create_spp_frc_data.R

@@ -147,6 +147,7 @@ create_spp_frc_data <- function(x,
                                 n_threads = 1,
                                 cache_limit = 1000,
                                 engine = "terra",
+                                rasterize_touches = FALSE,
                                 verbose = TRUE) {
   create_spp_data(
     x = x,
@@ -162,6 +163,7 @@ create_spp_frc_data <- function(x,
     n_threads = n_threads,
     cache_limit = cache_limit,
     engine = engine,
+    rasterize_touches = rasterize_touches,
     verbose = verbose
   )
 }

R/engine_spp_aoh_gdal.R

@@ -20,7 +20,8 @@ engine_spp_aoh_gdal <- function(range_data,
                                 extent,
                                 path,
                                 n_threads = 1,
-                                cache_limit = 1000) {
+                                cache_limit = 1000,
+                                rasterize_touches = FALSE) {
   # validate arguments
   assertthat::assert_that(
     inherits(range_data, "sf"),
@@ -37,7 +38,9 @@ engine_spp_aoh_gdal <- function(range_data,
     assertthat::is.string(path),
     assertthat::noNA(path),
     assertthat::noNA(n_threads),
-    assertthat::is.number(n_threads)
+    assertthat::is.number(n_threads),
+    assertthat::is.flag(rasterize_touches),
+    assertthat::noNA(rasterize_touches)
   )
 
   # create temporary directory for processing
@@ -115,6 +118,7 @@ engine_spp_aoh_gdal <- function(range_data,
     tiled = FALSE,
     nbits = 1,
     NAflag = 0,
+    touches = rasterize_touches,
     output_raster = FALSE,
     verbose = FALSE
   )

R/engine_spp_aoh_terra.R

@@ -24,6 +24,11 @@ NULL
 #' @param extent `terra::ext()` Object denoting the spatial extent for
 #'  data processing.
 #'
+#' @param touches `logical` Should raster cells that are overlap with any
+#'   part of a species' range be treated as covered?
+#'   Defaults to `FALSE`, such that only cells that have their centroid
+#'   covered by a species' range are treated as covered.
+#'
 #' @param path `character` File path to save resulting Area of Habitat data.
 #'
 #' @return An invisible `TRUE` indicating success.
@@ -36,6 +41,7 @@ engine_spp_aoh_terra <- function(range_data,
                                  lower_elevation,
                                  upper_elevation,
                                  extent,
+                                 rasterize_touches,
                                  path) {
   # validate arguments
   assertthat::assert_that(
@@ -50,6 +56,8 @@ engine_spp_aoh_terra <- function(range_data,
     assertthat::is.number(upper_elevation),
     assertthat::noNA(upper_elevation),
     inherits(extent, "SpatExtent"),
+    assertthat::is.flag(rasterize_touches),
+    assertthat::noNA(rasterize_touches),
     assertthat::is.string(path),
     assertthat::noNA(path)
   )
@@ -96,8 +104,10 @@ engine_spp_aoh_terra <- function(range_data,
   # mask data by species range
   spp_habitat_data <- terra::mask(
     x = spp_habitat_data,
-    mask =  terra_fasterize(
-      sf = range_data, raster = spp_habitat_data
+    mask = terra_fasterize(
+      sf = range_data,
+      raster = spp_habitat_data,
+      touches = rasterize_touches
     ),
     updatevalue = NA_integer_,
     filename = path,

R/misc_terra.R

@@ -69,14 +69,21 @@ terra_st_bbox <- function(x) {
 #' Quickly rasterize vector data
 #'
 #' This function converts a [sf::st_as_sf()] object to a [terra::rast()]
-#' object using [fasterize::fasterize()]. It is similar to
-#' [terra::rasterize()], except that it has greater performance.
+#' object
 #'
 #' @param sf [sf::st_sf()] Object.
 #'
 #' @param raster [terra::rast()] Object.
 #'
-#' @param ... Additional arguments passed to [fasterize::fasterize()].
+#' @param touches `logical` Should cells of `raster` that are overlap with any
+#'   part of `sf` be treated as covered by `sf`?
+#'   Defaults to `FALSE`, such that only cells that have their centroid
+#'   covered by `sf` are treated as covered.
+#'   Defaults to `FALSE`.
+#'
+#' @details
+#' If `touches = FALSE`, then [fasterize::fasterize()] is used to perform
+#' the processing. Otherwise, [terra::rasterize()] is used.
 #'
 #' @return A [terra::rast()] object.
 #'
@@ -94,38 +101,55 @@ terra_st_bbox <- function(x) {
 #' plot(nc_rast)
 #' }
 #' @noRd
-terra_fasterize <- function(sf, raster, ...) {
+terra_fasterize <- function(sf, raster, touches = FALSE) {
   # assert that arguments are valid
   assertthat::assert_that(
     inherits(sf, "sf"),
     inherits(raster, "SpatRaster")
   )
 
-  # convert raster to RasterLayer
-  raster <- withr::with_package(
-    "raster",
-    methods::as(raster[[1]], "Raster"),
-    verbose = FALSE
-  )
+  # run processing based on touches
+  if (isTRUE(touches)) {
+    ## create new column with values for rasterization
+    sf$value <- 1
 
-  # store raster filename
-  raster_filename <- raster::filename(raster)
+    ## rasterize data
+    out <- terra::rasterize(
+      x = terra::vect(sf),
+      y = raster,
+      field = "value",
+      touches = touches
+    )
 
-  # create result
-  out <- fasterize::fasterize(
-    sf = sf,
-    raster = raster,
-    ...
-  )
+  } else {
+    ## convert raster to RasterLayer
+    raster <- withr::with_package(
+      "raster",
+      methods::as(raster[[1]], "Raster"),
+      verbose = FALSE
+    )
+
+    ## store raster filename
+    raster_filename <- raster::filename(raster)
+
+    ## rasterize data
+    out <- fasterize::fasterize(
+      sf = sf,
+      raster = raster
+    )
+
+    ## ensure that result isn't affected by
+    ## https://github.com/ecohealthalliance/fasterize/issues/41
+    if (identical(raster::filename(out), raster_filename)) {
+      out <- raster::brick(out)
+    }
 
-  # ensure that result isn't affected by
-  # https://github.com/ecohealthalliance/fasterize/issues/41
-  if (identical(raster::filename(out), raster_filename)) {
-    out <- raster::brick(out)
+    ## convert result as terra object
+    out <- methods::as(out, "SpatRaster")
   }
 
-  # return result as terra object
-  methods::as(out, "SpatRaster")
+  # return result
+  out
 }
 
 #' On disk?

R/process_spp_aoh_on_local.R

@@ -23,7 +23,8 @@ process_spp_aoh_on_local <- function(x,
                                      path,
                                      engine = "terra",
                                      n_threads = 1,
-                                     cache_limit = 200) {
+                                     cache_limit = 200,
+                                     rasterize_touches = FALSE) {
   # normalize the file path
   path <- normalize_path(path, mustWork = FALSE)
 
@@ -37,6 +38,7 @@ process_spp_aoh_on_local <- function(x,
       lower_elevation = lower_elevation,
       upper_elevation = upper_elevation,
       extent = extent,
+      rasterize_touches = rasterize_touches,
       path = path
     )
   } else if (identical(engine, "gdal")) {
@@ -50,6 +52,7 @@ process_spp_aoh_on_local <- function(x,
       extent = extent,
       path = path,
       n_threads = n_threads,
+      rasterize_touches = rasterize_touches,
       cache_limit = cache_limit
     )
   } else if (identical(engine, "grass")) {