trends.Rmd

---
title: Piney Point trends synthesis
author: "MW Beck"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  html_document
---

```{r setup, echo = FALSE, message = F, warning = F, results = 'hide'}

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, fig.height = 6, fig.width = 9)

library(tidyverse)
library(lubridate)
library(sf)
library(ggmap)
library(ggspatial)
library(patchwork)
library(units)
library(grid)
library(scales)
library(here)
library(ggfx)
library(DescTools)

data(rswqdat)
data(bswqdat)
data(bsstatloc)
data(rsstatloc)
data(ppseg)
data(parms)
data(rsphydat)
data(rsphypts)
data(trnpts)
data(rstrndat)
data(rstrnpts)
data(bstransect)
data(rstrnwts)
data(rsallpts)
data(segmask)

source("R/funcs.R")

sf_use_s2(FALSE)

# segments
ppseg <- ppseg %>% 
  rename(area = Name) %>% 
  group_by(area) %>% 
  summarise() %>% 
  st_buffer(dist = set_units(0.0001, degree)) %>% 
  st_buffer(dist = set_units(-0.0001, degree)) %>% 
  mutate(
    area = factor(area)
  )

cols <- c("#E16A86", "#50A315", "#009ADE")
names(cols) <- levels(ppseg$area)

# nonbay
nonbay <- c('BH01', 'P Port 2', 'P Port 3', 'PM Out', '20120409-01', 'PPC41', 'P Port 4', 'PMB01', 'NGS-S Pond')
```

# {.tabset}

## Region delineation 

```{r, fig.height = 4, fig.width = 4, fig.margin = T}
##
# map

buffdist <- 0.01
northloc <- 'tr' 
scaleloc <- 'bl'

# layer extent as bbox plus buffer
dat_ext <- ppseg %>% 
  st_bbox %>% 
  st_as_sfc %>% 
  st_buffer(dist = set_units(buffdist, degree)) %>%
  st_bbox %>% 
  unname

# reference data for ggsn, MUST have geometry named column
ggsnref <- ppseg %>% 
  st_bbox %>% 
  st_as_sfc %>%
  st_buffer(dist = set_units(buffdist / 2, degree)) %>% 
  st_as_sf %>%
  st_cast('POINT') %>% 
  rename(geometry = x)

# stamen base map
bsmap1 <- get_stamenmap(bbox = dat_ext, maptype = 'toner-background', zoom = 10)

# change opacity of basemap
mapatt <- attributes(bsmap1)
bsmap1_transparent <- matrix(adjustcolor(bsmap1, 
                                         alpha.f = 0.2), 
                             nrow = nrow(bsmap1))
attributes(bsmap1_transparent) <- mapatt

# plot
map <- ggmap(bsmap1_transparent) +
  geom_sf(data = ppseg, aes(fill = area), color = 'black', inherit.aes = F, alpha = 0.8) +
  geom_sf_text(data = ppseg, aes(label = area), color = 'black', inherit.aes = F, alpha = 0.8, size = 6) +
  scale_fill_manual(values = cols, drop = F) +
  theme(
    legend.title = element_blank(), 
    panel.grid = element_blank(), 
    axis.title = element_blank(), 
    legend.justification = 'top',
    axis.text.y = element_text(size = 8), 
    axis.text.x = element_text(size = 8, angle = 30, hjust = 1),
    panel.background = element_rect(fill = 'white'),
    axis.ticks = element_line(colour = 'grey'),
    panel.border = element_rect(colour = 'grey', fill = NA), 
    legend.position = 'none'
  ) +
  annotation_scale(location = scaleloc) +
  annotation_north_arrow(location = northloc, which_north = "true", height = grid::unit(0.75, "cm"), 
                         width = grid::unit(0.75, "cm"))

map
```

## Sampling maps 

