inst/doc/reserve-design-rwanda.R

In multiscape: Multi-Action, Multi-Objective Spatial Planning

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "reserve-design-rwanda-",
  fig.width = 8,
  fig.height = 5.5,
  out.width = "100%",
  eval = FALSE,
  message = FALSE,
  warning = FALSE,
  error = FALSE
)


## ----load-data----------------------------------------------------------------
# library(multiscape)
# library(dplyr)
# library(ggplot2)
# library(sf)
# 
# Sys.setenv(
#   MAPP_RWANDA_DIRECTORY = "path/to/extracted/MaPP-Tutorial-Rwanda"
# )
# 
# source("data-raw/rwanda/prepare-data.R")
# rwanda_reserve <- readRDS(
#   "data-raw/rwanda/rwanda-reserve-inputs.rds"
# )
# 

## ----planning-context, fig.width=8, fig.height=5.5, fig.cap="Human-modification cost and existing protected planning units in the MaPP Rwanda training landscape."----
# planning_map <- rwanda_reserve$planning_units |>
#   mutate(
#     protection_status = if_else(
#       protected,
#       "Existing protected area",
#       "Available for selection"
#     )
#   )
# 
# ggplot(planning_map) +
#   geom_sf(aes(fill = cost), color = NA) +
#   geom_sf(
#     data = filter(planning_map, protected),
#     fill = NA,
#     color = "#1B7837",
#     linewidth = 0.25
#   ) +
#   scale_fill_viridis_c(
#     option = "C",
#     name = "Human\nmodification"
#   ) +
#   coord_sf(datum = NA) +
#   labs(
#     title = "Planning context",
#     subtitle = "Green outlines identify units with at least 75% protected-area coverage"
#   ) +
#   theme_minimal(base_size = 10) +
#   theme(
#     panel.grid = element_blank(),
#     plot.title = element_text(face = "bold")
#   )

## ----planning-context-static, echo=FALSE, eval=TRUE, fig.cap="Human-modification cost and existing protected planning units in the MaPP Rwanda training landscape."----
knitr::include_graphics(
  "figures/rwanda/reserve-design-rwanda-planning-context-1.png"
)

## ----feature-distributions, fig.width=9, fig.height=9-------------------------
# feature_map <- rwanda_reserve$feature_distribution |>
#   left_join(rwanda_reserve$features, by = c("feature" = "id")) |>
#   group_by(feature) |>
#   mutate(relative_amount = amount / max(amount)) |>
#   ungroup() |>
#   left_join(
#     select(rwanda_reserve$planning_units, id, geometry),
#     by = c("pu" = "id")
#   ) |>
#   st_as_sf()
# 
# ggplot(feature_map) +
#   geom_sf(aes(fill = relative_amount), color = NA) +
#   facet_wrap(~name, ncol = 4) +
#   scale_fill_viridis_c(
#     option = "C",
#     limits = c(0, 1),
#     name = "Relative\namount"
#   ) +
#   coord_sf(datum = NA) +
#   theme_void(base_size = 9) +
#   theme(
#     legend.position = "bottom",
#     strip.text = element_text(face = "bold", size = 8)
#   )

## ----feature-distributions-static, echo=FALSE, eval=TRUE, fig.cap="Spatial distributions of the thirteen conservation features. Values are normalized independently within each feature to emphasize spatial pattern rather than differences in total area."----
knitr::include_graphics(
  "figures/rwanda/reserve-design-rwanda-feature-distributions-1.png"
)

