inst/doc/Integrated_planning.R

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

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "integrated-ecosystem-services-",
  fig.width = 7,
  fig.height = 5,
  fig.retina = 2,
  out.width = "100%",
  message = FALSE,
  warning = FALSE,
  error = FALSE
)

objective_aliases <- c(
  "cost",
  "rodent_control",
  "carbon_stock",
  "water_yield",
  "sediment_retention"
)

## ----load-data----------------------------------------------------------------
library(multiscape)
library(dplyr)
library(ggplot2)
library(terra)

service_data <- load_ecosystem_services()
planning_units <- service_data$planning_units
ecosystem_services <- service_data$feature_raster

service_labels <- c(
  L1_control_inv = "Rodent and lagomorph control",
  L1_stock_inv = "Carbon stock",
  L1_rend_inv = "Seasonal water yield",
  L1_reten_inv = "Sediment retention"
)

nrow(planning_units)
names(ecosystem_services)

## ----plot-services, fig.width=8, fig.height=6---------------------------------
names(ecosystem_services) <- unname(service_labels[names(ecosystem_services)])
terra::plot(ecosystem_services)
names(ecosystem_services) <- names(service_labels)

## ----action-eligibility-------------------------------------------------------
planning_attributes <- planning_units |>
  sf::st_drop_geometry() |>
  select(id, cost, area, locked_in, locked_out)

actions <- data.frame(
  id = c("conservation", "restoration"),
  name = c("Conservation", "Restoration")
)

conservation_units <- planning_attributes |>
  filter(locked_in) |>
  transmute(pu = id, action = "conservation")

restoration_units <- planning_attributes |>
  filter(!locked_in, !locked_out) |>
  transmute(pu = id, action = "restoration")

feasible_actions <- bind_rows(
  conservation_units,
  restoration_units
)

action_costs <- feasible_actions |>
  left_join(
    planning_attributes |> select(id, cost),
    by = c("pu" = "id")
  )

feasible_actions |> count(action)

## ----create-problem, message=TRUE---------------------------------------------
problem <- create_problem(
  pu = planning_units,
  features = ecosystem_services,
  cost = "cost"
) |>
  add_actions(
    actions = actions,
    include_pairs = feasible_actions,
    cost = action_costs
  )

## ----constraints--------------------------------------------------------------
unit_area <- stats::median(planning_units$area)
stopifnot(all(is.finite(planning_units$area)),
          all(abs(planning_units$area - unit_area) < 1e-8))
total_commitment_units <- ceiling(0.20 * nrow(planning_units))
restoration_target_units <- total_commitment_units - nrow(conservation_units)
restoration_target_area <- restoration_target_units * unit_area

problem <- problem |>
  add_constraint_area(
    area = restoration_target_area,
    sense = "equal",
    tolerance = 0.25 * unit_area,
    actions = "restoration"
  ) |>
  add_constraint_locked_planning_units(
    locked_in = "locked_in",
    locked_out = "locked_out"
  )

## ----deficit-effects----------------------------------------------------------
restoration_ceilings <- problem$data$dist_features |>
  filter(pu %in% restoration_units$pu) |>
  group_by(feature) |>
  summarise(outcome = max(amount), .groups = "drop")

restoration_effects <- restoration_ceilings |>
  mutate(action = "restoration")

problem <- problem |>
  add_effects(restoration_effects)

restoration_ceilings

## ----objectives---------------------------------------------------------------
problem <- problem |>
  add_objective_min_cost(
    alias = "cost",
    include_pu_cost = FALSE,
    include_action_cost = TRUE,
    actions = "restoration"
  ) |>
  add_objective_max_effect(
    features = "L1_control_inv",
    actions = "restoration",
    alias = "rodent_control"
  ) |>
  add_objective_max_effect(
    features = "L1_stock_inv",
    actions = "restoration",
    alias = "carbon_stock"
  ) |>
  add_objective_max_effect(
    features = "L1_rend_inv",
    actions = "restoration",
    alias = "water_yield"
  ) |>
  add_objective_max_effect(
    features = "L1_reten_inv",
    actions = "restoration",
    alias = "sediment_retention"
  )

print(problem)

