Add dft_rast_consensus() for temporal mode smoothing

R/dft_rast_consensus.R

@@ -0,0 +1,101 @@
+#' Compute per-pixel mode across classified rasters
+#'
+#' Given multiple classified rasters of the same extent and resolution,
+#' compute the most frequent (mode) class at each pixel. Useful for
+#' temporal smoothing — averaging out single-year misclassification noise
+#' before running [dft_rast_transition()].
+#'
+#' @details
+#' Mode smoothing filters single-year misclassification but cannot
+#' distinguish noise from real change. A pixel that genuinely transitions
+#' mid-window may be voted back to its original class if the pre-change
+#' years outnumber the post-change years. See
+#' \href{https://github.com/NewGraphEnvironment/drift/issues/9}{drift#9}
+#' for discussion of weighted and breakpoint approaches.
+#'
+#' @param x A named list of classified `SpatRaster`s (e.g. from
+#'   [dft_rast_classify()]). Rasters with slightly different extents are
+#'   automatically resampled (nearest-neighbour) to the first raster's grid.
+#' @param confidence Logical. If `TRUE`, return a second layer with the
+#'   proportion of input rasters that agreed on the mode (e.g. 3/3 = 1.0,
+#'   2/3 = 0.67). Default `FALSE`.
+#'
+#' @return A `SpatRaster`. When `confidence = FALSE`, a single-layer factor
+#'   raster with the modal class. When `confidence = TRUE`, a two-layer
+#'   raster: `"consensus"` (factor) and `"confidence"` (numeric 0–1).
+#' @export
+#' @examples
+#' # Build 3 classified rasters from example data
+#' r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+#' r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+#' r23 <- terra::rast(system.file("extdata", "example_2023.tif", package = "drift"))
+#' classified <- dft_rast_classify(
+#'   list("2017" = r17, "2020" = r20, "2023" = r23), source = "io-lulc"
+#' )
+#'
+#' # Consensus raster (modal class)
+#' cons <- dft_rast_consensus(classified)
+#' terra::plot(cons)
+#'
+#' # With confidence layer
+#' cons2 <- dft_rast_consensus(classified, confidence = TRUE)
+#' terra::plot(cons2[["confidence"]])
+dft_rast_consensus <- function(x, confidence = FALSE) {
+  if (!is.list(x) || inherits(x, "SpatRaster")) {
+    stop("`x` must be a named list of SpatRasters, not a single SpatRaster.")
+  }
+  if (length(x) < 2) {
+    stop("`x` must contain at least 2 rasters.")
+  }
+
+  # Align all rasters to the first raster's grid
+  ref <- x[[1]]
+  x <- lapply(x, function(r) {
+    if (!terra::compareGeom(ref, r, stopOnError = FALSE)) {
+      terra::resample(r, ref, method = "near")
+    } else {
+      r
+    }
+  })
+
+  # Stack into matrix: rows = pixels, cols = rasters
+  vals <- do.call(cbind, lapply(x, function(r) terra::values(r)[, 1]))
+  n <- ncol(vals)
+
+  # Per-pixel mode and count
+  mode_result <- apply(vals, 1, function(row) {
+    row <- row[!is.na(row)]
+    if (length(row) == 0) return(c(NA_integer_, NA_real_))
+    tab <- tabulate(match(row, unique(row)))
+    uq <- unique(row)
+    idx <- which.max(tab)
+    c(uq[idx], tab[idx] / length(row))
+  })
+
+  mode_vals <- as.integer(mode_result[1, ])
+  conf_vals <- mode_result[2, ]
+
+  # Build output raster with factor levels from first input
+  r_out <- terra::rast(x[[1]])
+  terra::values(r_out) <- mode_vals
+
+  # Copy factor levels if present
+  if (terra::is.factor(x[[1]])) {
+    lvls <- terra::cats(x[[1]])[[1]]
+    # Only keep levels that appear in the consensus
+    present <- unique(mode_vals[!is.na(mode_vals)])
+    lvls_present <- lvls[lvls[[1]] %in% present, , drop = FALSE]
+    terra::set.cats(r_out, layer = 1, value = lvls_present)
+  }
+
+  names(r_out) <- "consensus"
+
+  if (confidence) {
+    r_conf <- terra::rast(x[[1]])
+    terra::values(r_conf) <- conf_vals
+    names(r_conf) <- "confidence"
+    r_out <- c(r_out, r_conf)
+  }
+
+  r_out
+}

