This file evaluates seagrass gains and losses in coverage from 2018 to 2020. The goal is to characterize the distribution of depth estimates for areas where changes have occurred. A hypothesis is that seagrass may be lost at deeper depths as compared to where it was gained, possibly related to changes in water clarity that can cause seagrass to grow in more shallow areas where they are less light limited. Source content is here.

The analysis is in several steps:

  1. A change analysis between 2018 and 2020 to identify areas where seagrass occurred in 2018 but not in 2020 (loss) and areas where seagrass did not occur in 2018 but did occur in 2020 (gain).
  2. A bathymetry layer is then used to estimate the average depth, as meters below MLLW, for each polygon identified as lost or gained.
  3. The summarized polygons are then intersected with bay segment and seagrass management areas to evaluate potential changes by major areas in the bay.

First, the relevant libraries and datasets are imported. All analyses are conducted using the NAD83(HARN) / Florida West (ftUS) projection.

knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)

library(sf)
library(tidyverse)
library(tbeptools)
library(mapview)
library(leaflet)
library(units)
library(stars)
library(raster)
library(here)
library(spatstat)

# https://github.com/NicolasWoloszko/stat_ecdf_weighted
source('R/stat_ecdf_weighted.R')

data(dem)
data(sgseg)
data(epcdata)
data(sgdat2018)
data(sgdat2020)
data(tnanndat)

# colors
cols <- c('green4', 'tomato1')
names(cols) <- c('gained', 'lost')

# for maps
colgrn <- c("#F7FCF5", "#E5F5E0", "#C7E9C0", "#A1D99B", "#74C476", "#41AB5D", 
            "#238B45", "#006D2C", "#00441B")
colred <- c("#FFF5F0", "#FEE0D2", "#FCBBA1", "#FC9272", "#FB6A4A", "#EF3B2C", 
            "#CB181D", "#A50F15", "#67000D")               

# this is the projstring for NAD83(HARN) / Florida West (ftUS)
# https://spatialreference.org/ref/epsg/2882/
prj4 <- '+proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs'

# reproject layers
tbseg <- tbseg %>% 
  st_transform(crs = prj4)
sgmanagement <- sgmanagement %>% 
  st_transform(crs = prj4)

# reproject dem, convert to stars 
demstr <- dem %>% 
  projectRaster(crs = prj4) %>% 
  st_as_stars()

# reproject seagrass layers
sgdat2018 <- sgdat2018 %>% 
  st_transform(crs = prj4)
stdat2020 <- sgdat2020 %>% 
  st_transform(crs = prj4)

This code chunk does the following:

  1. Filters continuous and patchy seagrass from the 2018 and 2020 layers, unions by features
  2. Conducts a true change analysis by taking the differences between the 2018 and 2020 seagrass layers, then their intersection
  3. Combines all results showing changes between seagrass categories (source as 2018 category, target as 2020 category)
  4. Saves the results because the change analysis takes a while.
# 2018 seagrass data
a <- sgdat2018 %>% 
  left_join(fluccs, by = 'FLUCCSCODE') %>%
  st_union(by_feature = TRUE) %>%
  mutate(
    Category = paste0(Category, ', 2018')
  )

# 2020 seagrass data
b <- sgdat2020 %>% 
  st_union(by_feature = TRUE) %>%
  left_join(fluccs,by = 'FLUCCSCODE') %>% 
  mutate(
    Category = paste0(Category, ', 2020')
  )
  
# so intersect doesnt complain about attributes
st_agr(a) = "constant"
st_agr(b) = "constant"

# union separate layers for faster intersect
aunion <- a %>% 
  st_union %>% 
  st_set_precision(1e5) %>% 
  st_make_valid() %>% 
  st_buffer(dist = 0)
bunion <- b %>% 
  st_union %>% 
  st_set_precision(1e5) %>% 
  st_make_valid() %>% 
  st_buffer(dist = 0)
  
# get full union
op1 <- st_difference(a, bunion)
op2 <- st_difference(b, aunion) %>%
  rename(Category.1 = Category)
op3 <- st_intersection(a, b)

