library(prioritizr)
# devtools::load_all()
# devtools::install_local(force=T)
# devtools::install_github("ecoquants/bbnj")
library(bbnj)
library(raster)
library(sf)
library(dplyr)
library(readr)
library(stringr)
library(glue)
library(here)
library(fs)
library(knitr)
library(formattable)
library(lwgeom)
library(purrr)
library(ggplot2)
library(formattable)
library(knitr)
library(ggplot2)
library(plotly)
area = raster::area
select = dplyr::select
Features
Source of features: s04a.biofish.alltime.mol50km.Rmd.
# borrowed from: ~/github/bbnj/inst/app/www/scenarios/
# s04a.biofish.alltime.mol50km.Rmd
prjres <- "_mol50km" # prjres in: View(projections_tbl)
P <- projections_tbl %>% filter(prjres == !!prjres)
# 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)
# 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")
Loop through Targets
rel_targets <- seq(0.1, 1, by=0.1)
for (rel_target in rel_targets[]){ # rel_target = rel_targets[1]
message(glue("rel_target: {rel_target} - {Sys.time()}"))
pfx <- glue("sensitivity_target/target_{rel_target}")
# problem ----
p <- problem(r_pu, s_features) %>%
add_min_set_objective() %>%
add_relative_targets(rel_target)
# solve ----
tif <- solve_log(p, pfx, redo=F)
}
Summarize Target Solutions by Percent Area of High Seas
rel_targets <- seq(0.1, 1, by=0.1)
d <- tibble(
target = seq(0.1, 1, by=0.1),
tif = glue("sensitivity_target/target_{target}_sol.tif"),
pct_hs = map_dbl(tif, function(x){ get_tif_area_stats(x)$pct_solution}))
d %>%
select(target, pct_hs) %>%
mutate(
target = formattable::percent(target, digits=0),
pct_hs = formattable::percent(pct_hs, digits=2)) %>%
kable()
| 10% |
7.16% |
| 20% |
15.15% |
| 30% |
23.70% |
| 40% |
32.85% |
| 50% |
42.57% |
| 60% |
52.83% |
| 70% |
63.60% |
| 80% |
74.87% |
| 90% |
86.35% |
| 100% |
99.33% |
g <- ggplot(data = d, aes(x=target, y=pct_hs)) +
#geom_line()
geom_area() + geom_point(color="red") +
scale_x_continuous(labels = scales::percent) +
scale_y_continuous(labels = scales::percent)
ggplotly(g)
rel_targets <- seq(0.1, 1, by=0.1)
d <- tibble(
target = seq(0.1, 1, by=0.1),
tif = glue("sensitivity_target/target_{target}_sol.tif"),
pct_hs = map_dbl(tif, function(x){ get_tif_area_stats(x)$pct_solution}))
d<- d %>%
select(target, pct_hs) %>%
mutate(
target = formattable::percent(target, digits=0),
pct_hs = formattable::percent(pct_hs, digits=2))
formattable(d,
align=c("c","c"),
col.names = c("Conservation Target", "Percentage of ABNJ in Solution"))
g <- ggplot(data = d, aes(x=target, y=pct_hs)) +
#geom_line()
geom_area() + geom_point(color="red") +
scale_x_continuous(labels = scales::percent) +
scale_y_continuous(labels = scales::percent)+
labs(x="Conservation Target", y="Percentage of ABNJ in Solution")
rel_targets <- seq(0.31, 0.39, by=0.01)
for (rel_target in rel_targets[]){ # rel_target = rel_targets[1]
message(glue("rel_target: {rel_target} - {Sys.time()}"))
pfx <- glue("sensitivity_target/target_{rel_target}")
# problem ----
p <- problem(r_pu, s_features) %>%
add_min_set_objective() %>%
add_relative_targets(rel_target)
# solve ----
tif <- solve_log(p, pfx, redo=T)
}
rel_targets <- seq(0.3, 0.4, by=0.01)
d2 <- tibble(
target = seq(0.3, 0.4, by=0.01),
tif = glue("sensitivity_target/target_{target}_sol.tif"),
pct_hs = map_dbl(tif, function(x){ get_tif_area_stats(x)$pct_solution}))
d2 %>%
select(target, pct_hs) %>%
mutate(
target = formattable::percent(target, digits=0),
pct_hs = formattable::percent(pct_hs, digits=2)) %>%
kable()
| 30% |
23.70% |
| 31% |
24.59% |
| 32% |
25.48% |
| 33% |
26.38% |
| 34% |
27.29% |
| 35% |
28.20% |
| 36% |
29.12% |
| 37% |
30.05% |
| 38% |
30.98% |
| 39% |
31.91% |
| 40% |
32.85% |
g2 <- ggplot(data = d2, aes(x=target, y=pct_hs)) +
#geom_line()
geom_area() + geom_point(color="red") +
scale_x_continuous(labels = scales::percent) +
scale_y_continuous(labels = scales::percent)
ggplotly(g2)