```{r, message = F, warning = F}
# piney point loc
pineypoint <- tibble(
  lon = -82.52469352586753, 
  lat = 27.629819505234703
) %>% 
  st_as_sf(coords = c('lon', 'lat'), crs = 4326)

# segments
areas <- ppseg %>% 
  st_intersection(segmask) %>% 
  st_cast('POLYGON') %>% 
  mutate(
    acres = st_area(.), 
    acres = set_units(acres, 'acres'), 
    acres = as.numeric(acres)
  ) %>% 
  dplyr::filter(acres > 1e4) %>% 
  st_buffer(dist = set_units(0.0005, degree)) %>% 
  st_buffer(dist = set_units(-0.0005, degree)) 

##
# map

buffdist <- 0.01
northloc <- 'tr' 
scaleloc <- 'bl'

# layer extent as bbox plus buffer
dat_ext <- areas %>% 
  st_bbox %>% 
  st_as_sfc %>% 
  st_buffer(dist = set_units(buffdist, degree)) %>%
  st_bbox %>% 
  unname

# reference data for ggsn, MUST have geometry named column
ggsnref <- areas %>% 
  st_bbox %>% 
  st_as_sfc %>%
  st_buffer(dist = set_units(buffdist / 2, degree)) %>% 
  st_as_sf %>%
  st_cast('POINT') %>% 
  rename(geometry = x)

# stamen base map
bsmap1 <- get_stamenmap(bbox = dat_ext, maptype = 'toner-background', zoom = 10)

# change opacity of basemap
mapatt <- attributes(bsmap1)
bsmap1_transparent <- matrix(adjustcolor(bsmap1, 
                                         alpha.f = 0.2), 
                             nrow = nrow(bsmap1))
attributes(bsmap1_transparent) <- mapatt

# samples by agency

tomap <- rsallpts %>% 
  filter(grepl('water quality', type)) %>% 
  mutate(
    source_lng = case_when(
      source_lng %in% c('City of St. Pete', 'ESA', 'Florida DEP', 'UF', 'USF', 'New College Fl.') ~ 'other',
      T ~ source_lng
    )
  ) %>% 
  .[areas, ]

tomap1 <- tomap %>% 
  filter(source_lng != 'other')
tomap2 <- tomap %>%
  filter(source_lng == 'other')

p1a <- ggmap(bsmap1_transparent) +
  geom_sf(data = tomap1, aes(color = source_lng), inherit.aes = F) +
  geom_sf(data = pineypoint, aes(fill = 'Piney Point', shape = 'Piney Point'), color = 'black', size = 3, inherit.aes= F) + 
  scale_shape_manual(values= c(24)) +
  scale_fill_manual(values = c('black')) +
  theme_bw() + 
  theme(
    panel.grid = element_blank(), 
    axis.title = element_blank(), 
    legend.justification = 'top',
    axis.text.y = element_text(size = 8), 
    axis.text.x = element_text(size = 8, angle = 30, hjust = 1),
    panel.background = element_rect(fill = 'white'),
    axis.ticks = element_line(colour = 'grey'),
    panel.border = element_rect(colour = 'grey', fill = NA), 
    legend.position = 'right',
    plot.caption = element_text(size = 6)
  ) +
  labs(
    color = 'Monitoring group',
    fill = NULL, 
    shape = NULL, 
    subtitle = "Water quality stations", 
    caption = ' '
  ) +
  guides(color = guide_legend(override.aes = list(size = 3))) + 
  annotation_north_arrow(location = northloc, which_north = "true", height = grid::unit(0.75, "cm"), 
                         width = grid::unit(0.75, "cm")) +
  annotation_scale(location = scaleloc)

p2a <- ggmap(bsmap1_transparent) +
  geom_sf(data = tomap2, aes(color = source_lng), inherit.aes = F) +
  geom_sf(data = pineypoint, aes(fill = 'Piney Point', shape = 'Piney Point'), color = 'black', size = 3, inherit.aes= F) + 
  scale_shape_manual(values = c(24)) +
  scale_fill_manual(values = c('black')) +
  theme_bw() + 
  theme(
    panel.grid = element_blank(), 
    axis.title = element_blank(), 
    legend.justification = 'top',
    axis.text.y = element_text(size = 8), 
    axis.text.x = element_text(size = 8, angle = 30, hjust = 1),
    panel.background = element_rect(fill = 'white'),
    axis.ticks = element_line(colour = 'grey'),
    panel.border = element_rect(colour = 'grey', fill = NA), 
    legend.position = 'right',
    plot.caption = element_text(size = 6)
  ) +
  labs(
    color = 'Monitoring group',
    fill = NULL, 
    shape = NULL, 
    subtitle = "Water quality stations", 
    caption = 'Other includes City of St. Pete, ESA, Florida DEP, UF, USF, New College Fl'
  ) +
  guides(color = guide_legend(override.aes = list(size = 3))) + 
  annotation_north_arrow(location = northloc, which_north = "true", height = grid::unit(0.75, "cm"), 
                         width = grid::unit(0.75, "cm")) +
  annotation_scale(location = scaleloc)

tomap <- rsallpts %>% 
  filter(grepl('seagrass and macroalgae', type)) %>% 
  .[areas, ]

p3a <- ggmap(bsmap1_transparent) +
  geom_sf(data = tomap, aes(color = source_lng), inherit.aes = F) +
  geom_sf(data = pineypoint, aes(fill = 'Piney Point', shape = 'Piney Point'), color = 'black', size = 3, inherit.aes= F) + 
  scale_shape_manual(values= c(24)) +
  scale_fill_manual(values = c('black')) +
  theme_bw() + 
  theme(
    panel.grid = element_blank(), 
    axis.title = element_blank(), 
    legend.justification = 'top',
    axis.text.y = element_text(size = 8), 
    axis.text.x = element_text(size = 8, angle = 30, hjust = 1),
    panel.background = element_rect(fill = 'white'),
    axis.ticks = element_line(colour = 'grey'),
    panel.border = element_rect(colour = 'grey', fill = NA), 
    legend.position = 'right', 
    plot.caption = element_text(size = 6)
  ) +
  labs(
    color = 'Monitoring group',
    fill = NULL, 
    shape = NULL, 
    subtitle = "Seagrass and macroalgae transects", 
    caption = ' '
  ) +
  guides(color = guide_legend(override.aes = list(size = 3))) + 
  annotation_north_arrow(location = northloc, which_north = "true", height = grid::unit(0.75, "cm"), 
                         width = grid::unit(0.75, "cm")) +
  annotation_scale(location = scaleloc)

tomap <- rsallpts %>% 
  filter(grepl('seagrass and macroalgae', type)) %>% 
  .[areas, ]

ppsegint <- ppseg %>% 
  st_intersection(segmask)

p4a <- ggmap(bsmap1_transparent) +
  geom_sf(data = ppsegint, aes(fill = area), inherit.aes = F, color = NA) +
  geom_sf(data = tomap, aes(shape = 'Transect locations'), color = 'black', inherit.aes = F) +
  geom_sf(data = pineypoint, aes(shape = 'Piney Point'), color = 'black', size = 3, inherit.aes= F) + 
  scale_shape_manual(values= c(17, 16)) +
  scale_fill_manual(values = c(cols)) +
  theme_bw() + 
  theme(
    panel.grid = element_blank(), 
    axis.title = element_blank(), 
    legend.justification = 'top',
    axis.text.y = element_text(size = 8), 
    axis.text.x = element_text(size = 8, angle = 30, hjust = 1),
    panel.background = element_rect(fill = 'white'),
    axis.ticks = element_line(colour = 'grey'),
    panel.border = element_rect(colour = 'grey', fill = NA), 
    legend.position = 'right', 
    plot.caption = element_text(size = 6)
  ) +
  labs(
    color = 'Monitoring group',
    fill = NULL, 
    shape = NULL, 
    subtitle = "Seagrass and macroalgae transects by area", 
    caption = ' '
  ) +
  guides(color = guide_legend(override.aes = list(size = 3))) + 
  annotation_north_arrow(location = northloc, which_north = "true", height = grid::unit(0.75, "cm"), 
                         width = grid::unit(0.75, "cm")) +
  annotation_scale(location = scaleloc)
```

```{r, message = F, warning = F, results = 'hide'}
jpeg(here('figure/mapwq.jpeg'), height = 4, width = 5, family = 'serif', units = 'in', res = 500)
print(p1a)
dev.off()
```

