In ecoquants/bbnj: Biodiversity Beyond National Jurisdiction

knitr::opts_chunk$set(echo=T, message=F)
library(prioritizr) # devtools::load_all("~/github/prioritizr")
#library(bbnj)       # 
devtools::load_all() # setwd(here()); devtools::install_local(force=T) 
# devtools::install_github("ecoquants/bbnj")
library(raster)
library(sf)
library(dplyr)
library(readr)
library(stringr)
library(glue)
library(here)
library(fs)
library(knitr)
library(formattable)
area   = raster::area
select = dplyr::select

#message("library done")

if (interactive()){
  wd <- file.path(here::here(), "inst/app/www/scenarios/")
  setwd(wd) # getwd()
  if (!exists("rmd")){
    #rmd <- file.path(getwd(), "s00a.bio.30pct.gl.mol50km.Rmd")
    rmd <- file.choose()
  }
} else {
  rmd <- knitr::current_input(dir = T)  
}

pfx <- rmd %>% path_ext_remove()
tif <- glue("{pfx}_sol.tif")

# variables ----
rel_target <- 0.3
prjres     <- "_mol50km" # prjres in: View(projections_tbl)
redo       <- T

# problem & solution ----
P <- projections_tbl %>% filter(prjres == !!prjres)
if (!file.exists(tif) | redo){

  # planning unit: ----
  r_pu_id <- get_d_prjres("r_pu_id", prjres) # plot(r_pu_id)
  r_pu <- setValues(r_pu_id, 1) %>% 
    mask(r_pu_id) # plot(r_pu)

  # biodiversity: now, not 2100 ----
  # devtools::load_all()
  #s_bio_gmbi <- get_d_prjres("s_bio_gmbi", prjres)
  #groups01: 1st taxonomic grouping
  s_bio_gmbi <- get_gmbi_grpsmdl_prjres("groups01", "_mol50km")
  lyrs_bio_now <- names(s_bio_gmbi) %>% 
    setdiff(str_subset(names(s_bio_gmbi), "rli")) %>% 
    setdiff(str_subset(names(s_bio_gmbi), "rls"))
  s_bio_now <- subset(s_bio_gmbi, lyrs_bio_now)

  #add red list sum for all species as separate layer
  rls_all <- get_gmbi_grpsmdl_prjres("groups00", prjres) %>%
    subset("groups00_rls_all")

  # features ----
  s_features <- stack(
    get_d_prjres("r_vgpm", prjres),
    s_bio_now,
    rls_all,
    #raster(s_fish_gfw, "mean_scaled_profits_with_subsidies") %>%
    #    gap_fill_raster() %>%
    #    rescale_raster(inverse=T),
    #raster(s_fish_ubc, "mcp_2004"),
    get_d_prjres("s_phys_seamounts",prjres),
    get_d_prjres("r_phys_vents",prjres),
    get_d_prjres("r_phys_scapes_hetero",prjres))

  names(s_features) <- c(
    "bio_vgpm",
    names(s_bio_now),
    "rls_all",
    #"fish_profit.subs"
    #"fish_mcp.2004",
    sprintf("phys_seamounts_%sm", c("0to200","gt200to800","gt800")),
    "phys_vents",
    "scapes_hetero")

  # problem ----
  # r_pu_0 <- r_pu
  #plot(r_pu)
  #r_pu <- setValues(r_pu, 0.0000000000000000001) %>% 
  r_pu <- setValues(r_pu, 1) %>% 
    mask(r_pu)

  p <- problem(r_pu, s_features) %>%
    add_min_set_objective() %>% 
    add_relative_targets(rel_target) %>% 
    add_boundary_penalties(penalty = 0.0000000000000000001)
  #                         works: = 0.00000000000000000001

  print(presolve_check(p))
  # range(r_pu): [0.0000000000000000001]
  #   [1] TRUE
  # range(r_pu): [0.0001,0.0001]
  #   Warning in presolve_check.OptimizationProblem(compile(x)) :
  #   planning units with (relatively) very high costs, note this may be a false positive
  # range(r_pu): [1,1]
  #   Warning in presolve_check.OptimizationProblem(compile(x)) :
  #   planning units with (relatively) very high costs, note this may be a false positive
  # [1] FALSE

  #presolve_check(p)
  # range(r_pu): [0.0000000000000000001]
  #   [1] TRUE
  # range(r_pu): [0.0001,0.0001]
  #   Warning in presolve_check.OptimizationProblem(compile(x)) :
  #   planning units with (relatively) very high costs, note this may be a false positive
  # range(r_pu): [1,1]
  #   Warning in presolve_check.OptimizationProblem(compile(x)) :
  #   planning units with (relatively) very high costs, note this may be a false positive
  # [1] FALSE


    # solve ----
  # redo=T
  tif <- solve_log(p, pfx, redo=redo, run_checks=F)
}

# debugging BLM w/ MV...
r_pu <- setValues(r_pu_id, 1) %>% 
  mask(r_pu_id)
r_sol <- raster(tif)
r1 <- r_pu
r2 <- r_sol
r_d  <- r2*10 - r1
d_subs <- tribble(
  ~old,   ~new,
  0, -9999,
  9,     0,
  -1,    -1,
  10,     1)
r_ds   <- subs(r_d, d_subs, "old", "new")
r_diff <- mask(r_ds, r_ds, maskvalue=-9999)

# why are all these 0?
# groups01_rls_crustaceans
# groups01_rls_euphausiids
# groups01_rls_gastropods
# groups01_rls_hydrozoans
# groups01_rls_seagrasses
# AH, do not want rls*

r <- raster(s_bio_now, "groups01_rls_seagrasses")
r

a_km2 <- prod(res(r_diff)) / (1000 * 1000)
table(values(r_diff))

tbl <- tibble(
  value = na.omit(values(rd))) %>%
  group_by(value) %>%
  summarize(
    ncells = n()) %>%
  mutate(
    label = recode(
      value,
      `-1` = "loss",
      `0`  = "same",
      `1`  = "gain"),
    area_km2 = ncells * a_km2,
    pct = area_km2 / sum(area_km2)) %>%
  dplyr::select(value, label, ncells, area_km2, pct)

map_r2png(rd, png)

`r basename(pfx)`

Projection: r P$name (r P$res resolution)
Planning unit cost: set to 1
Objective function:
- add_min_set_objective(): Minimize the cost of the solution whilst ensuring that all targets are met
Feature targets:
- r rel_target*100% of each input feature amount
- Biodiversity:
  - bio_vgpm: Vertically Generalized Production Model
  - nspp_*: number of species by 22 taxonomic group from AquaMaps [now]
  - rls_all: Red List sum of extinction risk from AquaMaps for all species [now]
- Physical:
  - phys_vents: hydrothermal vent count
  - phys_seamounts: seamounts count (3 depth classes)
  - phys_scapes_hetero: benthic heterogeneity
Regionalization:
- none, ie global
Boundary Penalty
- boundary penalty = 0.00000000000000000001