## ----build-problem------------------------------------------------------------
# problem <- create_problem(
#   pu = rwanda_reserve$planning_units,
#   features = rwanda_reserve$features,
#   dist_features = rwanda_reserve$feature_distribution,
#   cost = "cost",
#   pu_id_col = "id"
# ) |>
#   add_actions(
#     actions = rwanda_reserve$actions,
#     cost = 0
#   ) |>
#   add_effects(
#     effects = data.frame(
#       action = "reserve",
#       feature = rwanda_reserve$features$id,
#       effect = 0
#     )
#   ) |>
#   add_constraint_targets_relative(0.30) |>
#   add_constraint_locked_planning_units(
#     locked_in = rwanda_reserve$planning_units$protected
#   ) |>
#   add_spatial_boundary(
#     name = "boundary",
#     include_self = TRUE,
#     edge_factor = 0.5
#   ) |>
#   add_objective_min_cost(alias = "cost") |>
#   add_objective_min_fragmentation_planning_units(
#     relation_name = "boundary",
#     alias = "fragmentation"
#   )
# 
# problem

## ----problem-static, echo=FALSE, eval=TRUE, results='asis'--------------------
problem_output <- readLines(
  "results/rwanda/problem-summary.txt",
  warn = FALSE,
  encoding = "UTF-8"
)
cat("```text\n", paste(problem_output, collapse = "\n"), "\n```\n", sep = "")

## ----define-blm---------------------------------------------------------------
# blm_values <- c(0, 1e-5, 1e-4, 1e-3, 1e-2, 2e-2, 5e-2, 1e-1, 1)
# 
# blm_problem <- problem |>
#   set_method_weighted_sum(
#     aliases = c("cost", "fragmentation"),
#     runs = set_runs_manual(
#       data.frame(
#         weight_cost = 1,
#         weight_fragmentation = blm_values
#       )
#     ),
#     normalize_weights = FALSE,
#     objective_scaling = FALSE
#   ) |>
#   set_solver_gurobi(
#     gap_limit = 0,
#     verbose = TRUE
#   )

## ----obtain-solutions---------------------------------------------------------
# if (!requireNamespace("gurobi", quietly = TRUE)) {
#   stop("Reproducing this analysis requires Gurobi and a valid licence.")
# }
# 
# blm_solutions <- solve(blm_problem)

## ----gurobi-blm-log-static, echo=FALSE, eval=TRUE, results='asis'-------------
gurobi_blm_log <- readLines(
  "results/rwanda/gurobi-blm-log-excerpt.txt",
  warn = FALSE
)
cat("```text\n", paste(gurobi_blm_log, collapse = "\n"), "\n```\n", sep = "")

## ----inspect-blm-runs---------------------------------------------------------
# blm_run_summary <- get_runs(blm_solutions) |>
#   left_join(get_objectives(blm_solutions), by = "solution_id") |>
#   mutate(blm = blm_values[run_id]) |>
#   select(run_id, solution_id, blm, status, runtime, gap, cost, fragmentation)
# 
# blm_run_summary

## ----inspect-blm-runs-static, echo=FALSE, eval=TRUE---------------------------
blm_run_summary_static <- read.csv(
  "results/rwanda/blm-run-summary.csv",
  check.names = FALSE
)

knitr::kable(
  blm_run_summary_static,
  digits = c(0, 0, 5, NA, 3, 4, 0, 0),
  caption = "Stored run metadata and objective values for the BLM analysis."
)

## ----prepare-blm-results, include=FALSE---------------------------------------
# blm_results <- get_runs(blm_solutions) |>
#   left_join(get_objectives(blm_solutions), by = "solution_id") |>
#   mutate(blm = blm_values[run_id])
# 
# minimum_cost <- min(blm_results$cost, na.rm = TRUE)
# minimum_cost_fragmentation <- blm_results$fragmentation[
#   which.min(blm_results$cost)
# ]
# 
# blm_results <- blm_results |>
#   mutate(
#     additional_cost = 100 * (cost - minimum_cost) / minimum_cost,
#     boundary_km = fragmentation / 1000,
#     blm_label = paste0("BLM = ", format(blm, scientific = TRUE))
#   )

