SPARCLE Proposed Methods

This document walks through an example workflow for the CI-UCSB SPARCLE project. This project aims to support spatial planning by identifying where allocation of land use for conservation, agriculture, and energy can best support humans and biodiversity under climate change. The project launched in December 2021, and this workflow is meant to provide a starting point for discussion around the project methodology.

Click ‘code’ in the individual sections to view or hide the code used in the analysis. At the top right of this page, there is also the option to show all code.

library(raster)
library(sf)
library(terra)
library(Matrix)
library(rnaturalearth)
library(fasterize)
library(prioritizr)
library(tidyverse)
library(patchwork)
library(kableExtra)
library(calecopal)
source("methodsTrialx_functions.R")

Set up

Study area

Study area: California

Planning units: 5km pixels

sa <- terra::rasterize(
  x = shape,
  y = template) %>% 
  terra::crop(shape)

plot(sa, col = "navy", axes = FALSE)

ids <- terra::as.data.frame(
  sa, 
  cells = TRUE) %>% 
  dplyr::select(!layer) %>% 
  rename('id' = cell)

Cost

The cost of all land reflects the land area. This means that cost is consistent across the landscape. While in our California example using socio-economic cost estimates may make sense, this would be less applicable for identifying areas of global conservation importance.

The budge is percent of land conserved. In our example workflow, we solve separately for three budgets: 30% land conserved, 40% land conserved, and 50% land conserved. Note this is total land conserved, meaning existing protected areas will contribute to that total. In the case of California, approximately 22% of land is already conserved, so our 30% budget solution represents allocation of an additional 8%.

land_cost <- rep(25, length(ids$id)) # homogeneous - all pus cost 25
budgets <- c(0.3, 0.4, 0.5)

Conservation features

Conservation features: Species current & future distributions

Here we use species-level current and future distributions for over 17,000 terrestrial vertebrates. Future models are under a SSP3/RCP7.0 2050 scenario. Current and future ranges are weighted equally and used as individual features. We use sparse matrices for efficiency.

Current

We are using species ranges for over 17,000 species. Below is an example of one species range, as well as summaries of species richness for visualization. While we show summaries, please note that in the analysis we are using each species range rather than a summary of the biodiversity data.

taxas <- c("amphibians", "birds", "mammals", "reptiles")

for (t in seq_along(taxas)){
  taxa <- taxas[t]
  temp_mat <- readRDS(paste0("methodsTrial_data/interim/", 
                             roi_abb, "_", 
                             taxa, "_current_mat.tif"))
  if(t == 1){
    cur_mat <- temp_mat
  } else {
    cur_mat <- rbind(cur_mat, temp_mat)}
}

remove(temp_mat)

# view a single species
# custom function
# matrix_plot(
#   vals = cur_mat[114, ],
#   crop = sa,
#   main = rownames(cur_mat)[114])

# view richness
matrix_plot(
  vals = colSums(cur_mat),
  crop = sa,
  main = "Current species richness")

cur_species <- rownames(cur_mat)

cur_features_df <- data.frame(
  id = 1:length(cur_species),
  name = cur_species)

Future

for (t in seq_along(taxas)){
  taxa <- taxas[t]
  temp_mat <- readRDS(paste0("methodsTrial_data/interim/", 
                             roi_abb, "_", 
                             taxa, "_ssp37050_mat.tif"))
  if(t == 1){
    fut_mat <- temp_mat
  } else {
    fut_mat <- rbind(fut_mat, temp_mat)}
}

remove(temp_mat)

matrix_plot(
  vals = colSums(fut_mat),
  crop = sa,
  main = "Future species richness (2050)")

fut_species <- rownames(fut_mat)

fut_features_df <- data.frame(
  id = 1:length(fut_species),
  name = fut_species)

Targets

We use species-level relative targets (% of range) that are a function of range size (Rodrigues et al., 2004; Hansen et al., 2020).

# bring in global ranges
global_ranges <- read_csv("C:/Users/cbrock/OneDrive - Conservation International Foundation/Documents/Conservation International/Data/Sparse Matrices/current_and_projections_5km/summaries/allVerts_globalRanges.csv") %>% 
  dplyr::select(species, global_current_range_km2, global_future_range_km2) %>% 
  rename(name = species)

cur_ranges <- left_join(cur_features_df, global_ranges, by = "name") %>% 
  pluck("global_current_range_km2")
cur_targets <- calculate_targets(x = cur_ranges)

fut_ranges <- left_join(fut_features_df, global_ranges, by = "name") %>% 
  pluck("global_future_range_km2")
fut_targets <- calculate_targets(x = fut_ranges)
all_targets <- c(cur_targets, fut_targets)

ggplot() + 
  geom_histogram(aes(x = all_targets),
                 binwidth = 0.05) +
  scale_x_continuous(labels = scales::percent) + 
  theme_minimal() + 
  labs(x = "Relative target",
       y = "Number of species")