```{r, message = F, warning = F, results = 'hide'}
jpeg(here('figure/mapwqother.jpeg'), height = 4, width = 5, family = 'serif', units = 'in', res = 500)
print(p2a)
dev.off()
```

```{r, message = F, warning = F, results = 'hide'}
jpeg(here('figure/mapsg.jpeg'), height = 4, width = 5, family = 'serif', units = 'in', res = 500)
print(p3a)
dev.off()
```

```{r, message = F, warning = F, results = 'hide'}
jpeg(here('figure/mapsgarea.jpeg'), height = 4, width = 5, family = 'serif', units = 'in', res = 500)
print(p4a)
dev.off()
```

[download](figure/mapwq.jpeg)
```{r, out.width = '70%'}
knitr::include_graphics(here('figure/mapwq.jpeg'))
```

[download](figure/mapwqother.jpeg)
```{r, out.width = '70%'}
knitr::include_graphics(here('figure/mapwqother.jpeg'))
```

[download](figure/mapsg.jpeg)
```{r, out.width = '70%'}
knitr::include_graphics(here('figure/mapsg.jpeg'))
```

[download](figure/mapsgarea.jpeg)
```{r, out.width = '70%'}
knitr::include_graphics(here('figure/mapsgarea.jpeg'))
```

## Water Quality {.tabset .tabset-pills}

### TN

```{r}
wqplo_fun(rswqdat, bswqdat, ppseg, vr = 'tn', cols, logtr = TRUE)
```

### NH3

```{r}
wqplo_fun(rswqdat, bswqdat, ppseg, vr = 'nh34', cols, logtr = TRUE)
```

### TP

```{r}
wqplo_fun(rswqdat, bswqdat, ppseg, vr = 'tp', cols, logtr = TRUE)
```

### Chl

```{r}
wqplo_fun(rswqdat, bswqdat, ppseg, vr = 'chla', cols, logtr = TRUE)
```

### DO Sat.

```{r}
wqplo_fun(rswqdat, bswqdat, ppseg, vr = 'dosat', cols, logtr = FALSE)
```

### Aug/Sep 2021 v 2022

```{r, fig.height = 3, fig.width = 5}

nonbay <- c('BH01', 'P Port 2', 'P Port 3', 'PM Out', '20120409-01', 'PPC41', 'P Port 4', 'PMB01', 'NGS-S Pond')

# monitoring data
rswqtmp <- rswqdat %>% 
  filter(var %in% c('chla', 'tn', 'secchi')) %>% 
  filter(!station %in% nonbay) %>% 
  inner_join(rsstatloc, ., by = c('station', 'source')) %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  select(-qual, -bswqstation, -nrmrng, -source, -source_lng, -uni, -lbunis) %>% 
  filter(area == 'Area 1') %>% 
  filter(month(date) %in% c(8)) %>% 
  mutate(
    yr = year(date)
  ) %>% 
  group_by(yr, var) %>% 
  mutate(
    avev = case_when(
      var %in% c('chla', 'tn') ~ Gmean(val), 
      T ~ mean(val)
    ),
    valhi = case_when(
      var %in% c('chla', 'tn') ~ as.numeric(Gmean(val, conf.level = 0.95)[3]),
      T ~ as.numeric(t.test(val)$conf.int[2])
    ),
    vallo = case_when(
      var %in% c('chla', 'tn') ~ as.numeric(Gmean(val, conf.level = 0.95)[2]),
      T ~ as.numeric(t.test(val)$conf.int[1])
    )
  ) %>% 
  ungroup() %>% 
  select(yr, var, avev, vallo, valhi) %>% 
  unique() %>% 
  mutate(
    var = factor(var, levels = c('tn', 'chla', 'secchi'),
                 labels = c(
                   'Total~Nitrogen~(mg~L^{-1})',
                   'Chlorophyll-a~(ug~L^{-1})',
                   'Secchi~depth~(m)'
                   )
                 )
  )

ggplot(rswqtmp,  aes(x = factor(yr), y = avev, fill = factor(yr))) + 
  with_shadow(geom_bar(stat = 'identity'), sigma = 2.7, x_offset = 0, y_offset = 0) + 
  geom_errorbar(aes(ymin = vallo, ymax = valhi), width = 0) + 
  facet_wrap(~var, scales = 'free_y', labeller = label_parsed, strip.position = 'left') + 
  theme_minimal(base_size = 12) + 
  scale_fill_manual(values = c('#004F7E', '#5C4A42')) +
  theme(
    legend.position = 'top',
    strip.background = element_blank(), 
    strip.text = element_text(size = rel(1)), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_blank(), 
    panel.grid.minor = element_blank(),
    strip.placement = 'outside',
    panel.background = element_blank(), #element_rect(fill = 'grey95', color = 'white'), 
    panel.grid.major = element_blank(), 
    axis.text.x = element_blank()
  ) +
  labs(
    x = NULL,
    y = NULL, 
    fill = NULL,
    title = 'Water Quality in August'
  )
```

## Phytoplankton {.tabset .tabset-pills}

### Mixed algae

```{r}
phyplo_fun(rsphydat, rsphypts, ppseg, vr = 'mixed algae', cols)
```

### Centric diatoms

```{r}
phyplo_fun(rsphydat, rsphypts, ppseg, vr = 'Centric Diatoms', cols)
```

### Bacillariophyta

```{r}
phyplo_fun(rsphydat, rsphypts, ppseg, vr = 'Bacillariophyta', cols)
```

### Karenia brevis

```{r}
phyplo_fun(rsphydat, rsphypts, ppseg, vr = 'Karenia brevis', cols)
```

### Pseudo-nitzschia sp.

```{r}
phyplo_fun(rsphydat, rsphypts, ppseg, vr = 'Pseudo-nitzschia sp.', cols)
```

### Rhizosolenia setigera

```{r}
phyplo_fun(rsphydat, rsphypts, ppseg, vr = 'Rhizosolenia setigera', cols)
```

## Seagrass {.tabset .tabset-pills}

### Overall

Note that the "total" category applies if any seagrass or macroalgae were present along a transect and is not the sum of the individual frequency occurrences.  

