library(tidyverse)
library(lubridate)
library(readxl)
library(here)
here = here::here
library(glue)
library(fs)
library(sf)
library(mapview)
library(RColorBrewer)

library(dygraphs) # devtools::install_github("rstudio/dygraphs")
library(xts)
devtools::load_all(here("../infographiq"))

raw_csv   <- here("data/MARINe_raw_4c1e_9218_7d13.csv")
sites_csv <- here("data/MARINe_sites.csv")

# read_csv(raw_csv) %>% write_csv(here("data/MARINe_raw_4c1e_9218_7d13_head10k.csv"))

# https://www.eeb.ucsc.edu/pacificrockyintertidal/target/index.html
spp  <- tribble(
       ~sp,           ~sp_target,             ~sp_name,
  "MYTCAL",            "mytilus", "California Mussels",
  "CHTBAL", "chthamalus_balanus",    "Acorn Barnacles")

sanctuaries = c("cinms", "mbnms")

d_csv <- here("data/sanctuary_species_percentcover.csv")

# functions ----
get_nms_ply <- function(nms){
  # get polygon for National Marine Sanctuary
  
  nms_shp <- here(glue("data/shp/{nms}_py.shp"))
  
  if (!file.exists(nms_shp)){
    # download if needed
    
    # https://sanctuaries.noaa.gov/library/imast_gis.html
    nms_url <- glue("https://sanctuaries.noaa.gov/library/imast/{nms}_py2.zip")
    nms_zip <- here(glue("data/{nms}.zip"))
    shp_dir <- here("data/shp")
    
    download.file(nms_url, nms_zip)
    unzip(nms_zip, exdir = shp_dir)
    file_delete(nms_zip)
  }
  # read and convert to standard geographic projection
  read_sf(nms_shp) %>%
    st_transform(4326)
}

plot_intertidal_nms <- function(d_csv, NMS, spp, sp_name){

  # read in csv with fields site, date, pct_cover
  d <- read_csv(d_csv) %>%
    filter(nms==NMS, sp==spp) %>%
    select(-nms, -sp) %>%
    spread(site, pct_cover) # View(d_sites)
  
  # line colors
  ln_colors <- c(colorRampPalette(brewer.pal(11, "Set3"))(ncol(d)-2), "black")
  
  # convert to xts time object
  x <- select(d, -date) %>%
    as.xts(order.by=d$date)
  
  # plot dygraph
  #browser()
  dygraph(x, main=glue("{sp_name} in {NMS}")) %>%
    dyOptions(
      connectSeparatedPoints = TRUE,
      colors = ln_colors) %>%
    dySeries(NMS, strokeWidth = 3) %>%
    dyHighlight(highlightSeriesOpts = list(strokeWidth = 2)) %>%
    dyRangeSelector()
}

map_nms_sites <- function(nms){
  # nms <- "mbnms"
  NMS <- str_to_upper(nms)
  
  # get sites in nms
  sites_nms_shp <- here(glue("data/shp/{NMS}_sites.shp"))
  nms_ply <- get_nms_ply(nms)
  
  if (!file.exists(sites_nms_shp)){
    sites_nms_pts <- sites_pts %>%
      st_intersection(nms_ply)
    write_sf(sites_nms_pts, sites_nms_shp)
  }
  sites_nms_pts <- read_sf(sites_nms_shp)
  
  mapview(
    nms_ply, legend = TRUE, layer.name = "Sanctuary", zcol = "SANCTUARY") + 
    mapview(
      sites_nms_pts, legend = TRUE, layer.name = "Site",
      zcol = "site", col.regions = colorRampPalette(brewer.pal(11, "Set3")))
}
raw     <- read_csv(raw_csv)