# combine current + future
features_mat <- rbind(cur_mat, fut_mat)

cur_species <- paste0(rownames(cur_mat), " current")
fut_species <- paste0(rownames(fut_mat), " future")

species <- c(cur_species, fut_species)

cur_df <- data.frame(
  id = 1:length(cur_species),
  name = cur_species)

features_df <- data.frame(
  id = 1:length(species),
  name = species)

Constraints

Locked in

We lock in existing protected areas for all prioritizations using WDPA data (UNEP-WCMC and IUCN, 2021).

wdpa_rds <- readRDS(paste0("methodsTrial_data/interim/", roi_abb, "_wdpa.tif"))
# convert to indices
wdpa_bin <- wdpa_rds
wdpa_bin[wdpa_bin < 50] <- 0 
wdpa_bin[wdpa_bin >= 50] <- 1
wdpa <- which(wdpa_bin == 1)

matrix_plot(
  vals = wdpa_bin, 
  crop = sa,
  main = "Locked in: WDPA")

Locked-out

In all conservation prioritizations, we lock out land that is currently (2015) classified as agricultural or urban land under Chen et al 2022.

ag <- rast(paste0("methodsTrial_data/interim/", roi_abb, "_crops.tif"))
ag_mask <- ag
ag_mask[ag_mask < 1] <- NA
ag_vals <- cells(ag_mask)

urban <- rast(paste0("methodsTrial_data/interim/", roi_abb, "_urban.tif"))
urban_mask <- urban
urban_mask[urban_mask < 1] <- NA
urban_vals <- cells(urban_mask)

cur_ag_urban <- unique(c(ag_vals, urban_vals))
cur_ag_urban <- cur_ag_urban[which(cur_ag_urban %in% c(ids$id))]

# convert to indices which I don't love...
cur_ag_urban <- which(ids$id %in% cur_ag_urban)

plot(ag, 
     main = "Locked out agriculture",
     col = bin_cols,
     axes = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

plot(urban, 
     main = "Locked out urban",
     col = bin_cols,
     axes = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

# establish general problem function

# boundary data
boundary_data <- boundary_matrix(raster(sa))
boundary_data <- boundary_data[ids$id, ids$id] # crop to just land

# function for problem
my_problem <- function(x, 
                       features, 
                       rij_matrix, 
                       budget, 
                       targets,
                       locked_in,
                       locked_out) {
  p <- problem(
    x = x, 
    features = features,
    rij_matrix = rij_matrix) %>% 
    add_min_shortfall_objective(budget = budget) %>% 
    add_relative_targets(targets = targets) %>% 
    add_locked_out_constraints(locked_out = locked_out) %>% 
    add_locked_in_constraints(locked_in = locked_in) %>% 
    add_boundary_penalties(
      penalty = 0.00001,
      data = boundary_data) %>%
    add_gurobi_solver(
      gap = 0.05,
      time_limit = 172800,
      threads = 8)
}

Scenarios

In our example workflow, we’ve explored four scenarios: 1) Baseline Land Use + Current Species Targets; 2) Baseline Land Use + Climate Change Species Targets; 3) Future Land Use Expansion + Climate Change Species Targets; 4) Future Land Use Expansion that Avoids Biodiversity + Climate Change Species Targets.

In Scenario 1, we only lock out current land use for agriculture and urban areas and only include targets for current species distributions. In Scenario 2, we only lock out current land use and we include targets for both current and future species distributions. In Scenario 3, we lock out current land use for agriculture and urban areas and the results from Johnson et al., for areas of potential land use conversion for agriculture and renewable energy in a 2050 BAU scenario, and we use current and future species distributions. In Scenario 4, we first identify potential areas for land use conversion using conservation areas as a cost layer, Oakleaf et al., 2019 DPI as a features layer, and the summed overlap of Johnson et al. land use conversion with DPI as a target. We then lock out current land use for agriculture and urban areas and our new priorities for agriculture and energy expansion and include targets for current and future species distributions. See this information summarized in the table below.

Scenario 1

Baseline Land Use + Current Species Targets

Here we only consider current land use and include targets based on current species ranges.

Solve

for (i in seq_along(budgets)){
  
  rel_budget <- budgets[i] # % of total land 
  
  p <- my_problem(
    x = land_cost,
    features = features_df,
    rij_matrix = features_mat,
    budget = sum(land_cost) * rel_budget,
    targets = c(cur_targets, rep(0, length(fut_targets))),
    locked_in = wdpa,
    locked_out = cur_ag_urban)
  
  s <- solve(p)
  
  s_rast <- rasterize_soln(s)
  
  assign(paste0("s1_rast_", rel_budget),
         s_rast)
  
  # while we didn't include future targets in the problem, we will
  # bring in future targets now to evaluate target coverage
  tc_s <- eval_target_coverage_summary(p, s) %>% 
    mutate('budget' = rel_budget) %>% 
    mutate(relative_target = all_targets) %>% 
    rowwise() %>% 
    mutate(
      met = case_when(relative_held >= relative_target ~ TRUE,
           T ~ FALSE)) %>% 
    ungroup()
  
  assign(paste0("s1_tc_", rel_budget),
         tc_s)
  
  assign(paste0("s1_met_", rel_budget),
         sum(tc_s$met))
}

