Skip to contents

Background

The Tampa Bay Benthic Index (TBBI) [1,2]) is an assessment method that quantifies the ecological health of the benthic community in Tampa Bay. The index provides a complementary approach to evaluating environmental condition that is supported by other assessment methods currently available for the region (e.g., water quality report card, nekton index, etc.). The tbeptools package includes several functions described below to import data required for the index and plot the results to view trends over time. Each of the functions are described in detail below.

The TBBI uses data from the Tampa Bay Benthic Monitoring Program as part of the Environmental Protection Commission (EPC) of Hillsborough Country. The data are updated annually on a public site maintained by EPC, typically in December after Summer/Fall sampling. This is the same website that hosts water quality data used for the water quality report card. The required data for the TBBI are more extensive than the water quality report card and the data are made available as a zipped folder of csv files, available here. The process for downloading and working with the data are similar as for the other functions in tbeptools.

Data import and included datasets

Data for calculating TBBI scores can be imported into the current R session using the read_importbenthic() function. This function downloads the zipped folder of base tables used for the TBBI from the EPC site if the data have not already been downloaded to the location specified by the input arguments.

To download the data with tbeptools, first create a character path for the location where you want to download the zipped files. If one does not exist, specify a location and name for the downloaded file. Here, we name the folder benthic.zip and download it on our desktop.

path <- '~/Desktop/benthic.zip'
benthicdata <- read_importbenthic(path)

Running the above code will return the following error:

#> Error in read_importbenthic() : File at path does not exist, use download_latest = TRUE

We get an error message from the function indicating that the file is not found. This makes sense because the file does not exist yet, so we need to tell the function to download the latest file. This is done by changing the download_latest argument to TRUE (the default is FALSE).

benthicdata <- read_importbenthic(path, download_latest = T)
#> File ~/Desktop/benthic.zip does not exist, replacing with downloaded file...

#> trying URL 'https://epcbocc.sharepoint.com/:u:/s/Share/EdY7IaU75kFPnVhyoe1yITcBESvFZNAreBwabwyL8EfuqQ?e=1i8Rnf&download=1' length 37122877 bytes (35.4 MB)

Now we get an indication that the file on the server is being downloaded. When the download is complete, we’ll have the data downloaded and saved to the benthicdata object in the current R session.

If we try to run the function again after downloading the data, we get the following message. This check is done to make sure that the data are not unnecessarily downloaded if the current file matches the file on the server.

benthicdata <- read_importbenthic(path, download_latest = T)
#> File is current..

Every time that tbeptools is used to work with the benthic data, read_importbenthic() should be used to import the data. You will always receive the message File is current... if your local file matches the one on the server. However, data are periodically updated and posted on the server. If download_latest = TRUE and your local file is out of date, you will receive the following message:

#> Replacing local file with current...

Calculating TBBI scores

After the data are imported, you can view them from the assigned object. The data are provided as a nested tibble that includes three different datasets: station information, field sample data (salinity), and detailed taxa information.

benthicdata
#> # A tibble: 3 × 2
#>   name         value                 
#>   <chr>        <list>                
#> 1 stations     <tibble [4,723 × 9]>  
#> 2 fieldsamples <tibble [5,197 × 3]>  
#> 3 taxacounts   <tibble [145,171 × 7]>

The individual datasets can be viewed by extracting them from the parent object using the deframe() function from the tibble package.