# combine and make multipolygon to polygon, do by id otherwise it doesn't work
unidat <- bind_rows(op1, op2, op3) %>%
  mutate(
    yr = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category))),
    yr.1 = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category.1))),
    Category = ifelse(is.na(Category), paste0('other, ', yr), as.character(Category)),
    Category.1 = ifelse(is.na(Category.1), paste0('other, ', yr.1), as.character(Category.1)),
    idval = 1:nrow(.)
  ) %>%
  dplyr::select(idval, Category.1, Category) %>%
  dplyr::select(idval, source = Category, target = Category.1) %>% 
  group_by(idval, source, target) %>% 
  nest() %>% 
  mutate(
    data = purrr::map(data, st_cast, 'POLYGON')
  ) %>% 
  unnest('data') %>% 
  st_as_sf()

The unidat dataset is subset to only gains and losses (i.e., polygons can remain the same between years or change between seagrass categories). Areas of each polygon are calculated in acres. Polygons less than 0.023 acres are also removed because these are smaller than the cell size (30 by 30 feet, or 0.023 acres) of the bathymetry layer. The file is saved because it takes a long time to calculate.

# select only lost/gained
chgdat <- unidat %>% 
  ungroup() %>% 
  filter(target == 'other, 2020' | source == 'other, 2018') %>% 
  mutate(
    Acres = st_area(.), 
    Acres = set_units(Acres, 'acres'), 
    Acres = as.numeric(Acres), 
    var = case_when(
      target == 'other, 2020' ~ 'lost', 
      source == 'other, 2018' ~ 'gained'
    )
  ) %>% 
  dplyr::filter(Acres > 0.023) # remove slivers and those less than the pixel size (dem pixel size is about 30x30ft)

save(chgdat, file = here('data/chgdat.RData'))

Reload the data for the rendered file. The results show polygon types in 2018 and 2020 in the source and target columns. A polygon type of ‘other’ in the source column means seagrass was gained in 2020 and a polygon type of ‘other’ in the target column means seagrass was lost in 2020.

load(file = here('data/chgdat.RData'))
chgdat
## Simple feature collection with 4310 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 381530.6 ymin: 1150157 xmax: 531285.6 ymax: 1345319
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## # A tibble: 4,310 x 6
##    idval source       target                               geometry  Acres var  
##  * <int> <chr>        <chr>              <POLYGON [US_survey_foot]>  <dbl> <chr>
##  1     7 Patchy, 2018 other, 2020 ((411859.5 1250585, 411851.4 125~  0.841 lost 
##  2     8 Patchy, 2018 other, 2020 ((410664.3 1264177, 410671.8 126~  0.237 lost 
##  3     8 Patchy, 2018 other, 2020 ((411085.6 1262250, 411087.4 126~  0.973 lost 
##  4     8 Patchy, 2018 other, 2020 ((411220.3 1261246, 411225.4 126~  0.239 lost 
##  5    10 Patchy, 2018 other, 2020 ((391000.9 1273733, 390998.9 127~  0.869 lost 
##  6    15 Patchy, 2018 other, 2020 ((416303.4 1211170, 416301.8 121~ 15.0   lost 
##  7    16 Patchy, 2018 other, 2020 ((416374.9 1210817, 416377.2 121~  0.248 lost 
##  8    17 Patchy, 2018 other, 2020 ((415834.4 1209104, 415845.7 120~  0.902 lost 
##  9    22 Patchy, 2018 other, 2020 ((385843.3 1274880, 385809.3 127~  0.218 lost 
## 10    25 Patchy, 2018 other, 2020 ((382545.8 1284993, 382537.6 128~  0.109 lost 
## # ... with 4,300 more rows
m <- mapview(dem, layer.name = 'Depth (m)') + mapview(chgdat, zcol = 'var', layer.name = 'Seagrass', col.regions = cols)
m


Next, average depth for each polygon is esimated from the bathymetry layer. The chgdat layer (gained/lost) is aggregated by taking the average of the cells from the bathymetry layer that are within each polygon. Polygons are removed that have no mean depth estimate because they are on the boundary or outside of the bathymetry layer. This also takes a while and is not rendered in this file.

chgdatarea <- chgdat %>% 
  aggregate(demstr, ., mean, as_points = T) %>% 
  st_as_sf() %>% 
  rename(meandepth = GDAL.Band.Number.1) %>% 
  bind_cols(st_set_geometry(chgdat, NULL)) %>% 
  filter(!is.na(meandepth))

save(chgdatarea, file = here('data/chgdatarea.RData'))