```{r}
mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")
savsel <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme')

colpal <- colorRampPalette(RColorBrewer::brewer.pal(n = 8, name = 'Dark2'))
savlevs <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme', 'Ruppia maritima', 'Halophila engelmannii', 'Halophila decipiens')
savcol <- colpal(length(savlevs))
names(savcol) <- savlevs
savcol <- savcol[savsel]
mcrcol <- c('tomato1', 'lightgreen', 'burlywood3', 'lightblue')
names(mcrcol) <- mcrsel
mcrcol <- mcrcol[mcrsel]
cols <- c(mcrcol, savcol)

cols <- c(cols, Total = 'white')

# add area
trnsum <- rstrndat %>%
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>%
  dplyr::group_by(area, typ, date, station, taxa, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(area, typ, date, taxa) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  filter(taxa %in% c(mcrsel, savsel))

trnsumtots <- rstrndat %>%
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>%
  dplyr::group_by(area, typ, date, station, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(area, typ, date) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  mutate(taxa = 'Total')

toplo <- bind_rows(trnsum, trnsumtots)

dodge <- position_dodge(width=0) 

ggplot(toplo, aes(x = date, y = foest)) + 
  geom_line(aes(group = taxa), position = dodge) +
  geom_point(aes(fill = taxa, group = taxa), pch = 21, stat = 'identity', color = 'black', size = 4, position = dodge, stroke = 1.5) +
  facet_grid(typ ~ area) +
  theme_minimal(base_size = 14) + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b') + 
  scale_y_continuous(limits = c(0, 1))+
  scale_fill_manual(values = cols) +
  labs(
    y = 'Freq. occurrence'
  ) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_line(), 
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9), 
    # panel.spacing=unit(2, "lines"), 
    panel.background = element_rect(fill = 'grey95', color = 'white'), 
    panel.grid.major = element_line(color = 'grey90')
  )
```

The same plot as above but aggregated across all transects. 

```{r}
mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")
savsel <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme')

colpal <- colorRampPalette(RColorBrewer::brewer.pal(n = 8, name = 'Dark2'))
savlevs <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme', 'Ruppia maritima', 'Halophila engelmannii', 'Halophila decipiens')
savcol <- colpal(length(savlevs))
names(savcol) <- savlevs
savcol <- savcol[savsel]
mcrcol <- c('tomato1', 'lightgreen', 'burlywood3', 'lightblue')
names(mcrcol) <- mcrsel
mcrcol <- mcrcol[mcrsel]
cols <- c(mcrcol, savcol)

cols <- c(cols, Total = 'white')

# add area
trnsum <- rstrndat %>%
  dplyr::group_by(typ, date, station, taxa, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(typ, date, taxa) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  filter(taxa %in% c(mcrsel, savsel))

trnsumtots <- rstrndat %>%
  dplyr::group_by(typ, date, station, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(typ, date) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  mutate(taxa = 'Total')

toplo <- bind_rows(trnsum, trnsumtots)

dodge <- position_dodge(width=0) 

ggplot(toplo, aes(x = date, y = foest)) + 
  geom_line(aes(group = taxa), position = dodge) +
  geom_point(aes(fill = taxa, group = taxa), pch = 21, stat = 'identity', color = 'black', size = 4, position = dodge, stroke = 1.5) +
  facet_grid(typ ~ .) +
  theme_minimal(base_size = 14) + 
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_date(date_breaks = '1 month', date_labels = '%b') + 
  scale_fill_manual(values = cols) +
  labs(
    y = 'Freq. occurrence'
  ) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14), 
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_line(), 
    panel.grid.minor = element_blank(),
    panel.background = element_rect(fill = 'grey95', color = 'white'), 
    panel.grid.major = element_line(color = 'grey90')
  )
```

Macroalgae only. 

```{r, message = F, results = 'hide', warning = F}
mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")

colpal <- colorRampPalette(RColorBrewer::brewer.pal(n = 8, name = 'Dark2'))
mcrcol <- c('tomato1', 'lightgreen', 'burlywood3', 'lightblue')
names(mcrcol) <- mcrsel
mcrcol <- mcrcol[mcrsel]
cols <- c(mcrcol)

cols <- c(cols, Total = 'white')

# add area
trnsum <- rstrndat %>%
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>%
  dplyr::group_by(area, typ, date, station, taxa, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(area, typ, date, taxa) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  filter(taxa %in% c(mcrsel))

trnsumtots <- rstrndat %>%
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>%
  dplyr::group_by(area, typ, date, station, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(area, typ, date) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  mutate(taxa = 'Total') %>% 
  filter(typ == 'Macroalgae')

toplo <- bind_rows(trnsum, trnsumtots)

dodge <- position_dodge(width=0) 

p <- ggplot(toplo, aes(x = date, y = foest)) + 
  geom_line(aes(group = taxa), position = dodge) +
  geom_point(aes(fill = taxa, group = taxa), pch = 21, stat = 'identity', color = 'black', size = 4, position = dodge, stroke = 1.5) +
  facet_grid(~ area) +
  theme_minimal(base_size = 14) + 
  scale_x_date(date_breaks = '1 month', date_labels = '%b') + 
  scale_y_continuous(limits = c(0, 1))+
  scale_fill_manual(values = cols) +
  labs(
    y = 'Freq. occurrence'
  ) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_line(), 
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9), 
    panel.grid.minor = element_blank(),
    # panel.spacing=unit(2, "lines"), 
    panel.background = element_rect(fill = 'grey95', color = 'white'), 
    panel.grid.major = element_line(color = 'grey90')
  )

jpeg(here('figure/macrofo.jpeg'), height = 4, width = 9, units = 'in', res = 500, family = 'serif')
print(p)
dev.off()
```

```{r}
knitr::include_graphics(here('figure/macrofo.jpeg'))
```

### SAV vs macroalgae, by station