## Conservation Problem
##   planning units: matrix (17271 units)
##   cost:           min: 25, max: 25
##   features:       Afghanodon mustersi current, Ambystoma californiense current, Ambystoma gracile current, ... (1864 features)
##   objective:      Minimum shortfall objective [budget (215887.5)]
##   targets:        Relative targets [targets (min: 0, max: 0.960399690588772)]
##   decisions:      default
##   constraints:    <Locked in planning units [3766 locked units]
##                    Locked out planning units [2308 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (1e-05), zones]>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.05), node_file_start (-1), numeric_focus (0), presolve (2), threads (8), time_limit (172800), verbose (1)]

s1_rast <- sum(s1_rast_0.3, s1_rast_0.4, s1_rast_0.5, na.rm = TRUE)
plot_ranked(s1_rast,
            main = "Scenario 1")

s1_tcs <- bind_rows(s1_tc_0.3,
                    s1_tc_0.4,
                    s1_tc_0.5) %>% 
  mutate(budget = paste0(scales::percent(budget), " budget"))

s1_met <- s1_tcs %>% 
  group_by(budget) %>% 
  summarize(met = sum(as.numeric(met))) %>% 
  ungroup() %>% 
  mutate(pct_met = met/length(all_targets)) %>% 
  mutate(met = scales::comma(met),
         pct_met = scales::percent(pct_met))