Reload the data. The meandepth column shows the average depth within each polygon as meters below MLLW (negative values).

load(file = here('data/chgdatarea.RData'))
chgdatarea
## Simple feature collection with 3236 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 409944.4 ymin: 1150701 xmax: 531285.6 ymax: 1345319
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## First 10 features:
##      meandepth idval       source      target      Acres  var
## 1  -0.69418468    27  Cont., 2018 other, 2020 0.10473300 lost
## 2   0.43311293    31  Cont., 2018 other, 2020 0.04347834 lost
## 3  -1.19202473    49 Patchy, 2018 other, 2020 0.30698955 lost
## 4   0.21248606    49 Patchy, 2018 other, 2020 0.37404466 lost
## 5   0.08625199    49 Patchy, 2018 other, 2020 0.08304300 lost
## 6  -1.04209483    52 Patchy, 2018 other, 2020 0.15276761 lost
## 7  -1.14997934    52 Patchy, 2018 other, 2020 0.21398557 lost
## 8  -1.14677923    52 Patchy, 2018 other, 2020 0.37341319 lost
## 9  -1.26787841    53 Patchy, 2018 other, 2020 2.86281831 lost
## 10 -1.59664105    57 Patchy, 2018 other, 2020 2.40369859 lost
##                          geometry
## 1  POLYGON ((420115.7 1214353,...
## 2  POLYGON ((429474.7 1229396,...
## 3  POLYGON ((439513.2 1225601,...
## 4  POLYGON ((441174.7 1225210,...
## 5  POLYGON ((441836.9 1224821,...
## 6  POLYGON ((426139.4 1204520,...
## 7  POLYGON ((426142.6 1204596,...
## 8  POLYGON ((426455.1 1204737,...
## 9  POLYGON ((426557.5 1204604,...
## 10 POLYGON ((421091.5 1196938,...

Now the change layer with depth estimates is combined with segment and management polygons to summarize by areas of interest.

# intersect change with relevant boundaries
chgdatarea <- chgdatarea %>% 
  st_intersection(tbseg) %>% 
  st_intersection(sgmanagement) %>% 
  unite('allbnds', bay_segment, areas, sep = ', ', remove = F) %>% 
  mutate(
    bay_segment = factor(bay_segment, levels = c('OTB', 'HB', 'MTB', 'LTB'))
  )
chgdatarea
## Simple feature collection with 2620 features and 10 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 416881.4 ymin: 1157848 xmax: 528821.4 ymax: 1345074
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## First 10 features:
##        meandepth idval       source       target       Acres    var
## 114   0.11488018   629 Patchy, 2018  other, 2020 50.90067840   lost
## 235   0.04889106  1055 Patchy, 2018  other, 2020  9.64494620   lost
## 240  -0.17614305  1059 Patchy, 2018  other, 2020 96.57784211   lost
## 1287 -0.39206555  2729 Patchy, 2018  other, 2020  2.87894694   lost
## 1288 -0.43271106  2729 Patchy, 2018  other, 2020  0.24565248   lost
## 1289 -0.45060265  2729 Patchy, 2018  other, 2020  0.06399648   lost
## 1290 -0.40984065  2729 Patchy, 2018  other, 2020 13.52472502   lost
## 1291 -0.20258051  2729 Patchy, 2018  other, 2020  1.27926961   lost
## 1292 -0.11987968  2729 Patchy, 2018  other, 2020  1.18314739   lost
## 1836 -0.43993417  4222  other, 2018 Patchy, 2020  0.12147729 gained
##             long_name allbnds bay_segment areas                       geometry
## 114  Hillsborough Bay   HB, 2          HB     2 MULTIPOLYGON (((497592.4 13...
## 235  Hillsborough Bay   HB, 2          HB     2 POLYGON ((499438.1 1293914,...
## 240  Hillsborough Bay   HB, 2          HB     2 POLYGON ((503023.2 1306953,...
## 1287 Hillsborough Bay   HB, 3          HB     3 POLYGON ((521900.9 1295906,...
## 1288 Hillsborough Bay   HB, 3          HB     3 POLYGON ((522238.2 1295609,...
## 1289 Hillsborough Bay   HB, 3          HB     3 POLYGON ((522404.3 1294798,...
## 1290 Hillsborough Bay   HB, 3          HB     3 POLYGON ((525674.3 1281558,...
## 1291 Hillsborough Bay   HB, 3          HB     3 POLYGON ((526041.1 1281997,...
## 1292 Hillsborough Bay   HB, 3          HB     3 POLYGON ((526250.7 1283956,...
## 1836 Hillsborough Bay   HB, 3          HB     3 POLYGON ((522074.6 1295830,...