```{r}
mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")
savsel <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme')

# add area
savsum <- rstrndat %>%
  filter(typ == 'sav') %>% 
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  mutate(
    date = floor_date(date, unit = 'week')
  ) %>% 
  dplyr::group_by(area, date, station, taxa) %>%
  dplyr::summarise(
    nsites = length(unique(location)),
    bbest = sum(bb, na.rm = T) / nsites
  ) %>% 
  filter(taxa %in% savsel) %>% 
  group_by(area, date, station) %>% 
  summarize(
    savbbest = mean(bbest)
  )

mcrsum <- rstrndat %>%
  filter(typ == 'mcr') %>% 
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  mutate(
    date = floor_date(date, unit = 'week')
  ) %>% 
  dplyr::group_by(area, date, station, taxa) %>%
  dplyr::summarise(
    nsites = length(unique(location)),
    bbest = sum(bb, na.rm = T) / nsites
  ) %>% 
  filter(taxa %in% mcrsel) %>% 
  group_by(area, date, station) %>% 
  summarize(
    mcrbbest = mean(bbest)
  )

toplo <- full_join(mcrsum, savsum, by = c('area', 'station', 'date')) %>% 
  mutate(
    mo = month(date,label = T)
  )

ggplot(toplo, aes(x = mcrbbest, y = savbbest)) + 
  geom_point(aes(fill = mo), colour = 'grey', pch = 21, size = 3) +
  facet_grid(~ area, scales = 'free') +
  theme_minimal(base_size = 14) + 
  geom_smooth() +
  coord_cartesian(ylim = c(0, NA)) +
  labs(
    y = 'SAV mean abundance (bb)',
    x = 'Macroalgae mean abundance (bb)', 
    caption = 'Results retained by transect'
  ) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14)
  )
```

### Water quality vs macroalgae vs SAV

```{r}
# add area
savsum <- rstrndat %>%
  filter(typ == 'sav') %>%  
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  # mutate(
  #   date = floor_date(date, unit = 'week')
  # ) %>% 
  dplyr::group_by(area, date, station, taxa) %>%
  dplyr::summarise(
    nsites = length(unique(location)),
    bbest = sum(bb, na.rm = T) / nsites
  ) %>% 
  filter(taxa %in% savsel) %>% 
  dplyr::group_by(date, area) %>%
  dplyr::summarise(
    savbbest = mean(bbest, na.rm = T),
    .groups = 'drop'
  ) 

mcrsum <- rstrndat %>%
  filter(typ == 'mcr') %>% 
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  # mutate(
  #   date = floor_date(date, unit = 'week')
  # ) %>% 
  dplyr::group_by(area, date, station, taxa) %>%
  dplyr::summarise(
    nsites = length(unique(location)),
    bbest = sum(bb, na.rm = T) / nsites
  ) %>% 
  filter(taxa %in% mcrsel) %>% 
  dplyr::group_by(date, area) %>%
  dplyr::summarise(
    mcrbbest = mean(bbest, na.rm = T),
    .groups = 'drop'
  ) 

rswqtmp <- rswqdat %>% 
  filter(var %in%  c('tn', 'nh34', 'chla', 'tp')) %>% 
  filter(!station %in% nonbay) %>% 
  inner_join(rsstatloc, ., by = c('station', 'source')) %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  select(-qual, -bswqstation, -nrmrng, -source, -source_lng, -uni, -lbunis) %>% 
  # mutate(
  #   date = floor_date(date, unit = 'week')
  # ) %>% 
  group_by(date, area, var) %>% 
  summarise(
    val = mean(log(val)), 
    val = exp(val)
    ) %>% 
  spread(var, val)

toplo <- full_join(mcrsum, savsum, by = c('area', 'date')) %>% 
  inner_join(rswqtmp, by = c('area', 'date')) %>% 
  mutate(
    mo = month(date,label = T)
  )

bssz <- 14

mod1 <- summary(lm(chla ~ tn, toplo))
txt <- paste0('Slope ', round(coefficients(mod1)[2, 1], 2), ', R2 ', round(mod1$r.squared, 2), ', pval ', round(coefficients(mod1)[2, 4], 3))
p1 <- ggplot(toplo, aes(x = tn, y = chla)) + 
  geom_point(aes(fill = mo), colour = 'grey', size = 3, pch = 21) + 
  geom_smooth(method = 'lm') + 
  labs(
    x = 'Area averaged TN (mg/L)', 
    y = 'Area averaged Chla (ug/L)', 
    subtitle = txt, 
    title = 'Chla ~ TN'
  ) +
  theme_minimal(base_family = bssz)

mod2 <- summary(lm(mcrbbest ~ tn, toplo))
txt <- paste0('Slope ', round(coefficients(mod2)[2, 1], 2), ', R2 ', round(mod2$r.squared, 2), ', pval ', round(coefficients(mod2)[2, 4], 3))
p2 <- ggplot(toplo, aes(x = tn, y = mcrbbest)) + 
  geom_point(aes(fill = mo), colour = 'grey', size = 3, pch = 21) + 
  geom_smooth(method = 'lm') + 
  labs(
    x = 'Area averaged TN (mg/L)', 
    y = 'Area averaged macroalgae abundance', 
    subtitle = txt, 
    title = 'Macroalgae ~ TN'
  ) +
  theme_minimal(base_family = bssz)

mod3 <- summary(lm(savbbest ~ mcrbbest, toplo))
txt <- paste0('Slope ', round(coefficients(mod3)[2, 1], 2), ', R2 ', round(mod3$r.squared, 2), ', pval ', round(coefficients(mod3)[2, 4], 3))
p3 <- ggplot(toplo, aes(x = mcrbbest, y = savbbest)) +
  geom_point(aes(fill = mo), colour = 'grey', size = 3, pch = 21) +
  geom_smooth(method = 'lm') +
  labs(
    x = 'Area averaged macroalgae abundance',
    y = 'Area averaged SAV abundance',
    subtitle = txt,
    title = 'SAV ~ Macroalgae'
  ) +
  theme_minimal(base_family = bssz)


# mod4 <- summary(lm(chla ~ nh34, toplo))
# txt <- paste0('Slope ', round(coefficients(mod4)[2, 1], 2), ', R2 ', round(mod4$r.squared, 2), ', pval ', round(coefficients(mod1)[2, 4], 3))
# p4 <- ggplot(toplo, aes(x = nh34, y = chla)) + 
#   geom_point(aes(fill = mo), colour = 'grey', size = 3, pch = 21) + 
#   geom_smooth(method = 'lm') + 
#   labs(
#     x = 'Area averaged NH3 (mg/L)', 
#     y = 'Area averaged Chla (ug/L)', 
#     subtitle = txt, 
#     title = 'Chla ~ NH3'
#   ) +
#   theme_minimal(base_family = bssz)
# 
# mod5 <- summary(lm(mcrbbest ~ nh34, toplo))
# txt <- paste0('Slope ', round(coefficients(mod5)[2, 1], 2), ', R2 ', round(mod5$r.squared, 2), ', pval ', round(coefficients(mod5)[2, 4], 3))
# p5 <- ggplot(toplo, aes(x = nh34, y = mcrbbest)) + 
#   geom_point(aes(fill = mo), colour = 'grey', size = 3, pch = 21) + 
#   geom_smooth(method = 'lm') + 
#   labs(
#     x = 'Area averaged NH3 (mg/L)', 
#     y = 'Area averaged macroalgae abundance', 
#     subtitle = txt, 
#     title = 'Macroalgae ~ NH3'
#   ) +
#   theme_minimal(base_family = bssz)
# 
# mod6 <- summary(lm(savbbest ~ mcrbbest, toplo))
# txt <- paste0('Slope ', round(coefficients(mod6)[2, 1], 2), ', R2 ', round(mod6$r.squared, 2), ', pval ', round(coefficients(mod6)[2, 4], 3))
# p6 <- ggplot(toplo, aes(x = mcrbbest, y = savbbest)) + 
#   geom_point(aes(fill = mo), colour = 'grey', size = 3, pch = 21) + 
#   geom_smooth(method = 'lm') + 
#   labs(
#     x = 'Area averaged macroalgae abundance', 
#     y = 'Area averaged SAV abundance', 
#     subtitle = txt, 
#     title = 'SAV ~ Macroalgae'
#   ) +
#   theme_minimal(base_family = bssz)  

p1 + p2 + p3 + plot_layout(ncol = 3, guides = 'collect') & 
  theme(
    legend.position = 'top', 
    legend.title = element_blank()
  )
```