kable(s1_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,712	91.8%
40% budget	1,774	95.2%
50% budget	1,810	97.1%

ggplot() +
    geom_histogram(
      data = s1_tcs,
      aes(x = relative_held),
      binwidth = 0.05) + 
    scale_x_continuous(labels = scales::percent) + 
    labs(
      title = paste0("Feature representation by prioritization"),
      x = "Percent coverage of features",
      y = "Number of features") +
    theme_minimal() + 
  facet_wrap(~budget, nrow = 3)

Scenario 2

Baseline Land Use + Climate Change Species Targets

Here we only lock out current land use and we include targets for both current and future species distributions.

Solve

for (i in seq_along(budgets)){
  
  rel_budget <- budgets[i] # % of total land 
  
  p <- my_problem(
    x = land_cost,
    features = features_df,
    rij_matrix = features_mat,
    budget = sum(land_cost) * rel_budget,
    targets = all_targets,
    locked_in = wdpa,
    locked_out = cur_ag_urban)
  
  s <- solve(p)
  
  s_rast <- rasterize_soln(s)
  
  assign(paste0("s2_rast_", rel_budget),
         s_rast)
  
  tc_s <- eval_target_coverage_summary(p, s) %>% 
    mutate('budget' = rel_budget)

  assign(paste0("s2_tc_", rel_budget),
         tc_s)
  
  assign(paste0("s2_met_", rel_budget),
         sum(tc_s$met))
}

## Conservation Problem
##   planning units: matrix (17271 units)
##   cost:           min: 25, max: 25
##   features:       Afghanodon mustersi current, Ambystoma californiense current, Ambystoma gracile current, ... (1864 features)
##   objective:      Minimum shortfall objective [budget (215887.5)]
##   targets:        Relative targets [targets (min: 0.00991914163715929, max: 1)]
##   decisions:      default
##   constraints:    <Locked out planning units [2308 locked units]
##                    Locked in planning units [3766 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (1e-05), zones]>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.05), node_file_start (-1), numeric_focus (0), presolve (2), threads (8), time_limit (172800), verbose (1)]

s2_rast <- sum(s2_rast_0.3, s2_rast_0.4, s2_rast_0.5, na.rm = TRUE)
plot_ranked(s2_rast,
            main = "Scenario 2")

s2_tcs <- bind_rows(s2_tc_0.3,
                    s2_tc_0.4,
                    s2_tc_0.5) %>% 
  mutate(budget = paste0(scales::percent(budget), " budget"))

s2_met <- s2_tcs %>% 
  group_by(budget) %>% 
  summarize(met = sum(as.numeric(met))) %>% 
  ungroup() %>% 
  mutate(pct_met = met/length(all_targets)) %>% 
  mutate(met = scales::comma(met),
         pct_met = scales::percent(pct_met))

kable(s2_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,761	94.5%
40% budget	1,807	96.9%
50% budget	1,829	98.1%

ggplot() +
    geom_histogram(
      data = s2_tcs,
      aes(x = relative_held),
      binwidth = 0.05) + 
    scale_x_continuous(labels = scales::percent) + 
    labs(
      title = paste0("Feature representation by prioritization"),
      x = "Percent coverage of features",
      y = "Number of features") +
    theme_minimal() + 
  facet_wrap(~budget, nrow = 3)

Scenario 3

Future Land Use Expansion + Climate Change Species Targets

Here we lock out current land use for agriculture and urban areas and the results from Johnson et al., for areas of potential land use conversion for agriculture and renewable energy in a 2050 BAU scenario, and we use current and future species distributions. Renewable energy is a sum of land use conversion for concentrated solar power, photovoltaic solar power, wind power, and hydropower.

Constraints

Locked out

em_rast <- rast(paste0("methodsTrial_data/interim/", roi_abb, "_energy_matters.tif"))

em_ene <- sum(em_rast$CSP,
              em_rast$Hydro,
              em_rast$PV,
              em_rast$Wind)
names(em_ene) <- "ene"
em_ene[em_ene > 0] <- 1
em_en_mask <- em_ene
em_en_mask[em_en_mask == 0] <- NA

em_ag <- em_rast$Crops
names(em_ag) <- "ag"
em_ag[em_ag > 0] <- 1
em_ag[ag_vals] <- 0
em_ag_mask <- em_ag
em_ag_mask[em_ag_mask == 0] <- NA

em_en_cells <- cells(em_en_mask)
em_ag_cells <- cells(em_ag_mask)

plot(em_ag, 
     col = bin_cols,
     main = "Agriculture locked out",
     range = c(0, 1),
     axes = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

plot(em_ene, 
     col = bin_cols,
     main = "Energy locked out",
     range = c(0, 1),
     axes = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

p3_locked_out <- unique(c(em_en_cells,
                          em_ag_cells))
p3_locked_out <- p3_locked_out[which(p3_locked_out %in% c(ids$id))]
p3_locked_out <- which(ids$id %in% p3_locked_out)
p3_locked_out <- unique(c(cur_ag_urban, p3_locked_out))

Solve

for (i in seq_along(budgets)){
  
  rel_budget <- budgets[i] # % of total land 
  
  p <- my_problem(
    x = land_cost,
    features = features_df,
    rij_matrix = features_mat,
    budget = sum(land_cost) * rel_budget,
    targets = all_targets,
    locked_in = wdpa,
    locked_out = p3_locked_out)
  
  s <- solve(p)
  
  s_rast <- rasterize_soln(s)
  
  assign(paste0("s3_rast_", rel_budget),
         s_rast)
  
  tc_s <- eval_target_coverage_summary(p, s) %>% 
    mutate('budget' = rel_budget)

  assign(paste0("s3_tc_", rel_budget),
         tc_s)
  
    assign(paste0("s3_met_", rel_budget),
         tc_s$met)
}

## Conservation Problem
##   planning units: matrix (17271 units)
##   cost:           min: 25, max: 25
##   features:       Afghanodon mustersi current, Ambystoma californiense current, Ambystoma gracile current, ... (1864 features)
##   objective:      Minimum shortfall objective [budget (215887.5)]
##   targets:        Relative targets [targets (min: 0.00991914163715929, max: 1)]
##   decisions:      default
##   constraints:    <Locked out planning units [4134 locked units]
##                    Locked in planning units [3766 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (1e-05), zones]>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.05), node_file_start (-1), numeric_focus (0), presolve (2), threads (8), time_limit (172800), verbose (1)]

s3_rast <- sum(s3_rast_0.3, s3_rast_0.4, s3_rast_0.5, na.rm = TRUE)
plot_ranked(s3_rast,
            main = "Scenario 3")

s3_tcs <- bind_rows(s3_tc_0.3,
                    s3_tc_0.4,
                    s3_tc_0.5) %>% 
  mutate(budget = paste0(scales::percent(budget), " budget"))

s3_met <- s3_tcs %>% 
  group_by(budget) %>% 
  summarize(met = sum(as.numeric(met))) %>% 
  ungroup() %>% 
  mutate(pct_met = met/length(all_targets)) %>% 
  mutate(met = scales::comma(met),
         pct_met = scales::percent(pct_met))

kable(s3_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,749	93.83%
40% budget	1,796	96.35%
50% budget	1,814	97.32%

ggplot() +
    geom_histogram(
      data = s3_tcs,
      aes(x = relative_held),
      binwidth = 0.05) + 
    scale_x_continuous(labels = scales::percent) + 
    labs(
      title = paste0("Feature representation by prioritization"),
      x = "Percent coverage of features",
      y = "Number of features") +
    theme_minimal() + 
  facet_wrap(~budget, nrow = 3)

Scenario 4

Future Land Use Expansion that Avoids Biodiversity + Climate Change Species Targets

We first identify potential areas for land use conversion using conservation areas as a cost layer, Oakleaf et al., 2019 DPI as a features layer, and the summed overlap of Johnson et al. land use conversion with DPI as a target. We then lock out current land use for agriculture and urban areas and our new priorities for agriculture and energy expansion, and include targets for current and future species distributions.

Prioritize for ag + energy

When prioritizing units for agricultural and energy production, we modeled the ‘minimum set problem’ wherein the goal is to minimize the cost of the solution while ensuring that targets are met. This objective is similar to that used in Marxan (cite). We used the management zones functionality in prioritizr to allow for priority setting for agriculture and energy concurrently while using separate cost layers and targets. For the cost layer for energy, we used a combination of conservation priority areas and agriculture development potential, and for agriculture we used conservation priority areas and agriculture development potential. For priority areas for conservation, we used global conserved area expansion priorities for protecting vertebrate species under climate change from Roehrdanz et al (preprint).

Features

We will use the development potential index for cropland from Oakleaf et al., 2019.

dpi_rast <- rast(
  paste0("methodsTrial_data/interim/", roi_abb, "_dpi.tif"))

ag_dpi <- dpi_rast$crop/6 # 0-1
ag_features <- ag_dpi
plot(sa, col = "grey90", axes = FALSE,
     main = "Agriculture features", 
     legend = FALSE,
     mar = c(3.1, 3.1, 2.1, 7.1))
plot(ag_features, range = c(0, 1),
     col = blue_cols(40), add = TRUE,
     axes = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

# all renewables for now - in future may separate by sector
en_dpi <- sum(dpi_rast$csp,
              dpi_rast$hydro,
              dpi_rast$pv,
              dpi_rast$wind,
              na.rm = TRUE)/6 #0-1
names(en_dpi) <- "renewables"
en_dpi <- en_dpi/4 #0-1
en_features <- en_dpi
plot(sa, col = "grey90", axes = FALSE,
     main = "Renewable Energy Features", 
     legend = FALSE,
     mar = c(3.1, 3.1, 2.1, 7.1))
plot(en_features, range = c(0, 1),
     col = blue_cols(40), add = TRUE,
     axes = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

Targets

For now, we will use the sum of the dpi covered in the Johnson et al. land use conversion dataset as our targets. In the future, we aim to use yield (i.e., calories, megawatts) targets.

# targets
# will base on weighted sum of dpi overlapped by Johnson et al., data

em_ag_vals <- values(em_ag)
em_ag_vals <- em_ag_vals[em_ag_cells]
em_en_vals <- values(em_ene)
em_en_vals <- em_en_vals[em_en_cells]

ag_target <- sum((ag_dpi[em_ag_cells] * em_ag_vals), na.rm = TRUE) +
  sum(ag_dpi[ag_vals], na.rm = TRUE)

en_target <- sum((en_dpi[em_en_cells] * em_en_vals), na.rm = TRUE)

# create a matrix with the targets
# here each column corresponds to a different zone,
# each row corresponds to a different feature, and
# each cell value corresponds to the target
land_targets <- matrix(NA, ncol = 2, nrow = 1)
land_targets[, 1] <- ag_target
land_targets[, 2] <- en_target

Cost

For cost, we will use a layer for conservation priorities, as well as the DPI from the other sector (i.e., for agriculture we will use energy and vice versa). For conservation priorities, we will use global conserved area expansion priorities for protecting vertebrate species under climate change from Roehrdanz et al (preprint).

cons_priorities <- rast(paste0(
    "methodsTrial_data/interim/", 
    roi_abb, "_SPARCpriorities.tif")) ^ 2
# ^2 per PR note to better match extinction risk

plot(cons_priorities,
     col = blue_cols(40),
     axes = FALSE,
     range = c(0,1),
     main = "Conservation Priorities")
plot(shape, add = TRUE, lwd = 0.5)

ag_cost <- sum(en_dpi, cons_priorities, na.rm = TRUE)
ag_cost[ag_vals] <- 0
ag_cost <- raster(ag_cost)
names(ag_cost) <- "ag_cost"
en_cost <- sum(ag_dpi, cons_priorities, na.rm = TRUE) %>% raster()
names(en_cost) <- "en_cost"

land_costs <- stack(ag_cost, en_cost)

plot(rast(land_costs$ag_cost),
     col = blue_cols(40),
     axes = FALSE,
     box = FALSE,
     range = c(0, 2),
     main = "Costs for agriculture priorities")
plot(shape, add = TRUE, lwd = 0.5)

plot(rast(land_costs$en_cost),
     col = blue_cols(40),
     axes = FALSE,
     box = FALSE,
    range = c(0, 2),
     main = "Costs for energy priorities")
plot(shape, add = TRUE, lwd = 0.5)

Constraints

# locked in
# lock in ag for ag
no_en_lock <- ag
no_en_lock[] <- 0

p4a_locked_in <- stack(raster(ag), raster(no_en_lock))


# locked out
# lock out wdpa for both
wdpa_rast <- template
wdpa_rast[] <- 0
wdpa_rast[cells] <- wdpa_bin
wdpa_raster <- wdpa_rast %>% 
  prepare_raster(resample = FALSE) %>% 
  raster()

p4a_locked_out <- stack(wdpa_raster, wdpa_raster)

Solve

p4a <- problem(
  x = land_costs,
  features = zones(
    agriculture = raster(ag_features), 
    energy = raster(en_features))) %>% 
  add_min_set_objective() %>% 
  add_absolute_targets(land_targets) %>% 
  add_locked_in_constraints(p4a_locked_in) %>% 
  add_locked_out_constraints(p4a_locked_out) %>% 
  add_boundary_penalties(penalty = 0.000025)

s4a <- solve(p4a)
s4a[[1]][ag_vals] <- 0.50 # values that are already ag

p4a

## Conservation Problem
##   zones:          agriculture, energy (2 zones)
##   planning units: RasterStack (17063 units)
##   cost:           min: 0, max: 1.9996
##   features:       1 (1 features)
##   objective:      Minimum set objective 
##   targets:        Absolute targets [targets (min: 397.084897468487, max: 2184.57194999854)]
##   decisions:      default
##   constraints:    <Manually locked planning units [1740 locked units]
##                    Manually locked planning units [8248 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (2.5e-05), zones]>
##   portfolio:      default
##   solver:         default

plot(rast(s4a[[1]]), 
     main = "Agriculture priorities",
     col = c(bin_cols[1], blue_cols(3)[2], bin_cols[2]),
     axes = FALSE,
     box = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

plot(rast(s4a[[2]]), 
     main = "Energy priorities",
     col = bin_cols,
     axes = FALSE,
     box = FALSE)
plot(shape, add = TRUE, lwd = 0.5)

Light blue is already locked in (current land use), whereas dark blue are additional priorities.

Prioritize for conservation

Locked out

ag_pri <- rast(s4a[[1]])
ag_pri[ag_pri < 1] <- 0
ag_pri_mask <- ag_pri
ag_pri_mask[ag_pri_mask < 1] <- NA
ag_pri_vals <- cells(ag_pri_mask)

energy_pri <- rast(s4a[[2]])
energy_mask <- energy_pri
energy_mask[energy_mask < 1] <- NA
energy_vals <- cells(energy_mask)

p4_locked_out <- unique(c(ag_pri_vals, energy_vals))
p4_locked_out <- p4_locked_out[which(p4_locked_out %in% c(ids$id))]

# convert to indices which I don't love...
p4_locked_out <- which(ids$id %in% p4_locked_out)
p4_locked_out <- unique(c(cur_ag_urban, p4_locked_out)) # add current ag + urban

Solve

 for (i in seq_along(budgets)){
   
   rel_budget <- budgets[i] # % of total land 
   
   p <- my_problem(
     x = land_cost,
     features = features_df,
     rij_matrix = features_mat,
     budget = sum(land_cost) * rel_budget,
     targets = all_targets,
     locked_in = wdpa,
     locked_out = p4_locked_out)
   
   s <- solve(p)
   
   s_rast <- rasterize_soln(s)
   
   assign(paste0("s4_rast_", rel_budget),
          s_rast)
   
   tc_s <- eval_target_coverage_summary(p, s) %>% 
     mutate('budget' = rel_budget)
   
   assign(paste0("s4_tc_", rel_budget),
          tc_s)
   
   assign(paste0("s4_met_", rel_budget),
          tc_s$met)
 }

## Conservation Problem
##   planning units: matrix (17271 units)
##   cost:           min: 25, max: 25
##   features:       Afghanodon mustersi current, Ambystoma californiense current, Ambystoma gracile current, ... (1864 features)
##   objective:      Minimum shortfall objective [budget (215887.5)]
##   targets:        Relative targets [targets (min: 0.00991914163715929, max: 1)]
##   decisions:      default
##   constraints:    <Locked out planning units [3891 locked units]
##                    Locked in planning units [3766 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (1e-05), zones]>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.05), node_file_start (-1), numeric_focus (0), presolve (2), threads (8), time_limit (172800), verbose (1)]

s4_rast <- sum(s4_rast_0.3, s4_rast_0.4, s4_rast_0.5, na.rm = TRUE)
plot_ranked(s4_rast,
            main = "Scenario 4")

s4_tcs <- bind_rows(s4_tc_0.3,
                    s4_tc_0.4,
                    s4_tc_0.5) %>% 
  mutate(budget = paste0(scales::percent(budget), " budget"))

s4_met <- s4_tcs %>% 
  group_by(budget) %>% 
  summarize(met = sum(as.numeric(met))) %>% 
  ungroup() %>% 
  mutate(pct_met = met/length(all_targets)) %>% 
  mutate(met = scales::comma(met),
         pct_met = scales::percent(pct_met))

kable(s4_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,767	94.80%
40% budget	1,811	97.16%
50% budget	1,825	97.91%

ggplot() +
    geom_histogram(
      data = s4_tcs,
      aes(x = relative_held),
      binwidth = 0.05) + 
    scale_x_continuous(labels = scales::percent) + 
    labs(
      title = paste0("Feature representation by prioritization"),
      x = "Percent coverage of features",
      y = "Number of features") +
    theme_minimal() + 
  facet_wrap(~budget, nrow = 3)

plot_ag <- rast(s4a$agriculture)
plot_ag[plot_ag == 0.5] <- 40
plot_ag[plot_ag == 1] <- 50
plot_en <- rast(s4a$energy * 60)
plot_cons <- s4_rast
plot_urban <- urban * 70
plot_s4 <- sum(plot_cons, plot_ag, plot_en, plot_urban, na.rm = TRUE)

levels(plot_s4) <- NULL
levels(plot_s4) <- c(
  rep(" ", 3),
  "Current protected",
  rep(" ", 6),
  "Conservation priorities: 50% budget",
  rep(" ", 9),
  "Conservation priorities: 40% budget",
  rep(" ", 9),
  "Conservation priorities: 30% budget",
  rep(" ", 9),
  "Current agriculture",
  rep(" ", 9),
  "Agriculture priorities",
  rep(" ", 9),
  "Renewable energy priorities",
  rep(" ", 9),
  "Current urban")
plot(plot_s4,
     col = c("white", "#BEC4C8",
             blue_cols(25)[15],
             blue_cols(25)[20],
             blue_cols(25)[25],
             "#CEA69C",
             "#D4340C",
             "#5D8731",
             "#8D8DA7"),
     axes = FALSE,
     main = "Scenario 4 priorities for conservation, agriculture, and energy",
     mar = c(3.1, 3.1, 2.1, 14))
plot(shape, add = TRUE, lwd = 0.5)

# ## this is for plotting only 30p budget
# 
# plot_ag <- rast(s4a$agriculture)
# plot_ag[plot_ag == 0.5] <- 40
# plot_ag[plot_ag == 1] <- 50
# plot_en <- rast(s4a$energy * 60)
# plot_cons <- s4_rast * 3
# plot_urban <- urban * 70
# plot_s4 <- sum(plot_cons, plot_ag, plot_en, plot_urban, na.rm = TRUE)
# 
# levels(plot_s4) <- NULL
# levels(plot_s4) <- c(
#   rep(" ", 3),
#   "Current protected", # 3
#   rep(" ", 26),
#   # "Conservation priorities: 50% budget",
#   # rep(" ", 9),
#   # "Conservation priorities: 40% budget",
#   # rep(" ", 9),
#   "Conservation priorities: 30% budget", # 30
#   rep(" ", 9),
#   "Current agriculture",
#   rep(" ", 9),
#   "Agriculture priorities",
#   rep(" ", 9),
#   "Renewable energy priorities",
#   rep(" ", 9),
#   "Current urban")
# 
# png(paste0(roi_abb, "_priorities.png"),
#     width = 1000,
#     height = 600,
#     type = "cairo",
#     family = "ArialMT", 
#     pointsize = 18,)
# plot(plot_s4,
#      col = c("white", 
#              "#cccccc",
#              # blue_cols(25)[15],
#              # blue_cols(25)[20],
#              # blue_cols(25)[25],
#              "#5D8731",
#              "#F6A38E",
#              "#EE562F",
#              "#2C6CB5",
#              "#8D8DA7"),
#      axes = FALSE,
#      main = "Scenario 4 priorities for conservation, agriculture, and energy",
#      mar = c(3.1, 3.1, 2.1, 14))
# plot(shape, add = TRUE, lwd = 0.5)
# dev.off()

Summary

No Future Climate Change or LUC

plot_ranked(s1_rast,
            main = "Scenario 1")

kable(s1_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,712	91.8%
40% budget	1,774	95.2%
50% budget	1,810	97.1%

No Future LUC

plot_ranked(s2_rast,
            main = "Scenario 2")

kable(s2_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,761	94.5%
40% budget	1,807	96.9%
50% budget	1,829	98.1%

Economic-Driven LUC

plot_ranked(s3_rast,
            main = "Scenario 3")

kable(s3_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,749	93.83%
40% budget	1,796	96.35%
50% budget	1,814	97.32%

Economic + Biodiversity-Driven LUC

plot_ranked(s4_rast,
            main = "Scenario 4")

kable(s4_met,
      col.names = c("Budget",
                    "Targets met (n)",
                    "Targets met (%)")) %>% 
  kable_styling()

Budget	Targets met (n)	Targets met (%)
30% budget	1,767	94.80%
40% budget	1,811	97.16%
50% budget	1,825	97.91%

for (i in 1:4){
  met <- get(paste0("s", i, "_met")) %>% 
    mutate(scenario = paste0("Scenario ", i))
  assign(paste0("comp_", i),
         met)
}

compare <- bind_rows(comp_1, comp_2, comp_3, comp_4) %>% 
  mutate(met = as.numeric(str_remove(met, ",")),
         pct_met = as.numeric(str_remove(pct_met, "%"))/100)


ggplot(data = compare,
       aes(y = fct_rev(budget),
        x = pct_met)) + 
  geom_col(
    aes(fill = scenario,
        alpha = budget),
    position = "dodge",
    width = 0.5,
    show.legend = FALSE) +
  geom_text(
    aes(group = budget, 
        label = scales::percent(round(pct_met, 3))),
    position = position_dodge(width = 0.5),
    size = 4,
    hjust = 1) + 
  scale_alpha_discrete(range = c(0.5, 1)) + 
  scale_x_continuous(labels = scales::percent) + 
  scale_fill_manual(values = cal_palette("sierra1")) + 
  theme_minimal() + 
  coord_polar(clip = 'off') + 
  labs(x = NULL,
       y = NULL) + 
  facet_wrap(~scenario, nrow = 1)

Works cited

Chen, G., Li, X., & Liu, X. (2022). Global land projection based on plant functional types with a 1-km resolution under socio-climatic scenarios. Scientific Data, 9(1), 125. https://doi.org/10.1038/s41597-022-01208-6

Hanson, J. O., Rhodes, J. R., Butchart, S. H. M., Buchanan, G. M., Rondinini, C., Ficetola, G. F., & Fuller, R. A. (2020). Global conservation of species’ niches. Nature, 580(7802), 232–234. https://doi.org/10.1038/s41586-020-2138-7

Hanson JO, Schuster R, Morrell N, Strimas-Mackey M, Edwards BPM, Watts ME, Arcese P, Bennett J, Possingham HP (2022). prioritizr: Systematic Conservation Prioritization in R. R package version 7.1.1. Available at https://CRAN.R-project.org/package=prioritizr.

Johnson, J. A., Kennedy, C. M., Oakleaf, J. R., Baruch-Mordo, S., Polasky, S., & Kiesecker, J. (2021). Energy matters: Mitigating the impacts of future land expansion will require managing energy and extractive footprints. Ecological Economics, 187, 107106. https://doi.org/10.1016/j.ecolecon.2021.107106

Oakleaf, J. R., Kennedy, C. M., Baruch-Mordo, S., Gerber, J. S., West, P. C., Johnson, J. A., & Kiesecker, J. (2019). Mapping global development potential for renewable energy, fossil fuels, mining and agriculture sectors. Scientific Data, 6(1), 101. https://doi.org/10.1038/s41597-019-0084-8

Rodrigues, A. S. L., Andelman, S. J., Bakarr, M. I., Boitani, L., Brooks, T. M., Cowling, R. M., Fishpool, L. D. C., da Fonseca, G. A. B., Gaston, K. J., Hoffmann, M., Long, J. S., Marquet, P. A., Pilgrim, J. D., Pressey, R. L., Schipper, J., Sechrest, W., Stuart, S. N., Underhill, L. G., Waller, R. W., … Yan, X. (2004). Effectiveness of the global protected area network in representing species diversity. Nature, 428(6983), 640–643. https://doi.org/10.1038/nature02422

Roehrdanz, P., Hannah, L., Corcoran, D., Corlett, R., Enquist, B., Fajardo, J., Feng, X., Foden, W., Lovett, J., Maitner, B., Marquet, P., Merow, C., & Midgley, G. (preprint). Strategic Conservation of Global Vertebrates in Response to Climate Change. Available at SSRN: https://doi.org/10.2139/ssrn.3854499

UNEP-WCMC and IUCN (2021), Protected Planet: The World Database on Protected Areas (WDPA) [On-line], Cambridge, UK: UNEP-WCMC and IUCN. Available at: www.protectedplanet.net

SPARCLE Proposed Methods

Spatial Planning for Climate Change: Land Use for Conservation, Agriculture, and Energy

Cameryn Brock

6/10/2022

Set up

Study area

Cost

Conservation features

Current

Future

Targets

Constraints

Locked in

Locked-out

Scenarios

Scenario 1

Solve

Scenario 2

Solve

Scenario 3

Constraints

Locked out

Solve

Scenario 4

Prioritize for ag + energy

Features

Targets

Cost

Constraints

Solve

Prioritize for conservation

Locked out

Solve

Summary

No Future Climate Change or LUC

No Future LUC

Economic-Driven LUC

Economic + Biodiversity-Driven LUC

Works cited