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 66, 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/iwr2024_run66_2024-11-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))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 <- 2024
# 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)
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 13 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 406216.7 ymin: 2939554 xmax: 408028.2 ymax: 2941086
## Projected CRS: NAD83 / UTM zone 17N
## OBJECTID STATION_ID STATION_NAME WATERBODY_ID
## 1 1 1113S000120310 CAPE CORAL -1500FT FRM CANAL MOU 3240A1
## 2 2 1113S000120315 CAPE CORAL -1500FT FRM CANAL MOU 3240A
## 3 3 1113S000120320 CAPE CORAL -AT CANAL MOUTH 3240A1
## 4 4 1113S000120325 CAPE CORAL -4000FT FRM CANAL MOU 3240A1
## 5 5 1113S000120330 CAPE CORAL -2500FT FRM CANAL MOU 3240A1
## 6 6 1113S000120335 CAPE CORAL -3000FT FRM CANAL MOU 3240A1
## COMMENTS LATITUDE LONGITUDE RESULT_COUNT BEGIN_DATE END_DATE TMDL_RUN
## 1 submitted by PH 26.5875 -81.9319 150 1975 1975 Run 41
## 2 submitted by PH 26.5847 -81.9236 197 1975 1975 Run 41
## 3 submitted by PH 26.5847 -81.9278 208 1975 1975 Run 41
## 4 submitted by PH 26.5736 -81.9417 95 1975 1975 Run 41
## 5 submitted by PH 26.5764 -81.9417 93 1975 1975 Run 41
## 6 submitted by PH 26.5806 -81.9403 115 1975 1975 Run 41
## WATER_QUALITY_STATION BIOLOGY_STATION geometry
## 1 1 0 POINT (407203.9 2941086)
## 2 1 0 POINT (408028.2 2940770)
## 3 1 0 POINT (407610 2940773)
## 4 1 0 POINT (406216.7 2939554)
## 5 1 0 POINT (406219 2939864)
## 6 1 0 POINT (406361.8 2940328)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
## 1 1 Rock Creek CC01 1983B 2 6443.083
## 2 2 Rock Creek CC01 2052 3M 6443.083
## 3 3 Oyster Creek CC02 1983B 2 10175.167
## 4 4 Oyster Creek CC02 2067 3M 10175.167
## 5 5 Buck Creek CC03 1983B 2 2251.188
## 6 6 Buck Creek CC03 2068A 3M 2251.188
## geometry
## 1 LINESTRING (-82.33869 26.92...
## 2 MULTILINESTRING ((-82.33006...
## 3 MULTILINESTRING ((-82.33188...
## 4 MULTILINESTRING ((-82.3194 ...
## 5 LINESTRING (-82.31512 26.89...
## 6 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_ID, 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
## 1 21FLA 24010600 2052 CC01 1972 10 3 COLOR 15.0000 <NA> <NA>
## 2 21FLA 24010600 2052 CC01 1972 10 3 DO 4.0000 <NA> <NA>
## 3 21FLA 24010600 2052 CC01 1972 10 3 DOSAT 61.3873 $ $
## 4 21FLA 24010600 2052 CC01 1972 10 3 NO3O2 0.0400 + +
## 5 21FLA 24010600 2052 CC01 1972 10 3 ORGN 0.2800 <NA> <NA>
## 6 21FLA 24010600 2052 CC01 1972 10 3 SALIN 29.0000 <NA> <NA>
## class
## 1 3M
## 2 3M
## 3 3M
## 4 3M
## 5 3M
## 6 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 = 2024)
# get old estimates
tmpold <- anlz_tdlcrk(tidalcreeks, iwrraw, yr = 2023)
# 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 392 2 0 0 0
## Monitor 1 163 0 1 2
## Caution 0 1 16 0 0
## Investigate 0 0 1 18 1
## Prioritize 0 0 0 1 21