# see all
deframe(benthicdata)
#> $stations
#> # A tibble: 4,723 × 9
#>    StationID AreaAbbr Funding…¹ Progr…² Progr…³ Latit…⁴ Longi…⁵ date          yr
#>        <int> <chr>    <chr>       <int> <chr>     <dbl>   <dbl> <date>     <dbl>
#>  1       448 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  2       449 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#>  3       450 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  4       451 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  5       452 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#>  6       453 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  7       454 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  8       455 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  9       456 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#> 10       457 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#> # … with 4,713 more rows, and abbreviated variable names ¹​FundingProject,
#> #   ²​ProgramID, ³​ProgramName, ⁴​Latitude, ⁵​Longitude
#> 
#> $fieldsamples
#> # A tibble: 5,197 × 3
#>    StationID date       Salinity
#>        <int> <date>        <dbl>
#>  1       448 2002-05-21     30.5
#>  2       449 2002-05-20     30.8
#>  3       450 2002-05-21     30.5
#>  4       451 2002-05-21     30.6
#>  5       452 2002-05-20     30  
#>  6       453 2002-05-21     30.6
#>  7       454 2002-05-21     30.5
#>  8       455 2002-05-21     30.3
#>  9       456 2002-05-20     30.0
#> 10       457 2002-05-20     29.9
#> # … with 5,187 more rows
#> 
#> $taxacounts
#> # A tibble: 145,171 × 7
#>    StationID TaxaCountID TaxaListID FAMILY           NAME        TaxaC…¹ AdjCo…²
#>        <int>       <int>      <int> <chr>            <chr>         <dbl>   <dbl>
#>  1     11321      152055       1198 Anomiidae        Anomia sim…       1      25
#>  2     11321      152056       1224 Carditidae       Cardites f…       1      25
#>  3       584      233277        624 Cirratulidae     Chaetozone…       1      25
#>  4     11321      152057       1306 Veneridae        Chione ele…       3      75
#>  5      3132      239443        352 Pilargidae       Sigambra b…      14     350
#>  6     11321      152058       1088 Pyramidellidae   Turbonilla…       6     150
#>  7      2541      140530        258 NULL             NEMERTEA          8     200
#>  8      2542      140531        258 NULL             NEMERTEA          2      50
#>  9      2992      140631        262 NULL             Palaeoneme…       4     100
#> 10      3003      140635        285 Tetrastemmatidae Tetrastemm…       2      50
#> # … with 145,161 more rows, and abbreviated variable names ¹​TaxaCount,
#> #   ²​AdjCount

# get only station dat
deframe(benthicdata)[['stations']]
#> # A tibble: 4,723 × 9
#>    StationID AreaAbbr Funding…¹ Progr…² Progr…³ Latit…⁴ Longi…⁵ date          yr
#>        <int> <chr>    <chr>       <int> <chr>     <dbl>   <dbl> <date>     <dbl>
#>  1       448 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  2       449 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#>  3       450 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  4       451 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  5       452 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#>  6       453 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  7       454 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  8       455 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  9       456 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#> 10       457 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#> # … with 4,713 more rows, and abbreviated variable names ¹​FundingProject,
#> #   ²​ProgramID, ³​ProgramName, ⁴​Latitude, ⁵​Longitude

The anlz_tbbiscr() function uses the nested benthicdata to estimate the TBBI scores at each site. The TBBI scores typically range from 0 to 100 and are grouped into categories that describe the general condition of the benthic community. Scores less than 73 are considered “degraded”, scores between 73 and 87 are “intermediate”, and scores greater than 87 are “healthy”. Locations that were sampled but no organisms were found are assigned a score of zero and a category of “empty sample”. The total abundance (TotalAbundance, organisms/m2), species richness (SpeciesRichness) and bottom salinity (Salinity, psu) are also provided. Some metrics for the TBBI are corrected for salinity and bottom measurements taken at the time of sampling are required for accurate calculation of the TBBI.

tbbiscr <- anlz_tbbiscr(benthicdata)
tbbiscr
#> # A tibble: 4,723 × 14
#>    StationID AreaAbbr Funding…¹ Progr…² Progr…³ Latit…⁴ Longi…⁵ date          yr
#>        <int> <chr>    <chr>       <int> <chr>     <dbl>   <dbl> <date>     <dbl>
#>  1       448 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  2       449 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#>  3       450 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  4       451 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  5       452 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#>  6       453 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  7       454 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  8       455 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-21  2002
#>  9       456 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#> 10       457 MTB      Apollo B…       4 Benthi…    27.8   -82.4 2002-05-20  2002
#> # … with 4,713 more rows, 5 more variables: TotalAbundance <dbl>,
#> #   SpeciesRichness <dbl>, TBBI <dbl>, TBBICat <chr>, Salinity <dbl>, and
#> #   abbreviated variable names ¹​FundingProject, ²​ProgramID, ³​ProgramName,
#> #   ⁴​Latitude, ⁵​Longitude

Plotting results

The TBBI scores can be viewed as annual averages for each bay segment using the show_tbbimatrix() function. The show_tbbimatrix() plots the annual bay segment averages as categorical values in a conventional “stoplight” graphic. A baywide estimate is also returned, one based on all samples across all locations (“All”) and another weighted by the relative surface areas of each bay segment (“All (wt)”). The input to show_tbbimatrix() function is the output from the anlz_tbbiscr() function.

The matrix can also be produced as a plotly interactive plot by setting plotly = TRUE inside the function.

show_tbbimatrix(tbbiscr, plotly = T)

Additional sediment data

In addition to biological data, sediment contaminant concentrations are measured at sites within Tampa Bay. These include over 100 different constituents grouped broadly as metals, organics, physical, or other. The concentrations of these constituents can be compared relative to Threshold Effects Levels (TEL) or Potential Effects Levels (PEL), when available, as relative indications of the likelihood that the concentrations will have toxic effects on invertebrates that inhabit the sediments. The functions in tbeptools can be used to retrieve the sediment data and provide an indication of the concentrations relative to the TEL or PEL thresholds.

