get_ter_eez_wgcs_sf <- function(ter){
read_sf(eez_wgcs_geo) %>%
filter(territory==ter)
}
get_ter_depth_wgcs_r <- function(ter){
dir_depth <- glue("{dir_lyrs}/depth")
dir.create(dir_depth, showWarnings=F)
# ter = "Great Lakes"
ter_depth_wgcs_tif <- glue::glue("{dir_depth}/{ter}_depth_epsg4326.tif")
if (!file.exists(ter_depth_wgcs_tif)){
cat(" creating", ter_depth_wgcs_tif, "\n")
depth_wgcs_tif <- glue::glue("{dir_depth}/_depth_epsg4326.tif")
if (!file.exists(depth_wgcs_tif)){
cat(" creating", depth_wgcs_tif, "\n")
depth_gcs <- raster(depth_nc, layer = "elevation") # depth_gcs_0 <- depth_gcs
# burn in [Great Lakes Bathymetry | NCEI](https://www.ngdc.noaa.gov/mgg/greatlakes/greatlakes.html)
gl_tifs <- list.files(
"/Volumes/Best HD/nrel_data_big/greatlakes_ngdc",
".*\\.tif$", recursive = T, full.names = T)
r_gl <- do.call(raster::merge, lapply(gl_tifs, raster::raster)) # pl
r_gl_gcs <- raster::projectRaster(r_gl, depth_gcs) # plot(r_gl_gcs)
depth_gcs <- raster::merge(r_gl_gcs, depth_gcs)
depth_wgcs <- raster_wrap(depth_gcs)
depth_wgcs <- depth_wgcs * -1
writeRaster(depth_wgcs, depth_wgcs_tif, overwrite=T)
} else {
depth_wgcs <- raster(depth_wgcs_tif)
}
ter_eez_wgcs_sf <- get_ter_eez_wgcs_sf(ter) # plot(ter_eez_wgcs_sf["territory"])
ter_eez_wgcs_r <- fasterize(ter_eez_wgcs_sf, depth_wgcs) %>% trim() # plot(ter_eez_wgcs_sf["territory"]); plot(ter_eez_wgcs_r, add=T)
ter_depth_wgcs <- raster_clip(depth_wgcs, ter_eez_wgcs_r) # plot(ter_eez_wgcs_sf["territory"]); plot(ter_depth_wgcs, add=T)
ter_depth_wgcs[ter_depth_wgcs < 0] <- NA # plot(ter_eez_wgcs_sf["territory"]); plot(ter_depth_wgcs >= 0, add=T)
ter_depth_wgcs <- trim(ter_depth_wgcs) # plot(ter_eez_wgcs_sf["territory"]); plot(ter_depth_wgcs, add=T)
ter_depth_wgcs_tif <- sprintf("./depth/%s_depth_epsg4326.tif", ter)
writeRaster(ter_depth_wgcs, ter_depth_wgcs_tif, overwrite=T)
}
ter_depth_wgcs <- raster(ter_depth_wgcs_tif) # plot(ter_depth_wgcs)
ter_depth_wgcs # plot(ter_depth_wgcs)
}
ter_stack <- function(ter, lyr_params){
# ter = "West"
lyrs <- sort(lyr_params %>% filter(run==T) %>% .$key)
lyrs_in_ter <- setNames(rep(NA, length(lyrs)), lyrs) %>% as.list()
if (exists("s_ter", inherits=F)) rm(s_ter)
for (lyr in lyrs){ # lyr = "aquaculture" # "wind" # key <- "mpa" # keys[1] "efh" mpa" "wind"
lyr_info_yml <- glue::glue("{dir_lyrs}/{lyr}/_{lyr}.yml")
lyr_info <- yaml::read_yaml(lyr_info_yml)
# skip if layer not in territory
if (is.na(lyr_info$territories[[ter]])){
next
}
# get raster stack of layer components
msg(g(" {ter} - {lyr}"))
lyrs_in_ter[[lyr]] <- names(lyr_info$territories[[ter]]$components)
tifs <- purrr::map_chr(
lyr_info$territories[[ter]]$components,
~glue::glue("{dir_lyrs}/{lyr}/{.x[1]}"))
#if (lyr == "oceanuseatlas.hi" & ter == "Hawaii") browser("TODO: extend all to depth r")
#r_a <- fasterize::fasterize(x, y, field=field, fun="first", by=by)
s_lyr <- raster::stack(tifs)
if (length(lyrs_in_ter[[lyr]]) == 1 && lyr == lyrs_in_ter[[lyr]]){
names(s_lyr) <- lyr
} else {
names(s_lyr) <- glue::glue("{lyr}_{lyrs_in_ter[[lyr]]}")
}
# add layer stack to territory stack
if (!exists("s_ter", inherits=F)){
s_ter <- s_lyr
} else {
s_ter <- raster::stack(s_ter, s_lyr)
}
}
# add area to territory stack
s_ter <- raster::stack(s_ter, raster::area(s_ter))
names(s_ter)[raster::nlayers(s_ter)] <- "areakm2"
lyrs_in_ter$areakm2 <- "areakm2"
# add attribute of layers in territory stack
attr(s_ter, "lyrs_in_ter") <- lyrs_in_ter
s_ter
}
ter_stack_tbl <- function(s){
tbl_s <- raster::as.data.frame(s) %>%
tibble::as_tibble()
}
ter_bin <- function(tbl_s, bin_by = "depth"){
lims <- nrel_limits[[bin_by]]
tbl_b <- tbl_s %>%
dplyr::filter(!is.na(depth))
if (bin_by != "depth"){
tbl_b <- tbl_b %>%
dplyr::filter(
depth >= lims$depth$min,
depth <= lims$depth$max)
}
tbl_b$vals <- tbl_b[[bin_by]]
tbl_b <- tbl_b %>%
dplyr::mutate(
bin = cut(vals, lims$breaks, lims$break_labels, include.lowest = T)) %>%
dplyr::filter(
!is.na(bin))
tbl_b
}
ter_bin_lyr <- function(tbl_b){
tbl_b <- tbl_b %>%
dplyr::mutate(
bin = factor(as.character(bin), c(levels(bin), "Total"), ordered=T))
lyrs_rm <- c(names(nrel_limits), "vals")
tbl_b <- suppressWarnings(dplyr::select(tbl_b, -dplyr::one_of(lyrs_rm)))
tbl_b
}
ter_bin_lyr_smry <- function(tbl_b){
# totals
dt <- tbl_b %>%
#dplyr::mutate(
# bin = as.character(bin)) %>%
dplyr::group_by(bin) %>%
dplyr::summarise(
areakm2 = sum(areakm2))
dt <- dt %>%
tibble::add_row(
bin = "Total",
areakm2 = sum(dt$areakm2)) %>%
tidyr::spread(bin, areakm2)
dt <- tibble::tibble(
layer = " Total Area") %>%
dplyr::bind_cols(dt)
dt$rank <- 0
# dl: layer x bin = areakm2
#if (ter=="Atlantic Islands" & type == "wave" & element == "count")
#browser()
dl <- tbl_b %>%
tidyr::gather(key = layer, value = value, -bin, -areakm2) %>%
dplyr::filter(!is.na(value)) %>%
dplyr::group_by(layer, bin) %>%
dplyr::summarise(
areakm2 = sum(areakm2))
if (nrow(dl) == 0){
dl <- NULL
} else {
dl <- dl %>%
tidyr::spread(bin, areakm2) %>%
dplyr::ungroup()
dl$Total <- rowSums(dl[,-1], na.rm = TRUE)
# rank
dl <- dl %>%
dplyr::mutate(
rank = dplyr::dense_rank(1/Total))
}
# bind totals
tbl_bls <- dplyr::bind_rows(dt, dl)
tbl_bls
}
svg_bar <- function(pct, val){
pct <- ifelse(is.na(pct), 0, round(pct * 100, 3))
#cat("val before", val,"\n")
val <- ifelse(is.na(val), "0", format(round(val), big.mark = ","))
#cat("val after", val,"\n")
glue::glue('
<svg xmlns="http://www.w3.org/2000/svg" version="1.1" width="100%" height="14px" preserveAspectRatio="xMaxYMin meet">
<g>
<rect width="{pct}%" height="16" fill="DarkGray"></rect>
<text text-anchor="end" x="100%" y="12">{val}</text>
</g>
</svg>
')
# <text text-anchor="end" x="100%" y="11" font-family="Courier New" font-size="14px" color="black">{val}</text>
}
ter_bin_lyr_smry_plus <- function(ter, bin, lyr_params, lyr_categories){
ter_info <- get_ter_info(ter)
ter_bin_lyr_csv <- ter_info[[bin]]$count$layers_csv
tbl_bls <- suppressMessages(read_csv(ter_bin_lyr_csv, trim_ws=F))
bins <- tbl_bls %>%
select(-one_of(c("layer","rank"))) %>%
names()
dt <- tbl_bls %>%
filter(layer == " Total Area")
# percent, svg
#browser()
for (bin in bins){ # bin = bins[1]
col_pct <- glue::glue("{bin}_PCT")
col_svg <- glue::glue("{bin}_SVG")
tbl_bls[[col_pct]] <- tbl_bls[[bin]] / dt[[bin]]
tbl_bls[[col_svg]] <- purrr::map2_chr(tbl_bls[[col_pct]], tbl_bls[[bin]], svg_bar)
}
# key, component
tbl_bls <- tbl_bls %>%
dplyr::mutate(
key = stringr::str_replace(layer, "(.*)_(.*)", "\\1"),
component = ifelse(
key == layer, "",
stringr::str_replace(layer, "(.*)_(.*)", "\\2") %>% stringr::str_replace_all("\\.", " ")))
# title, category
tbl_blsp <- tbl_bls %>%
dplyr::left_join(
lyr_params %>%
dplyr::select(key, title) %>%
dplyr::bind_rows(tibble::tibble(
key = " Total Area",
title = " Total Area")),
by = "key") %>%
left_join(
lyr_categories %>%
select(key, use),
by="key")
tbl_blsp
}
print_ter_bin_lyr_smry_plus <- function(tbl_blsp){
# variables used, but not declared
# i_tbl (global, write)
# ter
# bin_html
# TODO: if html: sortable/searchable header
fmt <- get_fmt()
bins <- names(tbl_blsp)[2:which(names(tbl_blsp)=="Total")]
#browser()
th_bins <- htmltools::tagList(lapply(bins, htmltools::tags$th))
hdr = htmltools::withTags(table(
class = "display",
thead(
tr(
th(colspan=3, "Ocean Use"),
th(colspan=5, HTML(glue("Area (km<sup>2</sup>) by {bin_html}")))),
tr(
th("Category"),
th("Dataset"),
th("Layer"),
th_bins,
th("Rank")))))
i_tbl <<- i_tbl + 1
caption <- HTML(glue('Table {i_tbl}: Area (km<sup>2</sup>) of Ocean Use by {bin_html} bin for {ter} region.
Width of gray horizontal bars from left of cell indicate percent area occupied in given {bin_html} bin.'))
tbl_blsp %>%
select(-one_of(glue("{bins}"))) %>%
rename_at(
vars(one_of(glue("{bins}_SVG"))),
~stringr::str_replace(., "_SVG", "")) %>%
select(layer, one_of(glue("{bins}")))
d <- tbl_blsp[1,] %>%
bind_rows(
tbl_blsp[-1,] %>%
arrange(use, title, component)) %>%
rename(
Category = use,
Layer = title,
Component = component,
Rank = rank) %>%
select(-one_of(glue("{bins}"))) %>%
rename_at(
vars(one_of(glue("{bins}_SVG"))),
~stringr::str_replace(., "_SVG", "")) %>%
select(one_of(c("Category", "Layer", "Component", bins, "Rank")))
datatable(
d,
escape = F,
rownames = F,
container = hdr,
caption = caption,
options = list(
pageLength = nrow(tbl_blsp),
dom=c(
html_document = 'lftir',
pdf_document = 't',
word_document = 't')[[fmt]],
ordering=c(
html_document = TRUE,
pdf_document = FALSE,
word_document = FALSE)[[fmt]],
autoWidth = TRUE,
columnDefs = list(
list(width = '80px' , targets = 3:(2+length(bins))),
list(className = 'dt-right', targets = 3:(3+length(bins))))))
}
ter_bin_cnt <- function(tbl_b){
# table for cumulative count plot
cols_tally <- setdiff(names(tbl_b), c(names(nrel_limits), "bin", "vals", "areakm2"))
tbl_b$tally <- apply(tbl_b[, cols_tally], 1, function(x) sum(!is.na(x)))
tbl_b <- tbl_b %>%
dplyr::group_by(bin, tally) %>%
dplyr::summarise(
areakm2 = sum(areakm2)) #%>%
#filter(tally>0) # View(tbl_df2)
lyrs_rm <- c(names(nrel_limits), "vals")
tbl_bc <- suppressWarnings(dplyr::select(tbl_b, -dplyr::one_of(lyrs_rm)))
tbl_bc
}
plot_caption_ter_bin_cnt <- function(ter, bin_html){
glue::glue("Area of Cumulative Ocean Use by {bin_html} for {ter} region. Cumulative use are counted and summarized by {bin_html} and area.")
}
plot_ter_bin_cnt <- function(ter, bin){
fmt <- get_fmt()
ter_info <- get_ter_info(ter)
ter_bin_cnt_csv <- ter_info[[bin]]$count$cumulative_csv
tbl_bc <- suppressMessages(read_csv(ter_bin_cnt_csv, trim_ws=F))
tbl_bc$bin <- with(tbl_bc, factor(bin, unique(bin), ordered=T))
g <- ggplot(data = tbl_bc, aes(x = bin, y = areakm2, fill = tally)) +
geom_col() +
scale_fill_distiller(palette = "Spectral", name="Count")
if (fmt == "html_document"){
g <- g +
labs(x = bin_html, y = "Area (km<sup>2</sup>)")
g <- plotly::ggplotly(g)
} else {
g <- g +
labs(x = bin_html, y = expression(Area (km^2)))
}
g
}
ter_cnt_raster <- function(s){
lyrs_rm <- c(names(nrel_limits), "areakm2")
s_cnt <- raster::dropLayer(s, which(names(s) %in% lyrs_rm))
r_cnt <- sum(!is.na(s_cnt), na.rm=T) %>%
raster::mask(raster::raster(s, "depth"))
r <- raster_unwrap(r_cnt) # plot(r)
r_3857 <- raster_project_leaflet_nearest(r)
r_3857
}
ter_cnt_raster2 <- function(s){
# TODO: use yaml for epsg4326:*.tif \n epsg3857:\n big:*.tif\n small:*.tif
#ter_cnt_raster <- function(s){
lyrs_rm <- c(names(nrel_limits), "areakm2")
s_cnt <- raster::dropLayer(s, which(names(s) %in% lyrs_rm))
r_cnt <- sum(!is.na(s_cnt), na.rm=T) %>%
raster::mask(raster::raster(s, "depth"))
r <- raster_unwrap(r_cnt) # plot(r)
ext <- raster::extent(s)
ext@xmin < 180
r <- ter_cnt_raster(s)
r
}
map_ter_caption <- function(ter, type, element){
str_type <- ifelse(element == "count", "cumulative ocean use", glue::glue("viable {type}"))
glue::glue("Map of {str_type} in {ter} region.")
# TODO: if ( fmt==html_document & length(info$) ),
# pan to other side of dateline to see rest of map
}
#' Leaflet map of territory's ocean use count
#'
#' @param info
#'
#' @return
#' @export
#'
#' @examples
lmap_ter <- function(ter, type, element){
# ter='East'; type='wind'; element='limits'
ter_info <- get_ter_info(ter)
el <- ter_info[[type]][[element]]
b <- el$epsg3857_extent_init
r <- raster::raster(el$epsg3857_tif[1])
# TODO: fix creation vs here
if (b[1] > 180){
b[1:2] = b[1:2] - 360
}
vals <- getValues(r)
if (length(el$epsg3857_tif) == 2){
r2 <- raster::raster(el$epsg3857_tif[2])
vals <- c(vals, getValues(r2))
}
if (element == "limits"){
# TODO: for limits project to leaflet without interpolation so integer
brk_labels <<- nrel_limits[[type]]$break_labels
legend_title <- stringr::str_to_title(type)
vals <- round(vals)
r <- round(r)
if (length(el$epsg3857_tif) == 2){
r2 <- round(r2)
}
pal_col <- ifelse(type == "depth", "Blues", "Greens")
pal <- leaflet::colorFactor(
palette = pal_col,
domain = factor(round(vals), 1:length(brk_labels), brk_labels),
levels = 1:length(brk_labels),
na.color="#00000000", reverse=F)
} else {
legend_title <- "Count"
pal <- colorNumeric(
palette = 'Spectral', na.color="#00000000",
reverse=T, domain = vals)
}
m <- leaflet::leaflet(
options=leafletOptions(worldCopyJump=T)) %>%
#options=c(leafletOptions(), attributionControl=F, zoomControl=F, worldCopyJump=T)) %>%
leaflet::addProviderTiles(providers$Stamen.TonerLite, group = "B&W") %>%
leaflet::addProviderTiles(providers$Esri.OceanBasemap, group = "Ocean") %>%
leaflet::addRasterImage(r, project=F, colors=pal, opacity=0.8, group="Count")
if (length(el$epsg3857_tif) == 2){
m <- m %>%
leaflet::addRasterImage(r2, project=F, colors=pal, opacity=0.8, group="Count")
}
m <- m %>%
leaflet::addGraticule() %>% # interval=1
leaflet::addScaleBar("bottomleft") %>%
leaflet::fitBounds(b[1], b[3], b[2], b[4]) %>%
leaflet::addLayersControl(baseGroups = c("B&W", "Ocean"), overlayGroups=c("Count"))
if (element == "limits"){
m <- m %>%
leaflet::addLegend(
"bottomright",
colors = pal(1:length(brk_labels)),
labels = brk_labels,
opacity = 0.8, title = legend_title)
} else {
m <- m %>%
leaflet::addLegend("bottomright", pal, vals, opacity=0.8, title=legend_title)
}
m
}
#' Ggplot map of territory's ocean use count
#'
#' @param r
#'
#' @return
#' @export
#'
#' @examples
gmap_ter <- function(ter, type, element){
# ter <- "Alaska"; type <- "depth"; element <- "limits"
#gmap_ter("Alaska", "depth", "limits")
# ter <- ; type <- "depth"; element <-
ter_info <- get_ter_info(ter)
el <- ter_info[[type]][[element]]
# devtools::load_all("~/github/nrelutils")
r <- raster::raster(el$epsg4326_tif) %>% raster_trim()
b <- raster::extent(r) %>% as.vector()
bb <- extent(r) %>% as("SpatialPolygons") %>% fortify()
world2 <- ggplot2::map_data("world2")
p <- rasterVis::gplot(r) +
ggplot2::geom_polygon(
data = world2, aes(x=long, y=lat, group=group),
color = "black", fill = "darkgray", size=0.2) +
ggplot2::labs(x = "Longitude", y = "Latitude") +
ggsn::scalebar(
data = bb, location="bottomleft", dist=100, dd2km=T, model="WGS84", st.size=1) +
ggplot2::coord_fixed(xlim = b[1:2], ylim = b[3:4], ratio = 1) +
ggplot2::theme_bw() +
ggplot2::theme(
axis.title = element_text(size = ggplot2::rel(0.6)))
# TODO: add units to legend_title
#if (type == "depth" & element == "limits") browser()
if (element == "limits"){
brk_labels <<- nrel_limits[[type]]$break_labels
msg(g(" {type}: {paste(brk_labels, collapse=',')}"))
pal <- ifelse(type == "depth", "Blues", "Greens")
legend_title <- stringr::str_to_title(type)
p <- p +
ggplot2::geom_tile(aes(
fill=factor(value, levels=1:length(brk_labels), labels=brk_labels)), alpha=0.8) +
ggplot2::scale_fill_brewer(
palette = pal, name = legend_title, direction = 1, na.value = NA)
} else {
legend_title <- "Count"
p <- p +
ggplot2::geom_tile(aes(fill=value), alpha=0.8) +
ggplot2::scale_fill_distiller(
palette = "Spectral", name = legend_title, direction = -1, na.value = NA)
}
print(p)
}
get_fmt <- function(){
#rmd <- knitr::current_input()
pandoc_to <- opts_knit$get("rmarkdown.pandoc.to")
#if (!is.null(rmd)){
if (!is.null(pandoc_to)){
#fmt <- rmarkdown::default_output_format(rmd)$name
fmt <- switch(
pandoc_to,
html = "html_document",
latex = "pdf_document",
docx = "word_document")
} else {
fmt <- "html_document" # for interactive output
}
fmt
}
map_ter <- function(ter, type="all", element="count", redo_figs=F, return_null=F){
#msg(g(" {ter}_{type}_{element}_map.[png|pdf]"))
ter_info <- get_ter_info(ter)
if (is.na(ter_info[[type]][1])) return(NULL)
el <- ter_info[[type]][[element]]
map_pfx <- glue::glue("{dir_data}/territories/{ter}_{type}_{element}")
map_png <- glue::glue("{map_pfx}_map.png")
map_pdf <- glue::glue("{map_pfx}_map.pdf")
ter_info[[type]][[element]]$map_png <- map_png
ter_info[[type]][[element]]$map_pdf <- map_pdf
set_ter_info(ter_info)
fmt <- get_fmt()
b <- el$epsg4326_extent
dpi <- 300
w <- 6
h_max <- 8
h <- min(c(h_max, diff(b[3:4])/diff(b[1:2])*w))
msg(g(" map_ter({ter}, {type}, {element}); fmt={fmt}"))
if (!file.exists(map_png) | redo_figs){
message(glue::glue(" png({basename(map_png)})"))
png(filename=map_png, res=dpi, width=w*dpi, height=h*dpi)
gmap_ter(ter, type, element)
dev.off()
}
if (!file.exists(map_pdf) | redo_figs){
message(glue::glue(" pdf({basename(map_pdf)})"))
pdf(file=map_pdf, width=w, height=h)
gmap_ter(ter, type, element)
dev.off()
}
#browseURL(info$figure$png)
if (return_null) return(NULL)
res <- switch(
fmt,
html_document = lmap_ter(ter, type, element),
word_document = knitr::include_graphics(map_png),
pdf_document = knitr::include_graphics(map_pdf))
return(res)
}
get_lyr_info <- function(lyr){
dir_lyr <- glue::glue("{dir_lyrs}/{lyr}")
lyr_info_yml <- glue::glue("{dir_lyrs}/{lyr}/_{lyr}.yml")
if (!dir.exists(dir_lyr)){
dir.create(dir_lyr)
}
if (file.exists(lyr_info_yml)){
lyr_info <- yaml::read_yaml(lyr_info_yml)
} else {
lyr_info <- list(
key = lyr)
#browser()
lyr_info$territories <- setNames(
rep(list(NULL), length(eez_wgcs_sf$territory)),
eez_wgcs_sf$territory) %>% as.list()
yaml::write_yaml(lyr_info, lyr_info_yml)
}
lyr_info
}
set_lyr_info <- function(lyr_info){
lyr <- lyr_info$key
lyr_info_yml <- glue("{dir_lyrs}/{lyr}/_{lyr}.yml")
yaml::write_yaml(lyr_info, lyr_info_yml)
}
prep_lyr_ter <- function(lyr_p, lyr_ply, ter){
lyr <- lyr_p$key
lyr_info <- get_lyr_info(lyr)
lyr_ter_info <- lyr_info$territories[[ter]]
#browser()
if (is.null(lyr_ter_info) | lyr_p$redo){
# get eez
ter_eez_wgcs_sf <- get_ter_eez_wgcs_sf(ter) # ply_map(ter_eez_wgcs_sf, "territory")
# check if any intersection
ply <- sf_intersects(lyr_ply, ter_eez_wgcs_sf)
if (nrow(ply) == 0){
lyr_info$territories[[ter]] <- NA
set_lyr_info(lyr_info)
return(NULL)
}
# save ter lyr
lyr_ter_geo <- glue::glue("{dir_lyrs}/{lyr}/{ter}_{lyr}_epsg4326.geojson")
if (!file.exists(lyr_ter_geo)){
# ply <- sf_intersection(lyr_ply, ter_eez_wgcs_sf) # SLOW!
write_sf(ply, lyr_ter_geo)
} else {
ply <- read_sf(lyr_ter_geo)
}
# if (nrow(ply) == 0){
# lyr_info$territories[[ter]] <- NA
# set_lyr_info(lyr_info)
# } else {
# modify as needed
if (!is.na(lyr_p$mod_eval)){
ply <- eval(parse(text=lyr_p$mod_eval)) # ply_map(ply)
}
# get depth
ter_depth_wgcs_r <- get_ter_depth_wgcs_r(ter) # r_map(ter_depth_wgcs_r)
# convert to raster tif(s)
ply_to_tifs(ply, ter_depth_wgcs_r, ter, lyr, field=lyr_p$field, by=lyr_p$by)
}
#}
}
set_ter_info <- function(ter_info){
ter <- ter_info$territory
ter_info_yml <- glue("{dir_data}/territories/{ter}.yml")
yaml::write_yaml(ter_info, ter_info_yml)
}
get_ter_info <- function(ter){
ter_info_yml <- glue("{dir_data}/territories/{ter}.yml")
if (file.exists(ter_info_yml)){
ter_info <- yaml::read_yaml(ter_info_yml)
} else {
ter_info <- list(
territory = ter)
}
ter_info
}
make_ter_type_element <- function(ter, ter_s, ter_s_tbl, lyr_params, type, element){
# ter="Gulf of Mexico"; type="wind"; element="count"
msg(g(" make_ter_type_element({ter}, ..., {type}, {element})"))
ter_info <- get_ter_info(ter)
lyrs_not_count <- c(names(nrel_limits), "areakm2")
# skip if bin not in territory
if (type != "all"){
# TODO: make_ter_type_element(Alaska, ..., wave, limits) - 2018-01-04 16:32:25
# Quitting from lines 19-102 (nrel-uses.Rmd)
# Error in if (is.na(ter_info$layers[[type]])) { :
# argument is of length zero
# Calls: <Anonymous> ... eval -> eval -> make_ter_info -> make_ter_type_element
if (is.na(ter_info$layers[[type]])){
ter_info[[type]] <- NA
set_ter_info(ter_info)
msg(g(" skipping b/c {type} not in {ter}"))
return(ter_info)
}
}
# get raster of count or bin limits
if (element == "count"){
if (type == "all"){
msg(" calc all count")
ter_s_cnt <- raster::dropLayer(ter_s, which(names(ter_s) %in% lyrs_not_count))
# TODO: later, sum values after setting to 0 (but slow) vs assuming value of 1 if not NA
# ter_s_cnt[is.na(ter_s_cnt)] <- 0
r <- sum(!is.na(ter_s_cnt), na.rm=T) %>%
raster::mask(raster::raster(ter_s, "depth"))
} else {
# read raster from all for speedup. presumes make_ter_type_element() run for type='all' before others
r <- raster::raster(ter_info[["all"]][[element]]$epsg4326_tif)
}
} else {
# bin limits
r <- raster::raster(ter_s, type)
}
# filter to nrel_limits
msg(" filter to nrel_limits")
if (type != "all"){
r_depth <- raster::raster(ter_s, "depth")
# TODO: rerun with new is.na(r_bin) condition
if (type != "depth"){
r_bin <- raster::raster(ter_s, type)
r_na <- r_depth < nrel_limits[[type]]$depth$min |
r_depth > nrel_limits[[type]]$depth$max |
is.na(r_bin) |
r_bin < min(nrel_limits[[type]]$breaks)
r[r_na] <- NA
} else {
r[r_depth < min(nrel_limits$depth$breaks) |
r_depth > max(nrel_limits$depth$breaks)] <- NA
}
}
# assign breaks to bin limit rasters
msg(" assign breaks to bin limit rasters")
if (element == "limits"){
r <- raster::cut(r, breaks = nrel_limits[[type]]$breaks, include.lowest = T)
}
pfx <- glue::glue("{dir_data}/territories/{ter}_{type}_{element}")
# if all NA, set to NA, eg make_ter_type_element(Atlantic Islands, ..., tide, count)
if (is.na(raster::minValue(r)) & is.na(raster::maxValue(r))){
ter_info[[type]][[element]] = NA
set_ter_info(ter_info)
return(NULL)
}
msg(" trimming raster")
r <- raster_trim(r)
# register raster
msg(" register wrapped geographic raster")
tif_epsg4326 <- glue::glue("{pfx}_epsg4326.tif")
ter_info[[type]][[element]] <- list(
val_range = c(raster::minValue(r), raster::maxValue(r)),
epsg4326_tif = tif_epsg4326,
epsg4326_extent = raster::extent(r) %>% as.vector())
# write wrapped geographic raster
raster::writeRaster(r, tif_epsg4326, overwrite=T)
# handle mercator raster for leaflet
msg(" handle mercator raster for leaflet")
r_u <- raster_unwrap(r)
if (raster::extent(r_u)@xmin < -180){
# split into rasters around dateline and project for leaflet
tifs_epsg3857 <- glue::glue("{pfx}_epsg3857_{c('left','right')}.tif")
r_l <- raster::crop(r_u, raster::extent(-360,-180,-90,90), snap="in") %>% raster::shift(360) %>% raster_trim()
r_r <- raster::crop(r_u, raster::extent(-180, 180,-90,90), snap="in") %>% raster_trim()
r_l_epsg3857 <- raster_project_leaflet_nearest(r_l)
r_r_epsg3857 <- raster_project_leaflet_nearest(r_r)
# write rasters
raster::writeRaster(r_l_epsg3857, tifs_epsg3857[1], overwrite=T)
raster::writeRaster(r_r_epsg3857, tifs_epsg3857[2], overwrite=T)
# extent based on biggest longitudinal range for initial map display
e_l <- raster::extent(r_l)
e_r <- raster::extent(r_r)
x_l <- e_l %>% as.vector() %>% .[1:2] %>% diff()
x_r <- e_r %>% as.vector() %>% .[1:2] %>% diff()
b_epsg3857 <- c(e_l, e_r)[[which.max(c(x_l,x_r))]] %>% as.vector()
} else {
# project for leaflet
tifs_epsg3857 <- glue::glue("{pfx}_epsg3857.tif")
r_epsg3857 <- raster_project_leaflet_nearest(r_u)
# write raster
raster::writeRaster(r_epsg3857, tifs_epsg3857, overwrite=T)
# extent
b_epsg3857 <- raster::extent(r_u) %>% as.vector()
}
ter_info[[type]][[element]]$epsg3857_tif <- tifs_epsg3857
ter_info[[type]][[element]]$epsg3857_extent_init <- b_epsg3857
set_ter_info(ter_info)
# update map figures
msg(" update map figures")
map_ter(ter, type, element, redo_figs=T, return_null=F)
# table summaries
if (type != "all" & element == "count"){
msg(" table summaries")
bin <- type
ter_bin_lyr_csv <- glue("{pfx}_layers.csv")
ter_bin_cnt_csv <- glue("{pfx}_cumulative.csv")
ter_info[[type]][[element]]$layers_csv <- ter_bin_lyr_csv
ter_info[[type]][[element]]$cumulative_csv <- ter_bin_cnt_csv
tbl_b <- ter_bin(ter_s_tbl, bin)
tbl_bl <- ter_bin_lyr(tbl_b)
# devtools::load_all("../nrelutils")
tbl_bls <- ter_bin_lyr_smry(tbl_bl)
write_csv(tbl_bls, ter_bin_lyr_csv)
tbl_bc <- ter_bin_cnt(tbl_b)
write_csv(tbl_bc, ter_bin_cnt_csv)
}
set_ter_info(ter_info)
}
make_ter_info <- function(ter, lyr_params, bins){
ter_info <- get_ter_info(ter)
lyrs_not_count <- c(names(nrel_limits), "areakm2")
types <- c("all", bins)
# # quick fixes
# if ("raster" %in% names(ter_info)){
# ter_info$cnt <- ter_info$raster
# ter_info <- ter_info[-which(names(ter_info)=="raster")]
# names(ter_info$cnt)[which(names(ter_info$cnt)=="layers")] <- "layer_components"
# ter_info$cnt$layers_considered <- lyr_params$key[lyr_params$key != "aquaculture"]
# ter_info$cnt$layers_contained <- sort(unique(str_split(ter_info$cnt$layer_components, "_", simplify=T)[,1])) %>%
# setdiff("area")
# }
# if ("figure" %in% names(ter_info)){
# ter_info$cnt <- c(ter_info$cnt, ter_info$figure)
# ter_info <- ter_info[-which(names(ter_info)=="figure")]
# }
# yaml::write_yaml(ter_info, ter_info_yml)
# fetch and update territory if missing layers or bins
missing_layers <- setdiff(
lyr_params %>% filter(run==T) %>% .$key,
setdiff(names(ter_info$layers), lyr_params %>% filter(redo==T)))
# check that missing layers are available in
#browser()
for (lyr in missing_layers){ # lyr <- missing_layers[1]
lyr_in_ter <- !is.na(get_lyr_info(lyr)$territories[[ter]])
if (!lyr_in_ter){
msg(g(" missing or redo layer {lyr}, but not in {ter}, so ter_info$layers${lyr} = NA"))
ter_info$layers[[lyr]] <- NA
set_ter_info(ter_info)
ter_info <- get_ter_info(ter)
missing_layers <- setdiff(missing_layers, lyr)
}
}
missing_types <- setdiff(types, names(ter_info))
# TEMP: hard coded for new nrel bins
if (get_fmt()=="html_document") missing_types <- "wind"
# temp add map figures
# msg(" TEMP update map figures")
# map_ter(ter, "all", "count", redo_figs=T, return_null=T)
# for (bin in bins){
# map_ter(ter, bin, "count" , redo_figs=T, return_null=T)
# map_ter(ter, bin, "limits", redo_figs=T, return_null=T)
# }
if (length(missing_layers) > 0 | length(missing_types) > 0){
msg(g(" missing_layers: {paste(missing_layers, collapse=', ')}"))
msg(g(" missing_types : {paste(missing_types, collapse=', ')}"))
msg(" loading ter_s, ter_s_tbl")
ter_s <- ter_stack(ter, lyr_params)
ter_s_tbl <- ter_stack_tbl(ter_s)
# initiate layers
#if (!"layers" %in% names(ter_info)){
ter_info$layers <- attr(ter_s, "lyrs_in_ter")
ter_info$layers_not_count <- lyrs_not_count
set_ter_info(ter_info)
#}
if (length(missing_types) > 0 & length(missing_layers) == 0){
# if only missing type, redo needed type
for (type in missing_types){ # type = "wave"
make_ter_type_element(ter, ter_s, ter_s_tbl, lyr_params, type, "count")
if (type %in% bins){
make_ter_type_element(ter, ter_s, ter_s_tbl, lyr_params, type, "limits")
}
}
} else {
# otherwise if missing layer, redo all types
make_ter_type_element(ter, ter_s, ter_s_tbl, lyr_params, "all", "count")
for (bin in bins){
make_ter_type_element(ter, ter_s, ter_s_tbl, lyr_params, bin, "count")
}
}
# update layers into ter_info
ter_info <- get_ter_info(ter)
ter_info$layers <- attr(ter_s, "lyrs_in_ter")
ter_info$layers_not_count <- lyrs_not_count
set_ter_info(ter_info)
}
get_ter_info(ter)
}
fix_prep_lyr_ter <- function(){
# quick fix in prep_layers.R: from digest ({ter}_{lyr}_epsg4326.txt) to yaml (_{ter}.yml)
lyr_ter_txt <- glue("{dir_lyrs}/{lyr}/{ter}_{lyr}_epsg4326.txt")
if (file.exists(lyr_ter_txt)){
txt <- read_lines(lyr_ter_txt)
if (txt[1] == "0"){
lyr_info$territories[[ter]] <- NA
} else {
component_tif <- str_split(txt, ":", simplify=T)
component_tif <- setNames(component_tif[,2], component_tif[,1]) %>% as.list()
lyr_info$territories[[ter]]$components <- component_tif
}
# devtools::load_all("~/github/nrelutils")
set_lyr_info(lyr_info)
file.remove(lyr_ter_txt)
}
}
dt_lyrs_ter <- function(lyr_params){
fmt <- get_fmt()
lyrs_ter_tbl <- map_df(
1:nrow(lyr_params),
function(i){
lyr_p <- lyr_params[i,]
lyr <- lyr_p$key
lyrs_notna <- !is.na(get_lyr_info(lyr)$territories)
lyrs_checked <- ifelse(lyrs_notna, "✔", "")
tibble(
Source = lyr_p$source,
lyr = lyr,
Dataset = lyr_p$title) %>%
bind_cols(
lyrs_checked %>%
as.list() %>% as.tibble())
}) %>%
arrange(Source, Dataset) %>%
select(-lyr)
datatable(
lyrs_ter_tbl,
caption = "Existence of ocean use datasets across regions.",
options = list(
pageLength = nrow(lyrs_ter_tbl),
dom=c(
html_document = 'lftir',
pdf_document = 't',
word_document = 't')[[fmt]],
ordering=c(
html_document = TRUE,
pdf_document = FALSE,
word_document = FALSE)[[fmt]],
autoWidth = TRUE,
columnDefs = list(
list(className = 'dt-center', targets = 2:(ncol(lyrs_ter_tbl))))))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.