knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)
library(sf)
library(RODBC)
library(dplyr)
library(here)This workflow was developed as a test case for updating data used to estimate tidal creek assessment cateogires using FLDEP IWR run data. In this example, we use IWR Run 65, 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_run65_2024-01-30/"
# # 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 65 (url should not change betweem 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. It only takes a few minutes to run and does not eat up 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 <- 2023
# 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 run65
# 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: 182215.4 ymin: 3346961 xmax: 186005.4 ymax: 3349716
## Projected CRS: NAD83 / UTM zone 17N
## OBJECTID STATION_ID STATION_NAME WATERBODY_ID
## 1 1 21FLFSI WAK-5 WAK-5 1006
## 2 2 21FLFSI WAK-6 WAK-6 1006
## 3 3 21FLFSI WAK-7 WAK-7 1006
## 4 4 21FLFSI WAK-8 WAK-8 1006
## 5 5 21FLGW 22848 NW1-LR-2034 WAKULLA RIVER 1006
## 6 6 21FLGW 34879 WAKULLA RIVER 1006
## COMMENTS
## 1 New station assigned to 1006, Run 61 - P.Homann
## 2 New station assigned to 1006, Run 58 - R.Stuart
## 3 New station assigned to 1006, Run 58 - R.Stuart
## 4 New station assigned to 1006, Run 58 - R.Stuart
## 5 Joe Hand Run 22 9/13/2005
## 6 J Hand March 2009; Lat/Longs updated based on values in WIN, Run 58 - R.Stuart
## LATITUDE LONGITUDE RESULT_COUNT BEGIN_DATE END_DATE TMDL_RUN
## 1 30.23421 -84.29015 24 2019 2020 Run 61
## 2 30.23756 -84.29922 28 2019 2020 Run 58
## 3 30.23583 -84.30182 33 2019 2020 Run 58
## 4 30.23664 -84.30181 27 2019 2020 Run 58
## 5 30.22896 -84.28298 20 2004 2004 <NA>
## 6 30.21365 -84.26173 6591 2008 2022 <NA>
## WATER_QUALITY_STATION BIOLOGY_STATION geometry
## 1 1 0 POINT (183334 3349320)
## 2 1 0 POINT (182471.4 3349716)
## 3 1 0 POINT (182215.4 3349532)
## 4 1 0 POINT (182219 3349622)
## 5 1 1 POINT (184007.7 3348717)
## 6 1 1 POINT (186005.4 3346961)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
wBID.
# 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)
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 2068 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 uisng an inner join by station.
Finally, the WBID class is added. This is reproducing the
iwrraw file from the tbeptools 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'))
head(iwrraw)## sta wbid JEI year month day masterCode result rCode newComment
## 1 21FLSWFD764483 1983B CC01 2009 12 1 COND 52000.0 <NA> <NA>
## 2 21FLSWFD764483 1983B CC01 2009 12 1 DO 8.2 <NA> <NA>
## 3 21FLSWFD764483 1983B CC01 2009 12 1 DOSAT 114.7 + +
## 4 21FLSWFD764483 1983B CC01 2009 12 1 SALIN 34.3 <NA> <NA>
## 5 21FLSWFD764483 1983B CC01 2009 12 1 COND 52200.0 <NA> <NA>
## 6 21FLSWFD764483 1983B CC01 2009 12 1 DO 8.3 <NA> <NA>
## class
## 1 2
## 2 2
## 3 2
## 4 2
## 5 2
## 6 2Finally, 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 = 2023)
# get old estimates
tmpold <- anlz_tdlcrk(tidalcreeks, iwrraw, yr = 2022)
# 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 400 5 0 1 0
## Monitor 0 149 1 2 0
## Caution 0 0 16 1 0
## Investigate 0 0 2 14 0
## Prioritize 0 3 0 1 22