sites_pts <- raw %>%
  rename(
    site = marine_site_name) %>%
  group_by(site) %>%
  summarize(
    lat = first(`latitude (degrees_north)`),
    lon = first(`longitude (degrees_east)`)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

sites_pts %>%
  mutate(
    lon = st_coordinates(geometry)[,1],
    lat = st_coordinates(geometry)[,2]) %>%
  st_set_geometry(NULL) %>%
  write_csv(sites_csv)

Rocky Intertidal Scene

Rocky Intertidal Scene

CINMS

map_nms_sites("cinms")

Acorn Barnacles

plot_intertidal_nms(d_csv, "CINMS", "CHTBAL", "Acorn Barnacles")

California Mussels

plot_intertidal_nms(d_csv, "CINMS", "MYTCAL", "California Mussels")

MBNMS

map_nms_sites("mbnms")

Acorn Barnacles

plot_intertidal_nms(d_csv, "MBNMS", "CHTBAL", "Acorn Barnacles")

California Mussels

plot_intertidal_nms(d_csv, "MBNMS", "MYTCAL", "California Mussels")

Longitudinal Gradient

sites_tbl <- read_csv(sites_csv)

#plot_latitudinal_nms <- function(d_csv, NMS, spp, sp_name){
NMS="MBNMS"; spp="MYTCAL"; sp_name="California Mussels"

  # read in csv with fields site, date, pct_cover
  d <- read_csv(d_csv) %>%
    filter(nms==NMS, sp==spp) %>%
    select(-nms, -sp) %>%
    mutate(
      year = year(date)) %>%
    group_by(site, year) %>%
    summarize(
      pct_cover = mean(pct_cover)) %>%
    left_join(
      sites_tbl, by="site") %>%
    select(site, year, pct_cover, lon) %>%
    arrange(lon, year, pct_cover) # View(d)
  
  g <- ggplot(d, aes(x=lon, y=pct_cover, group=year, color=year)) +
    geom_line() + 
    scale_color_gradientn(
      colors=RColorBrewer::brewer.pal(11, "Spectral"),
      name = "Year") +
    coord_flip() +
    labs(
      title = glue("{sp_name} in {NMS}"),
      x = "Longitude", y = "% cover")
  
  plotly::ggplotly(g)
    #scale_x_log10()
  
#}
for (i in 1:length(sanctuaries)){
  
  # set sanctuary variables
  nms <- sanctuaries[i] # nms <- "mbnms"
  NMS <- str_to_upper(nms)
  
  # get sites in nms
  sites_nms_shp <- here(glue("data/shp/{NMS}_sites.shp"))
  if (!file.exists(sites_nms_shp)){
    nms_ply <- get_nms_ply(nms)
    sites_nms_pts <- sites_pts %>%
      st_intersection(nms_ply)
    write_sf(sites_nms_pts, sites_nms_shp)
  }
  sites_nms_pts <- read_sf(sites_nms_shp)
  
  # plot map of sanctuary and sites
  #m <- mapview(nms_ply) + sites_nms_pts
  #print(m)

  # iterate over species
  for (j in 1:nrow(spp)){
    
    # set species variables
    sp        <- spp$sp[j]
    sp_target <- spp$sp_target[j]
    sp_name   <- spp$sp_name[j]
    
    # average across plots for each site-species-date
    d_sites <- raw %>%
      rename(
        site = marine_site_name) %>%
      filter(
        site %in% sites_nms_pts$site,
        lumping_code      == sp,
        target_assemblage == sp_target) %>%
      mutate(
        date = ymd(`time (UTC)`)) %>%
      group_by(site, date) %>%
      summarize(
        pct_cover = mean(percent_cover)) # View(d_sites)
    
    # average across sites to nms-year
    d_nms <- d_sites %>%
      mutate(
        date = date(glue("{year(date)}-06-15"))) %>%
      group_by(date) %>%
      summarize(
        pct_cover = mean(pct_cover)) %>%
      ungroup() %>%
      mutate(
        site = NMS) # View(d_nms)
    
    # combine sites and sanctuary annual average
    d <- bind_rows(
      d_sites,
      d_nms) %>%
      mutate(
        nms = NMS,
        sp  = sp) # View(d)

    # write data to csv
    if (i == 1 & j == 1){
      write_csv(d, d_csv)
    } else {
      write_csv(d, d_csv, append=T)
    }
    
    # generate timeseries plot
    #p <- plot_intertidal_nms(d_csv, NMS, sp_name)
    #print(p)
  }
}
devtools::load_all(here("../infographiq"))
pwd <- here()
create_info_site(
  site_title       = "MARINe MBNMS",
  path_root        =  pwd,
  plot_function    = "plot_intertidal_nms",
  dir_svg          = 'svg',
  elements_csv     = 'svg_elements.csv',
  indicators_csv   = 'plot_indicators.csv',
  dir_rmd          = 'rmd',
  dir_web          = 'docs',
  #svg_paths        = list.files(file.path(path_root, dir_svg), '.*\\.svg$', full.names=T),
  #svg_names        = tools::file_path_sans_ext(basename(svg_paths)),
  site_yml_brew    = system.file('site_template/_site.yml.brew', package='infographiq'),
  index_md_brew    = system.file('site_template/index.md.brew', package='infographiq'),
  readme_md_brew   = system.file('site_template/README.md.brew', package='infographiq'),
  header           = system.file('site_template/_header.html', package='infographiq'),
  footer           = system.file('site_template/_footer.html', package='infographiq'),
  styles_css       = system.file('site_template/styles.css', package='infographiq'),
  render_modals    = T,
  preview_site     = T)