map layers

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)
library(units)
library(RColorBrewer)

area   = raster::area
select = dplyr::select

# set rainbow color palette

#RColorBrewer::display.brewer.all()
pal  <- colorRampPalette(brewer.pal(9, "YlOrRd"))
#pal  <- colorRampPalette(brewer.pal(11, "Reds"))
#cols <- rev(pal(255))
cols <- pal(255)

Features

Source of features: s04a.biofish.alltime.mol50km.Rmd.

# borrowed from: ~/github/bbnj/inst/app/www/scenarios/
#                s04a.biofish.alltime.mol50km.Rmd

data("projections_tbl")

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 by taxonomic group----
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)

# biodiversity nspp_all----
s_bio_now_all <- get_gmbi_grpsmdl_prjres("groups00", prjres) %>% 
  subset("groups00_nspp_all")

# 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)),
  gsub("groups02_2100_nspp","nspp_2100",names(s_bio_gmbi_future)),
  "rls_all_now",
  "rls_all_future",
  #"fish_profit.subs"
  #"fish_mcp.2004",
  lbls_seamounts,
  "phys_vents",
  "scapes_hetero")

dir_png <- here("map_layers")
redo    <- T
  
s04 <- raster(here("../bbnj/inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol.tif"))
s04 <- mask(s04, s04, maskvalue=0)
P

## # A tibble: 1 x 8
##   prj   default name    proj                      epsg res_num res   prjres
##   <chr> <lgl>   <chr>   <chr>                    <dbl>   <dbl> <chr> <chr> 
## 1 mol   FALSE   Mollwe… +proj=moll +lon_0=0 +x_… 54009   50000 50km  _mol5…

s04_sf <- read_sf(here("../bbnj/inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol_gcs.shp")) %>%
  st_transform(P$epsg)
  
s04 <- mask(s04, s04, maskvalue=0)

summary25pct_csv <- glue("{dir_png}/_summary25pct.csv")

cell_km2   <- prod(res(r_pu)) %>% as_units("m2") %>% set_units("km2")
hs_ncells  <- length(na.omit(values(r_pu_id)))
hs_areakm2 <- hs_ncells * cell_km2

countries <- rnaturalearth::ne_countries(returnclass = "sf") %>%
  st_transform(P$epsg)
graticules <- st_graticule(countries)