## ----plot-blm, fig.cap="Cost and spatial-configuration outcomes obtained by changing the boundary length modifier. Values are expressed relative to the minimum-cost solution (BLM = 0)."----
# ggplot(
#   blm_results,
#   aes(boundary_km, additional_cost)
# ) +
#   geom_path(color = "grey70", linewidth = 0.5) +
#   geom_point(color = "#2C7FB8", size = 2.8) +
#   ggrepel::geom_text_repel(
#     aes(label = blm_label),
#     size = 3,
#     min.segment.length = 0,
#     seed = 500
#   ) +
#   scale_y_continuous(labels = function(x) paste0(round(x, 1), "%")) +
#   labs(
#     x = "Exposed boundary (km)",
#     y = "Additional cost relative to the minimum-cost plan",
#     title = "BLM calibration",
#     subtitle = "Higher BLM values exchange additional cost for a more compact network"
#   ) +
#   theme_minimal(base_size = 11) +
#   theme(plot.title = element_text(face = "bold"))

## ----plot-blm-static, echo=FALSE, eval=TRUE, fig.cap="Cost and spatial-configuration outcomes obtained by changing the boundary length modifier. Boundary is shown in its original spatial units (converted to kilometres), while cost is expressed relative to the minimum-cost solution (BLM = 0)."----
knitr::include_graphics(
  "figures/rwanda/reserve-design-rwanda-plot-blm-1.png"
)

## ----define-epsilon-----------------------------------------------------------
# epsilon_limits <- seq(
#   min(blm_results$fragmentation),
#   max(blm_results$fragmentation),
#   length.out = 8
# )
# 
# epsilon_problem <- problem |>
#   set_method_epsilon_constraint(
#     primary = "cost",
#     aliases = c("cost", "fragmentation"),
#     runs = set_runs_manual(
#       data.frame(eps_fragmentation = epsilon_limits)
#     )
#   ) |>
#   set_solver_gurobi(
#     gap_limit = 0,
#     verbose = TRUE
#   )
# 
# epsilon_solutions <- solve(epsilon_problem)

## ----gurobi-epsilon-log-static, echo=FALSE, eval=TRUE, results='asis'---------
gurobi_log <- readLines(
  "results/rwanda/gurobi-epsilon-log-excerpt.txt",
  warn = FALSE
)
cat("```text\n", paste(gurobi_log, collapse = "\n"), "\n```\n", sep = "")

## ----prepare-epsilon-results, include=FALSE-----------------------------------
# epsilon_results <- get_runs(epsilon_solutions) |>
#   left_join(
#     get_objectives(epsilon_solutions),
#     by = "solution_id"
#   )
# 
# epsilon_results <- epsilon_results |>
#   filter(!is.na(solution_id), is.finite(cost), is.finite(fragmentation)) |>
#   mutate(
#     additional_cost = 100 * (cost - minimum_cost) / minimum_cost,
#     boundary_km = fragmentation / 1000
#   )

## ----inspect-epsilon-runs-----------------------------------------------------
# epsilon_run_summary <- get_runs(epsilon_solutions) |>
#   left_join(get_objectives(epsilon_solutions), by = "solution_id") |>
#   mutate(epsilon_fragmentation = epsilon_limits[run_id]) |>
#   select(
#     run_id,
#     solution_id,
#     epsilon_fragmentation,
#     status,
#     runtime,
#     gap,
#     cost,
#     fragmentation
#   )
# 
# epsilon_run_summary

## ----inspect-epsilon-runs-static, echo=FALSE, eval=TRUE-----------------------
epsilon_run_summary_static <- read.csv(
  "results/rwanda/epsilon-run-summary.csv",
  check.names = FALSE
)

knitr::kable(
  epsilon_run_summary_static,
  digits = c(0, 0, 0, NA, 3, 4, 0, 0),
  caption = paste(
    "Stored run metadata, epsilon limits, and objective values for the",
    "epsilon-constraint analysis."
  )
)

