minor fig reworking

R/figs.R

@@ -31,6 +31,50 @@ library(gridExtra)
 seglng <- c('Old Tampa Bay', 'Hillsborough Bay', 'Middle Tampa Bay', 'Lower Tampa Bay')
 segshr <- c('OTB', 'HB', 'MTB', 'LTB')
 
+# data used by more than one figure
+
+spp <- read_csv(here("data/tbm_combined_catch_env_factors.csv"))
+div <- read_csv(here("data/phy_tbni_sgrs.csv"))
+
+FLUC <- spp %>%
+  #filter(FLUCCSCODE %in% c(NA,9113,9116,9121)) 
+  mutate(
+    SAVcover = case_when(
+      FLUCCSCODE == 9113 ~ "patchy",
+      FLUCCSCODE == 9116 ~ "continuous",
+      FLUCCSCODE == 9121 ~ "algae",
+      TRUE ~ "none"  # Default case if none of the above
+      ),
+    Dominant = case_when(
+      DominantVeg == "Thalassia spp." ~ "Thalassia",
+      DominantVeg == "Halodule spp." ~ "Halodule",
+      DominantVeg == "Syringodium spp." ~ "Syringodium",
+      DominantVeg == "Ruppia spp." ~ "Ruppia",
+      #  DominantVeg == "Seagrasses: Mixed" ~ "Mixed",
+      DominantVeg == "None" ~ "None",
+      DominantVeg %in% c("Algae","Algae: Filamentous green","Algae: Filamentous red","Caulerpa spp.",
+                         "Gracillaria","Sargasum spp.") ~ "Algae",
+      TRUE ~ "Mixed/OtherSAV"  # Default case if none of the above 
+    ),
+    TBEP_seg = factor(TBEP_seg, levels = segshr),
+    SAVcover = factor(SAVcover, levels = c('continuous', 'patchy', 'algae', 'none'))
+  )
+
+BVcover <- FLUC %>% 
+  summarize(
+    mean_value = mean(BottomVegCover, na.rm = TRUE),
+    std_error = sd(BottomVegCover, na.rm = TRUE) / sqrt(n()),
+    Count = n(),
+    .by = c(TBEP_seg, SAVcover)
+  )
+
+Dominant <- FLUC %>% 
+  summarize(
+    mean_value = mean(BottomVegCover, na.rm = TRUE),
+    std_error = sd(BottomVegCover, na.rm = TRUE) / sqrt(n()),
+    Count = n(), 
+    .by = c(Dominant, TBEP_seg)
+  )
 
 # map -----------------------------------------------------------------------------------------
 
@@ -46,7 +90,12 @@ fl3 <- paste0(tempdir(), '/sgdat2022.RData')
 download.file('https://github.com/tbep-tech/hmpu-workflow/raw/master/data/sgdat2022.RData', destfile = fl3)
 load(file = fl3)
 
-data(fimsta)
+data(fimstations)
+fimsta <- fimstations %>% 
+  mutate(
+    yr = substr(Reference, 4, 7)
+  ) %>% 
+  filter(yr %in% c(1998:2021))
 
 sgdat99 <- sgdat1999 %>% 
   filter(FLUCCSCODE %in% c(9113, 9116)) %>% 
@@ -116,20 +165,21 @@ m2 <- ggplot() +
   ggspatial::annotation_map_tile(zoom = 11, type = 'cartolight', cachedir = system.file("rosm.cache", package = "ggspatial")) +
   geom_sf(data = sgdat16, fill = 'darkgreen', color = NA, inherit.aes = F) +
   #geom_sf(data = trnpts, color = 'black', inherit.aes = F) +
-  annotation_north_arrow(location = 'tl', style = north_arrow_orienteering(fill = c('black', 'black'), text_col = NA), 
-                         height = unit(0.5, "cm"), width = unit(0.5, "cm")) +
-  annotation_scale(location = 'br', text_cex = 1) +
+  # annotation_north_arrow(location = 'tl', style = north_arrow_orienteering(fill = c('black', 'black'), text_col = NA),
+  #                        height = unit(0.5, "cm"), width = unit(0.5, "cm")) +
+  # annotation_scale(location = 'br', text_cex = 1) +
   geom_sf(data = tbseg, fill = NA, color = NA, inherit.aes = F) +
   geom_sf(data = tbseglines, color = 'black', inherit.aes = F) +