The read_importsediment() function will retrieve all sediment data for Tampa Bay collected annually by the Environmental Protection Commission of Hillsborough County. The data are retrieved from the same location as the biological data used to calculate the TBBI.

path <- '~/Desktop/sediment.zip'
sedimentdata <- read_importsediment(path, download_latest = T)

After the data are imported, you can view them from the assigned object.

sedimentdata
#> # A tibble: 216,627 × 24
#>    Progr…¹ Progr…² Fundi…³    yr AreaA…⁴ Stati…⁵ Stati…⁶ Latit…⁷ Longi…⁸ Repli…⁹
#>      <int> <chr>   <chr>   <int> <chr>     <int> <chr>     <dbl>   <dbl> <chr>  
#>  1       4 Benthi… TBEP     1993 HB         2463 93HB15     27.8   -82.4 no     
#>  2       4 Benthi… TBEP     1993 HB         2464 93HB16     27.8   -82.5 no     
#>  3       4 Benthi… TBEP     1993 HB         2470 93HB23     27.9   -82.4 no     
#>  4       4 Benthi… TBEP     1993 LTB        2471 93LTB24    27.7   -82.6 no     
#>  5       4 Benthi… TBEP     1993 LTB        2472 93LTB25    27.6   -82.6 no     
#>  6       4 Benthi… TBEP     1993 LTB        2473 93LTB26    27.6   -82.7 no     
#>  7       4 Benthi… TBEP     1993 LTB        2474 93LTB27    27.6   -82.7 no     
#>  8       4 Benthi… TBEP     1993 LTB        2475 93LTB28    27.6   -82.7 no     
#>  9       4 Benthi… TBEP     1993 LTB        2476 93LTB29    27.6   -82.6 no     
#> 10       4 Benthi… TBEP     1993 LTB        2477 93LTB30    27.6   -82.6 no     
#> # … with 216,617 more rows, 14 more variables: SedResultsType <chr>,
#> #   Parameter <chr>, ValueAdjusted <dbl>, Units <chr>, Qualifier <chr>,
#> #   MDLnum <dbl>, PQLnum <dbl>, TEL <dbl>, PEL <dbl>, BetweenTELPEL <chr>,
#> #   ExceedsPEL <chr>, PELRatio <dbl>, PreparationDate <chr>,
#> #   AnalysisTimeMerge <chr>, and abbreviated variable names ¹​ProgramId,
#> #   ²​ProgramName, ³​FundingProject, ⁴​AreaAbbr, ⁵​StationID, ⁶​StationNumber,
#> #   ⁷​Latitude, ⁸​Longitude, ⁹​Replicate

The show_sedimentmap() function can be used to create maps of selected parameters relative to TEL and PEL values. Green points show concentrations below the TEL, yellow points show concentrations between the TEL and PEL, and red points show concentrations above the PEL. The applicable TEL and PEL values for the parameter are indicated in the legend. The selected stations are those that are sampled in the years within the yrrng argument.

show_sedimentmap(sedimentdata, param = 'Arsenic', yrrng = c(1993, 2021))

A single year of data can be shown as well.

show_sedimentmap(sedimentdata, param = 'Arsenic', yrrng = 2021)

A map showing only the concentrations is returned if TEL and PEL values are not available for a parameter.

show_sedimentmap(sedimentdata, param = 'Selenium', yrrng = c(1993, 2021))

The PEL ratio can also be used to assess relative sediment quality given the measured contaminants. The show_sedimentpelmap() function creates a map of average PEL ratios graded from A to F for benthic stations monitored in Tampa Bay. The PEL ratio is the contaminant concentration divided by the Potential Effects Levels (PEL) that applies to a contaminant, if available. Higher ratios and lower grades indicate sediment conditions that are likely unfavorable for invertebrates. The station average combines the PEL ratios across all contaminants measured at a station.

show_sedimentpelmap(sedimentdata, yrrng = c(1993, 2021))

The average PEL ratios and grades used to create the map can also be returned as a data frame using anlz_sedimentpel().