# TODO
for (i in nlayers(s_features)){  i = 1
#for (i in 1:5){ # i = 1
  
  r   <- raster(s_features, 1)
  lyr <- names(s_features)[1]
  
  r_tif      <- glue("{dir_png}/{lyr}.tif")
  map_png    <- glue("{dir_png}/{lyr}_map.png")
  map25_png  <- glue("{dir_png}/{lyr}_map25.png")
  hist25_png <- glue("{dir_png}/{lyr}_hist25.png")
  message(glue("{i}/{nlayers(s_features)}: {lyr}"))
  
  # write raster ----
  if (!file.exists(r_tif)){
    message("  writeRaster")
    writeRaster(r, r_tif)
  }
  
  P   <- get_tif_projection(r_tif)

  # values ----
  v <- values(r) %>% na.omit() %>% sort(decreasing = T)
  d <- tibble(v = v) %>% 
    filter(v > 0) %>% 
    mutate(
      v_cumsum    = cumsum(v),
      i           = 1:length(v),
      cumarea_km2 = i * cell_km2 %>% drop_units(),
      pct_totval  = v_cumsum/sum(v) * 100,
      pct_hsarea  = cumarea_km2/hs_areakm2 %>% drop_units() * 100)
  
  d_25 <- d %>% 
    filter(v_cumsum <= sum(v)*.25)
  
  if (nrow(d_25) == 0){
    d_25 <- d %>% slice(1)
  }
    
  v_at25     <- min(d_25$v)
  pcths_at25 <- max(d_25$pct_hsarea)
  title  <- glue("{lyr} — 25%: {round(pcths_at25, digits=2)}% hs >= {signif(v_at25)}")

  D_i <- tibble(
    lyr        = lyr,
    v_at25     = v_at25,
    pcths_at25 = pcths_at25)
  
  if (!file.exists(summary25pct_csv)){
    write_csv(D_i, summary25pct_csv)
  } else {
    read_csv(summary25pct_csv) %>% 
      filter(lyr != !!lyr) %>% 
      bind_rows(D_i) %>% 
      write_csv(summary25pct_csv)
  }
  
  # hist ----
  if (!file.exists(hist25_png) | redo){
    message("  hist25_png")
    #png(hist25_png, width=480*4, height = 480*4, res=300, units='px')
    p <- ggplot(d, aes(x=pct_hsarea, y=pct_totval)) + 
      geom_area(fill="gray70", color=NA) + 
      geom_area(data=d_25, aes(x=pct_hsarea, y=pct_totval), fill="gray10", color=NA) + 
      scale_fill_distiller("spectral") +
      theme_light() + 
      coord_cartesian(expand = F) +
      #scale_x_continuous(labels = scales::percent_format()) %>% 
      scale_x_continuous(labels = scales::percent_format(accuracy = 3)) %>% 
      scale_y_continuous(labels = scales::percent_format()) %>% 
      labs(
        title = title, 
        x="% high seas area", y="% total value")
    ggsave(
      hist25_png,
      p, 
      "png", width=6.4, height=6.4, dpi=300, units='in')
    #dev.off()
  }
  
  # map 25 ----
  if (!file.exists(map25_png)){
    message("  map25_png")

    r_v   <- mask(r, r >= v_at25, maskvalue=0)
    v_rng <- range(values(r), na.rm = T)
    pct_at25_4cols <- (v_rng[2] - v_at25) / (v_rng[2] - v_rng[1])
    
    png(map25_png, width=480*4, height = 480*4, res=300, units='px')
    
    cols_v <- cols[1:round(length(cols)*pct_at25_4cols)]

    # simple map
    plot(r, col = scales::alpha(cols, 0.2), legend=F, axes=F, box=F, 
         legend.mar = 3.1, legend.shrink = 0.3)
    plot(r_v, col = cols_v, legend=F, add=T)
    contour(r, add=T, levels=v_at25, col="red", lwd = 0.5)
    plot(s04_sf[1], col=scales::alpha("blue", 0.15), border=scales::alpha("blue", 0.25), lwd=0.5, add=T)
    plot(countries[1], col=gray(0.8), border=gray(0.7), lwd=0.5, add=T)
    plot(graticules[1], col=gray(0.6), lwd=0.5, add=T)
    plot(r, col = cols, legend.only=T, #legend.width = 1,
         legend.shrink = 0.3, legend.mar = 3.1,
         legend.args = list(text = lyr, side = 3, 
         font = 2, line = 0.2, cex = 0.5),
         axis.args = list(cex.axis = 0.5))

    dev.off()
    
    map25_png %>%
      magick::image_read() %>% magick::image_trim() %>%
      magick::image_write(map25_png)
  }
  
  # map ----
  if (!file.exists(map_png) | redo){
    message("  map_png")
    png(map_png, width=480*4, height = 480*4, res=300, units='px')
    
    #table(values(r))
    plot(r, col = cols, legend=T, axes=F, box=F, main=lyr)
    plot(s04, alpha=0.5, legend=F, axes=F, box=F)
    plot(st_geometry(countries), add=T, col=gray(0.8), border=gray(0.7), lwd=0.5)
    plot(st_geometry(graticules), add=T, col=gray(0.6), lwd=0.5)
    
    dev.off()
    map_png %>%
      magick::image_read() %>% magick::image_trim() %>%
      magick::image_write(map_png)
  }
}

# read in s04
s04 <- raster(here("../bbnj/inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol.tif"))
s04 <- mask(s04, s04, maskvalue=0)

s04_sf <- read_sf(here("../bbnj/inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol_gcs.shp")) %>%
  st_transform(P$epsg)
  
s04 <- mask(s04, s04, maskvalue=0)

highlights2<-read_sf(here("../bbnj-scripts/highlights/s04highlights_group2.shp"))%>%
  st_transform(P$epsg)

highlights1<-read_sf(here("../bbnj-scripts/highlights/s04highlights_group1.shp"))%>%
  st_transform(P$epsg)

# countries
countries <- rnaturalearth::ne_countries(returnclass = "sf") %>%
  st_transform(P$epsg)
