0.1 Gurobi setup

grbgetkey `cat ~/private/gurobi_key`
  • saved license to /Users/bbest/gurobi.lic

Testing: Gurobi Optimizer Quick Start Guide - Mac OS: Software Installation Guide: Solving a Simple Model - The Gurobi Command Line: Solving the model using the Gurobi command-line interface

gurobi_cl /Library/gurobi810/mac64/examples/data/coins.lp

Install R interface, per Solving the model using the Gurobi command-line interface: Gurobi Optimizer Example Tour: Gurobi Optimizer Quick Start Guide - Mac OS: R Interface: Installing the R Package:

tgz <- list.files('/Library/gurobi810/mac64/R', '.*tgz$', full.names = T)
install.packages(tgz, repos=NULL)

install.packages("slam", repos = "https://cloud.r-project.org")

Try out gurobi in R:

suppressPackageStartupMessages({
  library('gurobi')
})

model <- list()

model$A          <- matrix(c(1,2,3,1,1,0), nrow=2, ncol=3, byrow=T)
model$obj        <- c(1,1,2)
model$modelsense <- 'max'
model$rhs        <- c(4,1)
model$sense      <- c('<', '>')
model$vtype      <- 'B'

params <- list(OutputFlag=0)

#result <- gurobi(model, params)
result <- gurobi(model, list())
## Optimize a model with 2 rows, 3 columns and 5 nonzeros
## Variable types: 0 continuous, 3 integer (3 binary)
## Coefficient statistics:
##   Matrix range     [1e+00, 3e+00]
##   Objective range  [1e+00, 2e+00]
##   Bounds range     [0e+00, 0e+00]
##   RHS range        [1e+00, 4e+00]
## Found heuristic solution: objective 2.0000000
## Presolve removed 2 rows and 3 columns
## Presolve time: 0.00s
## Presolve: All rows and columns removed
## 
## Explored 0 nodes (0 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 12 available processors)
## 
## Solution count 2: 3 2 
## 
## Optimal solution found (tolerance 1.00e-04)
## Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0000%
print('Solution:')
## [1] "Solution:"
print(result$objval)
## [1] 3
print(result$x)
## [1] 1 0 1

0.2 Prioritizr

Solving a prioritzr problem with Gurobi

suppressPackageStartupMessages({
  library(prioritizr) # install.packages("prioritizr")
})

# formulate the problem
p <- problem(sim_pu_raster, sim_features) %>%
     add_min_set_objective() %>%
     add_relative_targets(0.1) %>%
     add_gurobi_solver()

# solve the problem
s <- solve(p)
## Optimize a model with 5 rows, 90 columns and 450 nonzeros
## Variable types: 0 continuous, 90 integer (90 binary)
## Coefficient statistics:
##   Matrix range     [2e-01, 9e-01]
##   Objective range  [2e+02, 2e+02]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [3e+00, 8e+00]
## Found heuristic solution: objective 2337.9617505
## Presolve time: 0.00s
## Presolved: 5 rows, 90 columns, 450 nonzeros
## Variable types: 0 continuous, 90 integer (90 binary)
## Presolved: 5 rows, 90 columns, 450 nonzeros
## 
## 
## Root relaxation: objective 1.931582e+03, 12 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0 1931.58191    0    4 2337.96175 1931.58191  17.4%     -    0s
## H    0     0                    1985.6818841 1931.58191  2.72%     -    0s
## 
## Explored 1 nodes (12 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 12 available processors)
## 
## Solution count 2: 1985.68 2337.96 
## 
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.985681884076e+03, best bound 1.931581908865e+03, gap 2.7245%
# plot solution
plot(s, col = c("grey90", "darkgreen"), main = "Solution",
     xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1))

0.3 Simple BBNJ

suppressPackageStartupMessages({
  library(tidyverse)
  library(glue)
  library(raster)
  library(sf)
  library(leaflet)
  library(RColorBrewer)
  library(rasterVis)
  library(prioritizr) # install.packages("prioritizr")
  
  #library(bbnj)
  devtools::load_all()
})


