In ecoquants/bbnj: Biodiversity Beyond National Jurisdiction

knitr::opts_chunk$set(echo=T, message=F, warning = 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)
library(lwgeom)
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
redo_map   <- 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)

  # countries <- rnaturalearth::ne_countries(returnclass = "sf") %>%
  #     st_transform(P$epsg)
  # graticules <- st_graticule(countries)
  # 
  # plot(r_pu_id, legend=F, axes=F, box=F)
  # plot(st_geometry(countries), add=T, col=gray(0.8), border=gray(0.7), lwd=0.5)
  # plot(st_geometry(graticules), add=T, col=gray(0.6), lwd=0.5)

  #mapview::mapview(r_pu_id_pr)

  # plot(r_pu_id)

  # r_pu <- setValues(r_pu_id, 1) %>%
  #   mask(r_pu_id) # plot(r_pu)

  # cost layer: fishing profitability ----
  r_fish_effort <- get_d_prjres("s_fish_gfw", prjres)%>%
    subset("fishing_KWH")

  # apply cost to planning units
  r_pu <- r_fish_effort %>% 
    reclassify(rcl=cbind(-Inf, 112774, NA), right=TRUE) %>% 
    gap_fill_raster(r_mask=r_pu_id)

  # biodiversity: now + 2100 ----
  s_bio_gmbi_now <- get_gmbi_grpsmdl_prjres("groups02", prjres)
  lyrs_bio_now <- names(s_bio_gmbi_now) %>% 
    setdiff(str_subset(names(s_bio_gmbi_now), "rli")) %>% 
    setdiff(str_subset(names(s_bio_gmbi_now), "rls"))
  s_bio_gmbi_now <- subset(s_bio_gmbi_now, lyrs_bio_now)

  s_bio_gmbi_future <- get_gmbi_grpsmdl_prjres("groups02_2100", prjres)
  lyrs_bio_future <- names(s_bio_gmbi_future) %>% 
    setdiff(str_subset(names(s_bio_gmbi_future), "rli")) %>% 
    setdiff(str_subset(names(s_bio_gmbi_future), "rls"))
  s_bio_gmbi_future <- subset(s_bio_gmbi_future, lyrs_bio_future)

  # rls now + 2100
  rls_all_now <- get_gmbi_grpsmdl_prjres("groups00", prjres) %>% 
    subset("groups00_rls_all")

  rls_all_future <- get_gmbi_grpsmdl_prjres("groups00_2100", prjres) %>% 
    subset("groups00_2100_rls_all")

  # features ----
  s_seamounts <- get_d_prjres("s_phys_seamounts",prjres)
  lu_seamounts <- c(lteq200m="0to200",gt200lteq800m="gt200to800",gt800m="gt800")
  lbls_seamounts <- sprintf("phys_seamounts_%sm", lu_seamounts[names(s_seamounts)])

  s_features <- stack(
    get_d_prjres("r_vgpm", prjres),
    s_bio_gmbi_now,
    s_bio_gmbi_future,
    rls_all_now,
    rls_all_future,
    #raster(s_fish_gfw, "mean_scaled_profits_with_subsidies") %>%
    #    gap_fill_raster(r_mask=r_pu_id) %>%
    #    rescale_raster(inverse=T),
    #raster(s_fish_ubc, "mcp_2004"),
    s_seamounts,
    get_d_prjres("r_phys_vents",prjres),
    get_d_prjres("r_phys_scapes_hetero",prjres))

  names(s_features) <- c(
    "bio_vgpm",
    gsub("^.*?_","",names(s_bio_gmbi_now)),
    gsub("^.*?_","",names(s_bio_gmbi_future)),
    "rls_all_now",
    "rls_all_future",
    #"fish_profit.subs"
    #"fish_mcp.2004",
    lbls_seamounts,
    "phys_vents",
    "scapes_hetero")

  tibble(

  )

  lyrs_csv <- glue("{fs::path_ext_remove(rmd)}_layers.csv")
  tibble(
    layer = names(s_features)) %>% 
    write.csv(lyrs_csv)

  # problem ----
  p <- problem(r_pu, s_features) %>%
    add_min_set_objective() %>% 
    add_relative_targets(rel_target)


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

`r basename(pfx)`

Projection: r P$name (r P$res resolution)
Planning unit cost:
- top quartile of fishing effort (KWH) (i.e. include all cells with fishing effort > 112,774 KWH)
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 23 taxonomic group from AquaMaps [now + future]
  - rls_all: Red List sum of extinction risk from AquaMaps for all species [now + future]
- Physical:
  - phys_vents: hydrothermal vent count
  - phys_seamounts: seamounts count (3 depth classes)
  - phys_scapes_hetero: benthic heterogeneity

# devtools::load_all()
report_solution(tif, redo=redo_map)

library(raster)
library(rasterVis)
library(here)

tif <- here("inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol.tif")
pdf <- here("inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol_levelplot.pdf")
r <- raster(tif)

#table(values(r))
#plot(r)
pdf(pdf)
rasterVis::levelplot(r, margin = list(FUN = 'mean'))
dev.off()