diff --git a/DESCRIPTION b/DESCRIPTION index eb2063a..b6fed45 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: drift Title: Detecting Riparian and Inland Floodplain Transitions -Version: 0.7.0 -Date: 2026-07-11 +Version: 0.8.0 +Date: 2026-07-17 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), diff --git a/NAMESPACE b/NAMESPACE index a9fa599..2c012de 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,4 +17,5 @@ export(dft_stac_classes) export(dft_stac_config) export(dft_stac_cube) export(dft_stac_fetch) +export(dft_transition_attribute) export(dft_transition_vectors) diff --git a/NEWS.md b/NEWS.md index a39c608..9350dc5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# drift 0.8.0 + +- `dft_transition_attribute()` tags change patches from `dft_transition_vectors()` with columns from any overlay polygon layer — fire perimeters, cutblocks, roads, tenures — so a driver can separate mapped transitions by cause without hand-rolling spatial joins (#42). Deliberately generic: drift carries no BC/domain knowledge; the caller supplies the overlay, the columns to carry (`cols`), and optionally a numeric temporal filter (`time_col` + `time_interval`, bounds inclusive) that keeps only overlay features whose time falls within the transition interval — a 2022 fire attributes a 2017→2023 loss, a 2012 fire does not. Two assignment modes for a patch that straddles multiple overlay features: `match_mode = "all"` (left join, one row per match) or `"largest"` (one row per patch by greatest overlap area). Largest-overlap assignment is intersection-based, so combining it with a custom `predicate` is an error rather than a silent mis-attribution; the overlay is reprojected to the patch CRS and run through `st_make_valid()` automatically, since real-world disturbance perimeters routinely fail validity checks. + # drift 0.7.0 - `dft_stac_cube()` gains `tile_size` (default `NULL`), the continuous-path twin of `dft_stac_fetch()`'s `tile_size` (#36): an opt-in that bounds the STAC *read* to the AOI footprint (#38). By default one gdalcubes cube is streamed over the whole AOI bounding box, so the COG streaming — the dominant cost (~10-30 min for a multi-year monthly Sentinel-2 fetch) — scales with the bbox, not the AOI; for a thin, diagonal floodplain corridor the bbox is largely empty. When `tile_size` (CRS units — metres for the default UTM CRS) is set, the bbox is split into a `res`-aligned grid and only tiles that intersect the AOI polygon are streamed — each carrying the full SCL mask, spectral index, and 2022 baseline-offset split — then mosaicked with `terra::merge()`, so a corridor reads close to its footprint. This is the `filter_geom`-independent path (the polygon clip that would do this in-cube segfaults on the pinned gdalcubes build). The cube always caches a `.tif` and a tiled read keys distinctly, so untiled caches are untouched and `tile_size = NULL` is byte-for-byte the previous behavior. Because the cube resamples with bilinear, a tiled cube faithfully reproduces the untiled cube (bilinear-aligned correlation ~0.997, per-layer means within ~1e-3, no tile seams) but lands on a bbox-anchored grid that is sub-pixel-offset from — not pixel-identical to — the untiled cube; the offset is immaterial to the per-pixel `dft_rast_break()` / `dft_rast_trend()` reducers. With `tile_size` set, `clip = FALSE` returns the AOI-intersecting tile union (with `NA` where empty tiles were skipped), not a gap-free bounding box. diff --git a/R/dft_transition_attribute.R b/R/dft_transition_attribute.R new file mode 100644 index 0000000..765f591 --- /dev/null +++ b/R/dft_transition_attribute.R @@ -0,0 +1,242 @@ +#' Attribute transition patches from an overlay polygon layer +#' +#' Tag change patches (from [dft_transition_vectors()]) with columns from any +#' overlay polygon layer — fire perimeters, cutblocks, roads, tenures — to help +#' separate mapped transitions by cause. Generic by design: drift carries no +#' domain knowledge; the caller supplies the overlay, the columns to carry, and +#' (optionally) the temporal filter. +#' +#' The overlay is transformed to the CRS of `patches` before joining, and run +#' through [sf::st_make_valid()] first — real-world disturbance perimeters +#' routinely fail GEOS validity checks, which would otherwise break the +#' largest-overlap computation. +#' +#' @param patches An `sf` object of change patches, typically from +#' [dft_transition_vectors()]. +#' @param overlay An `sf` polygon layer to attribute from (e.g. fire +#' perimeters, consolidated cutblocks). +#' @param cols Character vector of `overlay` column names to carry onto each +#' patch. Must not collide with existing `patches` column names. +#' @param predicate Spatial predicate function used to match patches to +#' overlay features, e.g. [sf::st_intersects()] (default) or +#' [sf::st_within()]. Only applies with `match_mode = "all"`: +#' largest-overlap assignment is inherently intersection-based +#' (`sf::st_join(largest = TRUE)` ignores the join predicate), so combining +#' a custom predicate with `match_mode = "largest"` is an error. +#' @param match_mode How a patch is assigned when it matches more than one +#' overlay feature: +#' - `"all"` (default) — plain left join; a patch straddling k overlay +#' features appears k times (`patch_id` repeats). +#' - `"largest"` — exactly one row per patch, assigned to the overlay +#' feature with the greatest intersection area (via +#' `sf::st_join(largest = TRUE)`; matching is by intersection, see +#' `predicate`). +#' @param time_col Character or `NULL`. Name of a **numeric** time column in +#' `overlay` used for temporal filtering (e.g. `FIRE_YEAR`, `HARVEST_YEAR`). +#' Must be supplied together with `time_interval`, and both must be on the +#' **same numeric scale** (see Details). `NULL` (default) skips the filter. +#' @param time_interval Length-2 numeric or `NULL`. The transition interval on +#' the same scale as `time_col`, e.g. `c(2017, 2023)` for calendar years. +#' Overlay features whose `time_col` falls outside the interval (both bounds +#' inclusive) are dropped before joining; features with an `NA` time are also +#' dropped. Patches from [dft_transition_vectors()] carry no time columns, so +#' the interval is always supplied explicitly. +#' +#' @details +#' # Temporal filter — how `time_col` and `time_interval` must be presented +#' +#' The filter is a plain numeric comparison +#' (`overlay[[time_col]] >= time_interval[1] & <= time_interval[2]`), so it is +#' scale-agnostic — the numbers may be calendar years, decimal years, months, +#' or epoch offsets — but **both arguments must be numeric and on the same +#' scale**. `time_col` must name a `numeric` (integer or double) column; +#' passing a `Date` or `POSIXct` column is a hard error, not a silent +#' mis-comparison. Values are not coerced or rounded: a decimal year like +#' `2018.5` is compared as-is. +#' +#' To filter on dates, convert to a numeric axis first, on **both** the column +#' and the interval: +#' +#' - Calendar year (simplest for annual disturbance data): +#' ``` +#' overlay$yr <- as.numeric(format(overlay$burn_date, "%Y")) +#' dft_transition_attribute(..., time_col = "yr", time_interval = c(2017, 2023)) +#' ``` +#' - Epoch days (`Date` stores days since 1970-01-01, so `as.numeric()` is the +#' coercion): +#' ``` +#' overlay$t <- as.numeric(overlay$burn_date) # days since epoch +#' dft_transition_attribute(..., time_col = "t", +#' time_interval = as.numeric(as.Date(c("2017-01-01", "2023-12-31")))) +#' ``` +#' (`POSIXct` coerces to *seconds* since epoch — keep the interval in seconds +#' to match.) +#' +#' Mixing scales (e.g. epoch-day column against a `c(2017, 2023)` year +#' interval) does not error — it silently matches nothing. Keep both on one +#' axis. +#' +#' @section Limitations: +#' - `overlay` is treated as a polygon layer. Passing point or line geometries +#' is not supported: `match_mode = "largest"` compares intersection *area*, +#' which is zero for non-polygon overlays. Use a point-in-polygon join +#' directly for such cases. +#' - Temporal filtering is numeric-only; `Date`/`POSIXct` columns must be +#' coerced by the caller (see Details). +#' +#' @return `patches` with the `cols` columns joined on (`NA` where a patch +#' matches no overlay feature). Under `match_mode = "all"` a patch matching +#' several overlay features is duplicated, one row per match; under +#' `"largest"` the result has exactly one row per input patch. +#' +#' @seealso [dft_transition_vectors()] for producing the input patches. +#' +#' @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") +#' patches <- dft_transition_vectors(result$raster, changes_only = TRUE) +#' +#' # synthetic disturbance overlay covering the western half of the AOI +#' bb <- sf::st_bbox(patches) +#' west <- sf::st_sf( +#' fire_year = 2018, +#' geometry = sf::st_as_sfc( +#' sf::st_bbox(c(bb["xmin"], bb["ymin"], +#' xmax = unname((bb["xmin"] + bb["xmax"]) / 2), bb["ymax"]), +#' crs = sf::st_crs(patches)) +#' ) +#' ) +#' +#' # tag each patch with the fire year where it overlaps (NA elsewhere) +#' tagged <- dft_transition_attribute(patches, west, cols = "fire_year", +#' match_mode = "largest") +#' table(tagged$fire_year, useNA = "ifany") +#' +#' # temporal filter: only overlay features within the transition interval +#' tagged_2017_2020 <- dft_transition_attribute( +#' patches, west, cols = "fire_year", match_mode = "largest", +#' time_col = "fire_year", time_interval = c(2017, 2020) +#' ) +dft_transition_attribute <- function(patches, + overlay, + cols, + predicate = sf::st_intersects, + match_mode = c("all", "largest"), + time_col = NULL, + time_interval = NULL) { + match_mode <- match.arg(match_mode) + + if (!inherits(patches, "sf")) { + cli::cli_abort(c( + "{.arg patches} must be an {.cls sf} object.", + "i" = "Use {.fn dft_transition_vectors} to create change patches." + )) + } + if (!inherits(overlay, "sf")) { + cli::cli_abort("{.arg overlay} must be an {.cls sf} object.") + } + if (!is.function(predicate)) { + cli::cli_abort( + "{.arg predicate} must be a function, e.g. {.fn sf::st_intersects}." + ) + } + if (identical(match_mode, "largest") && + !identical(predicate, sf::st_intersects)) { + cli::cli_abort(c( + "{.arg predicate} cannot be combined with {.code match_mode = \"largest\"}.", + "i" = "Largest-overlap assignment is intersection-based: + {.code sf::st_join(largest = TRUE)} ignores the join predicate." + )) + } + if (!is.character(cols) || length(cols) == 0 || anyNA(cols)) { + cli::cli_abort( + "{.arg cols} must be a non-empty character vector of {.arg overlay} column names." + ) + } + cols_missing <- setdiff(cols, names(overlay)) + if (length(cols_missing) > 0) { + cli::cli_abort( + "{.arg cols} not found in {.arg overlay}: {.val {cols_missing}}." + ) + } + # Only the active geometry column is guarded here; a secondary sfc column + # named in `cols` would slip through and leak a second geometry into the + # result. Left unguarded as a low-risk edge (overlays rarely carry two sfc + # columns); guard on names(overlay) types if that ever bites. + if (attr(overlay, "sf_column") %in% cols) { + cli::cli_abort( + "{.arg cols} must not include the {.arg overlay} geometry column + {.val {attr(overlay, 'sf_column')}}." + ) + } + cols_clash <- intersect(cols, names(patches)) + if (length(cols_clash) > 0) { + cli::cli_abort(c( + "{.arg cols} collide with existing {.arg patches} columns: {.val {cols_clash}}.", + "i" = "Rename the overlay columns before attributing." + )) + } + if (is.na(sf::st_crs(patches)) || is.na(sf::st_crs(overlay))) { + cli::cli_abort("Both {.arg patches} and {.arg overlay} must have a CRS.") + } + if (is.null(time_col) != is.null(time_interval)) { + cli::cli_abort( + "{.arg time_col} and {.arg time_interval} must be supplied together." + ) + } + + # Temporal filter: keep overlay features whose time falls within the + # transition interval (bounds inclusive); NA times are dropped + if (!is.null(time_col)) { + if (!is.character(time_col) || length(time_col) != 1 || + !time_col %in% names(overlay)) { + cli::cli_abort( + "{.arg time_col} must name a single column in {.arg overlay}." + ) + } + if (!is.numeric(time_interval) || length(time_interval) != 2 || + anyNA(time_interval)) { + cli::cli_abort( + "{.arg time_interval} must be a length-2 numeric vector with no NA." + ) + } + if (time_interval[1] > time_interval[2]) { + cli::cli_abort( + "{.arg time_interval} must be increasing, got {.val {time_interval[1]}} > {.val {time_interval[2]}}." + ) + } + if (!is.numeric(overlay[[time_col]])) { + cli::cli_abort(c( + "{.arg overlay} column {.val {time_col}} must be numeric.", + "i" = "Convert dates to a numeric axis first, e.g. + {.code as.numeric(format(date, \"%Y\"))} for calendar year or + {.code as.numeric(date)} for epoch days." + )) + } + tv <- overlay[[time_col]] + overlay <- overlay[!is.na(tv) & + tv >= time_interval[1] & tv <= time_interval[2], ] + } + + # Nothing to join: return patches with typed-NA cols so the schema matches + # the joined path (indexing with NA preserves each column's type) + if (nrow(patches) == 0 || nrow(overlay) == 0) { + for (col in cols) { + patches[[col]] <- overlay[[col]][rep(NA_integer_, nrow(patches))] + } + return(patches) + } + + overlay_sel <- sf::st_make_valid( + sf::st_transform(overlay[, cols], sf::st_crs(patches)) + ) + suppressWarnings(sf::st_join( + patches, overlay_sel, + join = predicate, + left = TRUE, + largest = identical(match_mode, "largest") + )) +} diff --git a/R/dft_transition_vectors.R b/R/dft_transition_vectors.R index a484b9f..3c84f6d 100644 --- a/R/dft_transition_vectors.R +++ b/R/dft_transition_vectors.R @@ -33,6 +33,9 @@ #' When `patch_area_min` or `changes_only` drop patches, `patch_id` is numbered #' over the surviving patches (dense `1..n`), not the pre-filter grid. #' +#' @seealso [dft_transition_attribute()] to tag the returned patches from an +#' overlay polygon layer (fire perimeters, cutblocks, ...). +#' #' @export #' @examples #' r17 <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift")) diff --git a/man/dft_transition_attribute.Rd b/man/dft_transition_attribute.Rd new file mode 100644 index 0000000..05eaff0 --- /dev/null +++ b/man/dft_transition_attribute.Rd @@ -0,0 +1,154 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_transition_attribute.R +\name{dft_transition_attribute} +\alias{dft_transition_attribute} +\title{Attribute transition patches from an overlay polygon layer} +\usage{ +dft_transition_attribute( + patches, + overlay, + cols, + predicate = sf::st_intersects, + match_mode = c("all", "largest"), + time_col = NULL, + time_interval = NULL +) +} +\arguments{ +\item{patches}{An \code{sf} object of change patches, typically from +\code{\link[=dft_transition_vectors]{dft_transition_vectors()}}.} + +\item{overlay}{An \code{sf} polygon layer to attribute from (e.g. fire +perimeters, consolidated cutblocks).} + +\item{cols}{Character vector of \code{overlay} column names to carry onto each +patch. Must not collide with existing \code{patches} column names.} + +\item{predicate}{Spatial predicate function used to match patches to +overlay features, e.g. \code{\link[sf:geos_binary_pred]{sf::st_intersects()}} (default) or +\code{\link[sf:geos_binary_pred]{sf::st_within()}}. Only applies with \code{match_mode = "all"}: +largest-overlap assignment is inherently intersection-based +(\code{sf::st_join(largest = TRUE)} ignores the join predicate), so combining +a custom predicate with \code{match_mode = "largest"} is an error.} + +\item{match_mode}{How a patch is assigned when it matches more than one +overlay feature: +\itemize{ +\item \code{"all"} (default) — plain left join; a patch straddling k overlay +features appears k times (\code{patch_id} repeats). +\item \code{"largest"} — exactly one row per patch, assigned to the overlay +feature with the greatest intersection area (via +\code{sf::st_join(largest = TRUE)}; matching is by intersection, see +\code{predicate}). +}} + +\item{time_col}{Character or \code{NULL}. Name of a \strong{numeric} time column in +\code{overlay} used for temporal filtering (e.g. \code{FIRE_YEAR}, \code{HARVEST_YEAR}). +Must be supplied together with \code{time_interval}, and both must be on the +\strong{same numeric scale} (see Details). \code{NULL} (default) skips the filter.} + +\item{time_interval}{Length-2 numeric or \code{NULL}. The transition interval on +the same scale as \code{time_col}, e.g. \code{c(2017, 2023)} for calendar years. +Overlay features whose \code{time_col} falls outside the interval (both bounds +inclusive) are dropped before joining; features with an \code{NA} time are also +dropped. Patches from \code{\link[=dft_transition_vectors]{dft_transition_vectors()}} carry no time columns, so +the interval is always supplied explicitly.} +} +\value{ +\code{patches} with the \code{cols} columns joined on (\code{NA} where a patch +matches no overlay feature). Under \code{match_mode = "all"} a patch matching +several overlay features is duplicated, one row per match; under +\code{"largest"} the result has exactly one row per input patch. +} +\description{ +Tag change patches (from \code{\link[=dft_transition_vectors]{dft_transition_vectors()}}) with columns from any +overlay polygon layer — fire perimeters, cutblocks, roads, tenures — to help +separate mapped transitions by cause. Generic by design: drift carries no +domain knowledge; the caller supplies the overlay, the columns to carry, and +(optionally) the temporal filter. +} +\details{ +The overlay is transformed to the CRS of \code{patches} before joining, and run +through \code{\link[sf:valid]{sf::st_make_valid()}} first — real-world disturbance perimeters +routinely fail GEOS validity checks, which would otherwise break the +largest-overlap computation. +} +\section{Temporal filter — how \code{time_col} and \code{time_interval} must be presented}{ +The filter is a plain numeric comparison +(\verb{overlay[[time_col]] >= time_interval[1] & <= time_interval[2]}), so it is +scale-agnostic — the numbers may be calendar years, decimal years, months, +or epoch offsets — but \strong{both arguments must be numeric and on the same +scale}. \code{time_col} must name a \code{numeric} (integer or double) column; +passing a \code{Date} or \code{POSIXct} column is a hard error, not a silent +mis-comparison. Values are not coerced or rounded: a decimal year like +\code{2018.5} is compared as-is. + +To filter on dates, convert to a numeric axis first, on \strong{both} the column +and the interval: +\itemize{ +\item Calendar year (simplest for annual disturbance data): + +\if{html}{\out{
}}\preformatted{overlay$yr <- as.numeric(format(overlay$burn_date, "\%Y")) +dft_transition_attribute(..., time_col = "yr", time_interval = c(2017, 2023)) +}\if{html}{\out{
}} +\item Epoch days (\code{Date} stores days since 1970-01-01, so \code{as.numeric()} is the +coercion): + +\if{html}{\out{
}}\preformatted{overlay$t <- as.numeric(overlay$burn_date) # days since epoch +dft_transition_attribute(..., time_col = "t", + time_interval = as.numeric(as.Date(c("2017-01-01", "2023-12-31")))) +}\if{html}{\out{
}} + +(\code{POSIXct} coerces to \emph{seconds} since epoch — keep the interval in seconds +to match.) +} + +Mixing scales (e.g. epoch-day column against a \code{c(2017, 2023)} year +interval) does not error — it silently matches nothing. Keep both on one +axis. +} + +\section{Limitations}{ + +\itemize{ +\item \code{overlay} is treated as a polygon layer. Passing point or line geometries +is not supported: \code{match_mode = "largest"} compares intersection \emph{area}, +which is zero for non-polygon overlays. Use a point-in-polygon join +directly for such cases. +\item Temporal filtering is numeric-only; \code{Date}/\code{POSIXct} columns must be +coerced by the caller (see Details). +} +} + +\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") +patches <- dft_transition_vectors(result$raster, changes_only = TRUE) + +# synthetic disturbance overlay covering the western half of the AOI +bb <- sf::st_bbox(patches) +west <- sf::st_sf( + fire_year = 2018, + geometry = sf::st_as_sfc( + sf::st_bbox(c(bb["xmin"], bb["ymin"], + xmax = unname((bb["xmin"] + bb["xmax"]) / 2), bb["ymax"]), + crs = sf::st_crs(patches)) + ) +) + +# tag each patch with the fire year where it overlaps (NA elsewhere) +tagged <- dft_transition_attribute(patches, west, cols = "fire_year", + match_mode = "largest") +table(tagged$fire_year, useNA = "ifany") + +# temporal filter: only overlay features within the transition interval +tagged_2017_2020 <- dft_transition_attribute( + patches, west, cols = "fire_year", match_mode = "largest", + time_col = "fire_year", time_interval = c(2017, 2020) +) +} +\seealso{ +\code{\link[=dft_transition_vectors]{dft_transition_vectors()}} for producing the input patches. +} diff --git a/man/dft_transition_vectors.Rd b/man/dft_transition_vectors.Rd index 77379ab..355d0cc 100644 --- a/man/dft_transition_vectors.Rd +++ b/man/dft_transition_vectors.Rd @@ -73,3 +73,7 @@ head(patches_large) changes <- dft_transition_vectors(result$raster, changes_only = TRUE) head(changes) } +\seealso{ +\code{\link[=dft_transition_attribute]{dft_transition_attribute()}} to tag the returned patches from an +overlay polygon layer (fire perimeters, cutblocks, ...). +} diff --git a/planning/archive/2026-07-issue-42-transition-attribute/.gitkeep b/planning/archive/2026-07-issue-42-transition-attribute/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/planning/archive/2026-07-issue-42-transition-attribute/README.md b/planning/archive/2026-07-issue-42-transition-attribute/README.md new file mode 100644 index 0000000..3ca0f1d --- /dev/null +++ b/planning/archive/2026-07-issue-42-transition-attribute/README.md @@ -0,0 +1,20 @@ +# Issue #42 — dft_transition_attribute(): tag transition patches from an overlay layer + +## Outcome + +Added `dft_transition_attribute()` (v0.8.0), a generic helper that tags change patches from +`dft_transition_vectors()` with columns from any overlay polygon layer (fire perimeters, +cutblocks, roads), with a spatial `predicate`, a `match_mode` of `"all"` (duplicating left +join) or `"largest"` (one row per patch by greatest overlap), and an optional inclusive +`year_col` + `interval` temporal filter. drift stays domain-free; the floodplains driver +supplies the tables and interval. `spacehakr::spk_join()` was evaluated as a dependency and +rejected — its reusable part is a thin veneer over `sf::st_join`, and the load-bearing pieces +(interval semantics, area-majority) didn't exist there. Key learnings: `sf::st_join(largest = +TRUE)` silently ignores the join predicate (matching is pure intersection), so the function +aborts on custom predicate + `"largest"` rather than mis-attributing; and `cols` naming the +overlay's geometry column (e.g. `geom` from a GeoPackage) passes naive name checks but breaks +the join contract — guard with `attr(overlay, "sf_column")`. Patches carry no interval columns, +so the transition interval is an explicit argument by design. + +Closed by: branch `42-dft-transition-attribute-tag-transition` (commits 84be3a8, 40553e4, +45fb59b, 7b9a30e) / PR pending diff --git a/planning/archive/2026-07-issue-42-transition-attribute/findings.md b/planning/archive/2026-07-issue-42-transition-attribute/findings.md new file mode 100644 index 0000000..cc51530 --- /dev/null +++ b/planning/archive/2026-07-issue-42-transition-attribute/findings.md @@ -0,0 +1,76 @@ +# Findings — dft_transition_attribute(): tag transition patches from an overlay layer (optional temporal filter) (#42) + +## Issue context + +### Problem + +`dft_transition_vectors()` produces per-patch change polygons (from_class -> to_class, area, +optional zone), but there is no way to attribute *why* a change happened. A driver that wants to +separate fire from harvest from classification noise has to hand-roll spatial joins against +external overlays (fire perimeters, cutblocks, roads) every time. + +### Proposed Solution + +A generic helper that tags transition patches from an overlay polygon layer — no domain context, +reusable for any disturbance/context source: + +```r +dft_transition_attribute(patches, overlay, cols = ..., predicate = st_intersects, + year_col = NULL, interval = NULL) +``` + +- `patches` — sf from `dft_transition_vectors()`. +- `overlay` — any polygon sf. +- `cols` — overlay columns to carry onto each patch (e.g. `fire_year`, `burn_severity`). +- `predicate` / area-majority option — how a patch is assigned when it straddles overlay polys. +- `year_col` + `interval` — optional temporal filter: keep only overlay features whose year + falls within the transition interval (e.g. a fire in 2022 attributes a 2017->2023 loss; a 2012 + fire does not). + +Returns `patches` with the overlay columns joined (NA where no overlap). + +### Why generic (not fire-specific) + +The same call attributes fire (`prot_historical_fire_polys` + `FIRE_YEAR`), harvest +(consolidated cutblocks + `HARVEST_YEAR`), roads, tenures, any context layer. drift stays free +of BC/domain knowledge; the driver supplies the tables + interval logic. Surfaced from the +floodplains driver (PCEA: 28% of loss is a 2022 burn; confined-valley groups: <1%). + +## spacehakr evaluation (rejected as dependency) + +- `spacehakr::spk_join()` covers ~60-70% of the need: predicate `st_join`, mask column + selection, left-join-NA contract. But: + - Its temporal filtering is `%in%`-membership on one column — cannot express an interval, + and cannot express per-row join conditions. + - Its one-to-many `collapse` option `toString()`-concatenates duplicate values into a comma + string (with type-change warning) — exactly wrong for attribution. No area-majority logic. + - The reusable part is a thin veneer over `sf::st_join` + dplyr, both already in drift Imports. +- Taking an Imports on an experimental, barely-used, GitHub-only package couples released drift + to spacehakr's uncertain fate. Decision: implement natively; steal the left-join-NA contract + and predicate-as-function *patterns* only. +- spacehakr fate (side note, user-level): pre-drift grab-bag; audit `grep -r "spk_"` across + repos, hoist live functions, likely archive. Not part of this issue. + +## Codebase exploration (2026-07-17) + +- `dft_transition_vectors()` (R/dft_transition_vectors.R) returns sf: `patch_id` (int dense + 1..n), `transition` (chr "Trees -> Rangeland"), `area_ha` (numeric, + `as.numeric(st_area)*1e-4`), geometry (projected CRS), optional zone col. **No interval/year + columns** — interval must be user-supplied in the new function. +- Zones pattern (dft_transition_vectors.R:163-166): `st_transform(zones[zone_col], st_crs(out))` + then `suppressWarnings(st_intersection(...))` — precedent for silent CRS transform + + suppressWarnings. +- No units pkg anywhere; areas are bare numerics with `switch(unit, "m2"=1, "ha"=1e-4, "km2"=1e-6)`. +- Validation style is transitional: new functions (dft_rast_break, dft_rast_trend, dft_stac_*) + use `cli::cli_abort(c(...))` with `{.fn}`/`{.val}` markup + `match.arg()` — follow that. +- No existing predicate-arg or area-majority code. `sf::st_join(..., largest = TRUE)` natively + implements area-majority — both match modes collapse to one st_join call. +- Roxygen: markdown on, `@seealso [fn()]` (no `@family`, no `@examplesIf`), runnable + `@examples` off `system.file("extdata", ...)` (example_2017.tif, example_2020.tif, + example_aoi.gpkg). +- Tests: fixture chain classify -> transition -> vectors on bundled rasters; synthetic sf from + `terra::as.polygons(terra::ext(...))`; exact-count regression assertions; no network skips. +- `_pkgdown.yml` has no `reference:` section — auto index, no edit needed. +- `match` rejected as param name (shadows base `match()`, per code-check base-name-shadowing + convention); user approved `match_mode` enum over sf-style `largest` logical for + extensibility. diff --git a/planning/archive/2026-07-issue-42-transition-attribute/progress.md b/planning/archive/2026-07-issue-42-transition-attribute/progress.md new file mode 100644 index 0000000..6e2e1cd --- /dev/null +++ b/planning/archive/2026-07-issue-42-transition-attribute/progress.md @@ -0,0 +1,13 @@ +# Progress — dft_transition_attribute(): tag transition patches from an overlay layer (optional temporal filter) (#42) + +## Session 2026-07-17 + +- Evaluated spacehakr::spk_join() as a dependency — rejected (see findings.md); implement natively +- Plan-mode exploration (Explore + Plan agents) — phases approved by user; `match_mode` enum chosen over `largest` logical via user question +- Created branch `42-dft-transition-attribute-tag-transition` off main +- Scaffolded PWF baseline (51d9a35) +- Phase 1: tests-first contract, 10 test blocks; confirmed fail on "could not find function"; code-check round 1 Clean (84be3a8) +- Phase 2: implemented `R/dft_transition_attribute.R`; code-check found 2 real bugs — `sf::st_join(largest = TRUE)` ignores the join predicate (now aborts on custom predicate + "largest") and `cols` naming the overlay geometry column slipped validation (now guarded via `attr(overlay, "sf_column")`); round 3 Clean; 52 assertions pass, full suite 402+ green (40553e4) +- Phase 3: roxygen + reciprocal @seealso, document(), lintr clean, `devtools::check()` 0 errors / 0 warnings / 1 unrelated timestamp NOTE (45fb59b) +- Phase 4: NEWS.md 0.8.0 entry + DESCRIPTION bump as final commit +- Next: /planning-archive, then /gh-pr-push diff --git a/planning/archive/2026-07-issue-42-transition-attribute/task_plan.md b/planning/archive/2026-07-issue-42-transition-attribute/task_plan.md new file mode 100644 index 0000000..c2eafd1 --- /dev/null +++ b/planning/archive/2026-07-issue-42-transition-attribute/task_plan.md @@ -0,0 +1,59 @@ +# Task: dft_transition_attribute(): tag transition patches from an overlay layer (optional temporal filter) (#42) + +`dft_transition_vectors()` produces per-patch change polygons (from_class -> to_class, area, +optional zone), but there is no way to attribute *why* a change happened. A driver that wants to +separate fire from harvest from classification noise has to hand-roll spatial joins against +external overlays (fire perimeters, cutblocks, roads) every time. + +Approved signature: + +```r +dft_transition_attribute(patches, overlay, cols, predicate = sf::st_intersects, + match_mode = c("all", "largest"), + year_col = NULL, interval = NULL) +``` + +## Phase 1: Tests first + +- [x] Create `tests/testthat/test-dft_transition_attribute.R` with shared fixture: bundled-raster chain `dft_rast_classify` -> `dft_rast_transition` -> `dft_transition_vectors(changes_only = TRUE)` (mirrors test-dft_transition_vectors.R) +- [x] Build synthetic overlay fixture: `sf::st_as_sf(terra::as.polygons(terra::ext(...)))` sub-extents with `fire_year` (numeric) + `cause` (character) cols; one adjacent-poly pair a known patch straddles unevenly +- [x] Test: returns sf; `cols` appended; `patch_id`/`transition`/`area_ha` preserved; NA cols where no overlap +- [x] Test: `match_mode = "largest"` -> exactly `nrow(patches)` rows; straddling patch gets larger-overlap poly's attribute +- [x] Test: `match_mode = "all"` duplicates straddling patch row (exact count) +- [x] Test: temporal filter — inclusive bounds at both ends; out-of-interval overlay -> all-NA cols, nrow preserved +- [x] Test: overlay in EPSG:4326 attributes identically to pre-projected overlay (silent transform) +- [x] Test: custom `predicate` (`sf::st_within`) honoured in `"all"` mode +- [x] Test: 0-row patches -> 0-row sf with correctly-typed `cols` +- [x] Test: validation errors — non-sf inputs; `cols` missing/colliding; `year_col` without `interval` (and vice versa); bad `interval` (length/type/reversed); non-numeric year column; bad `match_mode` +- [x] Confirm tests fail for the right reason (function not found) + +## Phase 2: Implementation + +- [x] Create `R/dft_transition_attribute.R` with approved signature +- [x] Validation block: `cli::cli_abort` with `{.fn}`/`{.val}`/`{.arg}` markup; `match.arg(match_mode)`; predicate-is-function; cols checks; paired year_col+interval; CRS present +- [x] Temporal filter (inclusive numeric interval, NA years dropped) +- [x] Short-circuit path (typed-NA cols) + join path (subset -> transform -> make_valid -> st_join) +- [x] Test file green; full `devtools::test()` no regressions + +## Phase 3: Docs and examples + +- [x] Roxygen: params incl. duplicate-row semantics of `"all"`, one-row guarantee of `"largest"`, inclusive bounds, numeric-years-only, silent CRS transform, `st_make_valid` in `@details`; `@seealso [dft_transition_vectors()]` + reciprocal in `R/dft_transition_vectors.R` +- [x] Runnable `@examples` from bundled data with inline synthetic overlay; with and without `interval`; no `\dontrun` +- [x] `devtools::document()`; NAMESPACE gains export +- [x] `devtools::check()` clean + +## Phase 4: Release bookkeeping + +- [x] NEWS.md entry (new unreleased heading): generic overlay attribution, both match modes, temporal filter, ref #42 +- [x] DESCRIPTION version bump 0.7.0 -> 0.8.0 as FINAL commit + +## Out of scope + +Date/POSIXct year columns; multi-value aggregation (concatenating fire years); any BC/fire/harvest domain helpers; units-package areas; spacehakr dependency. + +## Validation + +- [x] Tests pass +- [x] `/code-check` clean on each commit +- [x] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion diff --git a/tests/testthat/test-dft_transition_attribute.R b/tests/testthat/test-dft_transition_attribute.R new file mode 100644 index 0000000..b57c836 --- /dev/null +++ b/tests/testthat/test-dft_transition_attribute.R @@ -0,0 +1,305 @@ +# synthetic patches mimicking dft_transition_vectors() output schema: +# two 1-ha squares in a projected CRS, 200 m apart +attribute_test_patches <- function() { + sq <- function(x0) { + sf::st_polygon(list(rbind( + c(x0, 0), c(x0 + 100, 0), c(x0 + 100, 100), c(x0, 100), c(x0, 0) + ))) + } + sf::st_sf( + patch_id = 1:2, + transition = c("Trees -> Rangeland", "Trees -> Crops"), + area_ha = c(1, 1), + geometry = sf::st_sfc(sq(0), sq(300), crs = "EPSG:32609") + ) +} + +# overlay straddling patch 1 unevenly: poly with fire_year 2018 covers its +# left 30 m, poly with fire_year 2022 its right 70 m. Patch 2 is outside both. +# `ignore_me` exists to prove only `cols` are carried over. +attribute_test_overlay <- function() { + rect <- function(x0, x1) { + sf::st_polygon(list(rbind( + c(x0, -10), c(x1, -10), c(x1, 110), c(x0, 110), c(x0, -10) + ))) + } + sf::st_sf( + fire_year = c(2018, 2022), + cause = c("fire", "harvest"), + ignore_me = c("x", "y"), + geometry = sf::st_sfc(rect(-10, 30), rect(30, 150), crs = "EPSG:32609") + ) +} + +test_that("attributes patches from the bundled-raster fixture chain", { + 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, changes_only = TRUE) + + overlay <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), quiet = TRUE + ) + overlay$fire_year <- 2018 + + out <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "largest") + + expect_s3_class(out, "sf") + expect_equal(nrow(out), nrow(patches)) + expect_true(all(c("patch_id", "transition", "area_ha", "fire_year") %in% + names(out))) + # AOI covers the raster, so every patch is attributed + expect_false(any(is.na(out$fire_year))) +}) + +test_that("cols are appended, patch columns preserved, NA where no overlap", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + + out <- dft_transition_attribute(patches, overlay, + cols = c("fire_year", "cause"), + match_mode = "largest") + + expect_s3_class(out, "sf") + expect_equal(out$patch_id, patches$patch_id) + expect_equal(out$transition, patches$transition) + expect_equal(out$area_ha, patches$area_ha) + expect_true(all(c("fire_year", "cause") %in% names(out))) + # only requested cols carried over + expect_false("ignore_me" %in% names(out)) + # patch 2 touches no overlay feature + expect_true(is.na(out$fire_year[out$patch_id == 2])) + expect_true(is.na(out$cause[out$patch_id == 2])) +}) + +test_that("match_mode = 'largest' assigns straddling patch by greatest overlap", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + + out <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "largest") + + # one row per patch, no duplicates + expect_equal(nrow(out), nrow(patches)) + expect_equal(sort(out$patch_id), patches$patch_id) + # patch 1 straddles 2018 (30 m) and 2022 (70 m) -> larger overlap wins + expect_equal(out$fire_year[out$patch_id == 1], 2022) +}) + +test_that("match_mode = 'all' duplicates a straddling patch per match", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + + out <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "all") + + # patch 1 matches both overlay polys, patch 2 gets one NA row + expect_equal(nrow(out), 3L) + expect_equal(sum(out$patch_id == 1), 2L) + expect_setequal(out$fire_year[out$patch_id == 1], c(2018, 2022)) + expect_true(is.na(out$fire_year[out$patch_id == 2])) +}) + +test_that("temporal filter is inclusive on both time_interval bounds", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + + # inclusive: both 2018 and 2022 kept + out_all <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "all", + time_col = "fire_year", + time_interval = c(2018, 2022)) + expect_setequal(out_all$fire_year[out_all$patch_id == 1], c(2018, 2022)) + + # only the 2022 feature survives the filter + out_22 <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "largest", + time_col = "fire_year", + time_interval = c(2022, 2023)) + expect_equal(out_22$fire_year[out_22$patch_id == 1], 2022) + + # no feature in time_interval -> all-NA cols, nrow preserved + out_none <- dft_transition_attribute(patches, overlay, + cols = c("fire_year", "cause"), + match_mode = "largest", + time_col = "fire_year", + time_interval = c(2019, 2021)) + expect_equal(nrow(out_none), nrow(patches)) + expect_true(all(is.na(out_none$fire_year))) + expect_true(all(is.na(out_none$cause))) + expect_type(out_none$fire_year, "double") + expect_type(out_none$cause, "character") +}) + +test_that("overlay features with NA year are dropped by the temporal filter", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + overlay$fire_year[2] <- NA_real_ + + out <- dft_transition_attribute(patches, overlay, cols = "cause", + match_mode = "largest", + time_col = "fire_year", + time_interval = c(2000, 2030)) + # only the 2018 feature remains -> patch 1 attributed to it + expect_equal(out$cause[out$patch_id == 1], "fire") +}) + +test_that("overlay in a different CRS is transformed silently", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + overlay_4326 <- sf::st_transform(overlay, "EPSG:4326") + + out_proj <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "largest") + out_4326 <- dft_transition_attribute(patches, overlay_4326, + cols = "fire_year", + match_mode = "largest") + + expect_equal(out_4326$fire_year, out_proj$fire_year) + expect_equal(sf::st_crs(out_4326), sf::st_crs(patches)) +}) + +test_that("custom predicate is honoured", { + patches <- attribute_test_patches() + # one poly fully containing patch 1; patch 1 also *intersects* a second poly + # that it is not within, so st_within and st_intersects must differ + contains_p1 <- sf::st_polygon(list(rbind( + c(-10, -10), c(110, -10), c(110, 110), c(-10, 110), c(-10, -10) + ))) + clips_p1 <- sf::st_polygon(list(rbind( + c(90, -10), c(150, -10), c(150, 110), c(90, 110), c(90, -10) + ))) + overlay <- sf::st_sf( + cause = c("containing", "clipping"), + geometry = sf::st_sfc(contains_p1, clips_p1, crs = "EPSG:32609") + ) + + out <- dft_transition_attribute(patches, overlay, cols = "cause", + predicate = sf::st_within, + match_mode = "all") + + expect_equal(out$cause[out$patch_id == 1], "containing") + expect_true(is.na(out$cause[out$patch_id == 2])) +}) + +test_that("match_mode defaults to 'all'", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + + # no match_mode passed -> documented default is "all" + out <- dft_transition_attribute(patches, overlay, cols = "fire_year") + + # "all" duplicates the straddling patch 1 (two overlay matches) -> 3 rows + expect_equal(nrow(out), 3L) + expect_equal(sum(out$patch_id == 1), 2L) +}) + +test_that("invalid overlay geometry is repaired (st_make_valid), not fatal", { + patches <- attribute_test_patches() + + # self-intersecting bow-tie spanning patch 1 -- invalid until st_make_valid. + # This is the real-world case the make_valid call exists for (BC fire / + # cutblock perimeters routinely fail GEOS validity). "largest" runs + # st_intersection internally, so an unrepaired input would throw here. + bowtie <- sf::st_polygon(list(rbind( + c(-10, -10), c(110, 110), c(110, -10), c(-10, 110), c(-10, -10) + ))) + overlay <- sf::st_sf( + fire_year = 2020, + geometry = sf::st_sfc(bowtie, crs = "EPSG:32609") + ) + expect_false(all(sf::st_is_valid(overlay))) # fixture really is invalid + + out <- dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "largest") + + expect_equal(nrow(out), nrow(patches)) + expect_equal(out$fire_year[out$patch_id == 1], 2020) +}) + +test_that("0-row patches return a 0-row sf with correctly-typed cols", { + patches <- attribute_test_patches()[0, ] + overlay <- attribute_test_overlay() + + out <- dft_transition_attribute(patches, overlay, + cols = c("fire_year", "cause")) + + expect_s3_class(out, "sf") + expect_equal(nrow(out), 0L) + expect_true(all(c("fire_year", "cause") %in% names(out))) + expect_type(out$fire_year, "double") + expect_type(out$cause, "character") +}) + +test_that("validation errors are raised with informative messages", { + patches <- attribute_test_patches() + overlay <- attribute_test_overlay() + + expect_error(dft_transition_attribute(data.frame(x = 1), overlay, + cols = "fire_year"), + "sf") + expect_error(dft_transition_attribute(patches, data.frame(x = 1), + cols = "fire_year"), + "sf") + expect_error(dft_transition_attribute(patches, overlay, cols = character(0)), + "cols") + expect_error(dft_transition_attribute(patches, overlay, cols = "nope"), + "nope") + # collision with an existing patches column + overlay_clash <- overlay + overlay_clash$transition <- "boom" + expect_error(dft_transition_attribute(patches, overlay_clash, + cols = "transition"), + "transition") + # predicate must be a function + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + predicate = "st_intersects"), + "predicate") + # time_col and time_interval must come together + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_col = "fire_year"), + "time_interval") + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_interval = c(2018, 2022)), + "time_col") + # time_interval shape + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_col = "fire_year", + time_interval = 2018), + "time_interval") + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_col = "fire_year", + time_interval = c("a", "b")), + "time_interval") + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_col = "fire_year", + time_interval = c(2022, 2018)), + "time_interval") + # time_col must exist and be numeric + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_col = "nope", + time_interval = c(2018, 2022)), + "time_col") + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + time_col = "cause", + time_interval = c(2018, 2022)), + "numeric") + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + match_mode = "banana")) + # sf::st_join(largest = TRUE) ignores the join predicate, so a custom + # predicate with match_mode = "largest" would silently mis-attribute + expect_error(dft_transition_attribute(patches, overlay, cols = "fire_year", + predicate = sf::st_within, + match_mode = "largest"), + "largest") + # cols naming the overlay geometry column (when named differently from the + # patches geometry) must be rejected, not silently dropped by st_join + overlay_geom <- overlay + names(overlay_geom)[names(overlay_geom) == "geometry"] <- "geom" + sf::st_geometry(overlay_geom) <- "geom" + expect_error(dft_transition_attribute(patches, overlay_geom, + cols = c("fire_year", "geom")), + "geom") +})