scenarios
library(tidyverse)
library(glue)
library(here)
library(sf)
library(lwgeom)
library(raster)
library(leaflet)
library(mapview)
library(rmapshaper)
library(knitr)
library(DT)
library(bbnj)
#devtools::load_all()
select = dplyr::select
mapviewOptions(basemaps = c("Esri.OceanBasemap","Stamen.TonerLite"))
dir_data <- here("../bbnj/inst/data")
dir_gdata <- "/Volumes/GoogleDrive/My Drive/BOI- Collaboration Folders/High Seas MPAs/bbnj/data" # on MV's laptop
#wdpa_sfx <- "raw/WDPA/WDPA_Aug2019_marine-shapefile/WDPA_Aug2019_marine-shapefile-polygons.shp"
ebsarank_sfx <- "raw/EBSAs/Global_EBSAs_automated_final_1104_2016_WGS84_Rankings/EBSA_rankings_crit6.shp"
scenarios <- tribble(
~s, ~scenario,
# "s1", "s01a.bio.now.mol50km",
# "s2", "s02a.bio.alltime.mol50km",
# "s3", "s03a.biofish.now.mol50km",
"s4", "s04a.biofish.alltime.mol50km")
datatable(scenarios)
layers
#MV removed WDPA on 11/5/19 b/c taking very long and not essential for manuscript
layers <- tribble(
~layer, ~clip, ~color, ~fld_id, ~fld_name, ~raw_shp,
"ebsa", T, "purple", "EBSA_ID", "NAME", glue("{dir_data}/ebsa.shp"),
"ihor", F, "yellow", "seaid", "sea", glue("{dir_data}/ihor.shp"),
"mine", F, "red", "ID", "name", glue("{dir_data}/mine-claims.shp"),
"ebsa_rank",T, "pink", "EBSA_ID", "NAME", glue("{dir_gdata}/{ebsarank_sfx}"))
# "wdpa", T, "brown", "WDPA_PID", "NAME", glue("{dir_gdata}/{wdpa_sfx}"))
datatable(layers)
scenarios x layers
get_scenarios_over_layers <- function(
scenarios, layers,
dir_scenarios = here("../bbnj/inst/app/www/scenarios"),
dir_overlays = here("data_overlays"),
redo = T){ #,
#get_layers = F,
#get_scenarios = F){
# dir_scenarios = here("../bbnj/inst/app/www/scenarios"); dir_overlays = here("data_overlays"); redo = T
# prep layers ----
layers <- layers %>%
mutate(
hs_shp = ifelse(
clip,
file.path(dir_overlays, glue("{layer}_hs.shp")),
raw_shp),
exists_hs_shp = file.exists(hs_shp))
lyrs_clip <- layers %>% filter(!exists_hs_shp)
for (layer in lyrs_clip$layer){ layer = lyrs_clip$layer[1]
lyr <- layers %>% filter(layer==!!layer)
read_sf(lyr$raw_shp) %>%
st_buffer(dist=0) %>%
st_intersection(p_abnj) %>%
mutate(
area_km2 = st_area(geometry) %>%
units::set_units(km^2)) %>%
filter(
st_geometry_type(geometry) != "GEOMETRYCOLLECTION") %>% # to remove points
write_sf(lyr$hs_shp)
}
layers <- layers %>%
mutate(
sf_hs = map(hs_shp, function(x) read_sf(x) %>% st_make_valid() ))
# prep scenarios ----
scenarios <- scenarios %>%
mutate(
sol_shp = file.path(dir_scenarios, glue("{scenario}_sol_gcs.shp")),
sf = map(sol_shp, function(x) read_sf(x) %>% st_make_valid() ))
# iterate over layers x scenarios ----
# redo = T
#for (layer in layers$layer){ # layer = "ihor"
for (layer in layers$layer){ # layer = "ihor"
lyr <- layers %>% filter(layer == !!layer)
message(glue("layer: {layer}"))
p_l <- lyr$sf_hs[[1]]
# + area_km2 if missing
if (!"area_km2" %in% names(p_l)){
p_l <- p_l %>%
mutate(
area_km2 = st_area(geometry) %>% units::set_units(km^2))
}
# iterate over scenarios
for (s in scenarios$s){ # s = "s4"
scen <- scenarios %>% filter(s == !!s)
scenario <- scen$scenario
s_l_shp <- file.path(dir_overlays, glue("{scenario}_sol_{layer}.shp"))
message(glue(" s: {s}"))
# DEBUG
if (layer == "ihor" & s %in% c("s1","s2","s3")){
message(glue(" DEBUG skip ihor {s}"))
next()
}
if (!file.exists(s_l_shp) | redo){
message(glue(" {layer} x {s} intersection -- {Sys.time()}"))
# get poly of scenario, dissolved
p_s <- scen$sf[[1]] %>%
ms_dissolve() %>%
st_make_valid()
# intersect p_l x p_s
p_s_l <- p_l %>%
rename(l_hs_km2 = area_km2) %>%
st_intersection(p_s) %>%
st_make_valid()
# check for funky geometry type, eg GEOMETRYCOLLECTION
p_s_l_notpoly <- p_s_l %>%
filter(!st_geometry_type(geometry) %in% c("POLYGON","MULTIPOLYGON"))
if (nrow(p_s_l_notpoly) > 0){
# eg: ihor x s4: North Pacific Ocean
p_s_l_poly <- p_s_l %>%
filter(st_geometry_type(geometry) %in% c("POLYGON","MULTIPOLYGON"))
flds_fix <- names(p_s_l) %>% setdiff("geometry")
p_s_l_fixpoly <- p_s_l_notpoly %>%
st_collection_extract("POLYGON") %>%
group_by_at(flds_fix) %>%
summarize()
p_s_l <- rbind(
p_s_l_poly,
p_s_l_fixpoly)
#write_sf(p_s_l, s_l_shp)
}
p_s_l <- p_s_l %>%
mutate(
# area of layer-scenario (vs layer-highseas)
l_s_km2 = st_area(geometry) %>%
units::set_units(km^2),
# percent of layer-scenario in the highseas
l_s_pct = l_s_km2 / l_hs_km2)
sum_l_s_km2 <- sum(p_s_l$l_s_km2, na.rm=T)
message(glue(" sum_l_s_km2: {format(sum_l_s_km2, big.mark=',')}"))
write_sf(p_s_l, s_l_shp)
} # if (!file.exists(s_l_shp)
} # for (s in scenarios$s)
} # for (layer in layers$layer)
# summary ----
s_l <- expand.grid(
scenarios$s, layers$layer, stringsAsFactors=F) %>%
as_tibble() %>%
rename(s=Var1, layer=Var2) %>%
left_join(
layers %>%
select(layer, layer_fld_id = fld_id), by = "layer") %>%
left_join(
scenarios %>%
select(s, scenario), by = "s") %>%
mutate(
shp = file.path(dir_overlays, glue("{scenario}_sol_{layer}.shp")),
sf = map(shp, function(x) read_sf(x)))
smry_layers <- layers %>%
mutate(
lyr_hs_n = map_int(sf_hs, function(x) nrow(x)),
lyr_hs_km2 = map_dbl(sf_hs, function(x)
sum(st_area(x) %>% units::set_units(km^2), na.rm=T))) %>%
select(layer, lyr_hs_n, lyr_hs_km2)
smry <- s_l %>%
left_join(
smry_layers, by = "layer") %>%
mutate(
lyr_scen_n = map2_int(sf, layer_fld_id, function(x, fld_id){
#x %>% ms_dissolve(fld_id) %>% nrow() }),
nrow(x) }),
lyr_scen_km2 = map_dbl(sf, function(x) sum(x$l_s_km2, na.rm=T)),
lyr_scen_pct = lyr_scen_km2 / lyr_hs_km2) %>%
rename(
scen = s, lyr = layer) %>%
select(-shp, -sf)
# map ----
message(glue("map"))
m <- leaflet() %>%
addProviderTiles(providers$Esri.OceanBasemap, group = "ocean") %>%
addPolygons(data = p_abnj_s05, color="gray80", group = "abnj")
for (s in scenarios$s){ # s = scenarios$s[1] # s = "s1"
scen <- scenarios %>% filter(s == !!s)
m <- m %>%
addPolygons(
data=scen$sf[[1]], color="green", group = s)
}
for (l in layers$layer){ # l = layers$layer[1] # l = "ebsa"
lyr <- layers %>% filter(layer == !!l)
lyr_lbls <- lyr$sf_hs[[1]] %>%
rename(id=lyr$fld_id, name=lyr$fld_name) %>%
mutate(
lbl = glue("{name} ({id})"))
m <- m %>%
addPolygons(
data=lyr$sf_hs[[1]], color=lyr$color, group = l, label=lyr_lbls$lbl)
}
m <- m %>%
addLayersControl(
baseGroups = c("ocean"),
overlayGroups = c("abnj", scenarios$s, layers$layer),
options = layersControlOptions(collapsed = F)) %>%
hideGroup(c(
head(scenarios$s, -1),
head(layers$layer, -1)))
list(
summary = smry,
map = m,
layers = layers,
scenarios = scenarios,
scenario_layers = s_l)
}
tbl_scenario_layer <- function(s_l, s, layer){
# s = "s5"; layer = "wdpa"
lyr <- filter(s_l$layers, layer == !!layer)
s_l$scenario_layers %>%
filter(s == !!s, layer == !!layer) %>%
pull(sf) %>% .[[1]] %>%
st_drop_geometry() %>%
select(
id = str_sub(lyr$fld_id, 1, 7),
name = lyr$fld_name,
hs_km2 = l_hs_k2,
scenario_km2 = l_s_km2,
scenario_pct = l_s_pct) %>%
datatable() %>%
formatRound(c("hs_km2", "scenario_km2"), 1) %>%
formatPercentage("scenario_pct", 1)
}
s_l <- get_scenarios_over_layers(scenarios, layers, redo=F)
s_l$map