## ----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. ", ".\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. ", ".\n\n", "Olson, D. M., et al. 2001. Terrestrial Ecoregions of the World: A New Map ", "of Life on Earth. *BioScience* 51: 933-938. ", ".\n\n", "Rondinini, C., et al. 2011. Global habitat suitability models of ", "terrestrial mammals. *Philosophical Transactions of the Royal Society B* ", "366: 2633-2641. .\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. ", "." ))