anlz_sedimentpel(sedimentdata, yrrng = c(1993, 2021))
#> # A tibble: 2,189 × 7
#>       yr AreaAbbr StationNumber Latitude Longitude PELRatio Grade
#>    <int> <chr>    <chr>            <dbl>     <dbl>    <dbl> <fct>
#>  1  1993 HB       93HB15            27.8     -82.4  0.0157  B    
#>  2  1993 HB       93HB16            27.8     -82.5  0.0261  C    
#>  3  1993 HB       93HB23            27.9     -82.4  0.0174  B    
#>  4  1993 LTB      93LTB24           27.7     -82.6  0.0124  B    
#>  5  1993 LTB      93LTB25           27.6     -82.6  0.0189  B    
#>  6  1993 LTB      93LTB26           27.6     -82.7  0.00997 B    
#>  7  1993 LTB      93LTB27           27.6     -82.7  0.0125  B    
#>  8  1993 LTB      93LTB28           27.6     -82.7  0.0887  D    
#>  9  1993 LTB      93LTB29           27.6     -82.6  0.0350  C    
#> 10  1993 LTB      93LTB30           27.6     -82.6  0.0496  C    
#> # … with 2,179 more rows

Plots of bay segment averages of sediment concentrations for selected parameters can be created with show_sedimentave(). The plot includes appropriate lines for the TEL and PEL values, as well the grand mean across all segments. The former are omitted from the plot if unavailable for a selected parameter.

show_sedimentave(sedimentdata, param = 'Arsenic', yrrng = c(1993, 2021))

The same plot can be returned as an interactive plotly object using plotly = T.

show_sedimentave(sedimentdata, param = 'Arsenic', yrrng = c(1993, 2021), plotly = T)

The values used in the plot can be returned with anlz_sedimentave().

anlz_sedimentave(sedimentdata, param = 'Arsenic', yrrng = c(1993, 2021))
#> # A tibble: 7 × 8
#>   AreaAbbr   TEL   PEL Units   ave   lov   hiv grandave
#>   <fct>    <dbl> <dbl> <chr> <dbl> <dbl> <dbl>    <dbl>
#> 1 BCB       7.24  41.6 mg/kg  2.66 2.35   2.97     2.58
#> 2 HB        7.24  41.6 mg/kg  2.96 0.649  5.27     2.58
#> 3 LTB       7.24  41.6 mg/kg  2.92 2.45   3.39     2.58
#> 4 MR        7.24  41.6 mg/kg  1.66 1.36   1.96     2.58
#> 5 MTB       7.24  41.6 mg/kg  2.18 1.87   2.50     2.58
#> 6 OTB       7.24  41.6 mg/kg  3.29 0.704  5.87     2.58
#> 7 TCB       7.24  41.6 mg/kg  2.40 1.91   2.88     2.58

A similar plot of the bay segment averages for the average PEL ratios can be created with show_sedimentpelave(). The colors indicate the grades for A (green) to F (red).

show_sedimentpelave(sedimentdata, yrrng = c(1993, 2021))

The same plot can be returned as an interactive plotly object using plotly = T.

show_sedimentpelave(sedimentdata, yrrng = c(1993, 2021), plotly = T)

The values used in the plot can be returned with anlz_sedimentpelave().

anlz_sedimentpelave(sedimentdata, yrrng = c(1993, 2021))
#> # A tibble: 7 × 5
#>   AreaAbbr    ave    lov    hiv grandave
#>   <fct>     <dbl>  <dbl>  <dbl>    <dbl>
#> 1 BCB      0.0773 0.0590 0.0956   0.0691
#> 2 HB       0.133  0.0784 0.187    0.0691
#> 3 LTB      0.0342 0.0245 0.0438   0.0691
#> 4 MR       0.0629 0.0433 0.0825   0.0691
#> 5 MTB      0.0498 0.0375 0.0620   0.0691
#> 6 OTB      0.0690 0.0489 0.0890   0.0691
#> 7 TCB      0.0579 0.0285 0.0873   0.0691

References

[1] D.J. Karlen, T. Dix, B.K. Goetting, S.E. Markham, K.Campbell, J. Jernigan, J.Christian, K. Martinez, A. Chacour, Twenty-five-year Trends in the Benthic Community and Sediment Quality of Tampa Bay 1993-2017, Tampa Bay Estuary Program, St. Petersburg, Florida, 2020. https://drive.google.com/file/d/1XqwAbvhSRSDMHyGhVtrUTZpVmPqkFsWI/view?usp=drivesdk.

[2] D. Wade, K. Malloy, A. Janicki, R. Nijbroek, S. Grabe, Development of a Benthic Index to Establish Sediment Quality Targets for the Tampa Bay Estuary, Tampa Bay Estuary Program, St. Petersburg, Florida, 2006. https://drive.google.com/file/d/1jWYmbtBWACLXAM3F6xP9VqCD7qEbvOPF/view?usp=drivesdk.