## ----method-------------------------------------------------------------------
problem <- problem |>
  set_method_augmecon(
    primary = "cost",
    aliases = objective_aliases,
    runs = set_runs_grid(3),
    lexicographic = TRUE,
    augmentation = 1e-3,
    control = set_runs_control(
      stop_on_infeasible = FALSE,
      stop_on_no_solution = FALSE
    )
  ) |>
  set_solver_gurobi(gap_limit = 0)

## ----load-solutions-----------------------------------------------------------
solutions <- readRDS(file.path(
  "data",
  "integrated-ecosystem-services-solutions.rds"
))

## ----regenerate-solutions, eval=FALSE-----------------------------------------
# solutions <- solve(problem)
# 
# saveRDS(
#   solutions,
#   file.path("..", "data-raw", "restoration-full-solutions.rds"),
#   compress = "xz"
# )

## ----compact-regenerated-solutions, eval=FALSE--------------------------------
# source(file.path("..", "data-raw", "compact-restoration-cache.R"))
# saveRDS(
#   compact_restoration_solution_set(solutions),
#   file.path("data", "integrated-ecosystem-services-solutions.rds"),
#   compress = "xz"
# )

## ----inspect-results----------------------------------------------------------
efficient_solutions <- solutions |>
  solution_filter(feasible_only = TRUE, nondominated = TRUE) |>
  solution_unique(by = "decisions")

run_summary <- data.frame(
  requested_configurations = nrow(get_runs(solutions)),
  retained_efficient_plans = nrow(get_runs(efficient_solutions))
)
run_summary

get_runs(efficient_solutions)
get_objectives(efficient_solutions, format = "wide")

## ----prepare-objective-comparison, include=FALSE------------------------------
  objective_values <- get_objectives(
    efficient_solutions,
    format = "wide"
  )

  service_aliases <- objective_aliases[-1]
  objective_labels <- c(
    cost = "Cost",
    rodent_control = "Rodent control",
    carbon_stock = "Carbon stock",
    water_yield = "Water yield",
    sediment_retention = "Sediment retention"
  )

  minimum_cost_row <- which.min(objective_values$cost)
  minimum_cost_id <- objective_values$solution_id[minimum_cost_row]
  minimum_cost <- objective_values$cost[minimum_cost_row]

  compromise <- frontier_knee(
    efficient_solutions,
    objectives = objective_aliases,
    method = "ideal"
  )
  compromise_id <- compromise$solution_id

  service_best_ids <- setNames(
    vapply(
      service_aliases,
      function(service) {
        objective_values$solution_id[
          which.max(objective_values[[service]])
        ]
      },
      numeric(1)
    ),
    service_aliases
  )

  improvement_data <- do.call(
    rbind,
    lapply(service_aliases, function(service) {
      baseline <- objective_values[[service]][minimum_cost_row]
      ids <- objective_values$solution_id
      roles <- rep("Other efficient plan", length(ids))
      roles[ids == service_best_ids[[service]]] <- "Best for this service"
      roles[ids == compromise_id] <- "Five-objective compromise"
      roles[ids == minimum_cost_id] <- "Least cost"

      data.frame(
        solution_id = ids,
        service_alias = service,
        service = factor(
          objective_labels[[service]],
          levels = unname(objective_labels[service_aliases])
        ),
        additional_cost = 100 * (
          objective_values$cost / minimum_cost - 1
        ),
        benefit_increase = 100 * (
          objective_values[[service]] - baseline
        ) / baseline,
        role = roles,
        label = ifelse(
          roles == "Other efficient plan",
          "",
          paste("Solution", ids)
        )
      )
    })
  )