### Water quality vs macroalge vs SAV, correlations

This correlation matrix shows correlations between seagrasses, macroalgae, and water quality.  Seagrass and macroalgae are correlated at the transect scale, i.e., observations are compared at the same location over time.  The water quality comparisons to the other variables are at an aggregated spatial scale because of different sample locations for the transects and water quality.  

```{r, fig.height = 6, fig.width = 6}

mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")
savsel <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme')
wqsel <- c('chla', 'nh34', 'tn', 'tp')
wqlab <- c('Chl-a', 'NH3, NH4+', 'TN', 'TP')

# correlate transect data at transect scale -------------------------------

# add area
trnsum <- rstrndat %>%
  mutate(
    date = floor_date(date, unit = 'week')
  ) %>%
  dplyr::group_by(date, station, location, taxa) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%
  group_by(date, station, taxa) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  )

trncors <- trnsum %>% 
  ungroup %>% 
  filter(taxa %in% c(savsel, mcrsel)) %>% 
  spread(taxa, foest) %>% 
  # filter(area %in% 'Area 3') %>% 
  select(-station, -date)


trncrs <- crossing(var1 = names(trncors), var2 = names(trncors)) %>% 
  filter(var1 != var2) %>% 
  rownames_to_column() %>% 
  group_by(rowname) %>% 
  nest %>% 
  mutate(
    crs = map(data, function(x){
      
      # variables
      vr1 <- trncors[[x$var1]]
      vr2 <- trncors[[x$var2]]
      
      # pearson
      pr_ts <- cor.test(vr1, vr2, method = 'pearson')
      pr_cr <- round(pr_ts$estimate, 2)
      pr_pv <- p_ast(pr_ts$p.value)
      pr <- paste(pr_cr, pr_pv)
      
      out <- data.frame(pr = pr, stringsAsFactors = F)
      return(out)
      
    })
  ) %>% 
  unnest(c('data', 'crs')) %>% 
  ungroup %>% 
  select(-rowname)


# correlate water quality with transect data at aggregated scale ----------

# add area
trnsum <- rstrndat %>% 
  mutate(
    date = floor_date(date, unit = 'week')
  ) %>% 
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  dplyr::group_by(area, date, station, location, taxa) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%
  group_by(area, date, taxa) %>% 
  summarize(
    val = sum(pa) / length(pa) # freq occ.
  ) %>% 
  ungroup %>% 
  filter(taxa %in% c(savsel, mcrsel)) %>% 
  rename(var = taxa) %>%
  spread(var, val)

# water quality summary
rswqsum <- rswqdat %>% 
  filter(var %in%  c('tn', 'nh34', 'chla', 'tp')) %>% 
  inner_join(rsstatloc, ., by = c('source', 'station')) %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>% 
  filter(!station %in% nonbay) %>% 
  select(-qual, -bswqstation, -nrmrng, -source, -source_lng, -uni, -lbunis, -inrng, -lbs) %>% 
  mutate(
    date = floor_date(date, unit = 'week')
  ) %>%
  group_by(area, date, var) %>% 
  summarise(
    val = mean(log(val)), 
    val = exp(val)
  ) %>% 
  spread(var, val)

tocor <- full_join(trnsum, rswqsum, by = c('area', 'date'))

crs <- crossing(var1 = c(wqsel, mcrsel, savsel), var2 = c(wqsel, mcrsel, savsel)) %>% 
  filter(var1 != var2) %>% 
  rownames_to_column() %>% 
  group_by(rowname) %>% 
  nest %>% 
  mutate(
    crs = map(data, function(x){
      
      # variables
      vr1 <- tocor[[x$var1]]
      vr2 <- tocor[[x$var2]]
      
      # pearson
      pr_ts <- cor.test(vr1, vr2, method = 'pearson')
      pr_cr <- round(pr_ts$estimate, 2)
      pr_pv <- p_ast(pr_ts$p.value)
      pr <- paste(pr_cr, pr_pv)
      
      out <- data.frame(pr = pr, stringsAsFactors = F)
      return(out)
      
    })
  ) %>% 
  unnest(c('data', 'crs')) %>% 
  ungroup %>% 
  select(-rowname) %>% 
  filter(!(var1 %in% c(mcrsel, savsel) & var2 %in% c(mcrsel, savsel)))


# combine correlations and make plot --------------------------------------

prplo <- bind_rows(crs, trncrs) %>% 
  separate(pr, c('cor', 'sig'), sep = ' ') %>%  
  mutate(
    cor = as.numeric(cor), 
    var1 = factor(var1, levels = c(wqsel, mcrsel, savsel), labels =  c(wqlab, mcrsel, savsel)), 
    var2 = factor(var2, levels = c(wqsel, mcrsel, savsel), labels =  c(wqlab, mcrsel, savsel)), 
    sig = gsub('ns', '', sig)
  )

pbase <- theme(
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 8), 
  axis.text.y = element_text(size = 8),
  legend.position = c(0.5, 1.12),
  legend.direction = 'horizontal',
  plot.margin = unit(c(4,4,0,0), "lines"),
  strip.background = element_blank(), 
  strip.text.y = element_text(angle = 0, hjust = 0, vjust = 0.5), 
  panel.background = element_rect(fill = 'black')
) 

outlab <- data.frame(
  y = c(2.5, 6.5, 10), 
  lab = c('Water Quality', 'Macroalgae', 'Seagrasses')
)

p <- ggplot(prplo) + 
  geom_tile(aes(y = var1, x = var2, fill = cor), colour = 'black') + 
  geom_text(aes(y = var1, x = var2, label = sig)) +
  annotation_custom(grob = textGrob(label = outlab$lab[1], hjust = 0, gp = gpar(cex = 0.7)),
                    ymin = outlab$y[1], ymax = outlab$y[1], xmin = 12, xmax = 12) +
  annotation_custom(grob = textGrob(label = outlab$lab[2], hjust = 0, gp = gpar(cex = 0.7)),
                    ymin = outlab$y[2], ymax = outlab$y[2], xmin = 12, xmax = 12) +  
  annotation_custom(grob = textGrob(label = outlab$lab[3], hjust = 0, gp = gpar(cex = 0.7)),
                    ymin = outlab$y[3], ymax = outlab$y[3], xmin = 12, xmax = 12) +
  annotation_custom(grob = textGrob(label = outlab$lab[1], hjust = 0.5, gp = gpar(cex = 0.7)),
                    xmin = outlab$y[1], xmax = outlab$y[1], ymin = 12, ymax = 12) +
  annotation_custom(grob = textGrob(label = outlab$lab[2], hjust = 0.5, gp = gpar(cex = 0.7)),
                    xmin = outlab$y[2], xmax = outlab$y[2], ymin = 12, ymax = 12) +  
  annotation_custom(grob = textGrob(label = outlab$lab[3], hjust = 0.5, gp = gpar(cex = 0.7)),
                    xmin = outlab$y[3], xmax = outlab$y[3], ymin = 12, ymax = 12) +
  pbase +
  scale_y_discrete('', expand = c(0, 0)) + #, labels = parse(text = rev(labs))) + 
  scale_x_discrete('', expand = c(0, 0)) + #, labels = parse(text = rev(labs))) +
  scale_fill_gradientn('Correlation', colours = c(muted("blue"), "white", muted("red")), limits = c(-1, 1)) +
  guides(fill = guide_colourbar(barheight = 0.5, barwidth = 5, label.theme = element_text(size = 6, angle = 0))) +
  geom_hline(yintercept = 4.5, size = 1.5) +
  geom_hline(yintercept = 8.5, size = 1.5) +
  geom_vline(xintercept = 4.5, size = 1.5) +
  geom_vline(xintercept = 8.5, size = 1.5)

# Code to override clipping
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$clip[gt$layout$name == "panel"] <- "off"
grid.draw(gt)
```

