Assigning OTB loadings to sub-segments

This document provides a brief explanation of load assignments to sub-segments of Old Tampa Bay based on likely sub-basins. It address the problem of total loads assigned to Old Tampa Bay as a whole. An interim dataset provided by Janicki Env. that assigns monthly loading to gaged and ungaged portions of Old Tampa Bay is used for the assignments.

First, the required libraries and loading dataset are loaded.

library(haven)
library(dplyr)
library(sf)
library(tibble)
library(tidyr)
library(leaflet)
library(mapedit)

# OTB sub-segments
load(file = here::here('data/otbsub.RData'))

# basin loading data
# original here T:/03_BOARDS_COMMITTEES/05_TBNMC/TB_LOADS/2022_RA_Deliverables/2017-2021Annual&MonthlyLoadDatasets/NPS1721/Monthly
dat <- read_sas(here::here('data-raw/allnuts1721monthbasin20221027.sas7bdat')) |> 
  filter(BAY_SEG == 1)
head(dat)

# A tibble: 6 × 10
  BAY_SEG basin    source  year month tnload tpload tssload BODLOAD  h2oload
    <dbl> <chr>    <chr>  <dbl> <dbl>  <dbl>  <dbl>   <dbl>   <dbl>    <dbl>
1       1 02306647 NPS     2017     1  0.270 0.0538    8.32    1.53  177942.
2       1 02306647 NPS     2017     2  0.355 0.0717   11.1     2.04  237138.
3       1 02306647 NPS     2017     3  0.280 0.0572    8.85    1.63  189218.
4       1 02306647 NPS     2017     4  0.225 0.0435    6.73    1.24  143787.
5       1 02306647 NPS     2017     5  0.195 0.0358    5.54    1.02  118492.
6       1 02306647 NPS     2017     6  3.37  0.683   106.     19.5  2260524.

The basin column describes the gaged and ungaged basins for each loading estimate. The basin 206-1 is ungaged.

unique(dat$basin)

[1] "02306647" "02307000" "206-1"    "LTARPON"

We can view a map of these basins using the dbasin shapefile.

# dbasins shapefile
# original here T:/05_GIS/BOUNDARIES
shp <- st_read(here::here('data-raw/TBEP_dBasins_Correct_Projection.shp'), quiet = T) |>
  filter(BAYSEGNAME == 'Old Tampa Bay') |>
  st_transform(crs = 4326)

tomap <- shp |> 
  group_by(NEWGAGE) |> 
  summarise()

pal <- colorFactor(
  palette = 'viridis', 
  domain = tomap$NEWGAGE
)

leaflet() |> 
  addProviderTiles('CartoDB.Positron') |>
  addPolygons(data = tomap, weight =  0.5, fillColor = ~pal(NEWGAGE), fillOpacity = 0.9) |>
  addPolygons(data = otbsub, label = ~paste(subseg)) |> 
  addLegend(
    pal = pal,
    values = tomap$NEWGAGE,
    opacity = 0.9,
    title = "Basin",
    position = "bottomright"
  )

Loadings from the gaged portions (02307000, 02306647, 02307359, LTARPON) can be assigned to sub-segments based on drainage location (02307359, LTARPON drain to the NW subsegment; 02307000, 02306647 drain to the NE subsegment). The loadings for the ungaged portion (206-1) can be proportionately assigned based on the approximate area that likely drains to each sub-segment. The sub-basins in the dbasin area are interactively assigned to each sub-segment based on location in 206-1.

# this is all done interactively, saved file is loaded below

# gaged portion with sub-segment assignment
gagebas <- shp |> 
 mutate(
   subsegment = case_when(
     NEWGAGE %in% 'LTARPON' ~ 'NW', 
     NEWGAGE %in% '02307359' ~ 'NW', 
     NEWGAGE %in% c('02307000', '02306647') ~ 'NE', 
     T ~ NA_character_   
   )
 ) |> 
 filter(!is.na(subsegment))

# create leaflet polygon map with ungage
m <- leaflet() |> 
 addProviderTiles('CartoDB.Positron') |> 
 addPolygons(data = shp[shp$NEWGAGE == '206-1', ], weight =  0.5) |> 
 addPolygons(data = otbsub) |> 
 addDrawToolbar(
   targetGroup = 'draw',
   editOptions = editToolbarOptions()
 )

# use interactive selection
NW <- selectFeatures(ungage, map = m)
NE <- selectFeatures(ungage, map = m)
CW <- selectFeatures(ungage, map = m)
CE <- selectFeatures(ungage, map = m)
SW <- selectFeatures(ungage, map = m)
SE <- selectFeatures(ungage, map = m)