graticules <- st_graticule(countries)

#top quartile biodiversity now all
s_bio_now_all_v25<-quantile(s_bio_now_all, probs=0.75)

s_bio_now_all_25<- s_bio_now_all %>% 
  reclassify(rcl=cbind(-Inf,s_bio_now_all_v25,NA), right=TRUE)

  dir_png <- here("../bbnj-scripts/map_layers")
  map_png <- glue("{dir_png}/nspp_all_25.png")

png(map_png, width=480*6, height = 480*6, res=300, type="cairo", units='px', bg = "transparent")

plot(s_bio_now_all, col=rev( rainbow(30, start=0,end=0.7 ) ), breaks=seq(min(minValue(s_bio_now_all)),max(maxValue(s_bio_now_all)),length.out=100), alpha=0.2, legend=F, axes=F, box=F)
plot(s_bio_now_all_25, col=rev( rainbow(30, start=0,end=0.7 ) ), breaks=seq(min(minValue(s_bio_now_all)),max(maxValue(s_bio_now_all)),length.out=100), legend=F, axes=F, box=F, add=T)
plot(countries[1], col=gray(0.8), border=gray(0.7), lwd=0.5, add=T)
plot(graticules[1], col=gray(0.6), lwd=0.5, add=T)
plot(s04_sf[1], col=scales::alpha("black", 0.3), border=scales::alpha("black", 0.5), lwd=0.5, add=T)
plot(highlights1[1], col="orange",pch=19, cex=0.4,add=T)
plot(highlights2[1], col="orange",pch=19, cex=0.4,add=T)

dev.off()

D <- read_csv(summary30pct_csv)

if (exists("s_b")) rm(s_b)
r_sol <- r_pu - 1 # plot(R_b)
for (i in 1:nlayers(s_features)){ # i = 3
  
  r <- s_features[[i]]
  lyr <- names(s_features)[i]
  
  v_at30 <- D %>% 
    filter(lyr == !!lyr) %>% 
    pull(v_at30)
  
  r_sol[r >= v_at30] <- 1
  
  sol_areakm2  <- sum(values(r_sol), na.rm = T) * cell_km2
  pct_sol_hs <-  sol_areakm2 / hs_areakm2 * 100 # 84.47849 % high seas
  title <- glue("{i}/{nlayers(s_features)}: +{lyr}, {round(pct_sol_hs, digits=2)}% hs")
  
  Sys.sleep(1)
  plot(r_sol, main = title)
}
sol30_tif <- glue("{dir_png}/sol30.tif")
sol30_map_png <- glue("{dir_png}/sol30_map.png")

writeRaster(r_sol, sol30_tif)

png(sol30_map_png, width=480*4, height = 480*4, res=300, units='px')
#table(values(r))
plot(r_sol, legend=F, axes=F, box=F)
plot(st_geometry(countries), add=T, col=gray(0.8), border=gray(0.7), lwd=0.5)
plot(st_geometry(graticules), add=T, col=gray(0.6), lwd=0.5)
dev.off()
sol30_map_png %>%
  magick::image_read() %>% magick::image_trim() %>%
  magick::image_write(sol30_map_png)

  r <- r_fish_effort
  lyr <- names(r_fish_effort)
  #plot(r)
  dir_png <- here("inst/scripts/map_layers")
  r_tif   <- glue("{dir_png}/{lyr}.tif")
  map_png <- glue("{dir_png}/{lyr}_map.png")
  
  writeRaster(r, r_tif, overwrite=T)
  
  P <- get_tif_projection(r_tif)
  
  pal <- colorRampPalette(brewer.pal(7, "YlOrRd"))
  cols <- (pal(7))
  
  # plot ----
  countries <- rnaturalearth::ne_countries(returnclass = "sf") %>%
    st_transform(P$epsg)
  graticules <- st_graticule(countries)
  
  png(map_png, width=480*4, height = 480*4, res=300, type="cairo", units='px', bg = "transparent")
  
  plot(r, breaks=c(0,10000,30000,100000,300000,1000000,5000000,108528675),col=cols, legend=F, axes=F, box=F)
  plot(eez, add=T, col=gray(0.9), legend=F)
  plot(st_geometry(countries), add=T, col=gray(0.8), border=gray(0.7), lwd=0.5)
  plot(st_geometry(graticules), add=T, col=gray(0.6), lwd=0.5)
  
  dev.off()
  map_png %>%
    magick::image_read() %>% magick::image_trim() %>%
    magick::image_write(map_png)