### Biomass estimates {#biomass}

Simple regressions of weight vs abundance by major macroalgae group.  Regressions are fit through the origin so that zero weight is fixed at zero mass.

```{r, fig.height = 5, fig.width = 4}
rssub <- rstrndat %>% 
  filter(typ == 'mcr') %>% 
  select(station, location, date, taxa, bb)

bssub <- bstransect %>% 
  filter(var == 'Abundance') %>% 
  select(
    station = Transect, 
    location = Site, 
    date = Date,
    taxa = Savspecies,
    bb = aveval
  ) %>% 
  filter(date %in% rstrnwts$date & station %in% rstrnwts$station & location %in% rstrnwts$location) %>% 
  filter(grepl('^DA', taxa)) %>% 
  mutate(
    location = as.numeric(location), 
    taxa = case_when(
      grepl('Red', taxa) ~ 'Red', 
      grepl('Green', taxa) ~ 'Green', 
      grepl('Macroalgae', taxa) ~ 'Red',
      T ~ taxa# verified that the location/date/transect match with weights was red, only one instance
    )
  ) 
  
tojn <- bind_rows(rssub, bssub)

wtssub <- rstrnwts %>% 
  select(-genus) %>% 
  rename(taxa = group) %>% 
  filter(!grepl('and', taxa)) %>% 
  left_join(tojn, by = c('station', 'date', 'location', 'taxa')) %>% 
  mutate(
    weight = weight_g * 0.004 # g / 0.25m2 to kg / m2
  )

ymax <- wtssub %>% 
  pull(weight) %>% 
  max(na.rm = T)

ggplot(wtssub, aes(x = bb, y = weight)) + 
  geom_point() +
  facet_wrap(~taxa, ncol = 1) +
  stat_smooth_func_with_pval(geom = "text", method = "lm", hjust = 0, parse = T, xpos = 0, xpos2 = 0, ypos = ymax, ypos2 = 0.8 * ymax, formula = y ~ 0 + x) +
  geom_smooth(method = 'lm', se = F, formula = y ~ 0 + x) + 
  labs(
    x = 'Abundance (bb)', 
    y = 'Weight (kg / m2)'
  ) + 
  theme_minimal() + 
  theme(
    strip.background = element_blank(), 
    panel.grid.minor = element_blank()
  )
```