A tabular summary of acreage gained/lost by major bay segment:

chgdatarea %>% 
  st_set_geometry(NULL) %>% 
  group_by(bay_segment, var) %>%
  summarise(
    Acres = sum(Acres, na.rm = T), 
    .groups = 'drop'
  ) %>% 
  spread(var, Acres) %>% 
  mutate(
    change = gained - lost
  ) %>% 
  mutate_if(is.numeric, round, 0) %>% 
  knitr::kable()
bay_segment gained lost change
OTB 154 4474 -4320
HB 129 500 -371
MTB 102 543 -441
LTB 458 485 -27

The plot below shows the cumulative distribution of acres lost and gained by major bay segments.

ggplot(chgdatarea, aes(x = meandepth, color = var)) +
  stat_ecdf() +
  scale_color_manual(values = cols) + 
  facet_wrap(~bay_segment, ncol = 4) +
  theme_minimal() +
  labs(
    color = NULL,
    x = 'Mean depth below MLLW (m)', 
    y = 'ECDF'
  )

Because the loss/gain polygons vary in area and the average depth is based on the areal average, we can weight the cumulative distribution curves to favor polygons with more area.

ggplot(chgdatarea, aes(x = meandepth, color = var, weight = Acres)) +
  stat_ecdfwt() +
  scale_color_manual(values = cols) + 
  facet_wrap(~bay_segment, ncol = 4) +
  theme_minimal() +
  labs(
    color = NULL,
    x = 'Mean depth below MLLW (m)', 
    y = 'ECDF'
  )

We can show a similar result but mapped to the major polygon segment boundaries. In addition, we can show the difference between the weighted averages for polygons where seagrass was gained vs lost. This provides an overview of the net change in seagrass depth distributions, i.e., positive values show a net shift to shallower waters and negative values show a net shift to deeper waters. The net change in seagrass depth is also weighted by total acreage lost by polygon boundary (e.g., OTB, HB, etc.).

# color function
maxv <- 0.97
colfun <- colorNumeric(
  palette = c(rev(colgrn), colred),
  domain = c(-1 * maxv, maxv)
)

# summarize depth values by segment
# take difference to get net change
tomap <- chgdatarea %>%  
  st_set_geometry(NULL) %>% 
  group_by(bay_segment, var) %>% 
  mutate(
    wgt = sum(Acres, na.rm  = T)
  ) %>% 
  summarise(
    meandepth = weighted.mean(meandepth, w = Acres),
    wgt = unique(wgt),
    .groups = 'drop'
  ) %>%
  group_by(var) %>% 
  mutate(
    wgt = wgt / sum(wgt)
  ) %>% 
  ungroup() %>% 
  mutate(
    meandepth = meandepth * wgt
  ) %>% 
  dplyr::select(-wgt) %>% 
  spread(var, meandepth) %>% 
  mutate(
    sgdlt = gained - lost
  ) %>% 
  inner_join(tbseg, ., by = 'bay_segment') %>% 
  mutate(
    fillhx = colfun(sgdlt)
  ) %>% 
  st_transform(crs = 4326)

# map
mapin <- mapview(tomap, homebutton = F, popup = NULL, legend = F) %>% 
  .@map %>% 
  leafem::removeMouseCoordinates() %>% 
  clearShapes()
vls <- seq(-1 * maxv, maxv, length.out = 10)
mapin %>% 
  addPolygons(
    data = tomap,
    stroke = T,
    color = 'black',
    weight = 1,
    layerId = ~bay_segment,
    fillColor = ~fillhx,
    fillOpacity = 0.8
  ) %>% 
  addLegend("topright", pal = colfun, values = vls, opacity = 0.8, title = 'Net change (m)')


The plot below shows the cumulative distribution of acres lost and gained by major bay segments and seagrass management areas.

