Terra > 1.9.1 enables more efficient aggregation tables

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: mbg
 Title: Model-Based Geostatistics
-Version: 1.1.0
+Version: 1.1.1
 Authors@R: c(
     person("Nathaniel", "Henry", email = "nat@henryspatialanalysis.com",
         role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8150-4988")),
@@ -29,7 +29,7 @@ Imports:
     purrr,
     R6,
     sf,
-    terra,
+    terra (>= 1.9.1),
     tictoc
 Suggests:
     INLA,
@@ -40,7 +40,7 @@ Suggests:
 Additional_repositories: https://inla.r-inla-download.org/R/stable/
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE, r6 = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.3
 VignetteBuilder: knitr
 URL: https://henryspatialanalysis.github.io/mbg/
 BugReports: https://github.com/henryspatialanalysis/mbg/issues

NAMESPACE

@@ -5,7 +5,6 @@ export(aggregate_draws_to_polygons)
 export(aggregate_raster_to_polygons)
 export(build_aggregation_table)
 export(build_id_raster)
-export(calculate_pixel_fractions_single_polygon)
 export(dissolve_sf_by_attribute)
 export(fit_inla_model)
 export(generate_cell_draws_and_summarize)
@@ -28,12 +27,10 @@ importFrom(Matrix,rowMeans)
 importFrom(R6,R6Class)
 importFrom(assertthat,assert_that)
 importFrom(assertthat,has_name)
-importFrom(assertthat,is.scalar)
 importFrom(assertthat,noNA)
 importFrom(caret,train)
 importFrom(caret,trainControl)
 importFrom(data.table,as.data.table)
-importFrom(data.table,data.table)
 importFrom(glue,glue)
 importFrom(matrixStats,rowQuantiles)
 importFrom(purrr,map)

R/build_aggregation_table.R

@@ -1,73 +1,3 @@
-#' Calculate fractional pixels in a polygon
-#'
-#' @description Calculate the fraction of each pixel's area that falls within a single
-#'   polygon
-#'
-#' @details This is a helper function called by `build_aggregation_table()`.
-#'
-#' @param polygon [terra::SpatVector] object of length 1. The polygon to calculate
-#'   fractional areas across.
-#' @param id_raster [terra::SpatRaster] object. ID raster created for the set of all
-#'   polygons to be considered, created by `build_id_raster()`.
-#' @param polygon_id (optional). ID for this polygon. Must have length 1.
-#'
-#' @return data.table containing two or three columns:
-#'   * pixel_id: unique pixel ID from the ID raster
-#'   * area_fraction: fraction of the pixel area falling within this polygon
-#'   * polygon_id (optional): If `polygon_id` was defined, it is added to the table
-#'
-#' @examples
-#' \dontrun{
-#'   polygons <- sf::st_read(system.file('extdata/Benin_communes.gpkg', package = 'mbg'))
-#'   id_raster <- build_id_raster(polygons)
-#'   pixel_fractions <- calculate_pixel_fractions_single_polygon(
-#'     polygon = polygons[1, ], id_raster
-#'   )
-#'   head(pixel_fractions)
-#' }
-#'
-#' @concept aggregation
-#'
-#' @seealso build_aggregation_table
-#'
-#' @importFrom assertthat assert_that is.scalar
-#' @importFrom data.table data.table
-#' @importFrom terra vect same.crs extract
-#' @export
-calculate_pixel_fractions_single_polygon <- function(
-  polygon, id_raster, polygon_id = NULL
-){
-  # If polygon inherits sf, convert to SpatVector
-  if(inherits(polygon, 'sf')) polygon <- terra::vect(polygon)
-  # Input data checks
-  assertthat::assert_that(inherits(polygon, 'SpatVector'))
-  assertthat::assert_that(inherits(id_raster, 'SpatRaster'))
-  assertthat::assert_that(terra::same.crs(id_raster, polygon))
-  assertthat::assert_that(nrow(polygon) == 1)
-  if(!is.null(polygon_id)){
-    assertthat::assert_that(assertthat::is.scalar(polygon_id))
-  }
-
-  # Use terra::extract to get pixel area fractions falling within the polygon
-  # Drop the first column, which gives the polygon row ID
-  extract_wide <- terra::extract(
-    x = id_raster,
-    y = polygon,
-    fun = 'table',
-    exact = TRUE,
-    ID = FALSE
-  )
-  # Reshape long
-  extract_long <- data.table::data.table(
-    pixel_id = colnames(extract_wide) |> as.integer(),
-    area_fraction = unlist(extract_wide[1, ])
-  )
-  # Optionally add the polygon ID to the table
-  if(!is.null(polygon_id)) extract_long$polygon_id <- polygon_id
-
-  return(extract_long)
-}
-
 
 #' Validation: Build aggregation table
 #'
@@ -150,7 +80,7 @@ build_aggregation_table <- function(
   polygons, id_raster, polygon_id_field, verbose = FALSE
 ){
   # Overload some data.table variables to pass R CMD check
-  . <- pixel_id <- masked_pixel_id <- i.masked_pixel_id <- polygon_id  <-
+  . <- dummy_id <- pixel_id <- masked_pixel_id <- i.masked_pixel_id <- polygon_id  <-
     area_fraction <- NULL
 
   # If polygons inherit sf, convert to SpatVector
@@ -167,28 +97,25 @@ build_aggregation_table <- function(
 
   # Crop the polygons to the id raster
   polygons_cropped <- terra::crop(x = polygons, y = id_raster, ext = TRUE)
-  poly_ids <- polys_dt[[polygon_id_field]]
+  if(verbose){
+    dropped_rows <- nrow(polygons) - nrow(polygons_cropped)
+    if(dropped_rows > 0L) message(paste(
+      "Dropped", dropped_rows, "polygons not in the id_raster extent."
+    ))
+  }
+  poly_ids <- as.data.frame(polygons_cropped)[[polygon_id_field]]
 
-  # Build the aggregation table by calling zonal statistics for each polygon
-  agg_table <- lapply(poly_ids, function(poly_id){
-    # Create smaller spatial objects for calculating fractional zonal statistics
-    if(verbose) message('.', appendLF = FALSE)
-    one_poly <- polygons_cropped[poly_ids == poly_id, ]
-    id_raster_sub <- terra::crop(x = id_raster, y = one_poly, ext = TRUE, snap = 'out')
-    # Mask missing values with -1
-    missing_cells <- which(is.na(terra::values(id_raster_sub, mat = F)))
-    terra::values(id_raster_sub)[missing_cells] <- -1
-    # Get fractional pixel areas for the polygon
-    pixel_fractions <- calculate_pixel_fractions_single_polygon(
-      polygon = one_poly,
-      id_raster = id_raster_sub,
-      polygon_id = poly_id
-    )
-    # Drop -1 (masked) pixel IDs and return
-    return(pixel_fractions[pixel_id >= 0, ])
-  }) |> data.table::rbindlist()
-  # Add a newline to finish off the progress bar, if needed
-  if(verbose) message("")
+  agg_table <- terra::extract(
+    x = id_raster,
+    y = polygons_cropped,
+    fun = 'table',
+    exact = TRUE,
+    ID = TRUE
+  ) |>
+    data.table::as.data.table() |>
+    data.table::setnames(c('dummy_id', 'pixel_id', 'area_fraction')) |>
+    _[, polygon_id := sapply(dummy_id, function(x) poly_ids[x])] |>
+    _[, dummy_id := NULL]
 
   # Merge on 'masked_pixel_id'
   masked_pixel_table <- data.table::data.table(