Creek_select_tbeptools.Rmd

---
title: "Creek Select for tbeptools"
author: "M. W. Beck"
date:  "`r Sys.Date()`"
output: html_document
---

```{r setup, warning = F, message = F}
knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)

library(sf)
library(RODBC)
library(dplyr)
library(here)
```

```{r, include = F}
# get rdata from github
# try simple load, download if fail
rdataload <- function(flurl){
  
  fl <- basename(flurl)
  obj <- gsub('\\.RData$', '', fl)
  
  # try simple load
  ld <- try(load(url(flurl)), silent = T)
  
  # return x if load worked
  if(!inherits(ld, 'try-error')){
    out <- get(obj)
  }
  
  # download x if load failed
  if(inherits(ld, 'try-error')){
    
    fl <- paste(tempdir(), fl, sep = '/')
    download.file(flurl, destfile = fl, quiet = T)
    load(file = fl)
    out <- get(obj)
    suppressMessages(file.remove(fl))
    
  }
  
  return(out)
  
}
```

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:

* Install the R package dependencies above 
* Download/extract the `.accdb` files for a selected run from here: <https://publicfiles.dep.state.fl.us/dear/iwr/>
* Have access to the source tidal creek line file with creek ID (JEI) assignments created by MW (included in this [repository](https://tbep-tech.github.io/tidalcreek-stats/)).
* A WBID polygon shapefile for the appropriate run, which can be retrieved through web services (<https://geodata.dep.state.fl.us/datasets/FDEP::waterbody-ids-wbids/about>.  The Run66 WBID shapefiles are also available from the folder here: <http://publicfiles.dep.state.fl.us/dear/kristin/> (no need to download if using web services). 

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.

```{r}
# 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.

```{r}
# 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)
  
}

head(iwr)
```

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.

```{r}
# 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)
```

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.

```{r}
# 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)
```

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.

```{r}
# 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)
```

Finally, 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.

```{r}
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')])
```