ggplot(chgdatarea, aes(x = meandepth, color = var)) +
  ggplot2::stat_ecdf() +
  scale_color_manual(values = cols) + 
  facet_wrap(~allbnds) +
  theme_minimal() +
  labs(
    color = NULL,
    x = 'Mean depth below MLLW (m)', 
    y = 'ECDF'
  )

Now repeat but weight the cumulative distribution curve to favor loss/gain polygons with more area.

ggplot(chgdatarea, aes(x = meandepth, color = var, weight = Acres)) +
  stat_ecdfwt() +
  scale_color_manual(values = cols) + 
  facet_wrap(~allbnds) +
  theme_minimal() +
  labs(
    color = NULL,
    x = 'Mean depth below MLLW (m)', 
    y = 'ECDF'
  )

Again, create a similar map of net depth changes by bay segment but instead show the results by seagrass management areas.

# color function
maxv <- 0.4
colfun <- colorNumeric(
  palette = c(rev(colgrn), colred),
  domain = c(-1 * maxv, maxv)
)

# summarize depth values by segment
# take difference to get net change
tomap <- chgdatarea %>%  
  st_set_geometry(NULL) %>% 
  group_by(areas, var) %>% 
  mutate(
    wgt = sum(Acres, na.rm  = T)
  ) %>% 
  summarise(
    meandepth = weighted.mean(meandepth, w = Acres),
    wgt = unique(wgt),
    .groups = 'drop'
  ) %>%
  group_by(var) %>% 
  mutate(
    wgt = wgt / sum(wgt)
  ) %>% 
  ungroup() %>% 
  mutate(
    meandepth = meandepth * wgt
  ) %>% 
  dplyr::select(-wgt) %>% 
  spread(var, meandepth) %>% 
  mutate(
    sgdlt = gained - lost
  ) %>%   
  inner_join(sgmanagement, ., by = 'areas') %>% 
  mutate(
    fillhx = colfun(sgdlt)
  ) %>% 
  filter(!is.na(sgdlt)) %>% 
  st_transform(crs = 4326)

# map
mapin <- mapview(tomap, homebutton = F, popup = NULL, legend = F) %>% 
  .@map %>% 
  leafem::removeMouseCoordinates() %>% 
  clearShapes()
vls <- seq(-1 * maxv, maxv, length.out = 10)
mapin %>% 
  addPolygons(
    data = tomap,
    stroke = T,
    color = 'black',
    weight = 1,
    layerId = ~areas,
    fillColor = ~fillhx,
    fillOpacity = 0.8
  ) %>% 
  addLegend("topright", pal = colfun, values = vls, opacity = 0.8, title = 'Net change (m)')

The depth changes are super weird - most are opposite to our expectations that seagrass changes would be in more shallow areas if water quality conditions have changed between years. For context, water quality can also be evaluate. Changes in secchi depth (m), turbidity (NTU), and color 345 (Pt-Co units) are now evaluated, with a focus on years 2017 and 2019 since seagrass coverage years for 2018 and 2020 were taken from images captures in late fall/early winter of these years.

First, the water quality data are plotted as annual averages to get a sense of general trends.

# summarize water quality data
toplo <- epcdata %>% 
  dplyr::select(yr, bay_segment, sd = sd_raw_m, turb = `Turbidity_JTU-NTU`, col = Color_345_F45_PCU) %>% 
  gather('var', 'val', sd, turb, col) %>% 
  group_by(yr, bay_segment, var) %>% 
  filter(!grepl('[:alpha:]', val)) %>% 
  mutate(val = as.numeric(val)) %>% 
  filter(yr > 1990 & yr < 2020) %>% 
  na.omit %>% 
  summarise(
    avev = mean(val, na.rm = T), .groups = 'drop', 
    lov = t.test(val)$conf.int[1], 
    hiv = t.test(val)$conf.int[2]
  ) %>% 
  arrange(bay_segment, yr) %>% 
  group_by(bay_segment, var) %>% 
  mutate(
    bay_segment = factor(bay_segment, levels = c('OTB', 'HB', 'MTB', 'LTB')), 
    var = factor(var, levels = c('sd', 'turb', 'col'), labels = c('Secchi (m)', 'Turb (NTU)', 'Color (Pt-Co)')), 
    yrcat = case_when(
      yr %in% c(2017, 2019) ~ '2018, 2020 seagrass', 
      T ~ ''
    )
  ) %>% 
  ungroup()

