From notes | bbnj. Use same code as scenario_overlays? —> maybe no, instead try:

  1. take EEZ shapefile & turn into mollweide raster
  2. mask out HS from raster
  3. make vector from that raster
  4. calculate shared length of EEZ border with scenario 4
library(glue)
library(dplyr)
library(tidyr)
library(here)
library(glue)
library(bbnj)
# devtools::install_github("BenioffOceanInitiative/bbnj")
# devtools::install_local(here::here("../bbnj"))
library(fasterize)
library(sf)
library(lwgeom)
library(units)
library(tibble)
library(raster)
library(rmapshaper)
library(ggplot2)
library(readr)
library(DT)
library(scales)
library(rnaturalearth)
library(rnaturalearthdata)

dir_out       <- here("eez_scenario")
dir_bbnj      <- here("../bbnj")
dir_scenarios <- file.path(dir_bbnj, "inst/app/www/scenarios")
dir_data      <- file.path(dir_bbnj, "inst/data")
eez_s05_shp   <- file.path(dir_data, "eez_s05.shp")

#list.files(dir_scenarios)
scen      <- "s04a.biofish.alltime.mol50km"

scen_tif   <- glue("{dir_scenarios}/{scen}_sol.tif")
x_sol_shp  <- glue("{dir_out}/s04a_eez_bits.shp")
x_csv      <- glue("{dir_out}/s04a_eez.csv")
x_shp      <- glue("{dir_out}/s04a_eez.shp")
x_sov_csv  <- glue("{dir_out}/s04a_eez_sov.csv")
x_sov_shp  <- glue("{dir_out}/s04a_eez_sov.shp")
q_abnj_shp <- glue("{dir_out}/q_abnj.shp")
q_eez_shp  <- glue("{dir_out}/q_eez.shp")
q_sol_shp  <- glue("{dir_out}/q_s04a.shp")

r_scen <- raster(scen_tif)
plot(r_scen, main = "r_scen")

r_pu_id <- get_d_prjres("r_pu_id", "_mol50km") # plot(r_pu_id)
r_abnj  <- setValues(r_pu_id, 1) %>% mask(r_pu_id) # plot(r_abnj)

p_eez_s05     <- read_sf(eez_s05_shp)
p_eez_s05_mol <- p_eez_s05 %>% 
    filter(Territory1 != "Antarctica") %>% 
    filter(Pol_type == "200NM") %>% 
    st_transform(54009) # plot(p_eez_s05_mol["MRGID"])

if (!file.exists(x_csv)){
  
  r_eez_mrgid <- fasterize(p_eez_s05_mol, r_scen, field="MRGID")
  crs(r_eez_mrgid) <- crs(r_scen)
  # compareRaster(r_eez_mrgid, r_scen)
  # plot(r_eez_mrgid, main = "r_eez_mrgid")
  
  # q_*: polygons (p_*) that have been rasterized (r_*)
  if (!file.exists(q_eez_shp)){
    q_eez <- rasterToPolygons(r_eez_mrgid, dissolve=T) %>% 
      st_as_sf() %>%
      rename_at(vars(1), function(x) "eez_mrgid") %>%
      ms_dissolve(field="eez_mrgid") %>%
      st_buffer(0) %>% 
      st_make_valid() %>% 
      st_set_crs(54009)
    write_sf(q_eez , q_eez_shp)
  }
  if (!file.exists(q_sol_shp)){
    q_sol <- rasterToPolygons(r_scen, dissolve=T) %>% 
      st_as_sf() %>%
      rename_at(vars(1), function(x) "scen_sol") %>%
      ms_dissolve(field="scen_sol") %>%
      filter(scen_sol == 1) %>%
      st_buffer(0) %>% 
      st_make_valid() %>% 
      st_set_crs(54009)
    write_sf(q_sol , q_sol_shp)
  }
  if (!file.exists(q_abnj_shp)){
    q_abnj <- rasterToPolygons(r_abnj, dissolve=T) %>% 
      st_as_sf() %>%
      rename_at(vars(1), function(x) "abnj_1") %>%
      ms_dissolve(field="abnj_1") %>%
      st_buffer(0) %>% 
      st_make_valid() %>% 
      st_set_crs(54009)
    write_sf(q_abnj, q_abnj_shp)
  }

  q_abnj <- read_sf(q_abnj_shp)
  q_eez  <- read_sf(q_eez_shp) %>% 
    st_make_valid()
  q_sol  <- read_sf(q_sol_shp) %>% 
    st_make_valid()

  x_sol <- st_intersection(q_eez, q_sol) %>% 
    st_collection_extract("LINESTRING") %>% 
    group_by(eez_mrgid) %>% 
    summarize(
      sol_length_km = sum(st_length(geometry)) %>% set_units(km))
  
  x_abnj <- st_intersection(q_eez, q_abnj) %>% 
    st_collection_extract("LINESTRING") %>% 
    group_by(eez_mrgid) %>% 
    summarize(
      abnj_length_km = sum(st_length(geometry)) %>% set_units(km))
  
  x <- x_abnj %>% 
    left_join(
      x_sol %>% st_drop_geometry(),
      by = "eez_mrgid") %>% 
    mutate(
      sol_length_km = replace_na(sol_length_km, 0),
      pct_length    = sol_length_km / abnj_length_km) %>% 
    left_join(
      p_eez_s05 %>% st_drop_geometry(),
      by = c("eez_mrgid" = "MRGID"))
    
  x_sov <- x %>% 
    group_by(Sovereign1) %>% 
    summarise(
      sol_length_km  = sum(sol_length_km),
      abnj_length_km = sum(abnj_length_km)) %>% 
    mutate(
      pct_length = sol_length_km/abnj_length_km)

  write_csv(st_drop_geometry(x), x_csv)
  write_csv(st_drop_geometry(x_sov), x_sov_csv)
  write_sf(x_sov, x_sov_shp)
  write_sf(x, x_shp)
  write_sf(x_sol, x_sol_shp)
}
d      <- read_csv(x_csv)
d_sov  <- read_csv(x_sov_csv)
x      <- read_sf(x_shp)
#x_sov  <- read_sf(x_sov_shp)
q_abnj <- read_sf(q_abnj_shp)
q_eez  <- read_sf(q_eez_shp)
q_sol  <- read_sf(q_sol_shp)
x_sol  <- read_sf(x_sol_shp)

cntry <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  st_transform(54009)

theme_set(theme_bw())
ggplot() + 
  geom_sf(data = cntry , fill = "gray" , col = NA) +
  geom_sf(data = q_abnj, fill = "blue" , col = NA, alpha=0.2) + 
  geom_sf(data = q_sol , fill = "green", col = NA, alpha=0.3) + 
  geom_sf(
    data = q_eez %>% mutate(eez_mrgid = as.factor(eez_mrgid)),
    aes(fill = eez_mrgid), col = NA) +
  geom_sf(data = x_sol, col = "red") + 
  theme(legend.position="none")

d %>% 
  datatable(
    caption = "Intersection of Scenario 4a with individual EEZs.") %>% 
  formatPercentage(columns=c("pct_length"), digits=1) %>% 
  formatRound(columns=c("sol_length_km","abnj_length_km"), digits=1)

Download: s04a_eez.csv

d_sov %>% 
  datatable(
    caption = "Intersection of Scenario 4a with EEZs grouped by Sovereign1.") %>% 
  formatPercentage(columns=c("pct_length"), digits=1) %>% 
  formatRound(columns=c("sol_length_km","abnj_length_km"), digits=1)

Download: s04a_eez_sov.csv