## ----cost-service-improvement, echo=FALSE, fig.width=9, fig.height=6.5--------
  ggplot(
    improvement_data,
    aes(x = additional_cost, y = benefit_increase)
  ) +
    geom_hline(yintercept = 0, colour = "grey80", linewidth = 0.4) +
    geom_point(
      aes(colour = role),
      size = 2.6,
      alpha = 0.9
    ) +
    geom_text(
      aes(
        label = label,
        colour = role,
        hjust = dplyr::case_when(
          role == "Least cost" ~ 0,
          role == "Best for this service" ~ 1,
          TRUE ~ 0.5
        )
      ),
      vjust = -0.65,
      size = 3,
      show.legend = FALSE,
      check_overlap = FALSE
    ) +
    facet_wrap(~service, scales = "free_y", ncol = 2) +
    scale_colour_manual(
      values = c(
        "Other efficient plan" = "grey70",
        "Least cost" = "#D95F02",
        "Five-objective compromise" = "#7570B3",
        "Best for this service" = "#1B9E77"
      )
    ) +
    scale_x_continuous(
      labels = function(x) paste0(round(x), "%"),
      expand = expansion(mult = c(0.03, 0.16))
    ) +
    scale_y_continuous(
      labels = function(x) paste0(round(x, 1), "%"),
      expand = expansion(mult = c(0.05, 0.18))
    ) +
    labs(
      x = "Additional cost relative to the least-cost plan",
      y = "Increase in deficit benefit relative to the least-cost plan",
      colour = NULL
    ) +
    theme_minimal(base_size = 11) +
    theme(
      panel.grid.minor = element_blank(),
      strip.text = element_text(face = "bold"),
      legend.position = "bottom"
    )

## ----objective-heatmap, echo=FALSE, fig.width=9, fig.height=7-----------------
  objective_matrix <- as.matrix(
    objective_values[, objective_aliases, drop = FALSE]
  )
  objective_ranges <- apply(objective_matrix, 2, range)
  normalized_matrix <- sweep(
    objective_matrix,
    2,
    objective_ranges[1, ],
    "-"
  )
  normalized_matrix <- sweep(
    normalized_matrix,
    2,
    objective_ranges[2, ] - objective_ranges[1, ],
    "/"
  )
  # Express every column as regret: low cost is good; high benefits are good.
  normalized_matrix[, service_aliases] <- 1 - normalized_matrix[, service_aliases]
  normalized_matrix[!is.finite(normalized_matrix)] <- 0

  heatmap_data <- do.call(
    rbind,
    lapply(seq_along(objective_aliases), function(j) {
      data.frame(
        solution_id = objective_values$solution_id,
        objective = factor(
          objective_labels[[objective_aliases[j]]],
          levels = unname(objective_labels[objective_aliases])
        ),
        relative_performance = normalized_matrix[, j]
      )
    })
  )
  heatmap_data$is_best <- heatmap_data$relative_performance < 1e-10

  solution_order <- objective_values$solution_id[
    order(objective_values$cost)
  ]
  solution_labels <- as.character(solution_order)
  solution_labels[solution_order == minimum_cost_id] <- paste0(
    minimum_cost_id,
    " - least cost"
  )
  solution_labels[solution_order == compromise_id] <- paste0(
    compromise_id,
    " - compromise"
  )
  heatmap_data$solution <- factor(
    heatmap_data$solution_id,
    levels = solution_order,
    labels = solution_labels
  )

  ggplot(
    heatmap_data,
    aes(x = objective, y = solution, fill = relative_performance)
  ) +
    geom_tile(colour = "white", linewidth = 0.5) +
    geom_point(
      data = subset(heatmap_data, is_best),
      shape = 21,
      fill = "white",
      colour = "black",
      size = 2.2,
      stroke = 0.6
    ) +
    scale_fill_gradient2(
      low = "#2166AC",
      mid = "#F7F7F7",
      high = "#B2182B",
      midpoint = 0.5,
      limits = c(0, 1),
      breaks = c(0, 0.5, 1),
      labels = c("Best", "Intermediate", "Worst")
    ) +
    labs(
      x = NULL,
      y = "Solution (increasing cost from bottom to top)",
      fill = "Observed\nperformance"
    ) +
    theme_minimal(base_size = 11) +
    theme(
      panel.grid = element_blank(),
      axis.text.x = element_text(angle = 25, hjust = 1),
      legend.position = "right"
    )

## ----objective-extremes-------------------------------------------------------
  extremes <- frontier_extremes(
    efficient_solutions,
    objectives = objective_aliases,
    ties = "first"
  )
  extremes