dir_data <- "~/Gdrive Ecoquants/projects/bbnj/data/derived"
cell_res     <- 0.5  # n=  113,401
p_boundary_shp <- file.path(
  dir_data, glue("boundary/high_seas.shp"))
r_pu_tif <- file.path(
  dir_data, glue("boundary/high_seas_cellid_{cell_res}dd.tif"))

# get planning unit raster
r_pu <- raster(r_pu_tif)

pal <- colorRampPalette(brewer.pal(11, "Spectral"))
cols <- rev(pal(255))
plot(r_pu, col = cols, main="pu")

# get biodiversity features
bio_tifs <- list.files(
  file.path(dir_data, "biodiversity"), 
  "spp_.*_be_0\\.5dd\\.tif", full.names = T)

bio_vars <- tibble(
  tif = bio_tifs,
  var = str_replace(basename(tif), "spp_(.*)_be_0\\.5dd\\.tif", "\\1")) %>% 
  mutate(
    var = recode(var, `tunas.&.billfishes`="tunas.billfishes")) %>% 
  pull(var)

s_bio <- stack(bio_tifs)
names(s_bio) <- bio_vars

#plot(s_bio, col = cols, main="spp")
levelplot(s_bio, main="spp")

# formulate the problem
# TODO: check that s_bio being continuous is ok, and not binary
p <- problem(r_pu, s_bio) %>%
     add_min_set_objective() %>%
     add_relative_targets(0.1) %>%
     add_gurobi_solver()

# solve the problem
r_sol <- solve(p)

# plot solution
plot(r_sol, col = c("grey90", "darkgreen"), main = "Solution")

# devtools::load_all()
qmap_r(r_sol, "Solution", na0=T)

TODO: show pretty stack like

Plot many rasters on the same graph, one on top of the other, using a certain angle in R - Geographic Information Systems Stack Exchange

0.4 prioritizr vignette

prioritizr: Systematic Conservation Prioritization in R • prioritizr

0.4.1 Usage

# load prioritizr package
library(prioritizr)

0.4.2 Data

# load raster planning unit data
data(sim_pu_raster)

# print description of the data
print(sim_pu_raster)
## class       : RasterLayer 
## dimensions  : 10, 10, 100  (nrow, ncol, ncell)
## resolution  : 0.1, 0.1  (x, y)
## extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 190.1328, 215.8638  (min, max)
# plot the data
plot(sim_pu_raster)

# load polygon planning unit data
data(sim_pu_polygons)

# print first six rows of attribute table
head(sim_pu_polygons@data)
##       cost locked_in locked_out
## 1 215.8638     FALSE      FALSE
## 2 212.7823     FALSE      FALSE
## 3 207.4962     FALSE      FALSE
## 4 208.9322     FALSE       TRUE
## 5 214.0419     FALSE      FALSE
## 6 213.7636     FALSE      FALSE
# plot the planning units
spplot(sim_pu_polygons, zcol = "cost")

# specify file path for planning unit data
pu_path <- system.file("extdata/input/pu.dat", package = "prioritizr")

# load in the tabular planning unit data
# note that we use the data.table::fread function, as opposed to the read.csv
# function, because it is much faster
pu_dat <- data.table::fread(pu_path, data.table = FALSE)

# preview first six rows of the tabular planning unit data
# note that it has some extra columns other than id and cost as per the
# Marxan format
head(pu_dat)
##   id       cost status    xloc     yloc
## 1  3        0.0      0 1116623 -4493479
## 2 30   752727.5      3 1110623 -4496943
## 3 56  3734907.5      0 1092623 -4500408
## 4 58  1695902.1      0 1116623 -4500408
## 5 84  3422025.6      0 1098623 -4503872
## 6 85 17890758.4      0 1110623 -4503872
# load feature data
data(sim_features)