library(tidyr)
library(scales)

extract <- raster::extract

sol_shp <- here("inst/app/www/scenarios/s04a.biofish.alltime.mol50km_sol_gcs.shp")

sol_sf <- read_sf(sol_shp)

ply <- sol_sf %>% arrange(desc(are_km2)) %>% head(1)
# mapview::mapview(ply)

ply_mol_sp <- ply %>% st_transform(P$epsg) %>% as("Spatial")

d_ply <- extract(s_features, ply_mol_sp, df=T) %>% 
  as_tibble() %>% 
  select(-ID) %>% 
  gather("layer", "value") %>% 
  na.omit() %>% 
  mutate(
    source = "ply")

d_hs <- values(s_features) %>% 
  as_tibble() %>% 
  gather("layer", "value") %>% 
  na.omit() %>% 
  mutate(
    source = "hs")

d_vs <- bind_rows(
  d_hs,
  d_ply) %>% 
  mutate(
    score = rescale(value))

d_vs_avg <- d_vs %>% 
  group_by(source, layer) %>% 
  summarize(
    avg_score = mean(score))

d_vs_avg_w <- d_vs_avg %>% 
  ungroup() %>% 
  mutate(source = glue("avg_{source}")) %>% 
  spread(source, avg_score) %>% 
  mutate(
    pct_ply = avg_ply / avg_hs) %>% 
  arrange(desc(pct_ply))

# TODO: include s04 features
d_vs_avg_w # YES!

DT::datatable(d_vs_avg_w)
ggplot(d_vs, aes(x = score, y = layer, fill = source)) +
  #geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01) +
  geom_density_ridges(alpha = .5, color = "white") #+
  # geom_segment(
  #   aes(x = avg_pct, y = layer, xend = avg_pct, yend = layer + 1),
  #   data = d_vs_avg)

  #geom_vline(aes(xintercept = avg, y = layer), data = d_vs_avg) +
  #theme_ridges(font_size = 20, grid=TRUE, line_size=1, 
  #             center_axis_labels=TRUE)
  #scale_fill_viridis(name = "Temp. [F]", option = "C") #+
  #scale_fill_distiller("spectral", name = "value") #+
  #scale_fill_gradientn(colors = cols, name = "value") #+

png(hist_png, width=480*4, height = 480*4, res=300, units='px')
ggplot(d_features, aes(x = value, y = layer, fill = ..x..)) +
  #geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01) +
  geom_density_ridges_gradient() +
  
  #scale_fill_viridis(name = "Temp. [F]", option = "C") #+
  #scale_fill_distiller("spectral", name = "value") #+
  scale_fill_gradientn(colors = cols, name = "value") #+
#labs(title = 'Temperatures in Lincoln NE in 2016')
dev.off()


library(ggridges)
library(forcats)
Catalan_elections %>%
  mutate(YearFct = fct_rev(as.factor(Year))) %>%
  ggplot(aes(y = YearFct)) +
  geom_density_ridges(
    aes(x = Percent, fill = paste(YearFct, Option)), 
    alpha = .8, color = "white", from = 0, to = 100
  ) +
  labs(
    x = "Vote (%)",
    y = "Election Year",
    title = "Indy vs Unionist vote in Catalan elections",
    subtitle = "Analysis unit: municipalities (n = 949)",
    caption = "Marc Belzunces (@marcbeldata) | Source: Idescat"
  ) +
  scale_y_discrete(expand = c(0.01, 0)) +
  scale_x_continuous(expand = c(0.01, 0)) +
  scale_fill_cyclical(
    breaks = c("1980 Indy", "1980 Unionist"),
    labels = c(`1980 Indy` = "Indy", `1980 Unionist` = "Unionist"),
    values = c("#ff0000", "#0000ff", "#ff8080", "#8080ff"),
    name = "Option", guide = "legend"
  ) +
  theme_ridges(grid = FALSE)