# plot
ggplot(toplo, aes(x = yr, y = avev)) + 
  geom_ribbon(aes(ymin = lov, ymax = hiv), alpha = 0.3) +
  geom_line() + 
  geom_point(aes(size = yrcat, color = yrcat)) +
  facet_grid(var~bay_segment, scales = 'free_y') +
  theme_minimal(base_size = 14) + 
  theme(
    strip.background = element_blank(), 
    axis.title.x = element_blank(), 
    axis.text.x = element_text(size = 9, angle = 40, hjust = 1), 
    legend.position = 'top'
  ) + 
  scale_size_manual(values = c(-1, 2)) +
  scale_color_manual(values = c('black', 'tomato1')) +
  labs(
    y = 'Annual averages', 
    color = NULL, 
    size = NULL
  )

Just comparing the means and confidence intervals for 2017 and 2020:

# summarize water quality data
toplo <- epcdata %>% 
  dplyr::select(yr, bay_segment, sd = sd_raw_m, turb = `Turbidity_JTU-NTU`, col = Color_345_F45_PCU) %>% 
  gather('var', 'val', sd, turb, col) %>% 
  filter(!grepl('[:alpha:]', val)) %>% 
  filter(yr %in% c(2017, 2019)) %>% 
  mutate(
    val = as.numeric(val), 
    bay_segment = factor(bay_segment, levels = c('OTB', 'HB', 'MTB', 'LTB')), 
    var = factor(var, levels = c('sd', 'turb', 'col'), labels = c('Secchi (m)', 'Turb (NTU)', 'Color (Pt-Co)'))
  ) %>% 
  group_by(yr, var, bay_segment) %>% 
  summarise(
    lov = t.test(val)$conf.int[1],
    hiv = t.test(val)$conf.int[2],
    val = mean(val, na.rm = T), 
    .groups = 'drop'
  )

# plot
ggplot(toplo, aes(x = factor(yr), y = val, group = 1)) + 
  geom_errorbar(aes(ymin = lov, ymax = hiv), width = 0.5) +
  geom_line(linetype = 'dashed') + 
  geom_point(size = 3, color = 'tomato1') +
  facet_grid(var~bay_segment, scales = 'free_y') +
  theme_minimal(base_size = 14) + 
  theme(
    strip.background = element_blank(), 
    axis.title.x = element_blank(), 
    axis.text.x = element_text(size = 9, angle = 40, hjust = 1), 
    legend.position = 'top'
  ) +
  labs(
    y = 'Annual averages (+/- 95% CI)', 
    color = NULL, 
    size = NULL
  )

The water quality data can be aggregated by spatial units of interest. In this case, the same parameters are aggregated by seagrass management areas and the percent change between 2017 and 2019 is assessed.

wqdat <- epcdata %>% 
  dplyr::select(station = epchc_station, SampleTime, yr, mo, lat = Latitude, lng = Longitude, sd = sd_raw_m, turb = `Turbidity_JTU-NTU`, col = Color_345_F45_PCU) %>% 
  filter(yr %in% c(2019, 2017)) %>% 
  st_as_sf(coords = c('lng', 'lat'), crs = 4326) %>% 
  st_transform(crs = prj4) %>% 
  st_intersection(tbseg) %>% 
  st_intersection(sgmanagement) %>% 
  st_set_geometry(NULL) %>% 
  unite('allbnds', bay_segment, areas, sep = ', ', remove = T) %>% 
  dplyr::select(yr, allbnds, sd, turb, col) %>% 
  gather('var', 'val', sd, turb, col) %>% 
  filter(!grepl('[:alpha:]', val)) %>% 
  mutate(val = as.numeric(val)) %>% 
  group_by(yr, allbnds, var) %>% 
  summarise(val = mean(val, na.rm = T), .groups = 'drop') %>% 
  spread(yr, val) %>% 
  mutate(
    wqdlt = (`2019` - `2017`) / `2017`
  ) %>% 
  separate(allbnds, c('bay_segment', 'areas'), sep = ', ', remove = F)