Using the regression, a weight was assigned to all macroalgae samples where abundance (bb) was recorded.  The weights were then summed across all quadrats by transect, macroalgae group, month, and region ([here](https://tbep-tech.github.io/piney-point-analysis/trends#region-delineation)) to represent total sampled biomass. However, this is biased by relative sample effort (Area 1 has been sampled ~3x more than Area 3).  Below, the weights are standardized to unit effort.

```{r, fig.height = 3, fig.width = 5}
wtsmod <- wtssub %>% 
  group_by(taxa) %>% 
  nest() %>% 
  mutate(
    mod = purrr::map(data, lm, formula = weight ~ 0 + bb)
  ) %>% 
  select(taxa, mod)

wtsest <- rstrndat %>% 
  filter(taxa %in% c('Red', 'Green', 'Cyanobacteria')) %>% 
  group_by(taxa) %>% 
  nest %>% 
  left_join(wtsmod, by = 'taxa') %>% 
  mutate(
    weight_kgm2 = purrr::pmap(list(object = mod, newdata = data), predict)
  ) %>% 
  select(-mod) %>% 
  unnest(c('data', 'weight_kgm2'))
    
mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")
mcrcol <- c('tomato1', 'lightgreen', 'burlywood3', 'lightblue')
names(mcrcol) <- mcrsel
mcrcol <- mcrcol[unique(wtsest$taxa)]

# add area
wtsarea <- wtsest %>%
  inner_join(rstrnpts, ., by = 'station') %>% 
  st_intersection(ppseg) %>% 
  st_set_geometry(NULL) %>%
  dplyr::group_by(area, typ, date, station, taxa, location) %>%
  mutate(
    dateflr = floor_date(date, unit = 'month')
  ) 

wtssum <- wtsarea %>% 
  group_by(area, dateflr, taxa) %>% 
  summarize(
    kgest = sum(weight_kgm2) / 4 # divide by four to get total biomass across all quads (0.25 m2)
  )

ggplot(wtssum, aes(x = dateflr, y = kgest, fill = taxa)) + 
  geom_bar(pch = 16, stat = 'identity', color = 'grey') +
  facet_grid(~ area, scales = 'free') +
  theme_minimal(base_size = 14) + 
  scale_fill_manual(values = mcrcol) +
  labs(
    y = 'Total sampled weight (kg)'
  ) +
  coord_cartesian(ylim = c(0, NA)) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_line()
  )
```

Because of different sampling effort by month and location, weights were standardized by dividing by sample effort (number of unique transect visits in a month).  This can also be considered as kg / 1.25 m2 because five quadrats of 0.25 m2 are surveyed each transect.  So converting this to kg / m2, the weights are multiplied by 0.8 (1 / 1.25).

```{r, fig.height = 3, fig.width = 5}
wtssum <- wtsarea %>% 
  unite('datetrn', date, station, remove = F) %>% 
  group_by(area, dateflr, taxa) %>% 
  summarize(
    kgest = sum(weight_kgm2) / 4 / length(unique(datetrn)), 
    kgest = kgest * 0.8
  )

ggplot(wtssum, aes(x = dateflr, y = kgest, fill = taxa)) + 
  geom_bar(pch = 16, stat = 'identity', color = 'grey') +
  facet_grid(~ area, scales = 'free') +
  theme_minimal(base_size = 14) + 
  scale_fill_manual(values = mcrcol) +
  labs(
    y = 'Biomass (kg / m2)'
  ) +
  coord_cartesian(ylim = c(0, NA)) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_line()
  )

```

### August 2021 v 2022 {#Aug}

```{r, fig.height = 3, fig.width = 6}
mcrsel <- c("Red", "Green", "Brown", "Cyanobacteria")
savsel <- c('Thalassia testudinum', 'Halodule wrightii', 'Syringodium filiforme')

rstrndat2 <- rstrndat %>% 
  filter(station %in% c('S3T5b', 'S3T5a'))

# add area
trnsum <- rstrndat2 %>%
  dplyr::group_by(typ, date, station, taxa, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(typ, date, taxa) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  filter(taxa %in% c(mcrsel, savsel))

trnsumtots <- rstrndat2 %>%
  dplyr::group_by(typ, date, station, location) %>%
  dplyr::summarise(
    pa = as.numeric(any(bb > 0))
  ) %>%  
  mutate(
    date = floor_date(date, unit = 'month'), 
    typ = factor(typ, levels = c('mcr', 'sav'), labels = c('Macroalgae', 'Seagrass'))
  ) %>% 
  group_by(typ, date) %>% 
  summarize(
    foest = sum(pa) / length(pa)
  ) %>% 
  mutate(taxa = 'Total')

toplo <- bind_rows(trnsum, trnsumtots) %>% 
  mutate(taxa = factor(taxa)) %>% 
  filter(month(date) == 8) %>% 
  filter((taxa %in% 'Cyanobacteria') | (typ == 'Seagrass' & taxa == 'Total')) %>%
  # filter(taxa == 'Total') %>% 
  mutate(
    taxa = as.character(taxa),
    taxa = case_when(
      taxa == 'Total' ~ 'Seagrass', 
      T ~ taxa
    )
  )

ggplot(toplo, aes(x = date, y = foest)) + 
  with_shadow(geom_bar(aes(fill = taxa, group = taxa), stat = 'identity', color = 'black', position = position_dodge()), sigma = 2.7, x_offset = 0, y_offset = 0) +
  # facet_grid(typ ~ .) +
  theme_minimal(base_size = 14) + 
  scale_x_date(breaks = mdy(c('08-01-2021', '08-01-2022')), date_labels = '%b %Y') +
  scale_y_continuous(limits = c(0, 1))+
  scale_fill_manual(values = c(Cyanobacteria = 'lightblue', Seagrass = '#00806E')) +
  labs(
    y = 'Freq. occurrence',
    title = 'Transect S3T5 near Cockroach Bay'
  ) +
  theme(
    legend.position = 'top', 
    legend.title = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_text(size = 14), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_line(), 
    panel.grid.minor = element_blank(),
    # panel.spacing=unit(2, "lines"), 
    panel.background = element_blank(), 
    panel.grid.major = element_blank()
  )
```