# plot the distribution of suitable habitat for each feature
plot(sim_features, main = paste("Feature", seq_len(nlayers(sim_features))),
     nr = 2, box = FALSE, axes = FALSE)

0.4.3 Initialize a problem

  • TODO: check variety of options with Conservation planning problem — problem • prioritizr: “This function is used to specify the basic data used in a spatial prioritization problem: the spatial distribution of the planning units and their costs, as well as the features (e.g. species, ecosystems) that need to be conserved. After constructing this ConservationProblem-class object, it can be customized to meet specific goals using objectives, targets, constraints, and penalties. After building the problem, the solve function can be used to identify solutions.”

0.4.3.1 as raster

# create problem
p1 <- problem(sim_pu_raster, sim_features)

# print problem
print(p1)
## Conservation Problem
##   planning units: RasterLayer (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      none
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default
# print number of planning units
number_of_planning_units(p1)
## [1] 90
# print number of features
number_of_features(p1)
## [1] 5
  • “recommend initializing problems using raster data where possible”

0.4.3.2 as vector

# create problem with spatial vector data
# note that we have to specify which column in the attribute table contains
# the cost data
p2 <- problem(sim_pu_polygons, sim_features, cost_column = "cost")

# print problem
print(p2)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      none
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

0.4.3.3 as tabular

# set file path for feature data
spec_path <- system.file("extdata/input/spec.dat", package = "prioritizr")

# load in feature data
spec_dat <- data.table::fread(spec_path, data.table = FALSE)

# print first six rows of the data
# note that it contains extra columns
head(spec_dat)
##   id prop spf   name
## 1 10  0.3   1  bird1
## 2 11  0.3   1  nvis2
## 3 12  0.3   1  nvis8
## 4 13  0.3   1  nvis9
## 5 14  0.3   1 nvis14
## 6 15  0.3   1 nvis20
# set file path for planning unit vs. feature data
puvspr_path <- system.file("extdata/input/puvspr.dat", package = "prioritizr")

# load in planning unit vs feature data
puvspr_dat <- data.table::fread(puvspr_path, data.table = FALSE)

# print first six rows of the data
head(puvspr_dat)
##   species  pu     amount
## 1      26  56 1203448.84
## 2      26  58  451670.10
## 3      26  84  680473.75
## 4      26  85   97356.24
## 5      26  86   78034.76
## 6      26 111 4783274.17
# create problem
p3 <- problem(pu_dat, spec_dat, cost_column = "cost", rij = puvspr_dat)

# print problem
print(p3)
## Conservation Problem
##   planning units: data.frame (1751 units)
##   cost:           min: 0, max: 41569219.38232
##   features:       bird1, nvis2, nvis8, ... (17 features)
##   objective:      none
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

0.4.4 Add an objective

0.4.4.1 Minimum set objective

  • Minimum set objective: Minimize the cost of the solution whilst ensuring that all targets are met (Rodrigues et al. 2000). This objective is similar to that used in Marxan (Ball et al. 2009). For example, we can add a minimum set objective to a problem using the following code.
# create a new problem that has the minimum set objective
p3 <- problem(sim_pu_raster, sim_features) %>%
      add_min_set_objective()

# print the problem
print(p3)
## Conservation Problem
##   planning units: RasterLayer (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      Minimum set objective 
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

0.4.4.2 Maximum cover objective

  • Maximum cover objective: Represent at least one instance of as many features as possible within a given budget (Church et al. 1996).
# create a new problem that has the maximum coverage objective and a budget
# of 5000
p4 <- problem(sim_pu_raster, sim_features) %>%
      add_max_cover_objective(5000)

# print the problem
print(p4)
## Conservation Problem
##   planning units: RasterLayer (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      Maximum coverage objective [budget (5000)]
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