## ----compromise-solution------------------------------------------------------
  service_distances <- frontier_distances(
    efficient_solutions,
    objectives = service_aliases,
    reference = "ideal"
  )
  service_distances$cost <- objective_values$cost[
    match(service_distances$solution_id, objective_values$solution_id)
  ]
  service_distances$additional_cost <- 100 * (
    service_distances$cost / minimum_cost - 1
  )

  ecological_id <- service_distances$solution_id[
    which.min(service_distances$distance_to_ideal)
  ]
  specialist_id <- unname(service_best_ids[["sediment_retention"]])

  compromise[, c(
    "solution_id",
    objective_aliases,
    "knee_score",
    "knee_rank"
  )]

  data.frame(
    five_objective_compromise = compromise_id,
    closest_to_ecological_ideal = ecological_id,
    same_solution = compromise_id == ecological_id
  )

## ----cost-ecological-distance, echo=FALSE, fig.width=7.5, fig.height=5.5------
highlighted_ids <- unique(c(
  minimum_cost_id,
  compromise_id,
  unname(service_best_ids)
))

distance_plot_data <- service_distances
distance_plot_data$role <- "Other efficient plan"
distance_plot_data$role[
  distance_plot_data$solution_id %in% unname(service_best_ids)
] <- "Service-specific extreme"
distance_plot_data$role[
  distance_plot_data$solution_id == compromise_id
] <- "Five-objective compromise"
distance_plot_data$role[
  distance_plot_data$solution_id == minimum_cost_id
] <- "Least cost"

