diff --git a/NAMESPACE b/NAMESPACE index 7171146..7a39c4a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,3 +12,4 @@ export(dft_rast_transition) export(dft_stac_classes) export(dft_stac_config) export(dft_stac_fetch) +export(dft_transition_vectors) diff --git a/R/dft_transition_vectors.R b/R/dft_transition_vectors.R new file mode 100644 index 0000000..108a7c7 --- /dev/null +++ b/R/dft_transition_vectors.R @@ -0,0 +1,144 @@ +#' Vectorize transition raster into individual change patches +#' +#' Convert a transition `SpatRaster` (from [dft_rast_transition()]) into `sf` +#' polygons — one row per connected patch of pixels sharing the same +#' transition type. Useful for QA in GIS, spatial attribution to management +#' zones, and patch-level reporting. +#' +#' @param x A factor `SpatRaster` from [dft_rast_transition()] (the `$raster` +#' element). Must have a projected CRS. +#' @param zones Optional `sf` polygon layer for spatial attribution. Any +#' partitioning: sub-basins, parcels, climate regions, management units. +#' @param zone_col Character. Column name in `zones` identifying each zone. +#' Required when `zones` is supplied. +#' @param patch_area_min Numeric or `NULL`. Minimum patch area in m². Patches +#' smaller than this are dropped before returning. `NULL` (default) keeps all. +#' +#' @return An `sf` data frame (polygon geometry) with columns: +#' - `patch_id` (integer) — connected component ID +#' - `transition` (character) — transition label (e.g. "Trees -> Rangeland") +#' - `area_ha` (numeric) — patch area in hectares +#' - Zone column (if `zones` supplied) — from spatial intersection +#' +#' @export +#' @examples +#' 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") +#' result <- dft_rast_transition(classified, from = "2017", to = "2020") +#' +#' # Vectorize all transition patches +#' patches <- dft_transition_vectors(result$raster) +#' head(patches) +#' +#' # Filter to large patches only +#' patches_large <- dft_transition_vectors(result$raster, patch_area_min = 1000) +#' head(patches_large) +dft_transition_vectors <- function(x, + zones = NULL, + zone_col = NULL, + patch_area_min = NULL) { + if (!inherits(x, "SpatRaster")) { + stop("`x` must be a SpatRaster (e.g. from dft_rast_transition()$raster).", + call. = FALSE) + } + dft_check_crs(x, "dft_transition_vectors") + + if (!terra::is.factor(x)) { + stop("`x` must be a factor SpatRaster with transition labels.", call. = FALSE) + } + + if (!is.null(zones)) { + if (!inherits(zones, "sf")) { + stop("`zones` must be an sf object.", call. = FALSE) + } + if (is.null(zone_col) || !zone_col %in% names(zones)) { + stop("`zone_col` must name a column in `zones`.", call. = FALSE) + } + } + + if (!is.null(patch_area_min)) { + if (!is.numeric(patch_area_min) || length(patch_area_min) != 1 || + is.na(patch_area_min) || patch_area_min < 0) { + stop("`patch_area_min` must be a single non-negative number or NULL.", + call. = FALSE) + } + } + + lvls <- terra::cats(x)[[1]] + vals <- terra::values(x)[, 1] + + if (all(is.na(vals))) { + return(sf::st_sf( + patch_id = integer(0), transition = character(0), + area_ha = numeric(0), geometry = sf::st_sfc(crs = sf::st_crs(x)) + )) + } + + # Build unique patch IDs per transition type + patch_ids <- rep(NA_integer_, terra::ncell(x)) + patch_transition <- integer(0) + offset <- 0L + + for (i in seq_len(nrow(lvls))) { + code <- lvls$id[i] + mask <- which(vals == code) + if (length(mask) == 0) next + + r_mask <- terra::rast(x) + mask_vals <- rep(NA_integer_, terra::ncell(x)) + mask_vals[mask] <- 1L + terra::values(r_mask) <- mask_vals + p <- terra::patches(r_mask, directions = 8) + p_vals <- terra::values(p)[, 1] + + valid <- !is.na(p_vals) + unique_pids <- sort(unique(p_vals[valid])) + for (pid in unique_pids) { + offset <- offset + 1L + patch_ids[valid & p_vals == pid] <- offset + patch_transition <- c(patch_transition, code) + } + } + + if (offset == 0L) { + return(sf::st_sf( + patch_id = integer(0), transition = character(0), + area_ha = numeric(0), geometry = sf::st_sfc(crs = sf::st_crs(x)) + )) + } + + # Vectorize patch IDs + r_pid <- terra::rast(x) + names(r_pid) <- "pid" + terra::values(r_pid) <- patch_ids + polys <- terra::as.polygons(r_pid) + polys_sf <- sf::st_as_sf(polys) + + # Map patch IDs to transition labels + code_to_label <- stats::setNames(lvls$transition, lvls$id) + polys_sf$patch_id <- as.integer(polys_sf$pid) + polys_sf$transition <- as.character( + code_to_label[as.character(patch_transition[polys_sf$pid])] + ) + polys_sf$area_ha <- as.numeric(sf::st_area(polys_sf)) * 1e-4 + polys_sf$pid <- NULL + + # Filter by minimum area + if (!is.null(patch_area_min) && patch_area_min > 0) { + area_m2 <- as.numeric(sf::st_area(polys_sf)) + polys_sf <- polys_sf[area_m2 >= patch_area_min, ] + } + + # Select and order columns + out <- polys_sf[c("patch_id", "transition", "area_ha")] + + # Zone attribution + + if (!is.null(zones)) { + zones_proj <- sf::st_transform(zones[zone_col], sf::st_crs(out)) + out <- suppressWarnings(sf::st_intersection(out, zones_proj)) + } + + out +} diff --git a/man/dft_transition_vectors.Rd b/man/dft_transition_vectors.Rd new file mode 100644 index 0000000..95d568b --- /dev/null +++ b/man/dft_transition_vectors.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_transition_vectors.R +\name{dft_transition_vectors} +\alias{dft_transition_vectors} +\title{Vectorize transition raster into individual change patches} +\usage{ +dft_transition_vectors(x, zones = NULL, zone_col = NULL, patch_area_min = NULL) +} +\arguments{ +\item{x}{A factor \code{SpatRaster} from \code{\link[=dft_rast_transition]{dft_rast_transition()}} (the \verb{$raster} +element). Must have a projected CRS.} + +\item{zones}{Optional \code{sf} polygon layer for spatial attribution. Any +partitioning: sub-basins, parcels, climate regions, management units.} + +\item{zone_col}{Character. Column name in \code{zones} identifying each zone. +Required when \code{zones} is supplied.} + +\item{patch_area_min}{Numeric or \code{NULL}. Minimum patch area in m². Patches +smaller than this are dropped before returning. \code{NULL} (default) keeps all.} +} +\value{ +An \code{sf} data frame (polygon geometry) with columns: +\itemize{ +\item \code{patch_id} (integer) — connected component ID +\item \code{transition} (character) — transition label (e.g. "Trees -> Rangeland") +\item \code{area_ha} (numeric) — patch area in hectares +\item Zone column (if \code{zones} supplied) — from spatial intersection +} +} +\description{ +Convert a transition \code{SpatRaster} (from \code{\link[=dft_rast_transition]{dft_rast_transition()}}) into \code{sf} +polygons — one row per connected patch of pixels sharing the same +transition type. Useful for QA in GIS, spatial attribution to management +zones, and patch-level reporting. +} +\examples{ +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") +result <- dft_rast_transition(classified, from = "2017", to = "2020") + +# Vectorize all transition patches +patches <- dft_transition_vectors(result$raster) +head(patches) + +# Filter to large patches only +patches_large <- dft_transition_vectors(result$raster, patch_area_min = 1000) +head(patches_large) +} diff --git a/tests/testthat/test-dft_transition_vectors.R b/tests/testthat/test-dft_transition_vectors.R new file mode 100644 index 0000000..6bcde52 --- /dev/null +++ b/tests/testthat/test-dft_transition_vectors.R @@ -0,0 +1,109 @@ +test_that("dft_transition_vectors returns sf with expected columns", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + + patches <- dft_transition_vectors(result$raster) + + expect_s3_class(patches, "sf") + expect_true(all(c("patch_id", "transition", "area_ha") %in% names(patches))) + expect_true(nrow(patches) > 0) +}) + +test_that("total area matches raster summary", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + + patches <- dft_transition_vectors(result$raster) + + expect_equal(sum(patches$area_ha), sum(result$summary$area), tolerance = 0.1) +}) + +test_that("each patch has a valid transition label", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + + patches <- dft_transition_vectors(result$raster) + lvls <- terra::cats(result$raster)[[1]] + + expect_true(all(patches$transition %in% lvls$transition)) +}) + +test_that("patch_ids are unique", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + + patches <- dft_transition_vectors(result$raster) + + expect_equal(length(unique(patches$patch_id)), nrow(patches)) +}) + +test_that("patch_area_min filters small patches", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + + all_patches <- dft_transition_vectors(result$raster) + filtered <- dft_transition_vectors(result$raster, patch_area_min = 1000) + + expect_lt(nrow(filtered), nrow(all_patches)) + expect_true(all(as.numeric(sf::st_area(filtered)) >= 1000)) +}) + +test_that("zone attribution adds zone column", { + r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift")) + r20 <- terra::rast(system.file("extdata", "example_2020.tif", package = "drift")) + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), quiet = TRUE + ) + classified <- dft_rast_classify(list("2017" = r17, "2020" = r20), source = "io-lulc") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + + # Use AOI as a single zone + aoi$zone_name <- "test_zone" + patches <- dft_transition_vectors(result$raster, zones = aoi, + zone_col = "zone_name") + + expect_true("zone_name" %in% names(patches)) +}) + +test_that("errors on non-SpatRaster input", { + expect_error(dft_transition_vectors(data.frame(x = 1)), + "SpatRaster") +}) + +test_that("errors on non-factor SpatRaster", { + r <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift")) + expect_error(dft_transition_vectors(r), "factor") +}) + +test_that("errors on geographic CRS", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + r_geo <- terra::project(result$raster, "EPSG:4326", method = "near") + + expect_error(dft_transition_vectors(r_geo), "projected CRS") +}) + +test_that("errors when zones supplied without zone_col", { + 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") + result <- dft_rast_transition(classified, from = "2017", to = "2020") + aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), quiet = TRUE + ) + + expect_error(dft_transition_vectors(result$raster, zones = aoi), + "zone_col") +})