ungagebas <- list(
 NW = NW,  
 NE = NE,
 CW = CW,
 CE = CE,
 SW = SW,
 SE = SE
 ) |> 
 enframe(name = 'subsegment', value = 'data') |> 
 unnest('data') |> 
 st_sf()

allsubbas <- bind_rows(ungagebas, gagebas) 

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

The previously created file is loaded. This is how the assignment looks. The original gaged and ungaged areas are outlined in red.

load(file = here::here('data/allsubbas.RData'))

pal <- colorFactor(
  palette = 'viridis', 
  domain = allsubbas$subsegment
)

leaflet() |> 
  addProviderTiles('CartoDB.Positron') |>
  addPolygons(data = allsubbas, weight =  0.5, fillOpacity = 0.9, fillColor = ~pal(subsegment)) |>
  addPolygons(data = otbsub, label = ~paste(subseg)) |> 
  addPolygons(data = tomap, weight =  2, fillOpacity = 0, color = 'red') |>
  addLegend(
    pal = pal,
    values = allsubbas$subsegment,
    opacity = 0.9,
    title = "Sub-segment",
    position = "bottomright"
  )

Next, a data object that shows how to proportionally assign the loads to each subsegment based on the gaged and ungaged portions is created. The gaged portions are given a 1, meaning their loadings are used as is, and the ungaged portion is given a proportion based on the area attributed to each subsegment.

# get tn load proportion multipliers based on basin area attributed to ungaged portions
props <- allsubbas |> 
  group_by(NEWGAGE, subsegment) |> 
  summarise() |> 
  ungroup()
props$acres <- as.numeric(st_area(props) / 4047)
props <- props |> 
  mutate(
    prop = case_when(
      NEWGAGE != '206-1' ~ 1,
      TRUE ~ acres / sum(acres[NEWGAGE == '206-1'])
    )
  ) |> 
  st_set_geometry(NULL)
props

# A tibble: 10 × 4
   NEWGAGE  subsegment  acres   prop
 * <chr>    <chr>       <dbl>  <dbl>
 1 02306647 NE         12231. 1     
 2 02307000 NE         28852. 1     
 3 02307359 NW         23273. 1     
 4 206-1    CE          5563. 0.0642
 5 206-1    CW         23367. 0.270 
 6 206-1    NE         33685. 0.389 
 7 206-1    NW         14502. 0.167 
 8 206-1    SE          8147. 0.0940
 9 206-1    SW          1411. 0.0163
10 LTARPON  NW         13734. 1

The props object is joined to the original loading data using the gaged/ungaged column in each. The load is proportioned for each year/month/location by multiplying the original by the prop column in props.

datprop <- dat |>
  select(basin, source, year, month, tnload, tpload) |>
  left_join(props, by = c('basin' = 'NEWGAGE'), relationship = 'many-to-many') |>
  mutate(
    tnload = tnload * prop, 
    tpload = tpload * prop
  )
datprop

# A tibble: 1,380 × 9
   basin    source  year month tnload tpload subsegment  acres  prop
   <chr>    <chr>  <dbl> <dbl>  <dbl>  <dbl> <chr>       <dbl> <dbl>
 1 02306647 NPS     2017     1  0.270 0.0538 NE         12231.     1
 2 02306647 NPS     2017     2  0.355 0.0717 NE         12231.     1
 3 02306647 NPS     2017     3  0.280 0.0572 NE         12231.     1
 4 02306647 NPS     2017     4  0.225 0.0435 NE         12231.     1
 5 02306647 NPS     2017     5  0.195 0.0358 NE         12231.     1
 6 02306647 NPS     2017     6  3.37  0.683  NE         12231.     1
 7 02306647 NPS     2017     7  3.58  0.708  NE         12231.     1
 8 02306647 NPS     2017     8  7.54  1.53   NE         12231.     1
 9 02306647 NPS     2017     9 12.2   2.50   NE         12231.     1
10 02306647 NPS     2017    10  2.20  0.456  NE         12231.     1
# ℹ 1,370 more rows

The total loading is the same between the original dataset and the proportionally assigned dataset.

sum(datprop$tnload)

[1] 1991.495

sum(dat$tnload)

[1] 1991.495

sum(datprop$tpload)

[1] 361.6535

sum(dat$tpload)

[1] 361.6535

This example uses data only for the 2017 to 2021 RA period. For use with CASM, the monthly data must be interpolated to a daily time step and previous years must be added.