ggplot(
  distance_plot_data,
  aes(x = additional_cost, y = distance_to_ideal)
) +
  geom_point(
    aes(colour = role),
    size = 2.7,
    alpha = 0.9
  ) +
  ggrepel::geom_text_repel(
    data = subset(
      distance_plot_data,
      solution_id %in% highlighted_ids
    ),
    aes(
      label = paste("Solution", solution_id),
      colour = role
    ),
    seed = 123,
    box.padding = 0.35,
    point.padding = 0.25,
    min.segment.length = 0,
    max.overlaps = Inf,
    show.legend = FALSE,
    size = 3.1
  ) +
  scale_colour_manual(
    values = c(
      "Other efficient plan" = "grey70",
      "Least cost" = "#D95F02",
      "Five-objective compromise" = "#7570B3",
      "Service-specific extreme" = "#1B9E77"
    )
  ) +
  scale_x_continuous(
    labels = function(x) paste0(round(x), "%"),
    expand = expansion(mult = c(0.06, 0.18))
  ) +
  scale_y_continuous(
    expand = expansion(mult = c(0.08, 0.16))
  ) +
  labs(
    x = "Additional cost relative to the least-cost plan",
    y = "Joint ecosystem-service regret\n(0 = best observed joint service profile)",
    colour = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

## ----representative-plans, echo=FALSE-----------------------------------------
representative_plans <- data.frame(
  role = c(
    "Least-cost plan",
    "Five-objective compromise and water-yield best",
    "Rodent-control best",
    "Carbon-stock best",
    "Sediment-retention best"
  ),
  solution_id = c(
    minimum_cost_id,
    compromise_id,
    unname(service_best_ids[["rodent_control"]]),
    unname(service_best_ids[["carbon_stock"]]),
    unname(service_best_ids[["sediment_retention"]])
  )
)
representative_plans <- representative_plans[
  !duplicated(representative_plans$solution_id),
]

representative_summary <- representative_plans
representative_summary$additional_cost_percent <- round(
  100 * (
    objective_values$cost[
      match(
        representative_summary$solution_id,
        objective_values$solution_id
      )
    ] / minimum_cost - 1
  ),
  1
)

for (service in service_aliases) {
  baseline <- objective_values[[service]][minimum_cost_row]
  values <- objective_values[[service]][
    match(
      representative_summary$solution_id,
      objective_values$solution_id
    )
  ]
  representative_summary[[paste0(service, "_benefit_increase_percent")]] <- round(
    100 * (values - baseline) / baseline,
    1
  )
}

representative_summary

## ----map-selected-solutions, echo=FALSE, eval=FALSE, fig.keep='last', fig.width=11, fig.height=8----
# facet_labels <- setNames(
#   paste0(
#     "Solution ",
#     representative_plans$solution_id,
#     "\n",
#     representative_plans$role
#   ),
#   representative_plans$solution_id
# )
# 
# map_plot <- plot_spatial_actions(
#   efficient_solutions,
#   solutions = representative_plans$solution_id,
#   fill_values = c(
#     conservation = "#2E7D32",
#     restoration = "#E69F00"
#   ),
#   base_alpha = 0.08
# ) +
#   ggplot2::facet_wrap(
#     ~solution_id,
#     labeller = ggplot2::as_labeller(facet_labels)
#   )
# 
# solution_colours <- setNames(
#   rep("#1B9E77", nrow(representative_plans)),
#   representative_plans$solution_id
# )
# solution_colours[as.character(minimum_cost_id)] <- "#D95F02"
# solution_colours[as.character(compromise_id)] <- "#7570B3"
# 
# facet_order <- sort(representative_plans$solution_id)
# facet_colours <- unname(
#   solution_colours[as.character(facet_order)]
# )
# 
# style_facet_strip <- function(grob, fill) {
#   if (inherits(grob, "rect")) {
#     grob$gp$fill <- fill
#     grob$gp$col <- fill
#   }
#   if (inherits(grob, "text")) {
#     grob$gp$col <- "white"
#   }
#   if (!is.null(grob$children)) {
#     for (i in seq_along(grob$children)) {
#       grob$children[[i]] <- style_facet_strip(
#         grob$children[[i]],
#         fill
#       )
#     }
#   }
#   if (!is.null(grob$grobs)) {
#     for (i in seq_along(grob$grobs)) {
#       grob$grobs[[i]] <- style_facet_strip(
#         grob$grobs[[i]],
#         fill
#       )
#     }
#   }
#   grob
# }
# 
# map_grob <- ggplot2::ggplotGrob(map_plot)
# strip_index <- grep("^strip-t", map_grob$layout$name)
# strip_index <- strip_index[
#   order(
#     map_grob$layout$t[strip_index],
#     map_grob$layout$l[strip_index]
#   )
# ]
# 
# for (i in seq_along(strip_index)) {
#   map_grob$grobs[[strip_index[i]]] <- style_facet_strip(
#     map_grob$grobs[[strip_index[i]]],
#     facet_colours[i]
#   )
# }
# 
# grid::grid.newpage()
# grid::grid.draw(map_grob)

## ----map-selected-solutions-static, echo=FALSE, fig.width=11, fig.height=8----
knitr::include_graphics(
  "figures/ecosystem-services/integrated-ecosystem-services-map-1.png"
)

## ----selection-similarity, echo=FALSE, fig.width=8, fig.height=7--------------
  similarity_matrix <- selection_similarity(
    efficient_solutions,
    metric = "jaccard",
    format = "matrix"
  )

  similarity_data <- expand.grid(
    solution_x = rownames(similarity_matrix),
    solution_y = colnames(similarity_matrix),
    stringsAsFactors = FALSE
  )
  similarity_data$similarity <- as.vector(similarity_matrix)
  similarity_levels <- as.character(solution_order)
  similarity_data$solution_x <- factor(
    similarity_data$solution_x,
    levels = similarity_levels
  )
  similarity_data$solution_y <- factor(
    similarity_data$solution_y,
    levels = rev(similarity_levels)
  )

  ggplot(
    similarity_data,
    aes(x = solution_x, y = solution_y, fill = similarity)
  ) +
    geom_tile(colour = "white", linewidth = 0.25) +
    scale_fill_gradientn(
      colours = c(
        "#440154",
        "#31688E",
        "#35B779",
        "#FDE725"
      ),
      limits = range(similarity_data$similarity, na.rm = TRUE),
      breaks = c(0.6, 0.8, 1.0)
    ) +
    coord_fixed() +
    labs(
      x = "Solution",
      y = "Solution",
      fill = "Jaccard\nsimilarity"
    ) +
    theme_minimal(base_size = 10) +
    theme(
      panel.grid = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )

## ----selection-frequency------------------------------------------------------
  assignment_frequency <- selection_frequency(efficient_solutions)

  head(
    assignment_frequency[
      order(-assignment_frequency$frequency),
    ],
    15
  )