tests/testthat/test-dft_rast_consensus.R

@@ -0,0 +1,114 @@
+test_that("dft_rast_consensus returns factor raster with modal class", {
+  r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+  r23 <- terra::rast(system.file("extdata", "example_2023.tif", package = "drift"))
+  classified <- dft_rast_classify(
+    list("2017" = r17, "2020" = r20, "2023" = r23), source = "io-lulc"
+  )
+
+  cons <- dft_rast_consensus(classified)
+
+  expect_s4_class(cons, "SpatRaster")
+  expect_equal(terra::nlyr(cons), 1)
+  expect_true(terra::is.factor(cons))
+  expect_equal(names(cons), "consensus")
+})
+
+test_that("consensus raster has same dimensions as input", {
+  r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+  classified <- dft_rast_classify(
+    list("2017" = r17, "2020" = r20), source = "io-lulc"
+  )
+
+  cons <- dft_rast_consensus(classified)
+
+  expect_equal(terra::nrow(cons), terra::nrow(classified[["2017"]]))
+  expect_equal(terra::ncol(cons), terra::ncol(classified[["2017"]]))
+})
+
+test_that("confidence layer returned when requested", {
+  r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+  r23 <- terra::rast(system.file("extdata", "example_2023.tif", package = "drift"))
+  classified <- dft_rast_classify(
+    list("2017" = r17, "2020" = r20, "2023" = r23), source = "io-lulc"
+  )
+
+  cons <- dft_rast_consensus(classified, confidence = TRUE)
+
+  expect_equal(terra::nlyr(cons), 2)
+  expect_equal(names(cons), c("consensus", "confidence"))
+})
+
+test_that("confidence values are between 0 and 1", {
+  r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+  r23 <- terra::rast(system.file("extdata", "example_2023.tif", package = "drift"))
+  classified <- dft_rast_classify(
+    list("2017" = r17, "2020" = r20, "2023" = r23), source = "io-lulc"
+  )
+
+  cons <- dft_rast_consensus(classified, confidence = TRUE)
+  conf_vals <- terra::values(cons[["confidence"]])[, 1]
+  conf_vals <- conf_vals[!is.na(conf_vals)]
+
+  expect_true(all(conf_vals >= 0 & conf_vals <= 1))
+})
+
+test_that("unanimous pixels have confidence = 1", {
+  # Create 3 identical rasters — every pixel should have confidence 1.0
+  r <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  classified <- dft_rast_classify(list("a" = r, "b" = r, "c" = r), source = "io-lulc")
+
+  cons <- dft_rast_consensus(classified, confidence = TRUE)
+  conf_vals <- terra::values(cons[["confidence"]])[, 1]
+  conf_vals <- conf_vals[!is.na(conf_vals)]
+
+  expect_true(all(conf_vals == 1))
+})
+
+test_that("consensus of identical rasters equals the input", {
+  r <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  classified <- dft_rast_classify(list("a" = r, "b" = r), source = "io-lulc")
+
+  cons <- dft_rast_consensus(classified)
+  # Raw values should match
+  expect_equal(
+    terra::values(cons)[, 1],
+    terra::values(classified[["a"]])[, 1]
+  )
+})
+
+test_that("errors on single SpatRaster", {
+  r <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  expect_error(dft_rast_consensus(r), "named list")
+})
+
+test_that("errors on list with fewer than 2 rasters", {
+  r <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  classified <- dft_rast_classify(list("a" = r), source = "io-lulc")
+  expect_error(dft_rast_consensus(classified), "at least 2")
+})
+
+test_that("consensus works with dft_rast_transition", {
+  r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift"))
+  r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift"))
+  r23 <- terra::rast(system.file("extdata", "example_2023.tif", package = "drift"))
+  classified <- dft_rast_classify(
+    list("2017" = r17, "2020" = r20, "2023" = r23), source = "io-lulc"
+  )
+
+  # Use 2017+2020 as "early", 2023 alone as "late" (just to test integration)
+  early <- dft_rast_consensus(classified[c("2017", "2020")])
+  late <- classified[["2023"]]
+
+  result <- dft_rast_transition(
+    list(early = early, late = late),
+    from = "early", to = "late"
+  )
+
+  expect_type(result, "list")
+  expect_named(result, c("raster", "summary"))
+  expect_true(nrow(result$summary) > 0)
+})