diff --git a/NAMESPACE b/NAMESPACE index 7a90331..a3d45f2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ export(dft_class_table) export(dft_map_interactive) export(dft_rast_classify) export(dft_rast_summarize) +export(dft_rast_transition) export(dft_stac_classes) export(dft_stac_config) export(dft_stac_fetch) diff --git a/R/dft_map_interactive.R b/R/dft_map_interactive.R index 06460bb..5d40bf2 100644 --- a/R/dft_map_interactive.R +++ b/R/dft_map_interactive.R @@ -17,8 +17,9 @@ #' `"https://titiler.example.com"`). Only used when `x` contains COG URLs. #' Defaults to `getOption("drift.titiler_url")`. If `NULL` in COG mode, an #' error is raised prompting the user to set the option. -#' @param basemaps Named character vector of provider tile IDs. The first -#' element is the default basemap. Names become radio button labels. +#' @param basemaps Named character vector of provider tile IDs or tile URL +#' templates (starting with `http`). The first element is the default +#' basemap. Names become radio button labels. #' @param legend_position Legend placement passed to [leaflet::addLegend()]. #' Set to `NULL` to suppress the legend. #' @param zoom Initial zoom level. @@ -59,7 +60,8 @@ dft_map_interactive <- function(x, source = "io-lulc", titiler_url = getOption("drift.titiler_url"), basemaps = c("Light" = "CartoDB.Positron", - "Satellite" = "Esri.WorldImagery"), + "Esri Satellite" = "Esri.WorldImagery", + "Google Satellite" = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}"), legend_position = "bottomright", zoom = 14) { rlang::check_installed(c("leaflet", "leaflet.extras")) @@ -113,7 +115,13 @@ dft_map_interactive <- function(x, # Add basemaps — first is default for (i in seq_along(basemaps)) { - map <- leaflet::addProviderTiles(map, basemaps[[i]], group = names(basemaps)[[i]]) + bm <- basemaps[[i]] + nm <- names(basemaps)[[i]] + if (grepl("^https?://", bm)) { + map <- leaflet::addTiles(map, urlTemplate = bm, group = nm) + } else { + map <- leaflet::addProviderTiles(map, bm, group = nm) + } } # Add layers diff --git a/R/dft_rast_transition.R b/R/dft_rast_transition.R new file mode 100644 index 0000000..511297a --- /dev/null +++ b/R/dft_rast_transition.R @@ -0,0 +1,143 @@ +#' Detect land cover transitions between two time steps +#' +#' Compare two classified rasters cell-by-cell and return a transition raster +#' and summary table. Each pixel is encoded with its from→to class pair. +#' +#' @param x A named list of classified `SpatRaster`s (e.g. from +#' [dft_rast_classify()]). Names identify each time step (typically years). +#' @param from Character. Name of the "before" layer in `x`. +#' @param to Character. Name of the "after" layer in `x`. +#' @param class_table A tibble with columns `code`, `class_name`, `color`. +#' When `NULL`, loaded via [dft_class_table()] using `source`. +#' @param source Character. Used to load a shipped class table when +#' `class_table` is `NULL`. One of `"io-lulc"` or `"esa-worldcover"`. +#' @param from_class Character vector of class names to include as "from" +#' classes. When `NULL` (default), all classes are included. +#' @param to_class Character vector of class names to include as "to" +#' classes. When `NULL` (default), all classes are included. +#' @param unit Character. Area unit for the summary. One of `"ha"` +#' (default), `"km2"`, or `"m2"`. +#' +#' @return A list with two elements: +#' - `raster`: A `SpatRaster` with factor levels labelled +#' `"from_class -> to_class"`. Filtered-out transitions are set to `NA`. +#' - `summary`: A tibble with columns `from_class`, `to_class`, +#' `n_cells`, `area`, `pct`. +#' @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") +#' +#' # All transitions +#' result <- dft_rast_transition(classified, from = "2017", to = "2020") +#' result$summary +#' +#' # Only tree loss to agriculture +#' tree_loss <- dft_rast_transition(classified, from = "2017", to = "2020", +#' from_class = "Trees", +#' to_class = c("Crops", "Rangeland", "Bare Ground")) +#' tree_loss$summary +dft_rast_transition <- function(x, + from, + to, + class_table = NULL, + source = "io-lulc", + from_class = NULL, + to_class = NULL, + unit = "ha") { + unit <- match.arg(unit, c("ha", "km2", "m2")) + + if (!is.list(x) || inherits(x, "SpatRaster")) { + stop("`x` must be a named list of SpatRasters, not a single SpatRaster.") + } + if (!from %in% names(x)) stop("Layer '", from, "' not found in `x`.") + if (!to %in% names(x)) stop("Layer '", to, "' not found in `x`.") + + if (is.null(class_table)) { + class_table <- dft_class_table(source) + } + + r_from <- x[[from]] + r_to <- x[[to]] + + # Resolve class names to codes + code_lookup <- stats::setNames(class_table$class_name, class_table$code) + + # Get raw integer values (strip factor) + v_from <- terra::values(r_from)[, 1] + v_to <- terra::values(r_to)[, 1] + + # Map codes to class names + name_from <- code_lookup[as.character(v_from)] + name_to <- code_lookup[as.character(v_to)] + + # Encode transitions: from_code * 1000 + to_code (supports up to 999 classes) + trans_code <- v_from * 1000L + v_to + + # Build mask for filters + + keep <- rep(TRUE, length(trans_code)) + if (!is.null(from_class)) keep <- keep & (name_from %in% from_class) + if (!is.null(to_class)) keep <- keep & (name_to %in% to_class) + + # Also mask where either raster is NA + keep <- keep & !is.na(v_from) & !is.na(v_to) + + trans_code[!keep] <- NA_integer_ + + + # Build transition raster + r_trans <- terra::rast(r_from) + terra::values(r_trans) <- trans_code + + # Build factor table from observed transitions + valid <- !is.na(trans_code) + unique_codes <- sort(unique(trans_code[valid])) + + if (length(unique_codes) == 0) { + # No transitions found — return empty + terra::set.cats(r_trans, layer = 1, + value = data.frame(id = integer(0), transition = character(0))) + summary_tbl <- tibble::tibble( + from_class = character(0), to_class = character(0), + n_cells = integer(0), area = numeric(0), pct = numeric(0) + ) + return(list(raster = r_trans, summary = summary_tbl)) + } + + from_codes <- unique_codes %/% 1000L + to_codes <- unique_codes %% 1000L + labels <- paste0( + code_lookup[as.character(from_codes)], " -> ", + code_lookup[as.character(to_codes)] + ) + + lvl_df <- data.frame(id = unique_codes, transition = labels) + terra::set.cats(r_trans, layer = 1, value = lvl_df) + + # Summary table + m2_to_unit <- switch(unit, "m2" = 1, "ha" = 1e-4, "km2" = 1e-6) + res <- terra::res(r_from) + cell_area <- res[1] * res[2] * m2_to_unit + + freq_tbl <- terra::freq(r_trans) + freq_tbl <- freq_tbl[!is.na(freq_tbl$value), ] + total_valid <- sum(keep) + + # freq on a factor raster returns labels in $value — map back to codes + label_to_code <- stats::setNames(lvl_df$id, lvl_df$transition) + freq_codes <- label_to_code[as.character(freq_tbl$value)] + + summary_tbl <- tibble::tibble( + from_class = code_lookup[as.character(freq_codes %/% 1000L)], + to_class = code_lookup[as.character(freq_codes %% 1000L)], + n_cells = as.integer(freq_tbl$count), + area = freq_tbl$count * cell_area, + pct = round(freq_tbl$count / total_valid * 100, 2) + ) + + summary_tbl <- summary_tbl[order(summary_tbl$n_cells, decreasing = TRUE), ] + + list(raster = r_trans, summary = summary_tbl) +} diff --git a/man/dft_rast_transition.Rd b/man/dft_rast_transition.Rd new file mode 100644 index 0000000..7c5e72a --- /dev/null +++ b/man/dft_rast_transition.Rd @@ -0,0 +1,68 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_rast_transition.R +\name{dft_rast_transition} +\alias{dft_rast_transition} +\title{Detect land cover transitions between two time steps} +\usage{ +dft_rast_transition( + x, + from, + to, + class_table = NULL, + source = "io-lulc", + from_class = NULL, + to_class = NULL, + unit = "ha" +) +} +\arguments{ +\item{x}{A named list of classified \code{SpatRaster}s (e.g. from +\code{\link[=dft_rast_classify]{dft_rast_classify()}}). Names identify each time step (typically years).} + +\item{from}{Character. Name of the "before" layer in \code{x}.} + +\item{to}{Character. Name of the "after" layer in \code{x}.} + +\item{class_table}{A tibble with columns \code{code}, \code{class_name}, \code{color}. +When \code{NULL}, loaded via \code{\link[=dft_class_table]{dft_class_table()}} using \code{source}.} + +\item{source}{Character. Used to load a shipped class table when +\code{class_table} is \code{NULL}. One of \code{"io-lulc"} or \code{"esa-worldcover"}.} + +\item{from_class}{Character vector of class names to include as "from" +classes. When \code{NULL} (default), all classes are included.} + +\item{to_class}{Character vector of class names to include as "to" +classes. When \code{NULL} (default), all classes are included.} + +\item{unit}{Character. Area unit for the summary. One of \code{"ha"} +(default), \code{"km2"}, or \code{"m2"}.} +} +\value{ +A list with two elements: +\itemize{ +\item \code{raster}: A \code{SpatRaster} with factor levels labelled +\code{"from_class -> to_class"}. Filtered-out transitions are set to \code{NA}. +\item \code{summary}: A tibble with columns \code{from_class}, \code{to_class}, +\code{n_cells}, \code{area}, \code{pct}. +} +} +\description{ +Compare two classified rasters cell-by-cell and return a transition raster +and summary table. Each pixel is encoded with its from→to class pair. +} +\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") + +# All transitions +result <- dft_rast_transition(classified, from = "2017", to = "2020") +result$summary + +# Only tree loss to agriculture +tree_loss <- dft_rast_transition(classified, from = "2017", to = "2020", + from_class = "Trees", + to_class = c("Crops", "Rangeland", "Bare Ground")) +tree_loss$summary +} diff --git a/tests/testthat/test-dft_rast_transition.R b/tests/testthat/test-dft_rast_transition.R new file mode 100644 index 0000000..4e2a461 --- /dev/null +++ b/tests/testthat/test-dft_rast_transition.R @@ -0,0 +1,135 @@ +test_that("dft_rast_transition returns raster and 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") + + expect_type(result, "list") + expect_named(result, c("raster", "summary")) + expect_s4_class(result$raster, "SpatRaster") + expect_s3_class(result$summary, "tbl_df") + expect_true(terra::is.factor(result$raster)) +}) + +test_that("summary has expected columns and positive values", { + 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") + s <- result$summary + + expect_true(all(c("from_class", "to_class", "n_cells", "area", "pct") %in% names(s))) + expect_true(all(s$n_cells > 0)) + expect_true(all(s$area > 0)) + expect_true(all(s$pct > 0)) +}) + +test_that("pct sums to 100", { + 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") + expect_equal(sum(result$summary$pct), 100, tolerance = 0.1) +}) + +test_that("from_class filter works", { + 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", + from_class = "Trees") + + expect_true(all(result$summary$from_class == "Trees")) +}) + +test_that("to_class filter works", { + 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", + to_class = c("Crops", "Rangeland", "Bare Ground")) + + expect_true(all(result$summary$to_class %in% c("Crops", "Rangeland", "Bare Ground"))) +}) + +test_that("both from_class and to_class filter together", { + 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", + from_class = "Trees", + to_class = "Rangeland") + + if (nrow(result$summary) > 0) { + expect_true(all(result$summary$from_class == "Trees")) + expect_true(all(result$summary$to_class == "Rangeland")) + } +}) + +test_that("impossible filter returns empty 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", + from_class = "NONEXISTENT_CLASS") + + expect_equal(nrow(result$summary), 0) +}) + +test_that("unit parameter works", { + 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") + + ha <- dft_rast_transition(classified, from = "2017", to = "2020", unit = "ha") + km2 <- dft_rast_transition(classified, from = "2017", to = "2020", unit = "km2") + + # ha should be 100x km2 + expect_equal(ha$summary$area[1] / km2$summary$area[1], 100, tolerance = 0.01) +}) + +test_that("errors on single SpatRaster", { + r <- terra::rast(system.file("extdata", "example_2017.tif", package = "drift")) + expect_error(dft_rast_transition(r, from = "2017", to = "2020"), + "named list") +}) + +test_that("errors on missing layer name", { + 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") + + expect_error(dft_rast_transition(classified, from = "2017", to = "2099"), + "2099") +}) + +test_that("transition labels use arrow format", { + 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") + lvls <- terra::cats(result$raster)[[1]] + + expect_true(all(grepl(" -> ", lvls$transition))) +}) + +test_that("summary is sorted by n_cells descending", { + 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") + s <- result$summary + + if (nrow(s) > 1) { + expect_true(all(diff(s$n_cells) <= 0)) + } +}) diff --git a/vignettes/neexdzii-kwa.Rmd b/vignettes/neexdzii-kwa.Rmd index ef6a0d8..53d9037 100644 --- a/vignettes/neexdzii-kwa.Rmd +++ b/vignettes/neexdzii-kwa.Rmd @@ -116,6 +116,38 @@ summary_tbl |> theme_minimal() ``` +## Transition Detection + +Identify exactly which pixels changed from one class to another between +time steps. `dft_rast_transition()` compares two rasters cell-by-cell and +returns a transition raster plus a summary table. + +```{r transition-all} +result <- dft_rast_transition(classified, from = "2017", to = "2023") +knitr::kable(result$summary, digits = 2, caption = "All land cover transitions 2017–2023") +``` + +### Filter to Tree Loss + +Focus on the pixels where Trees in 2017 became agriculture-related classes +by 2023. At 10 m resolution, Crops, Rangeland, and Bare Ground can +represent different phases of the same land use depending on satellite +overpass timing. + +```{r transition-tree-loss} +tree_loss <- dft_rast_transition(classified, from = "2017", to = "2023", + from_class = "Trees", + to_class = c("Crops", "Rangeland", "Bare Ground")) +tree_loss$summary +``` + +### Plot Transition Raster + +```{r plot-transition, fig.height = 4, fig.width = 7} +terra::plot(result$raster, main = "Land cover transitions 2017–2023", + axes = FALSE, mar = c(1, 1, 2, 6)) +``` + ## Interactive Map Toggle between time periods to see how land cover changed across the floodplain. @@ -123,3 +155,68 @@ Toggle between time periods to see how land cover changed across the floodplain. ```{r map-interactive, fig.height = 5} dft_map_interactive(classified, aoi = aoi) ``` + +### Transition Map + +Each transition from Trees is shown as a separate toggleable layer, +making it easy to ground-truth specific change types against satellite +basemaps. + +```{r map-transition, fig.height = 5} +# All transitions from Trees (excluding stable Trees->Trees) +all_from_trees <- dft_rast_transition(classified, from = "2017", to = "2023", + from_class = "Trees") +to_classes <- all_from_trees$summary$to_class +to_classes <- to_classes[to_classes != "Trees"] + +# Color palette for transition types +trans_colors <- c("Rangeland" = "#e74c3c", "Crops" = "#e67e22", + "Bare Ground" = "#8e44ad", "Water" = "#3498db", + "Built Area" = "#2c3e50", "Flooded Vegetation" = "#1abc9c") + +map <- leaflet::leaflet() |> + leaflet::addProviderTiles("CartoDB.Positron", group = "Light") |> + leaflet::addProviderTiles("Esri.WorldImagery", group = "Esri Satellite") |> + leaflet::addTiles("https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}", + group = "Google Satellite") + +overlay_groups <- c() + +for (tc in to_classes) { + trans <- dft_rast_transition(classified, from = "2017", to = "2023", + from_class = "Trees", to_class = tc) + if (nrow(trans$summary) == 0) next + + label <- paste0("Trees -> ", tc) + col <- if (tc %in% names(trans_colors)) trans_colors[[tc]] else "#999999" + r_proj <- terra::project(trans$raster, "EPSG:4326") + map <- leaflet::addRasterImage(map, r_proj, group = label, + colors = col, project = FALSE) + overlay_groups <- c(overlay_groups, label) +} + +# AOI outline +map <- leaflet::addPolygons(map, data = sf::st_transform(aoi, 4326), + fill = FALSE, color = "red", weight = 2, group = "AOI") +overlay_groups <- c(overlay_groups, "AOI") + +# Legend +legend_colors <- vapply(to_classes, function(tc) { + if (tc %in% names(trans_colors)) trans_colors[[tc]] else "#999999" +}, character(1)) +map <- leaflet::addLegend(map, position = "bottomright", + colors = legend_colors, + labels = paste0("Trees -> ", to_classes), + title = "Tree Loss 2017-2023", opacity = 1) + +# Center and controls +ext <- terra::ext(terra::project(classified[["2017"]], "EPSG:4326")) +map <- leaflet::setView(map, lng = mean(c(ext[1], ext[2])), + lat = mean(c(ext[3], ext[4])), zoom = 14) +map <- leaflet::addLayersControl(map, + baseGroups = c("Light", "Esri Satellite", "Google Satellite"), + overlayGroups = overlay_groups, + options = leaflet::layersControlOptions(collapsed = FALSE)) +map <- leaflet.extras::addFullscreenControl(map) +map +```