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() # devtools::install_local(force=T)
library(raster)
library(sf)
library(dplyr)
library(stringr)
library(glue)
library(here)
library(fs)
library(knitr)
library(formattable)
area   = raster::area
select = dplyr::select

rmd <- knitr::current_input(dir = T)
#rmd <- file.path(getwd(), "s01a.bio.10pct.gl.now.Rmd")
pfx <- rmd %>% path_ext_remove()

r basename(pfx)

• Planning unit cost:

  • fishing profitability, ie minimize impacts to fishing

• Feature targets:

  • 30% of each input feature amount
  • **Biodiversity:
    • bio_vgpm: Vertically Generalized Production Model
    • nspp_*: number of species by taxonomic group from AquaMaps [now]
    • rls_*: Red List sum of extinction risk by taxonomic group from AquaMaps [now]
  • Physical:
    • phys_vents: hydrothermal vent count
    • phys_seamounts: seamounts count
    • phys_scape.*: benthic seascapes (11 types)

• Regionalization:

  • none, ie global

# variables: target 30% ----
rel_target <- 0.3
redo       <- F

# planning unit: inverse of fishing profitability ----
r_fish_profits <- raster(s_fish_gfw, "mean_scaled_profits_with_subsidies") # plot(r_fish_profits)
r_pu <- r_fish_profits %>%
  gap_fill_raster() # %>%
  # rescale_raster(inverse=T) # plot(r_pu)

# biodiversity: now, not 2100 ----
lyrs_bio_now <- names(s_bio_gmbi) %>% 
  setdiff(str_subset(names(s_bio_gmbi), "2100$")) %>% 
  # use taxonomic groups, so exclude all and NA layers
  setdiff(c("nspp_all", "nspp_na", "rls_all", "rls_na"))
s_bio_now <- subset(s_bio_gmbi, lyrs_bio_now)

# features ----
s_features <- stack(
  r_vgpm,
  s_bio_now,
  #raster(s_fish_gfw, "mean_scaled_profits_with_subsidies") %>%
  #    gap_fill_raster() %>%
  #    rescale_raster(inverse=T),
  #raster(s_fish_ubc, "mcp_2004"),
  r_phys_seamounts,
  r_phys_vents,
  s_phys_scapes)
names(s_features) <- c(
  "bio_vgpm",
  names(s_bio_now),
  #"fish_profit.subs"
  #"fish_mcp.2004",
  "phys_seamounts",
  "phys_vents",
  glue("phys_scape.{1:11}"))

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

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

r_pu_area <- area(r_pu_id) %>%  # in km2
  mask(r_pu_id)
#tif <- "/Users/bbest/github/bbnj/inst/app/www/scenarios/s03a.biofish.30pct.gl.now_sol.tif"
tif <- glue("{pfx}_sol.tif")
r_sol <- raster(tif)
area_km2 <- cellStats(r_pu_area * r_sol, "sum")

Area ($km^2$): r format(area_km2, big.mark=",")

tbl_target_representation(glue("{pfx}_rep.csv"))

ecoquants/bbnj documentation built on April 4, 2023, 9:08 p.m.