From 8f296a19494e15ccc448c6e6bd87b758b8bbd146 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Sat, 7 Mar 2026 12:26:28 -0800 Subject: [PATCH] Add dft_rast_consensus() for temporal mode smoothing Per-pixel mode across multiple classified rasters to filter single-year misclassification noise. Auto-aligns rasters with different extents via nearest-neighbour resampling. Optional confidence layer (proportion of inputs agreeing on mode). Fixes #8 Co-Authored-By: Claude Opus 4.6 --- R/dft_rast_consensus.R | 101 ++++++++++++++++++++ tests/testthat/test-dft_rast_consensus.R | 114 +++++++++++++++++++++++ 2 files changed, 215 insertions(+) create mode 100644 R/dft_rast_consensus.R create mode 100644 tests/testthat/test-dft_rast_consensus.R diff --git a/R/dft_rast_consensus.R b/R/dft_rast_consensus.R new file mode 100644 index 0000000..788df8b --- /dev/null +++ b/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 +} diff --git a/tests/testthat/test-dft_rast_consensus.R b/tests/testthat/test-dft_rast_consensus.R new file mode 100644 index 0000000..64bbc19 --- /dev/null +++ b/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) +})