-  geom_sf_text(data = segcent, aes(label = bay_segment), size = 4, color = 'black', inherit.aes = F) +
+  # geom_sf_text(data = segcent, aes(label = bay_segment), size = 4, color = 'black', inherit.aes = F) +
   # annotation_custom(ggplotGrob(minset), xmin = -9.185e6, xmax = -9.17e6, ymin = 3.22e6, ymax = 3.28e6) + 
-  annotation_custom(ggplotGrob(minset), xmin = bbox[3] - 0.1, xmax = bbox[3] + 0.015, ymin = bbox[4] - 0.1, ymax = bbox[4] + 0.06) + 
+  # annotation_custom(ggplotGrob(minset), xmin = bbox[3] - 0.1, xmax = bbox[3] + 0.015, ymin = bbox[4] - 0.1, ymax = bbox[4] + 0.06) + 
   coord_sf(xlim = bbox[c('xmin', 'xmax')], ylim = bbox[c('ymin', 'ymax')], crs = 4326) +
   labs(
     subtitle = '(b) Bay segments, seagrass 2016'
   ) +
   thm + 
-  theme(axis.text.x = element_blank()
+  theme(axis.text.x = element_blank(), 
+        axis.text.y = element_blank()
         )
 
 # xnrg <- ggplot_build(m2)$layout$panel_scales_x[[1]]$range$range
@@ -139,14 +189,14 @@ m3 <- ggplot() +
   ggspatial::annotation_map_tile(zoom = 11, type = 'cartolight', cachedir = system.file("rosm.cache", package = "ggspatial")) +
   geom_sf(data = sgdat22, fill = 'darkgreen', color = NA, inherit.aes = F) +
   #geom_sf(data = trnpts, color = 'black', inherit.aes = F) +
-  annotation_north_arrow(location = 'tl', style = north_arrow_orienteering(fill = c('black', 'black'), text_col = NA), 
-                         height = unit(0.5, "cm"), width = unit(0.5, "cm")) +
-  annotation_scale(location = 'br', text_cex = 1) +
+  # annotation_north_arrow(location = 'tl', style = north_arrow_orienteering(fill = c('black', 'black'), text_col = NA),
+  #                        height = unit(0.5, "cm"), width = unit(0.5, "cm")) +
+  # annotation_scale(location = 'br', text_cex = 1) +
   geom_sf(data = tbseg, fill = NA, color = NA, inherit.aes = F) +
   geom_sf(data = tbseglines, color = 'black', inherit.aes = F) +
-  geom_sf_text(data = segcent, aes(label = bay_segment), size = 4, color = 'black', inherit.aes = F) +
+  # geom_sf_text(data = segcent, aes(label = bay_segment), size = 4, color = 'black', inherit.aes = F) +
   # annotation_custom(ggplotGrob(minset), xmin = -9.185e6, xmax = -9.17e6, ymin = 3.22e6, ymax = 3.28e6) + 
-  annotation_custom(ggplotGrob(minset), xmin = bbox[3] - 0.1, xmax = bbox[3] + 0.015, ymin = bbox[4] - 0.1, ymax = bbox[4] + 0.06) + 
+  # annotation_custom(ggplotGrob(minset), xmin = bbox[3] - 0.1, xmax = bbox[3] + 0.015, ymin = bbox[4] - 0.1, ymax = bbox[4] + 0.06) + 
   coord_sf(xlim = bbox[c('xmin', 'xmax')], ylim = bbox[c('ymin', 'ymax')], crs = 4326) +
   labs(
     subtitle = '(c) Bay segments, seagrass 2022'
@@ -159,17 +209,18 @@ m3 <- ggplot() +
 m4 <- ggplot() + 
   ggspatial::annotation_map_tile(zoom = 11, type = 'cartolight', cachedir = system.file("rosm.cache", package = "ggspatial")) +
   # geom_sf(data = sgdat22, fill = 'darkgreen', color = NA, inherit.aes = F) +
-  annotation_north_arrow(location = 'tl', style = north_arrow_orienteering(fill = c('black', 'black'), text_col = NA), 
-                         height = unit(0.5, "cm"), width = unit(0.5, "cm")) +
-  annotation_scale(location = 'br', text_cex = 1) +
-  geom_sf(data = tbseg, fill = NA, color = NA, inherit.aes = F) +
+  # annotation_north_arrow(location = 'tl', style = north_arrow_orienteering(fill = c('black', 'black'), text_col = NA),
+  #                        height = unit(0.5, "cm"), width = unit(0.5, "cm")) +
+  # annotation_scale(location = 'br', text_cex = 1) +
+  geom_sf(data = tbseglines, color = 'black', inherit.aes = F) +
   geom_sf(data = fimsta, color = 'black', inherit.aes = F, size = 0.5, alpha = 0.5) + 
-  annotation_custom(ggplotGrob(minset), xmin = bbox[3] - 0.1, xmax = bbox[3] + 0.015, ymin = bbox[4] - 0.1, ymax = bbox[4] + 0.06) + 
+  # annotation_custom(ggplotGrob(minset), xmin = bbox[3] - 0.1, xmax = bbox[3] + 0.015, ymin = bbox[4] - 0.1, ymax = bbox[4] + 0.06) + 
   coord_sf(xlim = bbox[c('xmin', 'xmax')], ylim = bbox[c('ymin', 'ymax')], crs = 4326) +
   labs(
     subtitle = '(d) 21.3-m seines, 1998-2021'
   ) +
-  thm 
+  thm + 
+  theme(axis.text.y  = element_blank())
 
 m <- m1 + m2 + m3 + m4 + plot_layout(ncol = 2)
 
@@ -178,78 +229,57 @@ print(m)
 dev.off()
 
 #SAV summaries ---------------------------------------------------------------
-load(file = url('https://github.com/tbep-tech/tbep-os-presentations/raw/master/data/sgsegest.RData'))
-
-#import processed community data
-spp <- read_csv(here("data/tbm_combined_catch_env_factors.csv"))
-div <- read_csv(here("data/phy_tbni_sgrs.csv"))
-
-FLUC <- spp %>%
-  #filter(FLUCCSCODE %in% c(NA,9113,9116,9121)) 
-  mutate(SAVcover = case_when(
-    FLUCCSCODE == 9113 ~ "patchy",
-    FLUCCSCODE == 9116 ~ "continuous",
-    FLUCCSCODE == 9121 ~ "algae",
-    TRUE ~ "none"  # Default case if none of the above
-  )) 
-FLUC <- FLUC %>%
-  mutate(Dominant = case_when(
-    DominantVeg == "Thalassia spp." ~ "Thalassia",
-    DominantVeg == "Halodule spp." ~ "Halodule",
-    DominantVeg == "Syringodium spp." ~ "Syringodium",
-    DominantVeg == "Ruppia spp." ~ "Ruppia",
-    #  DominantVeg == "Seagrasses: Mixed" ~ "Mixed",
-    DominantVeg == "None" ~ "None",
-    DominantVeg %in% c("Algae","Algae: Filamentous green","Algae: Filamentous red","Caulerpa spp.",
-                       "Gracillaria","Sargasum spp.") ~ "Algae",
-    TRUE ~ "Mixed/OtherSAV"  # Default case if none of the above 
-  ))
-
-BVcover <- FLUC %>% 
-  group_by(TBEP_seg,SAVcover) %>%
-  summarize(mean_value = mean(BottomVegCover, na.rm = TRUE),
-            std_error = sd(BottomVegCover, na.rm = TRUE) / sqrt(n()),
-            Count=n())
-  BVcover <- BVcover %>%
-    mutate(segment = factor(TBEP_seg, levels = c('OTB', 'HB', 'MTB', 'LTB'))) 
 
-plot1 <- ggplot(BVcover, aes(x = reorder(SAVcover,Count), y = mean_value)) +
+toplo1 <- BVcover %>% 
+  mutate(
+    txtloc = ifelse(is.na(std_error), mean_value, mean_value + std_error)
+  )
+
+plot1 <- ggplot(toplo1, aes(x = reorder(SAVcover,Count), y = mean_value)) +
   geom_col(fill = "lightgreen") +  # Bar plot for means
   geom_errorbar(
     aes(ymin = mean_value - std_error, ymax = mean_value + std_error),
     width = 0.2
   ) +  # Error bars
-  facet_wrap(~segment,ncol=4)+
-  geom_text(aes(label = Count), colour ="black", size=2, nudge_y = 16) +
+  facet_wrap(~TBEP_seg,ncol=4)+
+  geom_text(aes(label = Count, y = txtloc), colour ="black", size=2, nudge_y = 16) +
   labs(
     x = "FLUCCSCODE",
     y = NULL
-  )  + scale_y_continuous(limits = c(0,120),breaks=breaks_extended(4)) + coord_flip()+
-theme_minimal() + theme(axis.title.x = element_blank(),
-                        axis.text.x = element_blank())
+  )  + 
+  scale_y_continuous(limits = c(0,120), breaks = seq(0, 100, by = 25)) + 
+  coord_flip() +
+  theme_minimal() +
+  theme(
+    axis.title.x = element_blank(),
+    axis.text.x = element_blank(),, 
+    panel.grid.minor = element_blank()
+  )
 
-Dominant <- FLUC %>% 
-  group_by(Dominant,TBEP_seg) %>%
-  summarize(mean_value = mean(BottomVegCover, na.rm = TRUE),
-            std_error = sd(BottomVegCover, na.rm = TRUE) / sqrt(n()),
-            Count=n())
-Dominant <- Dominant %>%
-  mutate(segment = factor(TBEP_seg, levels = c('OTB', 'HB', 'MTB', 'LTB'))
+toplo2 <- Dominant %>% 
+  mutate(
+    txtloc = ifelse(is.na(std_error), mean_value, mean_value + std_error)
   )
 
-plot2 <- ggplot(Dominant, aes(x = reorder(Dominant,Count), y = mean_value)) +
+plot2 <- ggplot(toplo2, aes(x = reorder(Dominant,Count), y = mean_value)) +
   geom_col(fill = "lightgreen") +  # Bar plot for means
   geom_errorbar(
     aes(ymin = mean_value - std_error, ymax = mean_value + std_error),
     width = 0.2
   ) +  # Error bars 
-  facet_wrap(~segment,ncol=4)+
-  geom_text(aes(label = Count), colour ="black", size=2, nudge_y = 16) + 
+  facet_wrap(~TBEP_seg,ncol=4)+
+  geom_text(aes(label = Count, y = txtloc), colour ="black", size=2, nudge_y = 16) + 
   labs(
     x = "Dominant SAV",
     y = "Percent SAV cover"
-  ) + scale_y_continuous(limits = c(0,120),breaks=breaks_extended(4)) + coord_flip()+
-  theme_minimal()+theme(strip.text.x=element_blank())
+  ) + 
+  scale_y_continuous(limits = c(0,120), breaks = seq(0, 100, by = 25)) + 
+  coord_flip() + 
+  theme_minimal() +
+  theme(
+    strip.text.x = element_blank(), 
+    panel.grid.minor = element_blank()
+  )
 
 sgsum <- plot1/plot2
 
@@ -259,14 +289,14 @@ dev.off()
 
 #SAV and TBNI by FLUCCSCODE------------------------------------------------------------------
 
-st1 <- ggplot(BVcover, aes(x = reorder(SAVcover,Count), y = mean_value)) +
+st1 <- ggplot(BVcover, aes(x = SAVcover, y = mean_value)) +
   geom_col(fill = "lightgreen") +  # Bar plot for means
   geom_errorbar(
     aes(ymin = mean_value - std_error, ymax = mean_value + std_error),
     width = 0.2
   ) +  # Error bars
-  facet_wrap(~segment,ncol=4)+
-  geom_text(aes(label = Count), colour ="black", size=3, nudge_y = 16) +
+  facet_wrap(~TBEP_seg,ncol=4)+
+  geom_text(aes(label = Count, y = mean_value + std_error), colour ="black", size=3, nudge_y = 10) +
   labs(
     x = "FLUCCSCODE",
     y = "Percent SAV cover",
@@ -276,7 +306,9 @@ st1 <- ggplot(BVcover, aes(x = reorder(SAVcover,Count), y = mean_value)) +
   theme_minimal() + 
   theme(
     axis.title.x = element_blank(),
-    axis.text.x = element_text(colour = 'black', angle = 60, size = 9, hjust = 1)
+    axis.text.x = element_text(colour = 'black', angle = 45, size = 9, hjust = 1),
+    panel.grid.minor = element_blank(), 
+    panel.grid.major.x = element_blank()
   )
 
 Dominantper <- FLUC %>% 
@@ -285,14 +317,13 @@ Dominantper <- FLUC %>%
     .by = c(Dominant, TBEP_seg, SAVcover)
   ) %>%
   mutate(
-    perc = (Count / 1475 * 100),
-    segment = factor(TBEP_seg, levels = c('OTB', 'HB', 'MTB', 'LTB'))
+    perc = (Count / 1475 * 100)
   ) %>% 
-  filter(!Dominant %in% c('Algae', 'None'))
+  filter(!Dominant %in% c('None'))
 
-st2 <- ggplot(Dominantper, aes(x = reorder(SAVcover,Count), y = perc, fill = Dominant)) +
+st2 <- ggplot(Dominantper, aes(x = SAVcover, y = perc, fill = Dominant)) +
   geom_bar (stat="identity") +  # Bar plot for percent
-  facet_wrap(~segment,ncol=4)+
+  facet_wrap(~TBEP_seg,ncol=4)+
  # geom_text(aes(label = Count), colour ="black", size=2, nudge_y = 16) + 
   labs(
     x = "FLUCCSCODE",
@@ -303,40 +334,42 @@ st2 <- ggplot(Dominantper, aes(x = reorder(SAVcover,Count), y = perc, fill = Dom
   theme_minimal() + 
   theme(
     axis.title.x = element_blank(),
-    axis.text.x = element_text(colour = 'black', angle = 60, size = 9, hjust = 1)
+    axis.text.x = element_text(colour = 'black', angle = 45, size = 9, hjust = 1), 
+    panel.grid.minor = element_blank(), 
+    panel.grid.major.x = element_blank()
   )
 
-TBNI <- div %>% 
-  mutate(SAVcover = case_when(
-    FLUCCSCODE == 9113 ~ "patchy",
-    FLUCCSCODE == 9116 ~ "continuous",
-    FLUCCSCODE == 9121 ~ "algae",
-    TRUE ~ "none"  # Default case if none of the above
-  ))
-
-TBNI_sav <- TBNI %>% 
-  group_by(SAVcover, TBEP_seg) %>%
-  summarize(mean_value = mean(TBNI_Score, na.rm = TRUE),
-            std_error = sd(TBNI_Score, na.rm = TRUE) / sqrt(n()),
-            Count=n())
-
-TBNI_sav <- TBNI_sav%>%
-  mutate(segment = factor(TBEP_seg, levels = c('OTB', 'HB', 'MTB', 'LTB'))
+TBNI_sav <- div %>% 
+  mutate(
+    SAVcover = case_when(
+      FLUCCSCODE == 9113 ~ "patchy",
+      FLUCCSCODE == 9116 ~ "continuous",
+      FLUCCSCODE == 9121 ~ "algae",
+      TRUE ~ "none"  # Default case if none of the above
+    ),
+    TBEP_seg = factor(TBEP_seg, levels = segshr),
+    SAVcover = factor(SAVcover, levels = c('continuous', 'patchy', 'algae', 'none'))
+  ) %>% 
+  summarize(
+    mean_value = mean(TBNI_Score, na.rm = TRUE),
+    std_error = sd(TBNI_Score, na.rm = TRUE) / sqrt(n()),
+    Count = n(), 
+    .by = c(SAVcover, TBEP_seg)
   )
 
-st3 <- ggplot(TBNI_sav, aes(x = reorder(SAVcover,Count), y = mean_value)) +
+st3 <- ggplot(TBNI_sav, aes(x = SAVcover, y = mean_value)) +
   geom_point() +
   geom_errorbar(
     aes(ymin = mean_value - std_error, ymax = mean_value + std_error),
     width = 0.2
   ) +  # Error bars 
-  facet_wrap(~segment, ncol = 4) +
-  theme_bw() + 
+  facet_wrap(~TBEP_seg, ncol = 4) +
+  theme_minimal() + 
   theme(panel.grid.minor =element_blank(),
         panel.grid.major.x =element_blank(),
         # plot.background = element_rect(fill = NA, color = NA),
         axis.text.y = element_text(colour = 'black', size = 9),
-        axis.text.x = element_text(colour = 'black', angle = 60, size = 9, hjust = 1),
+        axis.text.x = element_text(colour = 'black', angle = 45, size = 9, hjust = 1),
         strip.background = element_blank(),
         strip.text = element_text(size = 11),
         legend.position = 'none'
@@ -403,17 +436,17 @@ p1 <- ggplot(toplo1, aes(x = factor(year), y = acres)) +
   )
 
 TBNI <- div %>% 
-  group_by(sgyear, TBEP_seg) %>%
-  summarize(mean_value = mean(TBNI_Score, na.rm = TRUE),
-            std_error = sd(TBNI_Score, na.rm = TRUE) / sqrt(n()),
-            Count=n())
-
-TBNI <- TBNI%>%
-  mutate(segment = factor(TBEP_seg, levels = c('OTB', 'HB', 'MTB', 'LTB'))
-)
+  summarize(
+    mean_value = mean(TBNI_Score, na.rm = TRUE),
+    std_error = sd(TBNI_Score, na.rm = TRUE) / sqrt(n()),
+    Count = n(), 
+    .by = c(sgyear, TBEP_seg)
+  ) %>% 
+  mutate(
+    segment = factor(TBEP_seg, levels = segshr)
+  )
 
 p2 <- ggplot(TBNI, aes(x = factor(sgyear), y = mean_value)) +
-  geom_line() + 
   geom_point() +
   geom_errorbar(
     aes(ymin = mean_value - std_error, ymax = mean_value + std_error),