From notes | bbnj
. Use same code as scenario_overlays? —> maybe no, instead try:
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