0.4.4.3 Maximum features objective *

  • Maximum features objective: Fulfill as many targets as possible while ensuring that the cost of the solution does not exceed a budget (inspired by Cabeza & Moilanen 2001). This object is similar to the maximum cover objective except that we have the option of later specifying targets for each feature. In practice, this objective is more useful than the maximum cover objective because features often require a certain amount of area for them to persist and simply capturing a single instance of habitat for each feature is generally unlikely to enhance their long-term persistence.
# create a new problem that has the maximum features objective and a budget
# of 5000
p5 <- problem(sim_pu_raster, sim_features) %>%
      add_max_features_objective(budget = 5000)

# print the problem
print(p5)
## Conservation Problem
##   planning units: RasterLayer (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      Maximum representation objective [budget (5000)]
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

0.4.4.4 Maximum phylogenetic representation objective +

  • Maximum phylogenetic representation objective: Maximize the phylogenetic diversity of the features represented in the solution subject to a budget (inspired by Faith 1992; Rodrigues & Gaston 2002). This objective is similar to the maximum features objective except that emphasis is placed on phylogenetic representation rather than individual features. The prioritizr R package contains a simulated phylogeny that can be used with the simulated feature data (sim_phylogny).
library(ape) # ?plot.phylo
## 
## Attaching package: 'ape'
## The following objects are masked from 'package:raster':
## 
##     rotate, zoom
# load simulated phylogeny data
data(sim_phylogeny)
plot(sim_phylogeny)

# create a new problem that has the maximum phylogenetic representation
# objective and a budget of 5000
p6 <- problem(sim_pu_raster, sim_features) %>%
      add_max_phylo_objective(budget = 5000, tree = sim_phylogeny)

# print the problem
print(p6)
## Conservation Problem
##   planning units: RasterLayer (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      Phylogenetic representation objective [budget (5000)]
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

Could this be considered coarse to fine objectives? Habitats vs Species and within Species actual phylogenetic relationships?

0.4.4.5 Maximum utility objective (vs -)

  • Maximum utility objective: Secure as much of the features as possible without exceeding a budget. This objective is functionally equivalent to selecting the planning units with the greatest amounts of each feature (e.g. species richness). Generally, we don’t encourage the use of this objective because it will only rarely identify complementary solutions—solutions which adequately conserve a range of different features—except perhaps to explore trade-offs or provide a baseline solution with which to compare other solutions.
# create a new problem that has the maximum utility objective and a budget
# of 5000
p7 <- problem(sim_pu_raster, sim_features) %>%
      add_max_utility_objective(budget = 5000)

# print the problem
print(p7)
## Conservation Problem
##   planning units: RasterLayer (90 units)
##   cost:           min: 190.13276, max: 215.86384
##   features:       layer.1, layer.2, layer.3, ... (5 features)
##   objective:      Maximum utility objective [budget (5000)]
##   targets:        none
##   decisions:      default
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         default

0.4.5 Add targets

0.5 BBNJ Implementation Notes

0.5.1 Features

  • fisheries (UBC): projected change in catch potential of exploited species
    • ~ 1,000 species, ½ degree grid, up to 2060, starting with RCP 8.5 scenario aka “business as usual” at year 2050 compared with present
    • incorporates dynamic biogeochemistry and population dynames (larval / adult movement by year)
    • Vicky will send a few example outputs by species in landings (metric tons) for present and future (ie 2050 under RCP 8.5)
  • biodiversity (AquaMaps): projected change in distribution of species in the high seas across taxonomic groups
    • I came up with a web scraping function to extract the habitat parameters from the website (still need to extract FAO area restrictions)
    • I related these habitat parameters to higher resolution layers in the sdmpredictors R package now and into the future, particularly RCP 8.5 2050
    • Next, I will come up with a function to predict distribution by species now and in the future (ie 2050 under RCP 8.5), then difference this change and summarize across taxonomic groups

0.5.2 Zones

  1. surface (fishing)
  2. pelagic (bottom)