From b988b9694fdd87a148d4683ac506d32c4b8af044 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Mon, 9 Mar 2026 16:09:08 -0700 Subject: [PATCH] Restructure vignette with bookdown cross-references and transition filtering Rename vignette to land-cover-change.Rmd, convert to bookdown::html_document2 with numbered cross-references. Filter transition tables to >1% of total area, group tree loss classes as Agriculture, filter transition raster to significant changes only. Add bookdown to Suggests and VignetteBuilder. Relates to #11 Co-Authored-By: Claude Opus 4.6 --- CLAUDE.md | 2 +- DESCRIPTION | 3 +- vignettes/land-cover-change.Rmd | 228 ++++++++++++++++++++++++++++++++ vignettes/neexdzii-kwa.Rmd | 162 ----------------------- 4 files changed, 231 insertions(+), 164 deletions(-) create mode 100644 vignettes/land-cover-change.Rmd delete mode 100644 vignettes/neexdzii-kwa.Rmd diff --git a/CLAUDE.md b/CLAUDE.md index 025ca45..72eb95e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -505,7 +505,7 @@ drift::dft_map_interactive(classified, aoi = aoi) - Class colors come from drift's shipped class tables (IO LULC, ESA WorldCover) - For production COGs on S3, `dft_map_interactive()` serves tiles via titiler — set `options(drift.titiler_url = "...")` -- See the [drift vignette](https://www.newgraphenvironment.com/drift/articles/neexdzii-kwa.html) for a worked example (Neexdzii Kwa floodplain, 2017-2023) +- See the [drift vignette](https://www.newgraphenvironment.com/drift/articles/land-cover-change.html) for a worked example (Neexdzii Kwa floodplain, 2017-2023) # Code Check Conventions diff --git a/DESCRIPTION b/DESCRIPTION index f1a7fcf..76c210f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,6 +29,7 @@ Imports: terra, tibble Suggests: + bookdown, DT, flooded, gdalcubes, @@ -50,4 +51,4 @@ Remotes: NewGraphEnvironment/flooded, NewGraphEnvironment/gq Config/testthat/edition: 3 -VignetteBuilder: knitr +VignetteBuilder: knitr, bookdown diff --git a/vignettes/land-cover-change.Rmd b/vignettes/land-cover-change.Rmd new file mode 100644 index 0000000..68fda28 --- /dev/null +++ b/vignettes/land-cover-change.Rmd @@ -0,0 +1,228 @@ +--- +title: "Land Cover Change Detection for Floodplains" +output: bookdown::html_document2 +vignette: > + %\VignetteIndexEntry{Land Cover Change Detection for Floodplains} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 5 +) +``` + +This vignette demonstrates the drift pipeline using a small floodplain reach +on Neexdzii Kwa (Upper Bulkley River) in northern BC. We compare Esri IO LULC +land cover across 2017, 2020, and 2023 to track vegetation and land use change +in the riparian zone. + +The AOI polygon was delineated using the +[flooded](https://github.com/NewGraphEnvironment/flooded) package, which +identifies floodplain extents from DEMs and stream networks. + +The example data ships with the package — no STAC queries or database +connections needed. + +## Load Data + +```{r load-data} +library(drift) +library(terra) +library(sf) + +# AOI polygon (floodplain delineated via flooded package) +aoi <- sf::st_read( + system.file("extdata", "example_aoi.gpkg", package = "drift"), + quiet = TRUE +) + +# IO LULC rasters for 3 years +years <- c(2017, 2020, 2023) +rasters <- lapply(years, function(yr) { + terra::rast(system.file("extdata", paste0("example_", yr, ".tif"), + package = "drift")) +}) +names(rasters) <- years +``` + +## Classify + +Apply IO LULC class names and colors from the shipped class table. + +```{r classify} +classified <- dft_rast_classify(rasters, source = "io-lulc") + +# Check factor levels +terra::levels(classified[["2020"]])[[1]] +``` + +## Classified Rasters + +Figure \@ref(fig:plot-classified) shows the three classified time steps +side by side. + +```{r plot-classified, fig.cap = "Classified land cover for the Neexdzii Kwa floodplain reach across three time steps.", fig.height = 4, fig.width = 10} +stacked <- terra::rast(classified) +names(stacked) <- names(classified) +terra::plot(stacked, axes = FALSE, mar = c(1, 1, 2, 1)) +``` + +## Area Summary + +Table \@ref(tab:change-table) shows area by class for each year and the +net change between 2017 and 2023, sorted by magnitude of change. + +```{r summarize} +summary_tbl <- dft_rast_summarize(classified, source = "io-lulc", unit = "ha") +``` + +```{r change-table} +library(dplyr) +library(tidyr) + +change <- summary_tbl |> + dplyr::select(year, class_name, area) |> + tidyr::pivot_wider(names_from = year, values_from = area, values_fill = list(area = 0)) |> + dplyr::mutate( + change = `2023` - `2017`, + pct_change = round(change / `2017` * 100, 1) + ) |> + dplyr::arrange(dplyr::desc(abs(change))) + +knitr::kable(change, digits = 2, caption = "Net land cover change 2017--2023 (ha), sorted by absolute change.") +``` + +## Vegetation Change + +Trees and Rangeland show the clearest signal in Figure +\@ref(fig:plot-vegetation) — tree cover declining while rangeland expands. + +```{r plot-vegetation, fig.cap = "Dominant vegetation classes over time in the Neexdzii Kwa floodplain.", fig.height = 4, fig.width = 8} +library(ggplot2) + +summary_tbl |> + dplyr::filter(class_name %in% c("Trees", "Rangeland")) |> + ggplot(aes(x = year, y = area, fill = year)) + + geom_col() + + facet_wrap(~class_name, scales = "free_y") + + scale_fill_brewer(palette = "YlGnBu") + + labs(y = "Area (ha)", x = NULL, fill = "Year", + title = "Vegetation cover in Neexdzii Kwa floodplain") + + theme_minimal() +``` + +## Transition Detection + +`dft_rast_transition()` compares two rasters cell-by-cell and returns a +transition raster plus a summary table. Table \@ref(tab:transition-stable) +shows the area that remained in the same class (stable pixels), while Table +\@ref(tab:transition-change) shows pixels that changed class. Only +transitions representing more than 1% of the total area are shown. + +```{r transition-compute} +result <- dft_rast_transition(classified, from = "2017", to = "2023") + +stable <- result$summary |> + dplyr::filter(from_class == to_class) |> + dplyr::filter(pct >= 1) |> + dplyr::arrange(dplyr::desc(area)) + +changed <- result$summary |> + dplyr::filter(from_class != to_class) |> + dplyr::filter(pct >= 1) |> + dplyr::arrange(dplyr::desc(area)) +``` + +```{r transition-stable} +knitr::kable(stable, digits = 2, + caption = "Stable land cover 2017--2023 (only transitions >1% of total area shown).") +``` + +```{r transition-change} +knitr::kable(changed, digits = 2, + caption = "Land cover transitions 2017--2023 (only transitions >1% of total area shown).") +``` + +### Grouping Classes for Domain-Specific Analysis + +Fine-grained LULC classes can be grouped into categories relevant to +a specific analysis. Here we demonstrate grouping Crops, Rangeland, +and Bare Ground as "Agriculture" — at 10 m resolution these classes +can represent different phases of the same land use depending on +satellite overpass timing. + +Table \@ref(tab:tree-loss) shows the area of Trees in 2017 that +transitioned to agriculture-related classes by 2023. + +```{r tree-loss} +ag_classes <- c("Crops", "Rangeland", "Bare Ground") + +# All transitions from Trees to get total Trees-origin pixel count +all_from_trees <- dft_rast_transition(classified, from = "2017", to = "2023", + from_class = "Trees") +total_tree_cells <- sum(all_from_trees$summary$n_cells) + +# Filter to agriculture classes +tree_loss <- dft_rast_transition(classified, from = "2017", to = "2023", + from_class = "Trees", + to_class = ag_classes) + +# Relabel as Agriculture and compute pct of all Trees-origin pixels +tree_loss_tbl <- tree_loss$summary |> + dplyr::mutate(to_class = "Agriculture") |> + dplyr::group_by(from_class, to_class) |> + dplyr::summarize(n_cells = sum(n_cells), area = sum(area), .groups = "drop") |> + dplyr::mutate(pct_of_trees = round(n_cells / total_tree_cells * 100, 2)) + +knitr::kable(tree_loss_tbl, digits = 2, + caption = "Tree loss to agriculture (Crops + Rangeland + Bare Ground) 2017--2023. Percent is of all pixels classified as Trees in 2017.") +``` + +### Transition Raster + +Figure \@ref(fig:plot-transition) maps pixels that changed class +between 2017 and 2023. Only transitions representing more than 1% of +the total area are shown; minor transitions are masked. The AOI outline +is shown in red. + +```{r plot-transition, fig.cap = "Spatial distribution of land cover transitions 2017--2023 (only transitions >1% of total area shown).", fig.height = 4, fig.width = 7} +trans_vals <- terra::values(result$raster)[, 1] +lvls <- terra::cats(result$raster)[[1]] + +# Get codes for transitions >= 1% (excluding stable) +sig_labels <- changed$from_class # already filtered to >1% and from != to +sig_transitions <- paste0(changed$from_class, " -> ", changed$to_class) +sig_codes <- lvls$id[lvls$transition %in% sig_transitions] + +# Mask everything except significant transitions +change_vals <- rep(NA_integer_, length(trans_vals)) +change_vals[trans_vals %in% sig_codes] <- trans_vals[trans_vals %in% sig_codes] +r_change <- terra::rast(result$raster) +terra::values(r_change) <- change_vals + +# Keep only significant factor levels +change_lvls <- lvls[lvls$id %in% sig_codes, , drop = FALSE] +terra::set.cats(r_change, layer = 1, value = change_lvls) + +terra::plot(r_change, main = "Land cover transitions 2017\u20132023", + axes = FALSE, mar = c(1, 1, 2, 6)) +plot(sf::st_geometry(sf::st_transform(aoi, terra::crs(r_change))), + add = TRUE, border = "red", lwd = 2) +``` + +## Multi-Year Classification and Transition Overlay + +Toggle between classified time periods and overlay tree loss transition +layers to ground-truth change against multiple satellite basemaps +(Figure \@ref(fig:map-interactive)). + +```{r map-interactive, fig.cap = "Classified land cover by year (radio toggle) with tree loss transitions overlaid as toggleable layers.", fig.height = 5} +tree_trans <- dft_rast_transition(classified, from = "2017", to = "2023", + from_class = "Trees") +dft_map_interactive(classified, aoi = aoi, transition = tree_trans) +``` diff --git a/vignettes/neexdzii-kwa.Rmd b/vignettes/neexdzii-kwa.Rmd deleted file mode 100644 index de0193e..0000000 --- a/vignettes/neexdzii-kwa.Rmd +++ /dev/null @@ -1,162 +0,0 @@ ---- -title: "Land Cover Change — Neexdzii Kwa" -output: rmarkdown::html_vignette -vignette: > - %\VignetteIndexEntry{Land Cover Change — Neexdzii Kwa} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -```{r setup, include = FALSE} -knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>", - fig.width = 7, - fig.height = 5 -) -``` - -This vignette demonstrates the drift pipeline using a small floodplain reach -on Neexdzii Kwa (Upper Bulkley River) in northern BC. We compare Esri IO LULC -land cover across 2017, 2020, and 2023 to track vegetation and land use change -in the riparian zone. - -The AOI polygon was delineated using the -[flooded](https://github.com/NewGraphEnvironment/flooded) package, which -identifies floodplain extents from DEMs and stream networks. - -The example data ships with the package — no STAC queries or database -connections needed. - -## Load Data - -```{r load-data} -library(drift) -library(terra) -library(sf) - -# AOI polygon (floodplain delineated via flooded package) -aoi <- sf::st_read( - system.file("extdata", "example_aoi.gpkg", package = "drift"), - quiet = TRUE -) - -# IO LULC rasters for 3 years -years <- c(2017, 2020, 2023) -rasters <- lapply(years, function(yr) { - terra::rast(system.file("extdata", paste0("example_", yr, ".tif"), - package = "drift")) -}) -names(rasters) <- years -``` - -## Classify - -Apply IO LULC class names and colors from the shipped class table. - -```{r classify} -classified <- dft_rast_classify(rasters, source = "io-lulc") - -# Check factor levels -terra::levels(classified[["2020"]])[[1]] -``` - -## Plot Classified Rasters - -```{r plot-classified, fig.height = 4, fig.width = 10} -# Stack into a single multi-layer SpatRaster for panel plot -stacked <- terra::rast(classified) -names(stacked) <- names(classified) -terra::plot(stacked, axes = FALSE, mar = c(1, 1, 2, 1)) -``` - -## Summarize - -Compute area by class for each year. - -```{r summarize} -summary_tbl <- dft_rast_summarize(classified, source = "io-lulc", unit = "ha") -summary_tbl -``` - -## Area Change Table - -```{r change-table} -library(dplyr) -library(tidyr) - -change <- summary_tbl |> - dplyr::select(year, class_name, area) |> - tidyr::pivot_wider(names_from = year, values_from = area, values_fill = list(area = 0)) |> - dplyr::mutate( - change = `2023` - `2017`, - pct_change = round(change / `2017` * 100, 1) - ) |> - dplyr::arrange(dplyr::desc(abs(change))) - -knitr::kable(change, digits = 2, caption = "Land cover change 2017–2023 (ha)") -``` - -## Vegetation Change - -Bar plot of the dominant vegetation classes over time. Trees and Rangeland show -the clearest signal — tree cover declining while rangeland expands. - -```{r plot-vegetation, fig.height = 4, fig.width = 8} -library(ggplot2) - -summary_tbl |> - dplyr::filter(class_name %in% c("Trees", "Rangeland")) |> - ggplot(aes(x = year, y = area, fill = year)) + - geom_col() + - facet_wrap(~class_name, scales = "free_y") + - scale_fill_brewer(palette = "YlGnBu") + - labs(y = "Area (ha)", x = NULL, fill = "Year", - title = "Vegetation cover in Neexdzii Kwa floodplain") + - 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 classified time periods and overlay transition layers to -ground-truth change against multiple satellite basemaps. - -```{r map-interactive, fig.height = 5} -# All transitions from Trees -tree_trans <- dft_rast_transition(classified, from = "2017", to = "2023", - from_class = "Trees") - -dft_map_interactive(classified, aoi = aoi, transition = tree_trans) -```