Appendix - Land Use Land Cover Change
library(drift)
library(sf)
library(terra)
library(tmap)
library(maptiles)
library(ggplot2)
library(dplyr)
library(tidyr)
sf::sf_use_s2(FALSE)
lulc_dir <- here::here("data", "lulc")
# Pre-computed by scripts/lulc_classify.R
lulc_all <- readRDS(file.path(lulc_dir, "lulc_summary.rds"))
# Spatial inputs
subbasins <- sf::st_read(
file.path(lulc_dir, "subbasins.gpkg"), quiet = TRUE
) |> sf::st_transform(4326)
fp_file <- file.path(lulc_dir, "floodplain.gpkg")
floodplain <- sf::st_read(fp_file, layer = "co_ff04", quiet = TRUE) |> sf::st_transform(4326)
# name_basin comes from break_points.csv via fresh::frs_watershed_split5.1 Modelled Floodplain
Land cover change is assessed within the modelled floodplain — the area delineated by the Valley Confinement Algorithm [VCA; Nagel et al. (2014), flooded package] for coho rearing/spawning habitat on streams of 3rd order and greater, with lakes and wetlands included. The VCA uses a bankfull depth regression (Hall, Holzer, and Beechie 2007) scaled by flood_factor = 4 to approximate the functional floodplain — the area regularly influenced by high flows that sustains riparian vegetation, off-channel habitat, and large wood recruitment processes critical to salmon. The flood factor of 4 was selected as a conservative estimate between the empirically validated historical floodplain extent [flood factor = 3 on 10 m DEM; Hall, Holzer, and Beechie (2007)] and the recommended range of 5-7 for valley bottom delineation (Nagel et al. 2014), accounting for the coarser vertical resolution of the 25 m DEM used here. The modelled floodplain covers 148.2 km² along the Neexdzii Kwa (Bulkley River) and tributaries. All percentages and areas reported below refer to land cover within this modelled floodplain, not the full watershed.
# Use scenario-specific raster dir, fall back to legacy
fp_dir <- here::here("data", "lulc", "rasters", "co_ff04")
if (!dir.exists(fp_dir)) fp_dir <- here::here("data", "lulc", "rasters", "floodplain")
# Load classified rasters as named list
classified <- list(
"2017" = terra::rast(file.path(fp_dir, "classified_2017.tif")),
"2020" = terra::rast(file.path(fp_dir, "classified_2020.tif")),
"2023" = terra::rast(file.path(fp_dir, "classified_2023.tif"))
)
# Load transition raster, filter summary to agriculture classes only
ag_classes <- c("Crops", "Rangeland", "Bare Ground")
trans_rast <- terra::rast(file.path(fp_dir, "transition.tif"))
trans_cats <- terra::cats(trans_rast)[[1]]
trans_summary <- tibble::tibble(
id = trans_cats$id,
transition = trans_cats$transition,
from_class = "Trees",
to_class = sub("Trees -> ", "", trans_cats$transition)
) |>
dplyr::filter(to_class %in% ag_classes)
transition <- list(raster = trans_rast, summary = trans_summary)
# Build interactive map with drift
m <- drift::dft_map_interactive(
classified,
aoi = floodplain,
transition = transition,
source = "io-lulc"
)
# Transition groups for layer control
trans_groups <- trans_summary$transition
# High-contrast sub-basin palette for interactive map (Brewer Set3, 9 classes)
sb_pal <- RColorBrewer::brewer.pal(9, "Set3")
# Zoom to floodplain extent on load
bb <- sf::st_bbox(sf::st_transform(floodplain, 4326))
m <- m |>
leaflet::fitBounds(bb[["xmin"]], bb[["ymin"]], bb[["xmax"]], bb[["ymax"]]) |>
leaflet::addPolygons(
data = subbasins,
fillColor = sb_pal[seq_len(nrow(subbasins))],
fillOpacity = 0.3,
color = "#333333",
weight = 2, opacity = 0.8, group = "Sub-basins",
label = ~name_basin
) |>
leaflet::addLegend(
position = "bottomleft",
colors = sb_pal[seq_len(nrow(subbasins))],
labels = subbasins$name_basin,
title = "Sub-basin",
opacity = 0.8,
group = "Sub-basins"
) |>
leaflet::addProviderTiles("OpenTopoMap", group = "OpenTopoMap") |>
leaflet::addLayersControl(
baseGroups = c("Light", "Esri Satellite", "Google Satellite", "OpenTopoMap"),
overlayGroups = c(names(classified), trans_groups, "AOI", "Sub-basins"),
options = leaflet::layersControlOptions(collapsed = FALSE)
) |>
leaflet::hideGroup(setdiff(names(classified), names(classified)[1]))
mFigure 5.196: Land cover classification within the modelled floodplain (IO LULC 10 m). Use radio buttons to toggle between 2017, 2020, and 2023 classified layers. Transition overlays show where Trees converted to other classes between 2017 and 2023. Sub-basin boundaries are shown as a toggleable overlay.
cache <- here::here("data", "spatial")
streams_all <- readRDS(file.path(cache, "streams.rds"))
lakes <- readRDS(file.path(cache, "lakes.rds"))
roads <- readRDS(file.path(cache, "roads.rds"))
railway <- readRDS(file.path(cache, "railway.rds"))
towns <- sf::st_sf(
name = c("Houston", "Smithers", "Burns Lake", "Topley"),
geometry = sf::st_sfc(
sf::st_point(c(-126.648, 54.398)),
sf::st_point(c(-127.176, 54.779)),
sf::st_point(c(-125.764, 54.230)),
sf::st_point(c(-126.246, 54.566)),
crs = 4326
)
)
# Dissolved roads by route
roads_dissolved <- roads |>
dplyr::filter(!is.na(route)) |>
dplyr::group_by(route) |>
dplyr::summarise(geom = sf::st_union(geom), .groups = "drop")
# Bounding box from sub-basins with padding tuned for 8.5x11
bbox <- sf::st_bbox(subbasins)
bbox["ymax"] <- bbox["ymax"] + 0.06
bbox["ymin"] <- bbox["ymin"] - 0.06
bbox["xmin"] <- bbox["xmin"] - 0.04
bbox["xmax"] <- bbox["xmax"] + 0.04
# Clip streams and lakes to bbox
bbox_clip_sf <- sf::st_as_sfc(bbox) |> sf::st_set_crs(4326)
streams_clip <- streams_all |> sf::st_make_valid() |> sf::st_crop(bbox_clip_sf)
lakes_clip <- lakes |> sf::st_make_valid() |> sf::st_crop(bbox_clip_sf)
# Stream labels (5th order+ named)
stream_label_pts <- streams_clip |>
dplyr::filter(!is.na(gnis_name) & stream_order >= 5) |>
dplyr::group_by(gnis_name) |>
dplyr::summarise(geom = sf::st_union(geom), .groups = "drop") |>
dplyr::mutate(geometry = sf::st_point_on_surface(geom)) |>
sf::st_set_geometry("geometry") |>
sf::st_set_crs(4326)
# Lake labels (>1 km²)
lake_labels <- lakes_clip[!is.na(lakes_clip$name) & lakes_clip$area_km2 > 1, ]
# Hillshade basemap
bbox_sf <- sf::st_as_sfc(bbox) |> sf::st_set_crs(4326)
relief <- maptiles::get_tiles(bbox_sf, provider = "Esri.WorldShadedRelief",
zoom = 10, crop = TRUE)
basemap_stars <- stars::st_as_stars(relief)
# Earthy sub-basin palette (9 sub-basins)
earth_pal <- c("#A67C52", "#C4A882", "#8B9E6B", "#D4A76A",
"#7A8B6E", "#BFA07A", "#9B8E6E", "#C2B280", "#8FA87C")
tmap::tmap_mode("plot")
tmap::tm_shape(basemap_stars, bbox = bbox) +
tmap::tm_rgb() +
# Sub-basins
tmap::tm_shape(subbasins) +
tmap::tm_polygons(
fill = "name_basin",
fill.scale = tmap::tm_scale_categorical(values = earth_pal),
fill_alpha = 0.35,
col = "#5C4033", lwd = 1.5,
fill.legend = tmap::tm_legend(show = FALSE)
) +
tmap::tm_text("name_basin", size = 0.65, col = "#3B2716", fontface = "bold",
options = tmap::opt_tm_text(shadow = TRUE)) +
# Floodplain
tmap::tm_shape(floodplain) +
tmap::tm_polygons(fill = "#2d6a2e", fill_alpha = 0.30,
col = "#2d6a2e", lwd = 0.8) +
# Lakes
tmap::tm_shape(lakes_clip) +
tmap::tm_polygons(fill = "#c6ddf0", col = "#7ba7cc", lwd = 0.4, fill_alpha = 0.85) +
tmap::tm_shape(lake_labels) +
tmap::tm_text("name", size = 0.55, col = "#1a5276", fontface = "italic") +
# Streams
tmap::tm_shape(streams_clip |> dplyr::filter(stream_order >= 4)) +
tmap::tm_lines(col = "#7ba7cc", lwd = 0.4) +
tmap::tm_shape(streams_clip |> dplyr::filter(stream_order >= 5)) +
tmap::tm_lines(col = "#7ba7cc", lwd = 0.8) +
tmap::tm_shape(stream_label_pts) +
tmap::tm_text("gnis_name", size = 0.60, col = "#1a5276", fontface = "italic") +
# Railway
tmap::tm_shape(railway) +
tmap::tm_lines(col = "black", lwd = 1.2) +
tmap::tm_shape(railway) +
tmap::tm_lines(col = "white", lwd = 0.6, lty = "42") +
# Roads
tmap::tm_shape(roads_dissolved |> dplyr::filter(route == "16")) +
tmap::tm_lines(col = "#c0392b", lwd = 2.0) +
tmap::tm_shape(roads_dissolved |> dplyr::filter(route != "16")) +
tmap::tm_lines(col = "#e67e22", lwd = 1.2) +
# Towns
tmap::tm_shape(towns) +
tmap::tm_dots(fill = "black", size = 0.30) +
tmap::tm_text("name", size = 0.65, xmod = 0.8, ymod = -0.6,
col = "grey10", fontface = "bold") +
# Scale bar
tmap::tm_scalebar(breaks = c(0, 10, 20),
position = c("left", "bottom"),
text.size = 0.6) +
# Legend
tmap::tm_add_legend(
type = "polygons",
labels = "Modelled floodplain",
fill = "#2d6a2e",
col = "#2d6a2e",
lwd = 1
) +
tmap::tm_add_legend(
type = "lines",
labels = c("Highway 16", "Secondary highway", "Railway (CN)"),
col = c("#c0392b", "#e67e22", "black"),
lwd = c(2, 1.2, 1.2)
) +
tmap::tm_layout(
frame = TRUE,
inner.margins = c(0.005, 0.005, 0.005, 0.005),
outer.margins = c(0.002, 0.002, 0.002, 0.002),
legend.position = c("right", "top"),
legend.frame = TRUE,
legend.bg.color = "white",
legend.bg.alpha = 0.85
)5.2 Land Cover by Sub-Basin
IO Land Use Land Cover (10 m, Sentinel-2) is classified for 2017, 2020, and 2023 using the drift pipeline. Each sub-basin’s modelled floodplain is analyzed independently.
ag_classes <- c("Crops", "Rangeland", "Bare Ground")
# Compute change from 2017 to 2023 for each class and sub-basin (pct and area)
lulc_delta <- lulc_all |>
filter(year %in% c(2017, 2023)) |>
select(name_basin, year, class_name, pct, area) |>
pivot_wider(names_from = year, values_from = c(pct, area), names_sep = "_") |>
mutate(delta_pct = pct_2023 - pct_2017,
delta_ha = area_2023 - area_2017) |>
arrange(name_basin, class_name)
# Grouped agriculture delta
lulc_grouped <- lulc_all |>
mutate(
group = case_when(
class_name == "Trees" ~ "Trees",
class_name %in% ag_classes ~ "Agriculture",
TRUE ~ NA_character_
)
) |>
filter(!is.na(group)) |>
group_by(name_basin, year, group) |>
summarize(pct = sum(pct), area = sum(area), .groups = "drop")
lulc_grouped_delta <- lulc_grouped |>
filter(year %in% c(2017, 2023)) |>
pivot_wider(names_from = year, values_from = c(pct, area), names_sep = "_") |>
mutate(delta_pct = pct_2023 - pct_2017,
delta_ha = area_2023 - area_2017)major <- lulc_all |>
group_by(class_name) |>
summarize(max_pct = max(pct), .groups = "drop") |>
filter(max_pct >= 5) |>
pull(class_name)
lulc_all |>
filter(class_name %in% major) |>
mutate(year = factor(year)) |>
ggplot(aes(x = year, y = pct, fill = class_name)) +
geom_col(position = "dodge") +
facet_wrap(~name_basin, scales = "free_y") +
scale_fill_brewer(palette = "Set2", name = "Land Cover") +
theme_minimal() +
labs(
y = "% of Modelled Floodplain",
x = NULL
)
Figure 5.197: Land cover composition of modelled floodplain by sub-basin (IO LULC 10 m, 2017-2023).
5.3 Agriculture Definition
For this analysis, “Agriculture” groups three IO LULC classes that can represent different phases of the same land use at 10 m resolution: Rangeland (fallow, grazed pasture), Crops (actively growing), and Bare Ground (tilled or harvested). A single field can be classified differently depending on the satellite overpass date relative to the growing season.
lulc_all |>
filter(class_name %in% ag_classes) |>
mutate(year = factor(year),
class_name = factor(class_name, levels = c("Bare Ground", "Crops", "Rangeland"))) |>
ggplot(aes(x = year, y = pct, fill = class_name)) +
geom_col() +
facet_wrap(~name_basin, nrow = 3) +
scale_fill_manual(
values = c("Rangeland" = "#FFBB22", "Crops" = "#E49635", "Bare Ground" = "#A59B8F"),
name = NULL
) +
theme_minimal() +
labs(
y = "% of Modelled Floodplain",
x = NULL
)
Figure 5.198: Composition of the Agriculture superclass within each sub-basin’s modelled floodplain. Rangeland dominates; Crops and Bare Ground are minor components.
5.4 Trees vs Agriculture
lulc_grouped |>
mutate(year = factor(year)) |>
ggplot(aes(x = year, y = pct, color = group, group = group)) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
facet_wrap(~name_basin, nrow = 2) +
scale_color_manual(values = c("Trees" = "#006400", "Agriculture" = "#c0392b"), name = NULL) +
theme_minimal() +
labs(
y = "% of Modelled Floodplain",
x = NULL
)
Figure 5.199: Trees and Agriculture (Crops + Rangeland + Bare Ground) as percentage of modelled floodplain by sub-basin, 2017-2023.
5.5 Sub-Basin Ranking
Tree loss and agriculture gain do not sum to zero because other IO LULC classes (Built Area, Water, Flooded Vegetation) also gain or lose area between years. The other_ha column captures this residual — positive values indicate net area flowing into non-tree, non-agriculture classes (e.g., Built Area expansion, water level changes). Classification noise at 10 m pixel boundaries and differences in Sentinel-2 overpass timing between years also contribute.
# One-row-per-sub-basin summary with key metrics
fp_areas <- lulc_all |>
filter(year == 2017) |>
group_by(name_basin) |>
summarize(fp_ha = round(sum(area), 1), .groups = "drop")
summary_tab <- lulc_grouped_delta |>
select(name_basin, group, delta_pct, delta_ha) |>
pivot_wider(names_from = group, values_from = c(delta_pct, delta_ha)) |>
left_join(fp_areas, by = "name_basin") |>
mutate(
tree_loss_pct = round(delta_pct_Trees, 1),
tree_loss_ha = round(delta_ha_Trees, 1),
ag_gain_pct = round(delta_pct_Agriculture, 1),
ag_gain_ha = round(delta_ha_Agriculture, 1),
other_ha = round(-tree_loss_ha - ag_gain_ha, 1)
) |>
select(name_basin, fp_ha, tree_loss_pct, tree_loss_ha, ag_gain_pct, ag_gain_ha, other_ha) |>
arrange(tree_loss_ha)my_caption <- "Tree cover loss and agriculture gain within modelled floodplain by sub-basin (2017-2023). Negative values = loss. Agriculture = Crops + Rangeland + Bare Ground. other_ha = residual area change attributed to remaining classes (Built Area, Water, Flooded Vegetation, etc.)."
my_tab_caption()knitr::kable(summary_tab, digits = 1,
caption = "Tree cover loss and agriculture gain within modelled floodplain by sub-basin (2017-2023). other_ha = residual attributed to remaining classes.")ranking_ha <- lulc_grouped_delta |>
select(name_basin, group, delta_ha) |>
pivot_wider(names_from = group, values_from = delta_ha) |>
mutate(net_shift = Agriculture - Trees) |>
arrange(Trees)
ranking_ha |>
pivot_longer(cols = c(Trees, Agriculture), names_to = "class", values_to = "delta") |>
mutate(name_basin = factor(name_basin, levels = ranking_ha$name_basin)) |>
ggplot(aes(x = name_basin, y = delta, fill = class)) +
geom_col(position = "dodge") +
scale_fill_manual(values = c("Trees" = "#006400", "Agriculture" = "#c0392b"), name = NULL) +
geom_hline(yintercept = 0, linewidth = 0.5) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(
y = "Change in Area (ha, 2017-2023)",
x = NULL
)
Figure 5.200: Change in tree cover and agriculture area (ha) within modelled floodplain by sub-basin, 2017-2023. Agriculture = Crops + Rangeland + Bare Ground.
5.6 Land Cover Composition
my_caption <- "Land cover composition (% of modelled floodplain) by sub-basin and year. fp_ha = total modelled floodplain area."
my_tab_caption()wide <- lulc_all |>
left_join(fp_areas, by = "name_basin") |>
select(name_basin, fp_ha, year, class_name, pct) |>
pivot_wider(names_from = year, values_from = pct) |>
arrange(name_basin, class_name)
wide |>
my_dt_table(cols_freeze_left = 3, page_length = 25)wide <- lulc_all |>
left_join(fp_areas, by = "name_basin") |>
select(name_basin, fp_ha, year, class_name, pct) |>
pivot_wider(names_from = year, values_from = pct) |>
arrange(name_basin, class_name)
knitr::kable(wide, digits = 1,
caption = "Land cover composition (% of modelled floodplain) by sub-basin and year. fp_ha = total modelled floodplain area.")