## ----compare-methods, fig.cap="Solutions generated using fixed BLM coefficients and explicit epsilon limits. Epsilon-constraint reveals the attainable cost of progressively reducing exposed boundary."----
# comparison_data <- bind_rows(
#   blm_results |>
#     transmute(
#       solution_id,
#       additional_cost,
#       boundary_km,
#       method = "BLM weighted sum"
#     ),
#   epsilon_results |>
#     transmute(
#       solution_id,
#       additional_cost,
#       boundary_km,
#       method = "Epsilon-constraint"
#     )
# )
# 
# ggplot(
#   comparison_data,
#   aes(
#     boundary_km,
#     additional_cost,
#     color = method,
#     shape = method
#   )
# ) +
#   geom_point(size = 3, alpha = 0.9) +
#   scale_color_manual(
#     values = c(
#       "BLM weighted sum" = "#D95F02",
#       "Epsilon-constraint" = "#1B9E77"
#     )
#   ) +
#   scale_x_continuous(labels = function(x) paste0(round(x, 1), "%")) +
#   scale_y_continuous(labels = function(x) paste0(round(x, 1), "%")) +
#   labs(
#     x = "Exposed boundary (km)",
#     y = "Additional cost relative to the minimum-cost plan",
#     color = NULL,
#     shape = NULL,
#     title = "Two ways to express the same planning tension"
#   ) +
#   theme_minimal(base_size = 11) +
#   theme(
#     legend.position = "bottom",
#     plot.title = element_text(face = "bold")
#   )

## ----compare-methods-static, echo=FALSE, eval=TRUE, fig.cap="Solutions generated using fixed BLM coefficients and explicit epsilon limits. Moving toward a lower exposed boundary requires progressively greater cost."----
knitr::include_graphics(
  "figures/rwanda/reserve-design-rwanda-compare-methods-1.png"
)

## ----select-plans, include=FALSE----------------------------------------------
# cost_range <- diff(range(epsilon_results$cost))
# fragmentation_range <- diff(range(epsilon_results$fragmentation))
# 
# epsilon_ranked <- epsilon_results |>
#   mutate(
#     normalized_cost = if (cost_range > 0) {
#       (cost - min(cost)) / cost_range
#     } else {
#       0
#     },
#     normalized_fragmentation = if (fragmentation_range > 0) {
#       (fragmentation - min(fragmentation)) / fragmentation_range
#     } else {
#       0
#     },
#     distance_to_ideal = sqrt(
#       normalized_cost^2 + normalized_fragmentation^2
#     )
#   )
# 
# representative_ids <- c(
#   "Minimum cost" = epsilon_ranked$solution_id[which.min(epsilon_ranked$cost)],
#   "Observed compromise" = epsilon_ranked$solution_id[
#     which.min(epsilon_ranked$distance_to_ideal)
#   ],
#   "Minimum fragmentation" = epsilon_ranked$solution_id[
#     which.min(epsilon_ranked$fragmentation)
#   ]
# )
# 
# representative_ids <- representative_ids[!duplicated(representative_ids)]
# 
# action_results <- get_actions(epsilon_solutions)
# 
# selected_units <- action_results |>
#   filter(
#     solution_id %in% unname(representative_ids),
#     selected == 1
#   ) |>
#   mutate(
#     plan = names(representative_ids)[
#       match(solution_id, representative_ids)
#     ]
#   ) |>
#   select(pu, solution_id, plan)
# 
# representative_maps <- bind_rows(lapply(
#   names(representative_ids),
#   function(plan_name) {
#     rwanda_reserve$planning_units |>
#       select(id, protected) |>
#       mutate(plan = plan_name)
#   }
# )) |>
#   left_join(selected_units, by = c("id" = "pu", "plan" = "plan")) |>
#   mutate(
#     status = case_when(
#       protected ~ "Existing protected area",
#       !is.na(solution_id) ~ "Additional selected unit",
#       TRUE ~ "Not selected"
#     ),
#     plan = factor(plan, levels = names(representative_ids))
#   )

