knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)
library(sf)
library(RODBC)
library(dplyr)
library(here)This workflow was developed to update data used to estimate tidal creek assessment categories using FLDEP IWR run data. In this example, we use IWR Run 67, but in theory it should work with any IWR run. The only requirements are:
.accdb files for a selected run
from here: https://publicfiles.dep.state.fl.us/dear/iwr/Below we set a path to the extracted location for the IWR database, load the tidal creek line file, and load the wbid polygon file.
# set database path
pth <- "~/Desktop/iwr2026_Run67_2025-01-15/"
# # to explore
# con <- odbcConnectAccess2007('C:/Users/mbeck/Desktop/iwr2021_run61_2021-05-27/iwr_run61_05272021.accdb')
# sqlTables(con)
# sqlFetch(con, 'station list')
# source creek line file
tdlcrk <- st_read(here('shapes/TidalCreek_ALL_Line.shp'), quiet = T)
# wbid data, run 66 (url should not change between runs)
# wbid <- st_read('https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/WBIDS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson', quiet = T) %>%
# st_transform(crs = st_crs(tdlcrk))
wbid <- st_read('~/Desktop/WBIDs_Run67/WBIDs_Run67.shp', quiet = T) |> # used for 67 for manual download, online not yet updated, google found it
st_transform(crs = st_crs(tdlcrk))Now we import the IWR run data for the entire state using an ODBC
connection to the extracted folder. The following code just iterates
through the .accdb files in the extracted folders,
retrieves the data within a ten year window from the current run, and
pulls out parameters of interest. Note that 11 years are subset based on
the current year. This is done to accommodate IWR runs where incomplete
data are present for the current year. This is done in a loop using a
SQL request. It only takes a few minutes to run and does not use a lot
of memory.
# this retrieves most recent ten years of iwr data, including the current year and ten years prior
# filters only parameters in mcode below
# used ODBC connect to access file
# entire state
# .accdb files to query
dbs <- list.files(pth, pattern = 'rawDataDB.*\\.accdb$', full.names = T)
# the current year of interest, used to build the query for ten prior years
curyr <- 2025
# query building tools
mcode <- c("CHLAC","COLOR","COND","DO","DOSAT","NO3O2","ORGN","SALIN","TKN","TN","TP","TSS","TURB") %>%
paste0("masterCode = '", ., "'") %>%
paste0(., collapse = ' or ')
qry_fun <- function(tab, curyr){
out <- paste0("select sta, year, month, day, masterCode, result, rCode from ", tab, " where year <= ", curyr,
" and (", mcode, ")"
)
return(out)
}
# loop through the accdb files
iwr <- NULL
for(db in dbs){
cat(basename(db), '\n')
# setup the connection
con <- odbcConnectAccess2007(db)
# build the query
qry <- basename(db) %>%
gsub('DB|\\.accdb', '', .) %>%
qry_fun(curyr)
# retrieve the data
res <- sqlQuery(con, qry)
# add it to the output
iwr <- rbind(iwr, res)
}## rawDataDB1.accdb
## rawDataDB2.accdb
## rawDataDB3.accdb
## rawDataDB4.accdbhead(iwr)## sta year month day masterCode result rCode
## 1 112WRD 02231400 2021 9 1 COND 201.00000 <NA>
## 2 112WRD 02231400 2021 9 1 DO 5.30000 <NA>
## 3 112WRD 02231400 2021 9 1 ORGN 0.27500 U
## 4 112WRD 02231400 2021 9 1 TN 0.17425 +
## 5 112WRD 02234600 2015 8 24 COND 361.00000 <NA>
## 6 112WRD 02234600 2015 8 24 DO 6.00000 <NA>Now we need to get station lat/lon data for intersection with the tidal creek line layer. This is for the entire state and is retrieved from an ODBC connection to the extracted folder.
# from geojson for run66
# https://geodata.dep.state.fl.us/datasets/FDEP::impaired-waters-rule-iwr-stations/about
# staraw <- st_read('https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/IMPAIRED_WATERS/MapServer/1/query?outFields=*&where=1%3D1&f=geojson', quiet = T)
staraw <- st_read('~/Desktop/Stations_Run67/IWR_stations_Run67.shp', quiet = T) # used for 67 for manual download, online not yet updated, google found it
stas <- staraw %>%
filter(STATION_ID %in% iwr$sta) %>%
st_transform(crs = st_crs(tdlcrk))
# pull from access, currently doesn't work
# stas <- list.files(pth, pattern = '^iwr\\_run', full.names = T) %>%
# odbcConnectAccess2007 %>%
# sqlFetch('station list') %>%
# filter(STA %in% unique(iwr$sta)) %>%
# st_as_sf(coords = c('LONG', 'LAT'), crs = 4326) %>%
# st_transform(crs = st_crs(tdlcrk))
head(stas)## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -95894.59 ymin: 3362583 xmax: 434397.4 ymax: 3442532
## Projected CRS: NAD83 / UTM zone 17N
## STATION_ID WATERBODY_ STATION_NA
## 1 1110NET 110057 375H APALACHICOLA R. AT CHATTAHOOCHEE
## 2 1110NET 110062 10C ESCAMBIA RIVER AT CENTURY
## 3 1113MONA120062 10C ESCAMBIA RIVER AT CENTURY FLA.
## 4 1113MONA123929 60A CHATTAHOOCHEE R AT NEALS LANDING
## 5 1113MONO030126 2224A RIBAULT RIVER AT LEM TURNER RD BRIDGE
## 6 1113MONO030128 2224A RIBAULT RIVER AT FERNANDO RD BRIDGE
## COMMENTS LATITUDE LONGITUDE BEGIN_DATE END_DATE
## 1 <NA> 30.7089 -84.8622 <NA> <NA>
## 2 <NA> 30.9658 -87.2339 <NA> <NA>
## 3 <NA> 30.9658 -87.2339 <NA> <NA>
## 4 <NA> 30.9772 -85.0061 <NA> <NA>
## 5 Reassigned from 2224 to 2224A, Run 50 - P.Homann 30.3949 -81.6829 <NA> <NA>
## 6 Reassigned from 2224 to 2224A, Run 50 - P.Homann 30.3933 -81.7135 <NA> <NA>
## RESULT_COU FIRST_RUN WATER_QUAL BIOLOGY_ST geometry
## 1 <NA> <NA> 1 0 POINT (130040.7 3403715)
## 2 <NA> <NA> 1 0 POINT (-95894.59 3442532)
## 3 <NA> <NA> 1 0 POINT (-95894.59 3442532)
## 4 <NA> <NA> 1 0 POINT (117314.8 3433968)
## 5 <NA> <NA> 1 0 POINT (434397.4 3362743)
## 6 <NA> <NA> 1 0 POINT (431456.7 3362583)The next step is to create a polyline file that is similar to the
source file, but includes an intersection with the WBID layer. It also
includes creek ID (JEI), class, and creek length (total across WBIDs).
It is specific to southwest Florida. This reproduces the
tidalcreeks sf data object in tbeptools for the current
wBIDs and must be saved to the data folder for the
package.
# create tidalcreeks polyline file with wbid, class, jei info -------------
# intersect creek lines with wbids
tidalcreeks <- st_intersection(tdlcrk, wbid) %>%
st_transform(4326) %>%
arrange(CreekID) %>%
mutate(
id = 1:nrow(.)
) %>%
select(id, name = Name, JEI = CreekID, wbid = WBID, class = CLASS, Creek_Length_m = Total_m)
# fix rownames
row.names(tidalcreeks) <- 1:nrow(tidalcreeks)
# save to Desktop for tbeptools data folder
save(tidalcreeks, file = '~/Desktop/tidalcreeks.RData', compress = 'xz')
head(tidalcreeks)## Simple feature collection with 6 features and 6 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -82.34154 ymin: 26.89042 xmax: -82.29655 ymax: 26.95527
## Geodetic CRS: WGS 84
## id name JEI wbid class Creek_Length_m geometry
## 1 1 Rock Creek CC01 1983B 2 6443.083 LINESTRING (-82.33869 26.92...
## 2 2 Rock Creek CC01 2052 3M 6443.083 MULTILINESTRING ((-82.33006...
## 3 3 Oyster Creek CC02 1983B 2 10175.167 MULTILINESTRING ((-82.33188...
## 4 4 Oyster Creek CC02 2067 3M 10175.167 MULTILINESTRING ((-82.3194 ...
## 5 5 Buck Creek CC03 1983B 2 2251.188 LINESTRING (-82.31512 26.89...
## 6 6 Buck Creek CC03 2068A 3M 2251.188 LINESTRING (-82.30026 26.89...The final step is to extract the IWR station data at the state level
to tidal creeks in southwest Florida. This is done by buffering the
original creek line file (200m, flat ends), intersecting the buffered
file with the station lat/lon file, and using an inner join by station.
Finally, the WBID class is added. This is reproducing the
iwrraw file for tbeptools and must be saved in the
data folder for the package.
# wbid classes from wbid, to join
classadd <- wbid %>%
select(wbid = WBID, class= CLASS) %>%
st_set_geometry(NULL) %>%
unique
# assigns stations to creek id basedon buffer, wbid already in dataset
# pull station data with the intersect, add class column
iwrraw <- st_buffer(tdlcrk, dist = 200, endCapStyle = 'FLAT') %>%
st_intersection(stas, .) %>%
select(sta = STATION_ID, wbid = WATERBODY_, JEI = CreekID) %>%
st_set_geometry(NULL) %>%
inner_join(., iwr, by = 'sta', relationship = 'many-to-many') %>%
mutate(
newComment = rCode
) %>%
left_join(classadd, by = c('wbid'))
# save to Desktop for tbeptools data folder
save(iwrraw, file = '~/Desktop/iwrraw.RData', compress = 'xz')
head(iwrraw)## sta wbid JEI year month day masterCode result rCode newComment class
## 1 21FLA 24010600 2052 CC01 1972 10 3 COLOR 15.0000 <NA> <NA> 3M
## 2 21FLA 24010600 2052 CC01 1972 10 3 DO 4.0000 <NA> <NA> 3M
## 3 21FLA 24010600 2052 CC01 1972 10 3 DOSAT 61.3873 $ $ 3M
## 4 21FLA 24010600 2052 CC01 1972 10 3 NO3O2 0.0400 + + 3M
## 5 21FLA 24010600 2052 CC01 1972 10 3 ORGN 0.2800 <NA> <NA> 3M
## 6 21FLA 24010600 2052 CC01 1972 10 3 SALIN 29.0000 <NA> <NA> 3MFinally, we want to compare the results from the previous year’s run with the new results to see how the scores have changed. We use functions from the tbeptools package to estimate creek scores, then compare in a matrix.
library(tbeptools)
# get new estimates
tmpnew <- anlz_tdlcrk(tidalcreeks, iwrraw, yr = 2025)
# get old estimates
tmpold <- anlz_tdlcrk(tidalcreeks, iwrraw, yr = 2024)
# compare
tmpnewcmp <- tmpnew %>%
select(wbid, JEI, class, scorenew = score)
tmpoldcmp <- tmpold %>%
select(wbid, JEI, class, scoreold = score)
levs <- c('No Data', 'Monitor', 'Caution', 'Investigate', 'Prioritize')
cmp <- full_join(tmpnewcmp, tmpoldcmp, by = c('wbid', 'JEI', 'class')) |>
mutate(
scorenew = factor(scorenew, levels = levs),
scoreold = factor(scoreold, levels = levs)
)
table(cmp[, c('scorenew', 'scoreold')])## scoreold
## scorenew No Data Monitor Caution Investigate Prioritize
## No Data 390 1 0 2 0
## Monitor 1 159 1 0 2
## Caution 0 2 21 1 0
## Investigate 0 2 0 17 3
## Prioritize 1 0 0 0 20