wqdat
## # A tibble: 39 x 7
##    allbnds bay_segment areas var   `2017` `2019`   wqdlt
##    <chr>   <chr>       <chr> <chr>  <dbl>  <dbl>   <dbl>
##  1 HB, 2   HB          2     col    14.2   19.8   0.389 
##  2 HB, 2   HB          2     sd      1.21   1.34  0.106 
##  3 HB, 2   HB          2     turb    4.23   2.40 -0.432 
##  4 HB, 3   HB          3     col    15.8   17.3   0.0982
##  5 HB, 3   HB          3     sd      1.12   1.35  0.209 
##  6 HB, 3   HB          3     turb    4.82   2.58 -0.464 
##  7 HB, 4   HB          4     col    13.6   16.3   0.197 
##  8 HB, 4   HB          4     sd      1.25   1.6   0.28  
##  9 HB, 4   HB          4     turb    4.69   2.71 -0.423 
## 10 HB, 6   HB          6     col     8.53  13.8   0.621 
## # ... with 29 more rows

Maps of the above percent changes by water quality parameter and seagrass management areas.

Secchi:

var <- 'sd'
# color function
maxv <- wqdat %>% 
  filter(var == !!var) %>% 
  pull(wqdlt) %>% 
  abs %>% 
  max
colfun <- colorNumeric(
  palette = c(rev(colred), colgrn),
  domain = c(-1 * maxv, maxv)
)

# summarize depth values by segment
# take difference to get net change
tomap <- wqdat %>%  
  dplyr::select(areas, var, wqdlt) %>% 
  filter(var == !!var) %>% 
  mutate(
    fillhx = colfun(wqdlt), 
    areas = as.numeric(areas)
  ) %>% 
  inner_join(sgmanagement, ., by = 'areas') %>%
  st_transform(crs = 4326)

# map
mapin <- mapview(tomap, homebutton = F, popup = NULL, legend = F) %>% 
  .@map %>% 
  leafem::removeMouseCoordinates() %>% 
  clearShapes()
vls <- seq(-1 * maxv, maxv, length.out = 10)
mapin %>% 
  addPolygons(
    data = tomap,
    stroke = T,
    color = 'black',
    weight = 1,
    layerId = ~areas,
    fillColor = ~fillhx,
    fillOpacity = 0.8
  ) %>% 
  addLegend("topright", pal = colfun, values = vls, opacity = 0.8, title = 'Secchi % change')

Turbidity:

var <- 'turb'
# color function
maxv <- wqdat %>% 
  filter(var == !!var) %>% 
  pull(wqdlt) %>% 
  abs %>% 
  max
colfun <- colorNumeric(
  palette = c(rev(colgrn), colred),
  domain = c(-1 * maxv, maxv)
)

# summarize depth values by segment
# take difference to get net change
tomap <- wqdat %>%  
  dplyr::select(areas, var, wqdlt) %>% 
  filter(var == !!var) %>% 
  mutate(
    fillhx = colfun(wqdlt), 
    areas = as.numeric(areas)
  ) %>% 
  inner_join(sgmanagement, ., by = 'areas') %>%
  st_transform(crs = 4326)

# map
mapin <- mapview(tomap, homebutton = F, popup = NULL, legend = F) %>% 
  .@map %>% 
  leafem::removeMouseCoordinates() %>% 
  clearShapes()
vls <- seq(-1 * maxv, maxv, length.out = 10)
mapin %>% 
  addPolygons(
    data = tomap,
    stroke = T,
    color = 'black',
    weight = 1,
    layerId = ~areas,
    fillColor = ~fillhx,
    fillOpacity = 0.8
  ) %>% 
  addLegend("topright", pal = colfun, values = vls, opacity = 0.8, title = 'Turb. % change')

Color:

var <- 'col'
# color function
maxv <- wqdat %>% 
  filter(var == !!var) %>% 
  pull(wqdlt) %>% 
  abs %>% 
  max
colfun <- colorNumeric(
  palette = c(rev(colgrn), colred),
  domain = c(-1 * maxv, maxv)
)

# summarize depth values by segment
# take difference to get net change
tomap <- wqdat %>%  
  dplyr::select(areas, var, wqdlt) %>% 
  filter(var == !!var) %>% 
  mutate(
    fillhx = colfun(wqdlt), 
    areas = as.numeric(areas)
  ) %>% 
  inner_join(sgmanagement, ., by = 'areas') %>%
  st_transform(crs = 4326)

# map
mapin <- mapview(tomap, homebutton = F, popup = NULL, legend = F) %>% 
  .@map %>% 
  leafem::removeMouseCoordinates() %>% 
  clearShapes()