## ----representative-maps, fig.width=9, fig.height=4.8, fig.cap="Representative reserve designs from the epsilon-constraint frontier. Green units were already protected; orange units are additional selections made under different spatial limits."----
# action_plot <- plot_spatial_actions(
#   epsilon_solutions,
#   solutions = unname(representative_ids),
#   layout = "single",
#   ncol = 2,
#   fill_values = c(reserve = "#F39C12"),
#   base_fill = "#F2F2F2",
#   base_alpha = 1
# )
# 
# action_plot +
#   facet_wrap(
#     ~solution_id,
#     ncol = 2,
#     labeller = as_labeller(
#       setNames(names(representative_ids), representative_ids)
#     )
#   ) +
#   geom_sf(
#     data = filter(rwanda_reserve$planning_units, protected),
#     fill = "#238B45",
#     color = NA
#   ) +
#   labs(title = "Alternative reserve designs")

## ----representative-maps-static, echo=FALSE, eval=TRUE, fig.cap="Representative reserve designs from the epsilon-constraint frontier. Existing protected units occur in every plan; orange units are additions selected to meet the 30% targets under different spatial limits."----
knitr::include_graphics(
  "figures/rwanda/reserve-design-rwanda-representative-maps-1.png"
)

## ----selection-frequency, fig.width=7.5, fig.height=5.2, fig.cap="Selection frequency across the epsilon-constraint solutions. Existing protected units are shown separately because they were locked into every plan."----
# frequency_data <- action_results |>
#   filter(selected == 1) |>
#   count(pu, name = "selected_runs") |>
#   mutate(frequency = selected_runs / nrow(epsilon_results))
# 
# frequency_map <- rwanda_reserve$planning_units |>
#   left_join(frequency_data, by = c("id" = "pu")) |>
#   mutate(
#     frequency = coalesce(frequency, 0),
#     frequency_available = if_else(protected, NA_real_, frequency)
#   )
# 
# ggplot(frequency_map) +
#   geom_sf(aes(fill = frequency_available), color = NA) +
#   geom_sf(
#     data = filter(frequency_map, protected),
#     fill = "#238B45",
#     color = NA
#   ) +
#   scale_fill_viridis_c(
#     option = "B",
#     limits = c(0, 1),
#     na.value = "#238B45",
#     name = "Selection\nfrequency"
#   ) +
#   coord_sf(datum = NA) +
#   labs(
#     title = "Spatial stability across the frontier",
#     subtitle = "Green units were already protected and locked into every solution"
#   ) +
#   theme_minimal(base_size = 10) +
#   theme(
#     panel.grid = element_blank(),
#     plot.title = element_text(face = "bold")
#   )

## ----selection-frequency-static, echo=FALSE, eval=TRUE, fig.cap="Selection frequency across the epsilon-constraint solutions. Existing protected units are shown separately because they were locked into every plan."----
knitr::include_graphics(
  "figures/rwanda/reserve-design-rwanda-selection-frequency-1.png"
)

## ----references, results='asis', echo=FALSE, eval=TRUE------------------------
cat(paste0(
  "TNC. 2024. *The Marxan Planning Platform: Tutorial (Rwanda).* ",
  "Global Science, The Nature Conservancy, Arlington, Virginia, USA. ",
  "<https://marxanplanning.org>.\n\n",
  "Kennedy, C. M., et al. 2019. Managing the middle: a shift in conservation ",
  "priorities based on the global human modification gradient. ",
  "*Global Change Biology* 25: 811-826. ",
  "<https://doi.org/10.1111/gcb.14549>.\n\n",
  "Olson, D. M., et al. 2001. Terrestrial Ecoregions of the World: A New Map ",
  "of Life on Earth. *BioScience* 51: 933-938. ",
  "<https://doi.org/10.1641/0006-3568(2001)051%5B0933:TEOTWA%5D2.0.CO;2>.\n\n",
  "Rondinini, C., et al. 2011. Global habitat suitability models of ",
  "terrestrial mammals. *Philosophical Transactions of the Royal Society B* ",
  "366: 2633-2641. <https://doi.org/10.1098/rstb.2011.0113>.\n\n",
  "Theobald, D. M., et al. 2020. Earth transformed: detailed mapping of ",
  "global human modification from 1990 to 2017. ",
  "*Earth System Science Data* 12: 1953-1972. ",
  "<https://doi.org/10.5194/essd-12-1953-2020>."
))