Solutions 4 - rescale
2022-06-08
Cumulative sum plots of layers for determining target per layer and allowing weights to have differential impact on scenario solution.
Get planning units as hexagons
# remotes::install_github("BenioffOceanInitiative/bbnj")
# devtools::install_local("../bbnj")
librarian::shelf(
BenioffOceanInitiative/bbnj, dplyr, DT,
exactextractr, fs, ggplot2, glue, here, mapview,
prioritizr,
purrr, raster, rlang, scales, sf, stringr,
tibble, tidyr, units)
source(here("libs/plots.R"))
select <- dplyr::select
mapviewOptions(fgb = FALSE)
mapviewOptions(
vector.palette = mapviewPalette("mapviewSpectralColors"))
# devtools::load_all(here("~/Github/bbest/prioritizr"))
# devtools::install_local(here("~/Github/bbest/prioritizr"))
#install.packages("/Library/gurobi951/macos_universal2/R/gurobi_9.5-1_R_4.1.1.tgz")
# copied into /Library/Frameworks/R.framework/Versions/4.2/Resources/library/gurobi
hex_geo <- "data/abnj_hex_res2.geojson"
# hex_geo <- "data/abnj_hex_res4.geojson"
pu_rds <- glue("{path_ext_remove(hex_geo)}_area.rds")
if (!file.exists(pu_rds)){
pu <- read_sf(hex_geo)
# calculate area
pu <- pu %>%
# st_make_valid() %>% # dateline issues
filter(st_is_valid(.)) %>% # 8 hexagons dropped
mutate(
area_km2 = st_area(geometry) %>%
set_units(km^2) %>%
drop_units(),
area_pct = area_km2 / sum(area_km2))
# sum(pu$area_pct) # check adds to 12
saveRDS(pu, file = pu_rds)
} else {
pu <- readRDS(pu_rds)
# paste(names(pu), collapse = ", ")
# hexid, on_dtln, lon, lat, FID, geometry, area_km2, area_pct
}
# mapview(pu, zcol="area_pct")
# plot(pu)Get layers
lyrs7_grd <- here("data/lyrs7_mol.grd")
if (!file.exists(lyrs7_grd)){
tifs1 <- list.files(here("../bbnj-scripts/map_layers"), "tif$", full.names = T)
tifs2 <- list.files(here("../bbnj-scripts/presentation_maps"), "tif$", full.names = T)[-1] # drop redundant bio_vgpm.tif
lyrs <- stack(c(tifs1, tifs2))
sort(names(lyrs))
s_nspp <- subset(lyrs, sort(names(lyrs))[26:48])
r_nspp <- calc(s_nspp, sum, na.rm = T)
s_mounts <- subset(lyrs, sort(names(lyrs))[49:51])
r_mounts <- calc(s_mounts, sum, na.rm = T)
r_vents <- raster(lyrs, layer = "phys_vents")
r_scapes <- raster(lyrs, layer = "scapes_hetero")
r_vgpm <- raster(lyrs, layer = "bio_vgpm")
r_rls <- raster(lyrs, layer = "rls_all")
r_fishing <- raster(lyrs, layer = "fishing_KWH")
# TODO: consider log/log10/normalize transforms of nspp, fish effort
#plot(r_vgpm)
s_lyrs <- stack(
r_fishing,
r_vgpm,
r_nspp,
r_rls,
r_vents,
r_mounts,
r_scapes)
names(s_lyrs) <- c(
"fishing",
"vgpm",
"nspp",
"rls",
"vents",
"mounts",
"scapes")
terra::writeRaster(s_lyrs, lyrs7_grd)
}
s_lyrs7 <- terra::rast(lyrs7_grd)
# plot(s_lyrs7)Extract layers to planning units
# extract feature values
f_0 <- exact_extract(
s_lyrs7, pu, fun = "mean", append_cols = "hexid", progress = F) %>%
tibble() %>%
rename_with(~str_replace(., "mean.", ""), everything())
# extra calculations
f_1 <- f_0 %>%
mutate(
fishing = 1/fishing, # TODO: log(1/fishing)?
spp = rls * nspp,
benthic = vents + mounts) %>% # TODO: rescale(vents) + rescale(mounts)?
mutate(
across(where(is.numeric), ~replace_na(., 0)))
summary(f_1)## hexid fishing vgpm nspp
## Length:3120 Min. :0.0000000 Min. : 0.0 Min. : 0.0
## Class :character 1st Qu.:0.0000000 1st Qu.: 191.3 1st Qu.: 179.4
## Mode :character Median :0.0000098 Median : 268.5 Median : 330.7
## Mean :0.0003546 Mean : 302.7 Mean : 295.5
## 3rd Qu.:0.0000438 3rd Qu.: 386.0 3rd Qu.: 397.2
## Max. :0.9072177 Max. :1568.3 Max. :1337.7
## rls vents mounts scapes
## Min. : 0.00 Min. :0.000000 Min. :0.00000 Min. : 0.000
## 1st Qu.:14.42 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.: 3.686
## Median :26.96 Median :0.000000 Median :0.09812 Median : 4.721
## Mean :25.43 Mean :0.003204 Mean :0.18451 Mean : 4.677
## 3rd Qu.:36.04 3rd Qu.:0.000000 3rd Qu.:0.26481 3rd Qu.: 5.846
## Max. :83.72 Max. :0.637052 Max. :2.98653 Max. :11.000
## spp benthic
## Min. : 0 Min. :0.00000
## 1st Qu.: 2900 1st Qu.:0.00000
## Median : 8692 Median :0.09886
## Mean : 9128 Mean :0.18771
## 3rd Qu.: 13266 3rd Qu.:0.26988
## Max. :111985 Max. :2.98653# rescale
f <- f_1 %>%
mutate(
across(where(is.numeric), rescale))
summary(f)## hexid fishing vgpm nspp
## Length:3120 Min. :0.0000000 Min. :0.0000 Min. :0.0000
## Class :character 1st Qu.:0.0000000 1st Qu.:0.1220 1st Qu.:0.1341
## Mode :character Median :0.0000108 Median :0.1712 Median :0.2472
## Mean :0.0003908 Mean :0.1930 Mean :0.2209
## 3rd Qu.:0.0000482 3rd Qu.:0.2461 3rd Qu.:0.2969
## Max. :1.0000000 Max. :1.0000 Max. :1.0000
## rls vents mounts scapes
## Min. :0.0000 Min. :0.000000 Min. :0.00000 Min. :0.0000
## 1st Qu.:0.1722 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.3351
## Median :0.3220 Median :0.000000 Median :0.03286 Median :0.4291
## Mean :0.3038 Mean :0.005029 Mean :0.06178 Mean :0.4252
## 3rd Qu.:0.4305 3rd Qu.:0.000000 3rd Qu.:0.08867 3rd Qu.:0.5314
## Max. :1.0000 Max. :1.000000 Max. :1.00000 Max. :1.0000
## spp benthic
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.02589 1st Qu.:0.00000
## Median :0.07761 Median :0.03310
## Mean :0.08151 Mean :0.06285
## 3rd Qu.:0.11846 3rd Qu.:0.09037
## Max. :1.00000 Max. :1.00000pu_f <- pu %>% # 129,906 × 8
left_join(f, by = "hexid") # 129,906 × 10
get_pu_v_ks <- function(pu_f, v, ks){
var <- as_label(enquo(v))
d_v <- pu_f %>%
st_drop_geometry() %>%
arrange(desc({{v}})) %>%
mutate(
var = var,
km2 = cumsum(area_km2),
val = cumsum({{v}})) %>%
select(var, km2, val)
val_sum <- d_v %>% pull(val) %>% last()
d_k <- lapply(
ks,
function(k){
i <- max(which(d_v$km2 < d_v$km2[nrow(d_v)]*k))
tibble(
var = var,
k = k,
km2 = d_v$km2[i],
val = d_v$val[i],
pct_val = val/val_sum) }) %>%
bind_rows()
p <- d_v %>%
ggplot(aes(x = km2, y = val)) +
geom_area(fill = "gray", color=NA)
add_k_paths <- function(k, d_k){
# https://colorbrewer2.org/#type=sequential&scheme=Blues&n=3
col <- c(
`0.3` = "#9ecae1",
`0.5` = "#3182bd")[as.character(k)]
x <- d_k %>% filter(k == !!k) %>% pull(km2)
y <- d_k %>% filter(k == !!k) %>% pull(val)
y_pct <- d_k %>% filter(k == !!k) %>% pull(pct_val)
z <- tribble(
~x, ~y,
x, 0,
x, y,
0, y)
list(
geom_path(data = z, aes(x,y), color=col),
annotate(
"text", x = x, y = y,
label = glue("{k*100}%: {comma(x, 0.1)} ({percent(y_pct, 0.1)})")))
}
p <- p +
lapply(
ks,
add_k_paths,
d_k) +
labs(
x = "area_km2",
y = var) +
theme_minimal()
list(
k=d_k, p=p)
}
# kp1 <- get_pu_v_ks(pu_f, fishing, c(0.3, 0.5))
# names(k_nspp)
# k_nspp$k
# kp1$p
k_p <- sapply(
c("fishing", "benthic", "scapes", "vgpm", "spp"),
function(var){
get_pu_v_ks(pu_f, !!sym(var), c(0.3, 0.5))
},
simplify = F)
p <- map(k_p, "p")
k <- map(k_p, "k") %>%
bind_rows() %>%
select(-km2) %>%
rename(pct_area = k)
datatable(k) %>%
formatPercentage(c("pct_area", "pct_val"), 1)walk(p, print)
s_name <- "area: 30% ; targets: max; weights: 1"Solution - area: 30% ; targets: max; weights: 1
# ocean area
a <- 0.3
# feature weights
f_w <- c(
benthic = 1,
fishing = 1,
scapes = 1,
spp = 1,
vgpm = 1)
f_w <- f_w[sort(names(f_w))]
# feature targets (max proportion given ocean area)
f_t <- k %>%
filter(pct_area == !!a) %>%
select(var, pct_val) %>%
arrange(var) %>%
deframe()
# problem ----
p <- pu_f %>%
problem(features = names(f_w), cost_column = "area_pct") %>%
add_min_shortfall_objective(a) %>%
add_relative_targets(f_t) %>%
add_feature_weights(f_w)
# solve ----
s <- solve(p, force = T)## Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (mac64[x86])
## Thread count: 6 physical cores, 12 logical processors, using up to 1 threads
## Optimize a model with 6 rows, 3125 columns and 16654 nonzeros
## Model fingerprint: 0x18810aa2
## Variable types: 5 continuous, 3120 integer (3120 binary)
## Coefficient statistics:
## Matrix range [3e-10, 1e+00]
## Objective range [2e-03, 8e-01]
## Bounds range [1e+00, 1e+00]
## RHS range [3e-01, 6e+02]
## Found heuristic solution: objective 5.0000001
## Found heuristic solution: objective 5.0000000
## Presolve removed 0 rows and 43 columns
## Presolve time: 0.02s
## Presolved: 6 rows, 3082 columns, 16606 nonzeros
## Found heuristic solution: objective 4.9964437
## Variable types: 5 continuous, 3077 integer (3077 binary)
## Root relaxation presolve removed 0 rows and 3 columns
## Root relaxation presolved: 6 rows, 3079 columns, 16599 nonzeros
##
##
## Root relaxation: objective 3.273993e-01, 22 iterations, 0.01 seconds (0.01 work units)
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 0.32740 0 2 4.99644 0.32740 93.4% - 0s
## H 0 0 0.3278579 0.32740 0.14% - 0s
##
## Explored 1 nodes (22 simplex iterations) in 0.03 seconds (0.02 work units)
## Thread count was 1 (of 12 available processors)
##
## Solution count 3: 0.327858 4.99644 5
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 3.278578334443e-01, best bound 3.273992762236e-01, gap 0.1399%# mapview(s %>% filter(solution_1==1), layer.name = s_name)
# gmap(s, "solution_1", s_name)
s %>% filter(solution_1==1) %>% pull(area_pct) %>% sum()## [1] 0.2999407# calculate feature representation statistics based on the prioritization
s_cov <- eval_target_coverage_summary(p, s[, "solution_1"])
datatable(s_cov)So scapes target met, whereas benthic has
highest shortfall of 16.9% shortfall.
s_name <- "area: 30% ; targets: max; weights: ∆ benthic"Solution - area: 30% ; targets: max; weights: ∆ benthic
Let’s ensure similar results.
# ocean area
a <- 0.3
# feature weights
f_w <- c(
benthic = 0.1,
fishing = 1,
scapes = 1,
spp = 1,
vgpm = 1)
f_w <- f_w[sort(names(f_w))]
# problem ----
p <- pu_f %>%
problem(features = names(f_w), cost_column = "area_pct") %>%
add_min_shortfall_objective(a) %>%
add_relative_targets(f_t) %>%
add_feature_weights(f_w)
# solve ----
s <- solve(p, force = T)## Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (mac64[x86])
## Thread count: 6 physical cores, 12 logical processors, using up to 1 threads
## Optimize a model with 6 rows, 3125 columns and 16654 nonzeros
## Model fingerprint: 0xe52cecef
## Variable types: 5 continuous, 3120 integer (3120 binary)
## Coefficient statistics:
## Matrix range [3e-10, 1e+00]
## Objective range [7e-04, 8e-01]
## Bounds range [1e+00, 1e+00]
## RHS range [3e-01, 6e+02]
## Found heuristic solution: objective 4.1000001
## Found heuristic solution: objective 4.1000000
## Presolve removed 0 rows and 43 columns
## Presolve time: 0.02s
## Presolved: 6 rows, 3082 columns, 16606 nonzeros
## Found heuristic solution: objective 4.0964437
## Variable types: 5 continuous, 3077 integer (3077 binary)
## Root relaxation presolve removed 0 rows and 3 columns
## Root relaxation presolved: 6 rows, 3079 columns, 16599 nonzeros
##
##
## Root relaxation: objective 6.642393e-02, 38 iterations, 0.01 seconds (0.01 work units)
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 0.06642 0 4 4.09644 0.06642 98.4% - 0s
## H 0 0 0.0709670 0.06642 6.40% - 0s
##
## Cleanup yields a better solution
##
## Explored 1 nodes (38 simplex iterations) in 0.03 seconds (0.02 work units)
## Thread count was 1 (of 12 available processors)
##
## Solution count 4: 0.0707813 0.070967 4.09644 4.1
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 7.078130146066e-02, best bound 6.642392785458e-02, gap 6.1561%#mapview(s %>% filter(solution_1==1), layer.name = "solution")
# gmap(s, "solution_1", s_name)
# s %>% filter(solution_1==1) %>% pull(area_pct) %>% sum()
# calculate feature representation statistics based on the prioritization
s_cov <- eval_target_coverage_summary(p, s[, "solution_1"])
datatable(s_cov)f_w.benthic: shortfall
- 0.1: 29.4%
- 1: 16.9%
- 2: 7.8%
- 4: 2.1%
- 10: 0.0453%
- 100: 0.0453%
- 1000: 0%
# turn off subsequent chunks
knitr::opts_chunk$set(eval = F, echo = F)