vls <- seq(-1 * maxv, maxv, length.out = 10)
mapin %>% 
  addPolygons(
    data = tomap,
    stroke = T,
    color = 'black',
    weight = 1,
    layerId = ~areas,
    fillColor = ~fillhx,
    fillOpacity = 0.8
  ) %>% 
  addLegend("topright", pal = colfun, values = vls, opacity = 0.8, title = 'Color % change')

The changes in water quality measurements can be compared to seagrass depth change estimates. The same analysis as above to evaluate net depth changes by management areas is repeated, then joined with the water quality percent change estimates by seagrass management areas where water quality data are available.

# summarize seagrass net change
sgdat <- chgdatarea %>% 
  st_set_geometry(NULL) %>% 
  group_by(allbnds, var) %>% 
  mutate(
    wgt = sum(Acres, na.rm  = T)
  ) %>% 
  summarise(
    meandepth = weighted.mean(meandepth, w = Acres),
    wgt = unique(wgt),
    .groups = 'drop'
  ) %>%
  group_by(allbnds) %>% 
  mutate(
    wgt = wgt / sum(wgt)
  ) %>% 
  ungroup() %>% 
  mutate(
    meandepth = meandepth * wgt
  ) %>% 
  dplyr::select(-wgt) %>% 
  spread(var, meandepth) %>% 
  mutate(
    sgdlt = gained - lost
  ) %>% 
  na.omit() %>% 
  separate(allbnds, c('bay_segment', 'areas'), sep = ', ', remove = F)

# join with water quality 
toplo <- inner_join(wqdat, sgdat, by = c('allbnds', 'bay_segment', 'areas')) %>% 
  dplyr::select(allbnds, bay_segment, var, sgdlt, wqdlt) %>% 
  mutate(
    bay_segment = factor(bay_segment, levels = c('OTB', 'HB', 'MTB', 'LTB')), 
    var = factor(var, levels = c('sd', 'turb', 'col'), labels = c('Secchi (m)', 'Turb (NTU)', 'Color (Pt-Co)'))
  )

# plot
ggplot(toplo, aes(x = wqdlt, y = sgdlt)) + 
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  geom_point(aes(color = bay_segment), size = 3) + 
  facet_wrap(~var, scales = 'free_x') +
  geom_smooth(method = 'lm', color = 'black') +
  theme_minimal(base_size = 16) + 
  theme(
    strip.background = element_blank(), 
    axis.text.x = element_text(size = 9, angle = 40, hjust = 1), 
    legend.position = 'top', 
    panel.grid.minor = element_blank()#, 
    # panel.grid = element_blank()
  ) +
  labs(
    y = 'Seagrass net depth change (m)', 
    color = 'Bay segment', 
    x = '% change'
  )

Just for fun, let’s plot TN load estimates by source in Old Tampa Bay.

levs <- c('NPS', 'IPS', 'GWS', 'DPS', 'AD')
names(levs) <- c('Non-point source', 'Industrial', 'Groundwater', 'Domestic', 'Atmospheric')

toplo <- tnanndat %>% 
  filter(bay_segment %in% c('Old Tampa Bay')) %>%
  mutate(
    source = factor(source,
      levels = c('NPS', 'IPS', 'GWS', 'DPS', 'AD'),
      labels = c('Non-point source', 'Industrial', 'Groundwater', 'Domestic', 'Atmospheric')
    ), 
    source = reorder(source, -tn_load),
    reg = source == 'Non-point source'
  )

p <- ggplot(toplo, aes(x = year, y = tn_load, color = source)) + 
  geom_line(size = 1) + 
  geom_smooth(data = toplo[toplo$source == 'Non-point source', ], method = 'lm', se = F, color = 'black', linetype = 'dashed') +
  # scale_color_brewer(palette = 'Dark2') +
  scale_y_continuous(expand = c(0, 0)) + 
  labs(
    y = 'TN (tons / yr)', 
    x = NULL, 
    color = NULL
  ) + 
  theme_minimal() + 
  theme(
    panel.grid.major.x = element_blank(), 
    panel.grid.minor.x = element_blank(), 
    panel.grid.minor.y = element_blank(), 
    axis.ticks.x = element_line(), 
    legend.position = 'top'
  )
p