04_StreamCat_Covariates_and_SSN_Object.Rmd

---
title:  "James EDA StreamCat Covariates"
author: "Michael McManus & Travis Linscome-Hatfield"
date: "12/03/2024"
output:
  html_document:
    fig_caption: yes
    number_sections: false
    toc: yes
    toc_float:  yes
    code_folding: hide
    self_contained: yes
    theme: lumen
editor_options: 
  chunk_output_type: console
  markdown: 
    wrap: sentence
---

## Outline

We now turn our attention to the covariates, or predictor variables, we want to use to model and predict VSCI.
The SSN object we are using was built using the STARS Toolbox in ArcGIS.
Now the building of an SSN object can be done in R using <https://cran.r-project.org/package=SSNbler>.
This script shows how additional covariates from EPA's Stream Catchment dataset, or StreamCat, can be incorporated in the SSN object.

Log Keeping log below for now when reach final version will delete.
Text above the log I plan to keep to provide background for users.

11/19/2024 have code running end to end.

On 10/15/2024 added waterbody data for Travis to evaluate on github.

On 09/19/2024 Ellen provided update waterbody datasets as csv files.

09/11/2024 Ellen and I worked out how to get NHD waterbody associated with 199 monitoring observations.

08/28/2024 explored STATSGO_Set1 and STATSGO_Set2 variables from StreamCat.
Also, brought in waterbodies data created by Ellen D'Amico from Pegasus.
Ellen used comid from edges, but I am not sure they all match to the featureid of the points.

On 08/15/2024 added data, figures_sfs, and ssn_object folders and trying to run code pointing to those folders.

On 08/14/2024 pushing to github repository.

On 07/31/2024 bringing in WWTP_VA.vcs to do a left merge with DFobs so I can compare stations with WWTPs to stations without WWTPs in the VSCI response.

On 07/11/2024 bring the new ssn into R.
On 07/10/2024 downloaded James_071024.ssn that had PRISM climate data, including normals and deviations from normals.

On 05/14/2024 running code to get SSN output for SSN models for the 4 geographies and because 05/08/2024 updated R, RTools, and RStudio.

I am going to bring in two SSNs, the 040424 version and the latest 041924 version.
I will take the PRISM climate variables out of 041924 and put into 040424.
Downloaded James_041924.ssn_NEWEST.zip and it is stored in James_041924.ssn with file modified on 04/22/2024 at 7:34am.
Ellen removed netid and pid columns and reran through SSN compilation.
Bestglm still gave error

On 04/19/2024 imported James041924.ssn, still could not get bestglm to run.

On 04/17/2024 imported James041724.ssn that now include annual maximum temperature, annual average temperature, and annual average precipitation.
Also, I will simplify coding of wetlands to be absence/presence across all 4 geographies: watershed, watershed-riparian, catchment, and catchment-riparian.
For grass cover at catchment and catchment-riparian I will also use binary code for absence/presence.

On 04/16/2024 imported James041924.ssn, which now has PRISM climate data from StreamCat.

This code will only use vsci as the response variable.
Also, I will examine if impervious surface, forest, hay, and grass at 0 values show dramatically different vsci y-intercept, when x = 0, compared to vsci y-intercepts when x \> 0.

On 04/09/2024, I had a call with Mike Dumelle, and he suggested 1) use untransformed response variable, 2) include binary absence/presence along with continuous landscape percentages, and 3) try including year.

On 04/04/2024 now using James040424.ssn as Ellen D'Amico from Pegasus updated Preds_2021_2022 points.
I entered WQ data as described in section 1.10 Missing Water Chemistry Data below.

On 01/09/2023 read in Pegasus' SSN object from 01/09/2024, which has additional StreamCat variables.
I am interested in the shrub and wetland metrics.
In previous SSN, Note that VDEQ switched fonts between Wadeable_ProbMon_2001-2018_Final_Final.xlsx and Wadeable_ProbMon_2001-2020.xlsx.
I used janitor package to put all variable names in lowercase.

On 12/26/2023 read in Pegasus' SSN object from 11/16/2023, which as fields of COMID, AreaSqkM, TotDASqKM, h2oAreaKm2, and rcaAreaKM2.
Using this new SSN object means I re-entered missing water chemistry data for two stations in ArcGIS.

Now using SSN2 from CRAN on 11/16/2023.

Note on SSNs I receive from Pegasus, I have to hand enter some of the water chemistry values for stations 2-WLS023.10 and 2AMLC000.84.
This is because Emma Jones at VDEQ sent FinalFieldAnalytes_MissingJamesSites.csv by email on 09/29/2023.
That csv file has water chemistry entries for those two stations, which in the original data file were NAs.
I edit those entries in ArcGIS before I do import into R.
Once the SSN is in R I can't export back to ArcGIS.

On 11/09/2023 unzipped.
James.ssn_110923 and saved to James.ssn.
This James.ssn_11092023 has additonal covariates produced by Ellen D'Amico, GIS analyst with Pegasus.
The James.ssn_102323 had stream order added to edges.
On 11/08/2023 unzipped James.ssn_102323 into ssn_objects folder and renamed it as James.ssn so I would not need to change R code.
Can I use SSN2 package for some of ESDA of James River?
Yes.

# Load Libraries

```{r setup, collapse=TRUE}

library(ggplot2)
library(gridExtra)
library(car)
library(bestglm)
library(sf)
library(dplyr)
library(gpairs)
library(moments)
library(fitdistrplus)
library(openxlsx)
library(kableExtra)
library(scales)
library(gstat) # for semivariogram cloud
library(lattice) # for random 
library(tmap)
library(mapview)
library(leaflet)
library(leafpop) # for popups in mapview
library(leafsync) # to sync obs and prediction maps
library(units)
library(GGally)
library(readxl)
library(tidyr)
library(spmodel)
library(purrr)
library(SSN2)
library(dummy)
library(performance)
library(see)
library(effects)
library(ggeffects)
library(janitor)
library(nngeo) # nearest neighbor distances
library(plotly)

knitr::opts_chunk$set(message=FALSE, warning=FALSE,collapse = T)

```

# 1.0 SSN and Initial Covariates

Bring in the SSN object.
Pulling out DFobs gives an sf and data frame object in which variables can be added, edited, modified, transformed, etc. to the obs and preds and then put back into the SSN object.

Also, get some summary stats on the VDEQ stations.
```{r load data}
j_ssn1a <- SSN2::ssn_import("ssn_object/James_071024_pluspreds.ssn", predpts = "sites")


names(j_ssn1a)
summary(j_ssn1a)

# pull set of observation data out from ssn and clean names
DFobs <- SSN2::ssn_get_data(j_ssn1a) %>% clean_names(.)

# now Torgegram function gives all 3 distances & plots them
# perhaps I had typo somewhere
# ztg1 <- SSN2::Torgegram(VSCIVCPMI ~ 1, j_ssn1a, type = c("flowcon", "flowuncon", "euclid"))
# 
# names(ztg1)
# plot(ztg1, main = "VSCI")
# torg <- ztg1[["flowcon"]] # returns flowcon & see small np
# torg
# names(torg)
# class(torg)
# ggplot(torg, aes(x=dist, y=gamma,size=np)) + geom_point()
# 



```

Steps are 1) select station_id from DFobs with geometry, 3) use unique.data.frame to get unique stations, and then 4) use st_nn
```{r dist_knn1}
# See binding distances to join result from
# https://michaeldorman.github.io/nngeo/articles/intro.html

z <- dplyr::select(DFobs, station_id) #199 obs
head(z)
z1 <- unique.data.frame(z) #181 unique stations
z1_nn <- st_nn(z1, z1, k =2, returnDist = TRUE, progress = False)
str(z1_nn)

head(z1_nn[[1]]) # list station obs & nn obs
head(z1_nn[[2]]) # list distance to self & knn1

station_nn <- as.data.frame(z1_nn$nn)
# go to 2nd of two list and grab first entry by row, which is 0
dists = sapply(z1_nn[[2]], "[", 1)
# go to 2nd of two list and grab second entry by row, which is knn1 distance
dists_knn1 = sapply(z1_nn[[2]], "[", 2)
head(dists_knn1)
summary(dists_knn1)
# Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# 0.161  2014.809  4765.854  5034.192  7320.395 16006.036 

# based on VDEQ stream order (SO) entries, compare to SSN sites and edges
vahusb <- as.data.frame(DFobs) %>%
  group_by(vahusb, order) %>%
  summarize(
            count = dplyr::n()) %>% 
  print(n = Inf)

# remove datasets not used downstream
rm(z, z1, z1, z1_nn, station_nn, dists, dists_knn1,vahusb)
```


# 2.0 Add New StreamCat Covariate

## 2.00 Explore WWTP Data

Having the COMID/FeatureID for each of the obs site and prediction points is what allows us to work with StreamCat Covariates.
StreamCat covariates have a COMID/FeatureID associated with them, with those identifiers coming from the reaches in NHDPlusV2.

A typical StreamCat workflow is to bring a csv file for a specific variable for all of the catchments in the state.
The left join of the DFobs from the SSN to the csv using `r by = join_by(FeatureID == COMID))`

returns only the observations satisfying that match.
Evaluating if wastewater treatment plant (wwtp) data from StreamCat should be considered as a predictor.

```{r wwtp}
wwtp_ds1 <- read.csv("data/WWTP_VA.csv") %>% 
  clean_names(.)

DFobs <- left_join(DFobs, wwtp_ds1, by = join_by(feature_id == comid))

DFobs <- DFobs %>% 
  mutate(
    wwtp = case_when(
      wwtp_all_dens_ws > 0 ~ "yes",
      wwtp_all_dens_ws == 0 ~ "no"
    )
  )

DFobs$wwtp <- factor(DFobs$wwtp, levels = c("yes","no"))
summary(DFobs$wwtp)

# remove file if not needed downstream
rm(wwtp_ds1)

# initially created zDFobs for EDA
# zDFobs <- left_join(DFobs, wwtp_ds1, by = join_by(FeatureID == COMID))
# 
# names(zDFobs)
# 
# zDFobs <- zDFobs %>% 
#   mutate(
#     wwtp = case_when(
#       WWTPAllDensWs > 0 ~ "yes",
#       WWTPAllDensWs == 0 ~ "no"
#     )
#   )
# 
# names(zDFobs)
# 
# summary(zDFobs$WWTPAllDensCat)
# summary(zDFobs$wwtp)
# # in VIEW filter wwtp = yes returns 68 entries
# # in VIEW filter wwtp = no returns 131 entries
# 
# ggplot(zDFobs, aes(x = as.factor(wwtp), y = log(SpCond))) + geom_boxplot()
# 
# ggplot(zDFobs, aes(x = as.factor(wwtp), y = log(TP))) + geom_boxplot()
# 
# ggplot(zDFobs, aes(x = as.factor(wwtp), y = log(Turb))) + geom_boxplot()
# 
# ggplot(zDFobs, aes(x = as.factor(wwtp), y = pH)) + geom_boxplot()
# 
# ggplot(zDFobs, aes(x = as.factor(wwtp), y = PctImp_W)) + geom_boxplot()
# 
# zDFobs %>%
#   group_by(wwtp) %>%
#   summarize(min = fivenum(PctImp_W)[1],
#             low = fivenum(PctImp_W)[2],
#             median = fivenum(PctImp_W)[3],
#             mean = mean(PctImp_W, na.rm = TRUE),
#             upper = fivenum(PctImp_W)[4],
#             max = fivenum(PctImp_W)[5],
#             count = dplyr::n()) %>% 
#   print(n = Inf)
# 
# ggplot(zDFobs, aes(x = as.factor(wwtp), y = VSCIVCPMI)) + geom_boxplot()
# 
# # Trend station 2-TYE008.77 has wwtp upstream.
# trend1 <- dplyr::filter(zDFobs, (StationID != StationID_) & wwtp == "yes")
# 
# impwwtp <- ggplot(zDFobs, aes( x = PctImp_W, y = WWTPAllDensWs)) + geom_point()
# ggplotly(impwwtp)
# 
# wwtp_3highest <- dplyr::filter(zDFobs, WWTPAllDensWs > 0.05)
# 
# ggplot(zDFobs, aes( x = PctImp_W, y = VSCIVCPMI, color = as.factor(wwtp))) + geom_point()
# 
# rm(zDFobs)

```

## 2.01 Explore STATSGO_Set1 Variables

Based on Schmid et al. 2011 and Hill et al. 2017, consider SandWs as a possible covariate.
Four observations at \~ 55% or greater really stand out and are all from JL.
Median SandWs and IQR for JL much greater than other for subbasins.
JU shows positive slope to VSCI, 3 other subbasins show negative slope, and JL is flat.
Based on EDA decided not include ws_sand as covariate.

```{r statsgo_set1}
# statsgo_set1 <- read.csv("data/STATSGO_Set1_VA.csv") %>% clean_names(.)
# 
# zDFobs <- left_join(DFobs, statsgo_set1, by = join_by(feature_id == comid))
# 
# names(zDFobs)
# summary(zDFobs$sand_ws)
# 
# ggplot(zDFobs, aes(x = sand_ws, y = log(sp_cond))) + geom_point()
# 
# ggplot(zDFobs, aes(x = sand_ws, y = log(tp))) + geom_point()
# 
# ggplot(zDFobs, aes(x = sand_ws, y = log(turb))) + geom_point()
# 
# ggplot(zDFobs, aes(x = sand_ws, y = p_h)) + geom_point()
# 
# ggplot(zDFobs, aes(x = sand_ws, y = PctImp_W)) + geom_point()
# 
# mapview(zDFobs, zcol = "sand_ws", cex = "sand_ws", alpha.regions = .8, legend = TRUE, layer.name = 'Mean Percent Sand Ws', popup = popupTable(zDFobs, zcol = c("sand_ws", "station_id", "vahusb", "year")))
# 
# ggplot(zDFobs, aes(x = as.factor(vahusb), y = sand_ws)) + geom_boxplot()
# 
# filter(zDFobs, vahusb == "JL") %>% select(sand_ws)
# 
# ggplot(zDFobs, aes( x = sand_ws, y = vscivcpmi, color = as.factor(vahusb))) + geom_point() + geom_smooth(method = lm)
# 
# ggplot(zDFobs, aes( x = sand_ws, y = vscivcpmi)) + geom_point() + geom_smooth(method = lm)
# 
# rm(zDFobs)
```

## 2.03 Explore STATSGO_Set2 Variables
Mean seasonal water table depth (cm) of soils within within a wateshed is abbreviated as wt_dep_ws.
wt_dep_ws values decrease going from west to east.
JU, JM, and JR all have about same median, and much lower medians in JR and JL.
Mean organic matter content (% by weight) of soils within a watershed is abbreviated as om_ws.
om_ws shows west to east gradient on map, but lot of overlap among boxplots by vAHUSB.
Mean permeability (cm/hour) of soils within a watershed is abbreviated as perm_ws.
perm_ws spatially variable as highest medians in JU and JL and lower medians in JM, JR, and JA.
Mean depth (cm) to bedrock of soils within watershed is abbreviated as rck_dep_ws.
rck_dep_ws increases as go from west to east.
rck_dep_ws has the largest correlation, -0.5, with VSCICPMI so likely only use that as covariate from STATSGO data.
JL and JA have largest rck_dep_ws but with very little variation.
rck_dep_ws is highly, negatively correlated with elev_ws.

```{r statsgo_set2}
statsgo_set2 <- read.csv("E:/R_vdeq_nhdplus/STATSGO_Set2_VA.csv") %>% clean_names(.)

names(statsgo_set2)

DFobs <- left_join(DFobs, statsgo_set2, by = join_by(feature_id == comid))
names(DFobs)

rm(statsgo_set2)
# EDA code commented out
# mapview(zzDFobs, zcol = "RckDepWs", cex = "RckDepWs", alpha.regions = .8, legend = TRUE, layer.name = 'Mean Depth to Bedrock (cm)', popup = popupTable(zzDFobs, zcol = c("RckDepWs", "StationID_", "VAHUSB", "Year_")))
# 
# ggplot(zzDFobs, aes(x = as.factor(VAHUSB), y = RckDepWs)) + geom_boxplot()
# 
# 
# as.data.frame(dplyr::select(zzDFobs,c(VSCIVCPMI,SandWs,WtDepWs,OmWs, PermWs, RckDepWs))) %>%
# ggcorr(., method=c("pairwise", "spearman"), label = TRUE)
# 
# st_set_geometry(zzDFobs, NULL) %>% dplyr::select(.,c(VSCIVCPMI,SandWs,WtDepWs,OmWs, PermWs, RckDepWs,ElevWs)) %>% ggpairs(.)
# 
# ggplot(zzDFobs, aes( x = RckDepWs, y = VSCIVCPMI, color = as.factor(VAHUSB))) + geom_point() + geom_smooth(method = lm)
# 
# ggplot(zzDFobs, aes( x = RckDepWs, y = VSCIVCPMI)) + geom_point() + geom_smooth(method = lm)
# 
# filter(zzDFobs, VAHUSB == "JU") %>% select(c(SandWs, RckDepWs,ElevWs)) %>% arrange(ElevWs) %>% print(n =Inf)
# 
# zzDFobs %>%
#   group_by(VAHUSB) %>%
#   summarize(min = fivenum(RckDepWs)[1],
#             low = fivenum(RckDepWs)[2],
#             median = fivenum(RckDepWs)[3],
#             mean = mean(RckDepWs, na.rm = TRUE),
#             upper = fivenum(RckDepWs)[4],
#             max = fivenum(RckDepWs)[5],
#             count = dplyr::n()) %>% 
#   print(n = Inf)

```

## 2.04 Explore Waterbodies Data

Ellen D'Amico calculated waterbody variables.
Have Travis evaluate coding to evaluate an interaction between presence/absence of waterbody and distance (km) to nearest upstream waterbody.

```{r waterbodies}
wb_ds1 <- read.csv("data/ObservationPoints_DistancesUpstream_100424.csv") %>% clean_names(.)
class(wb_ds1)
names(wb_ds1)
head(wb_ds1)
str(wb_ds1)
wb_metadata <- read.csv("data/waterbody_variables.csv") %>% clean_names(.)

summary(wb_ds1$num_waterbody_up, na.rm = TRUE)
summary(wb_ds1$dist_nearest_up_wb_km, na.rm = TRUE)
wb_ds1$dist_nearest_up_wb_km <- round(wb_ds1$dist_nearest_up_wb_km, 1)

summary(wb_ds1$dist_nearest_up_wb_km, na.rm = TRUE)
# file below provides a link to waterbody, and subsequently to DFobs data.
StationIDs_UniqLocIDs <- read.csv("data/StationIDs_UniqLocIDs.csv") %>% clean_names(.)

names(StationIDs_UniqLocIDs)

wb_ds2 <- left_join(StationIDs_UniqLocIDs, wb_ds1, by=join_by(uniq_loc_id))
# produces wb_ds2 with 199 obs

DFobs <- left_join(DFobs, wb_ds2, by = join_by(station_id_2==station_id))
names(DFobs)

DFobs <- DFobs %>% 
  mutate(
    wbc = case_when(
      incl_nhdwb == "yes" ~ "present",
      incl_nhdwb != "yes" ~ "absent"
    )
  )

class(DFobs$wbc)

DFobs$wbc <- factor(DFobs$wbc, levels = c("present", "absent"))
summary(DFobs$wbc)
names(DFobs)

rm(wb_ds1, wb_ds2, wb_metadata, StationIDs_UniqLocIDs)

# EDA code below commented out
# DFobs_wb <- left_join(DFobs, wb_ds2, by = join_by(StationID_))
# 
# DFobs_wb2 <- DFobs_wb %>% 
#   mutate(
#     wbc = case_when(
#       inclNHDWB == "yes" ~ "present",
#       inclNHDWB != "yes" ~ "absent"
#     )
#   )
# 
# class(DFobs_wb2$wbc)
# 
# DFobs_wb2$wbc <- factor(DFobs_wb2$wbc, levels = c("present", "absent"))
# summary(DFobs_wb2$wbc)
# 
# DFobs_wb2$DistNearestUpWb_km[is.na(DFobs_wb2$DistNearestUpWb_km)] <- 0
# 
# DFobs_wb2$wbc_nrstupwb <- (DFobs_wb2$wbc*DFobs_wb2$DistNearestUpWb_km)
# 
# DFobs_wb2$wbc_nrstupwb == DFobs_wb2$DistNearestUpWb_km
# # Does the equality above mean I could skip having to test main effect A (wbc), main effect B (DistNearestUpWb_km) and the A x B interaction (wbc_nrstupwb)? Could I just test a single covariate of DistNearestUpWb_km where 0 is absence of upstream waterbodies and values > 0 are km to nearest upstream waterbody?
# 
# as.data.frame(DFobs_wb2) %>%
#   group_by(wbc) %>%
#   summarize(min = fivenum(VSCIVCPMI)[1],
#             q1 = fivenum(VSCIVCPMI)[2],
#             median = fivenum(VSCIVCPMI)[3],
#             mean = mean(VSCIVCPMI, na.rm = TRUE),
#             q3 = fivenum(VSCIVCPMI)[4],
#             max = fivenum(VSCIVCPMI)[5],
#             count = dplyr::n()) %>% 
#   print(n = Inf)
# 
# ggplot(DFobs_wb2, aes(x = wbc, y = VSCIVCPMI)) + geom_boxplot()
# 
# ggplot(DFobs_wb2, aes(x = DistNearestUpWb_km, y = VSCIVCPMI)) + geom_point() + geom_smooth(method = "lm", se = TRUE)
# 
# ggplot(DFobs_wb2, aes(x = DistNearestUpWb_km, y = VSCIVCPMI)) + geom_point() + facet_wrap(vars(wbc))  + geom_smooth(method = "lm", se =TRUE)
# 
# summary(DFobs_wb$NumWaterbodyUp, na.rm = TRUE)
# summary(DFobs_wb$DistNearestUpWb_km, na.rm = TRUE)
# 
# ggplot(DFobs_wb |> filter(inclNHDWB == "yes"), aes(x = NumWaterbodyUp, y = VSCIVCPMI)) + geom_point() 
# 
# ggplot(DFobs_wb |> filter(inclNHDWB == "yes"), aes(x = DistNearestUpWb_km, y = VSCIVCPMI)) + geom_point() + facet_wrap(vars(as.factor(VAHUSB)))  + geom_smooth(method = "lm", se =TRUE)
# 
# ggplot(DFobs_wb |> filter(inclNHDWB == "yes"), aes(x = wbArea_up_SqKM, y = VSCIVCPMI)) + geom_point()  + geom_smooth(method = "lm", se =TRUE)
# 
# ggplot(DFobs_wb |> filter(inclNHDWB == "yes"), aes(x = clanl, y = VSCIVCPMI)) + geom_point()  + geom_smooth(method = "lm", se =TRUE)
# 
# filter(DFobs_wb, inclNHDWB == "yes") %>% 
#   select(., StationID_,VSCIVCPMI, inclNHDWB, DistNearestUpWb_km, NumWaterbodyUp, wbArea_up_SqKM, clanl) %>% arrange(VSCIVCPMI) %>% 
#   print(n = Inf)
# 
# filter(DFobs_wb, inclNHDWB == "yes") %>% 
#   select(., StationID_,VSCIVCPMI, inclNHDWB, DistNearestUpWb_km, NumWaterbodyUp, wbArea_up_SqKM, clanl) %>% arrange(StationID_) %>% 
#   print(n = Inf)

# wrote shapefile to explore in GIS
# st_write(DFobs_wb, "outputs/DFobs_wb_100424.shp")
```

# 3.0 Entering, Modifying, and Transforming Covariates

## 3.10 Missing Water Chemistry Data

I emailed Emma Jones, VDEQ, on September 28, 2023 regarding two stations, 2AMLC000.84 and 2-WLS023.10, in the James in 2011 that had missing water chemistry.
She replied the next day saying the NA entries may have occurred because the internal data querying system, called Logi,may not have pulled the data for those sites.
Emma attached a spreadsheet, FinalFieldAnalyteMetals_MissingJamesSites.csv, that contained numeric entries for those two stations.
I entered those values in ArcMap so that the data would come in with the James SSN object.
Also, that way when I look at the James SSN object in ArcMap all the fields have the correct entries and not zeroes.
That csv file is pulled in below so it can be checked against the SSN object by comparing the values for those stations in the csv file to those in DFobs.

```{r missing_waterchem}
# nachem_ds1 <- read.csv("data/FinalFieldAnalyteMetals_MissingJamesSites.csv")

```

## 3.20 Total Habitat (RBP) score

Bring Total Habitat Score (TotHab) in from Wadeable_ProbMon_2001-2018_Final_Final.xslx spreadsheet so it can be joined to DFobs and then j_ssn1a.
These 2 stations: 2-JKS070.97 and 2-DDY000.75_2017 do not have tothab as Emma confirmed in her 11/24/2023 email.
Both sites are in Central Appalachian Ridges and Valleys.
Both have high VSCI scores of 73.8 and 84.5 (the latter is max VSCI), respectively.

I imputed their tothab scores by taking averages.
For trend station 2-DDY000.75_2017, I could average the scores of 172.5, 173.5, and 178 from 2011, 2013, and 2015, respectively.
For 2-JKS070.97 on 3rd order, could I average nearby 2-JKS0 sites?
2-JKS076.16, has tothab of 162.0 on 3rd order, is upstream of 2-JKS070.97, about 16 km apart.
Perhaps also use Back Creek site, 2-BCC001.90 (has tothab of 189.0 on 2nd order), that flows parallel to Jackson River, where the 2 sites are, and Back Creek site is near confluence to Jackson River.
2-BCC001.90 is about 9 stream km from 2-JKS070.97.
The new variable tothab is created.

```{r tothab}

tothab_ds1 <- read_xlsx("data/Wadeable_ProbMon_2001-2018_Final_Final.xlsx", range = "Wadeable_ProbMon_2001-2018!D1:BK814")
                        
tothab_ds2 <- tothab_ds1 |>
  filter(SubBasin == "James") |>
  dplyr::select(StationID_Trend, TotHab)

tothab_ds2 <- tothab_ds2 %>% mutate_at(c('TotHab'), as.numeric)

summary(tothab_ds2$TotHab)

# 2-DDY000.75_2017 is on Daddy Run headwater of Calfpasture River
# https://stackoverflow.com/questions/32829358/dplyr-filter-with-sql-like-wildcard
dr_na <- filter(tothab_ds2, grepl("2-DDY000.75", StationID_Trend, fixed = TRUE))
summary(dr_na$TotHab)

tothab_ds3 <- tothab_ds2|>
  mutate(
    TotHab = case_when(StationID_Trend == "2-DDY000.75_2017" ~ 176.4,
       TRUE ~ TotHab))


# 2-JKS070.97 is on Jackson River
jr_na <- filter(tothab_ds3, StationID_Trend == "2-JKS076.16"| StationID_Trend == "2-BCC001.90")

summary(jr_na$TotHab)

tothab_ds4 <- tothab_ds3|>
  mutate(
    TotHab = case_when(StationID_Trend == "2-JKS070.97" ~ 175.5,
       TRUE ~ TotHab))

summary(tothab_ds4$TotHab)

names(tothab_ds4)

# ggplot(tothab_ds4, aes(TotHab)) + geom_histogram()
# 
# ggplot(tothab_ds4, aes(TotHab)) + geom_boxplot()
# ggplot(tothab_ds4, aes(x=TotHab)) + geom_density()

tothab_ds4 <- rename(tothab_ds4,c(tothab = TotHab, st_id_tren = StationID_Trend ))
head(tothab_ds4)


DFobs <- full_join(DFobs, tothab_ds4, by = join_by(station_id_2==st_id_tren))

# remove datasets not needed downstream
rm(tothab_ds1, tothab_ds2, dr_na, tothab_ds3, jr_na, tothab_ds4)
```

## 3.25 PRISM data
Get better metadata description of the variables.
Correlations of PRISM climate data.

```{r prism}
# p_df1 <- SSN2::ssn_get_data(prism)
# p_df1 <- dplyr::select(p_df1, c(StationID_, precip_mm, Tmean, Tmax)) %>% 
#   dplyr::rename(., st_id_tren = StationID_) %>% 
#   clean_names(.)
# head(p_df1)
# p_dfz1 <- st_set_geometry(p_df1, NULL)
# DFobs <- left_join(DFobs, p_dfz1, by = join_by(st_id_tren))


as.data.frame(dplyr::select(DFobs,c(vscivcpmi,precip_mm,tmean,tmax, d_precip_m, d_t_mean, d_t_max))) %>%
ggcorr(., method=c("pairwise", "spearman"), label = TRUE)
# vsci & precip_mm correlation -0.1
# vsci & tmean/tmax correlations -0.5
```

## 3.30 Factors: Ecoregion, Bioregion, VAHUSB, Year, Status and Trend Stations

Based on parallel coordinate plots I made in 06_SSN Model_Predict.Rmd I made ju as yes/no factor to use in SSN partition factor argument.
Emma Jones on 02/05/2024 sent shapefiles for VAHUSB polygons, Virginia hydrologic unit subbasins.
These 5 polygons from West to East are: JU James River Upper, (Mountain) n = 71; JM James River Middle, (Piedmont) n = 79; JR James River - Rivanna River n = 12; JA James River - Appomattox River n = 29; and JL James River Lower, (Tidal) n = 8.
I created a station type (st_type) factor to indicate a status or trend station.
Trend stations have repeated visits.
Jason Hill's email of 01/03/2024 mentioned station 2-JKS028.69_2004 began as trend, wadeable station, but was difficult to sample.
It is now a trend, boatable station.
Because I am only working with wadeable stations I assigned the one observation of 2-JKS028.69_2004 as a status station.
The factors station type and year (year_f) are considered random effects in the analysis.

```{r ecoregions_vahusb_trend}
# DFobs$HUC8 <- as.factor(DFobs$HUC8)
DFobs$station_id <- as.factor(DFobs$station_id)
DFobs$eco_region <- as.factor(DFobs$eco_region)
DFobs$bio_region <- factor(DFobs$bio_region, levels = c("Mountain", "Piedmont", "Coast"))
DFobs$year_f <- as.factor(DFobs$year)
DFobs$vahusb <- factor(DFobs$vahusb, levels = c("JU", "JM", "JR", "JA", "JL"))
summary(DFobs$vahusb)

DFobs <- DFobs %>% 
  mutate(
    ju = case_when(
      vahusb == "JU" ~ "yes",
      .default = "no"
    )
  )
DFobs$ju <- factor(DFobs$ju, levels = c("yes","no"))
summary(DFobs$ju)

# ggplot(DFobs, aes(x=vahusb, y = vscivcpmi)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x = tothab, y = vscivcpmi, colour = vahusb)) + geom_point() + geom_smooth(method = "lm")
# 
# ggplot(DFobs, aes(x = tothab, y = vscivcpmi)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(vars(vahusb))

# use code to evaluate vsci by year for trend sites if a random effect for trend stations is needed
DFobs <- DFobs %>%
  mutate(
    st_type = case_when(
      station_id_2 == "2-JKS028.69_2004" ~ "status",
      station_id_2 != "2-JKS028.69_2004" & station_id_2 != station_id ~ "trend",
      station_id_2 != "2-JKS028.69_2004" & station_id_2 == station_id ~ "status",
    )
  )
summary(as.factor(DFobs$st_type))

trend <- dplyr::filter(DFobs, st_type == "trend")
# returned 38 obs
# trend stations suggestive of including random effect
# ggplot(trend, aes(x=year_f, y = vscivcpmi)) + geom_point() + facet_wrap(vars(vahusb_nam, station_id))

# ggplot(DFobs, aes(x=st_type, y = vscivcpmi)) + geom_boxplot()

# ggplot(DFobs, aes(x=st_type, y = tothab)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x=st_type, y = do)) + geom_boxplot()

# remove data frames not needed downstream
rm(trend)
```

## 3.3.1 SFS Trend Stations

```{r sfs_trend_stations}
# trend_stations <- ggplot(trend, aes(x=year_f, y = vscivcpmi)) + geom_point() + facet_wrap(vars(vahusb_nam, station_id), nrow = 2, ncol=4) +labs(x = "Survey Year", y = "VSCI")
# 
# png(file="figures_sfs/trend_stations.png",width=12,height=3,units="in",res=150)
# trend_stations
# dev.off()

```

## 3.3.2 Factors: Wetlands, Impervious Surface, Grass & TP

On January 26, 2024, Jason Hill at VDEQ sent an email describing total phosphorus (tp) concentrations that impacted VSCI based on a relative risk analysis.
At 0.02 mg/L tp and less see little impact, at 0.05 mg/L and above start seeing some impacts, and at 0.1 mg/L and above more likely to have impacts.
I believe this is from a statewide relative risk analysis.
I will use these levels, instead of the ones I had made based on quartertiles.

```{r wetlands_impervious_grass_tp}

# Watershed
summary(DFobs$pct_wet_w)

DFobs <- DFobs %>%
  mutate(
    pct_wet_w_f = case_when(
      pct_wet_w == 0 ~ "none",
      pct_wet_w > 0 & pct_wet_w <= 0.040 ~ "low",
      pct_wet_w > 0.040 & pct_wet_w <= 1.885 ~ "medium",
      pct_wet_w > 1.885 & pct_wet_w <= 16.820 ~ "high",
    )
  )

DFobs$pct_wet_w_f <- factor(DFobs$pct_wet_w_f, levels = c("none", "low", "medium", "high"))

summary(as.factor(DFobs$pct_wet_w_f))

## binary absence/presence of watershed wetlands
DFobs <- DFobs %>%
  mutate(
    bin_wet_w = case_when(
      pct_wet_w == 0 ~ 0,
      pct_wet_w > 0 ~ 1
    )
  )
# DFobs$bin_wet_w <- factor(DFobs$bin_wet_w, levels = c("absence", "presence"))
summary(as.factor(DFobs$bin_wet_w))

# ggplot(DFobs, aes(x=bin_wet_w, y = vscivcpmi)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x=pct_wet_w_f, y = tothab)) + geom_boxplot()
# 
# ggplot(DFobsz, aes(x=pct_wet_w_f, y = pct_imp_w_emplog)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x = pct_wet_w_f, y = vscivcpmi)) + geom_point() + facet_wrap(vars(eco_region))

# Watershed-riparian wetlands
summary(DFobs$pct_wet_wr)
DFobs <- DFobs %>%
  mutate(
    pct_wet_wr_f = case_when(
      pct_wet_wr == 0 ~ "none",
      pct_wet_wr > 0 & pct_wet_wr <= 0.170 ~ "low",
      pct_wet_wr > 0.170 & pct_wet_wr <= 8.530 ~ "medium",
      pct_wet_wr > 8.530 & pct_wet_wr <= 52.570 ~ "high",
    )
  )

DFobs$pct_wet_wr_f <- factor(DFobs$pct_wet_wr_f, levels = c("none", "low", "medium", "high"))

summary(as.factor(DFobs$pct_wet_wr_f))

## binary absence/presence of watershed riparian wetlands
DFobs <- DFobs %>%
  mutate(
    bin_wet_wr = case_when(
      pct_wet_wr == 0 ~ 0,
      pct_wet_wr > 0 ~ 1
    )
  )
# DFobs$bin_wet_wr <- factor(DFobs$bin_wet_wr, levels = c("absence", "presence"))
summary(as.factor(DFobs$bin_wet_wr))


# ggplot(DFobs, aes(x=bin_wet_wr, y = vscivcpmi)) + geom_boxplot()

# ggplot(DFobs, aes(x=pct_wet_wr_f, y = tothab)) + geom_boxplot()
# 
# ggplot(DFobsz, aes(x=pct_wet_wr_f, y = pct_imp_w_emplog)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x = pct_wet_wr_f, y = vscivcpmi)) + geom_point() + facet_wrap(vars(eco_region))

# Catchment wetlands
summary(DFobs$pct_wet_c)

DFobs <- DFobs %>%
  mutate(
    pct_wet_c_f = case_when(
      pct_wet_c == 0 ~ "none",
      pct_wet_c > 0 & pct_wet_c <= 0.040 ~ "low",
      pct_wet_c > 0.040 & pct_wet_c <= 2.940 ~ "medium",
      pct_wet_c > 2.940 & pct_wet_c <= 62.90 ~ "high",
    )
  )

DFobs$pct_wet_c_f <- factor(DFobs$pct_wet_c_f, levels = c("none", "low", "medium", "high"))

summary(as.factor(DFobs$pct_wet_c_f))

summary(as.factor(DFobs$pct_wet_wr_f))

## binary absence/presence of catchment wetlands
DFobs <- DFobs %>%
  mutate(
    bin_wet_c = case_when(
      pct_wet_c == 0 ~ 0,
      pct_wet_c > 0 ~ 1
    )
  )
# DFobs$bin_wet_c <- factor(DFobs$bin_wet_c, levels = c("absence", "presence"))
summary(as.factor(DFobs$bin_wet_c))


# ggplot(DFobs, aes(x=bin_wet_c, y = vscivcpmi)) + geom_boxplot()
 
# ggplot(DFobs, aes(x=pct_wet_c_f, y = tothab)) + geom_boxplot()
# 
# ggplot(DFobsz, aes(x=pct_wet_c_f, y = pct_imp_w_emplog)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x = pct_wet_c_f, y = vscivcpmi)) + geom_point() + facet_wrap(vars(eco_region))

# Catchment-Riparian wetlands
# Change to none, low, or high
summary(DFobs$pct_wet_cr)

DFobs <- DFobs %>%
  mutate(
    pct_wet_cr_f = case_when(
      pct_wet_cr == 0 ~ "none",
      pct_wet_cr > 0 & pct_wet_cr <= 11.85 ~ "low",
      pct_wet_cr > 11.85 ~ "high",
    )
  )

DFobs$pct_wet_cr_f <- factor(DFobs$pct_wet_cr_f, levels = c("none", "low", "high"))

summary(as.factor(DFobs$pct_wet_cr_f))

## binary absence/presence catchment riparian wetlands
DFobs <- DFobs %>%
  mutate(
    bin_wet_cr = case_when(
      pct_wet_cr == 0 ~ 0,
      pct_wet_cr > 0 ~ 1
    )
  )
# DFobs$bin_wet_cr <- factor(DFobs$bin_wet_cr, levels = c("absence", "presence"))
summary(as.factor(DFobs$bin_wet_cr))


# ggplot(DFobs, aes(x=bin_wet_cr, y = vscivcpmi)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x=pct_wet_cr_f, y = tothab)) + geom_boxplot()
# 
# ggplot(DFobsz, aes(x=pct_wet_cr_f, y = pct_imp_w_emplog)) + geom_boxplot()
# 
# ggplot(DFobs, aes(x = pct_wet_cr_f, y = vscivcpmi)) + geom_point() + facet_wrap(vars(eco_region))

# Catchment-Riparian percent impervious
# pct_im_rp_c
summary(DFobs$pct_imp_rp_c)

DFobs <- DFobs %>%
  mutate(
    pct_imp_rp_c_f = case_when(
      pct_imp_rp_c == 0 ~ "none",
      pct_imp_rp_c > 0 & pct_imp_rp_c <= 1.015 ~ "low",
      pct_imp_rp_c > 1.015 ~ "high"
    )
  )

DFobs$pct_imp_rp_c_f <- factor(DFobs$pct_imp_rp_c_f, levels = c("none", "low", "high"))

summary(as.factor(DFobs$pct_imp_rp_c_f))

# ggplot(DFobs, aes(x=pct_imp_rp_c_f, y = vscivcpmi)) + geom_boxplot()


# Catchment percent grass cover
## binary absence/presence
DFobs <- DFobs %>%
  mutate(
    bin_grs_c = case_when(
      pct_grs_c == 0 ~ 0,
      pct_grs_c > 0 ~ 1
    )
  )
# DFobs$bin_grs_c <- factor(DFobs$bin_grs_c, levels = c("absence", "presence"))
summary(as.factor(DFobs$bin_grs_c))

# ggplot(DFobs, aes(x=bin_grs_c, y = vscivcpmi)) + geom_boxplot()

# Catchment-riparian percent grass cover
DFobs <- DFobs %>%
  mutate(
    bin_grs_cr = case_when(
      pct_grs_cr == 0 ~ 0,
      pct_grs_cr > 0 ~ 1
    )
  )
# DFobs$bin_grs_cr <- factor(DFobs$bin_grs_cr, levels = c("absence", "presence"))
summary(as.factor(DFobs$bin_grs_cr))

# ggplot(DFobs, aes(x=bin_grs_cr, y = vscivcpmi)) + geom_boxplot()

# tp set to low, medium, and high
summary(DFobs$tp)

DFobs <- DFobs %>%
  mutate(
    tp_f = case_when(
      tp <= 0.02 ~ "low",
      tp > 0.02 & tp <= 0.05 ~ "medium",
      tp > 0.05 ~ "high"
    )
  )

DFobs$tp_f <- factor(DFobs$tp_f, levels = c("low", "medium", "high"))

summary(as.factor(DFobs$tp_f))

# ggplot(DFobs, aes(x=tp_f, y = vscivcpmi)) + geom_boxplot()
```

## 3.40 Transform Onsite and Area Covariates

Turb and no3 had NA entries.
Emma Jones, VDEQ, sent email with csv file, FinalFieldAnalyte_MissingJamesSites.csv, that had analyte entries for stations 2AMLC000.84 and 2-WLS023.10.
The entries I made to James.ssn_101623 appear to have come in as character.
The columns in DFobs of nh4,no3,tkn, ortho_p, turb, and tss are charcter.
Convert them to numeric, log transform some of the covariates and put DFobs2 back into SSN.

Not considering total suspended solids, total phosphorus, total Kjeldahl nitrogen, and ammonia for analysis as logged scatter plots of those variables showed vertical bars of points.
That pattern might arise from detection limits.
Need to speak to VDEQ as to how those variables might be used.

```{r transform_onsite}

glimpse(DFobs) # shows pct_for_c and pct_for_w as chr
# In ArcGIS these 2 fields are not numeric so have to mutate
DFobs2 <- DFobs %>% mutate_at(c('pct_for_c', 'pct_for_w'), as.numeric)

DFobs2$l_no3 <- log(DFobs2$no3)
DFobs2$l_tn <- log(DFobs2$tn)
DFobs2$l_tp <- log(DFobs2$tp)
# DFobs2$l_tss <- log(DFobs2$tss)
DFobs2$l_turb <- log(DFobs2$turb)
DFobs2$l_tds <- log(DFobs2$tds)
DFobs2$l_spc <- log(DFobs2$sp_cond)
DFobs2$vsci <- round(DFobs2$vscivcpmi,1)
# DFobs2$l_catkm2 <- log(DFobs2$AreaSqKM)
# DFobs2$l_wskm2 <- log(DFobs2$TotDASqKM)
DFobs2$l_ws_area <- log(DFobs2$ws_area)
DFobs2$l_ws_area_r <- log(DFobs2$ws_area_r)
DFobs2$l_cat_area <- log(DFobs2$cat_area)
DFobs2$l_cat_area_r <- log(DFobs2$cat_area_r)

# commented out summary stats and EDA code below
# summary(DFobs2$pct_for_w)
# summary(DFobs2$pct_for_wr)
# summary(DFobs2$pct_for_c)
# summary(DFobs2$pct_for_cr)
# 
# summary(DFobs2$pct_imp_w)
# summary(DFobs2$pct_imp_rp_w)
# summary(DFobs2$pct_imp_c)
# summary(DFobs2$pct_imp_rp_c)
# 
# summary(DFobs2$pct_grs_w)
# summary(DFobs2$pct_grs_wr)
# summary(DFobs2$pct_grs_c)
# summary(DFobs2$pct_grs_cr)
# 
# summary(DFobs2$pct_hay_w)
# summary(DFobs2$pct_hay_wr)
# summary(DFobs2$pct_hay_c)
# summary(DFobs2$pct_hay_cr)
# 
# summary(DFobs2$l_tp)
# ggplot(DFobs2, aes(tp)) + stat_ecdf()
# 
# str(DFobs2)
# 
# str(DFobs2$turb)
# str(DFobs2$no3)
# 
# #across years
# ggplot(DFobs2, aes(x=as.factor(year), y=vsci)) + geom_boxplot(varwidth = TRUE)
# ggplot(DFobs2, aes(x=as.factor(year), y=tothab)) + geom_boxplot()
# #across stream orders
# ggplot(DFobs2, aes(x=as.factor(order), y=tothab)) + geom_boxplot()
# 
# #across different polygon geographies
# ggplot(DFobs2, aes(x=eco_region, y=tothab)) + geom_boxplot()
# ggplot(DFobs2, aes(x=eco_region, y=vsci)) + geom_boxplot()
# # clear downward trend from mountains to coast
# ggplot(DFobs2, aes(x=no3)) + geom_histogram()
# ggplot(DFobs2, aes(x=no3)) + geom_boxplot()
# ggplot(DFobs2, aes(x=no3)) + geom_density()
# 
# ggplot(DFobs2, aes(x = l_tp, y = vsci)) + geom_point() + geom_smooth()
# 
# ggplot(DFobs2, aes(x = log(tss), y = vsci)) + geom_point()
# 
# ggplot(DFobs2, aes(x = tothab, y = vsci)) + geom_point() + facet_wrap(vars(year))
# ggplot(DFobs2, aes(x = tn, y = no3)) + geom_point()
# 
# 
# summary(DFobs2$turb)
# 
# summary(DFobs2$no3)
# 
# ggplot(DFobs2, aes(x= tothab)) + geom_boxplot() + facet_wrap(vars(year))
# 
# ggplot(DFobs2, aes(x=turb, y = as.factor(year))) + geom_boxplot()
# 
# #64.  Tidyverse five number summary
# tothab_by_year <- DFobs2 %>%
#   group_by(year) %>%
#   summarize(min = fivenum(tothab)[1],
#             low = fivenum(tothab)[2],
#             median = fivenum(tothab)[3],
#             mean = mean(tothab, na.rm = TRUE),
#             upper = fivenum(tothab)[4],
#             max = fivenum(tothab)[5],
#             count = dplyr::n()) %>% 
#   print(n = Inf)
# 
# pct_imp_w_by_vahusb <- DFobs2 %>%
#   group_by(vahusb) %>%
#   summarize(min = fivenum(pct_imp_w)[1],
#             low = fivenum(pct_imp_w)[2],
#             median = fivenum(pct_imp_w)[3],
#             mean = mean(pct_imp_w, na.rm = TRUE),
#             upper = fivenum(pct_imp_w)[4],
#             max = fivenum(pct_imp_w)[5],
#             count = dplyr::n()) %>% 
#   print(n = Inf)
# 
# vsci_by_vahusb <- DFobs2 %>%
#   group_by(vahusb) %>%
#   summarize(min = fivenum(vsci)[1],
#             low = fivenum(vsci)[2],
#             median = fivenum(vsci)[3],
#             mean = mean(vsci, na.rm = TRUE),
#             upper = fivenum(vsci)[4],
#             max = fivenum(vsci)[5],
#             count = dplyr::n()) %>% 
#   print(n = Inf)
# 
# ggplot(DFobs2, aes(x = l_tn, y = vsci, colour = tp_f)) + geom_point() + geom_smooth()
# 
# ggplot(DFobs2, aes(x = l_tn, y = vsci, colour = tp_f)) + geom_point() + geom_smooth(method = "lm", se = TRUE)

```

## 3.41 New Variables into SSN
All the variable manipulations have been done on DFobs. We now need to put those changes into the SSN object using the SSN2::ssn_put_data function.
```{r put DFobs2 into SSN}
names(DFobs2)

class(DFobs2)
# str(DFobs2$ju)
# note ssn_put_data requires sf object and SSN2 object
j_ssn2 <-  SSN2::ssn_put_data(DFobs2,j_ssn1a)
# just doing this assignment so not have to rename objects
j_ssn3 <- j_ssn2

```

## 3.42 SFS Log_tn and tothab Plots
Scatter plots used in SFS 2024 talk.
```{r sfs_logtn_plot}
# scat_tn1 <- ggplot(DFobs2, aes(x=l_tn, y=vsci)) + geom_point() + geom_smooth(method = "lm", se = FALSE) + labs(x="Total Nitrogen (mg/L, logged)", y = "VSCI")
# 
# png(file="figures_sfs/scat_tn1.png",width=6,height=3,units="in",res=150)
# scat_tn1
# dev.off()
# 
# scat_tothab1 <- ggplot(DFobs2, aes(x=tothab, y=vsci)) + geom_point() + geom_smooth(method = "lm", se = FALSE) + labs(x="Total Habitat", y = "VSCI")
# 
# png(file="figures_sfs/scat_tothab1.png",width=6,height=3,units="in",res=150)
# scat_tothab1
# dev.off()

```

# 4.0 Empirical Logit Transformation

## 4.10 Empirical Logit

This code came from Travis Linscome-Hatfield at Neptune. Applying empirical logit only to obs.

Another transformation explored is the empirical logit for the covariates "pct_for_W", "pct_imp_w", "pct_crop_w", "pct_hay_w", "pct_grs_w", etc. The empirical logit code needs to be applied to a data frame without any geometry attached. This means using with st_set_geometry(DFobs2, NULL) to make that data frame. After empirical logit transformation is done, then join data frame to sf object and put the sf object into the SSN.

$$ 
\begin{eqnarray}
P_i \, & \text{are} & \, \text{a set of percentages} \\
f_i & = & P_i/100 \, \text{are a set of frequencies} \\
F_\text{int} \, & \text{is} & \, \text{the set of} \, f_i \text{with 0s and 1s removed} \\
\delta & = & \min \left\{ F_\text{int} \, \cup \, 1-F_\text{int} \right\} \\
 & = & \min_{F_{int}} \left\{  f_i, 1-f_i   \right\}  \\
f^*_i & = & \left\{   
\begin{array}{ll}
\delta / 2, & f_i \leq \delta \\
1 - \delta / 2, & f_i \geq 1 - \delta \\
f_i, & \text{otherwise}
\end{array}
\right. \\
l_i & = &\log\left( \frac{f^*_i}{1-f^*_i} \right)
\end{eqnarray}
$$

```{r emplogit_eda}

################################################################
################################################################
# DFobs2 %>% dplyr::select(ends_with(("_c")))
# DFobs2 %>% dplyr::select(starts_with(("pct_"))) %>% distinct(.)

# VARIABLE ADJUSTMENT ZONE 4
### Variables to apply empirical logit transformation
emplog_vars <- c("pct_for_w","pct_imp_w","pct_crop_w","pct_hay_w","pct_grs_w","pct_shrb_w","pct_for_wr","pct_imp_rp_w","pct_crop_wr","pct_hay_wr","pct_grs_wr","pct_shrb_wr","pct_for_c","pct_imp_c","pct_crop_c","pct_hay_c","pct_grs_c", "pct_shrb_c","pct_for_cr","pct_imp_rp_c","pct_crop_cr","pct_hay_cr","pct_grs_cr","pct_shrb_cr")

# remove geometry so empirical logit can be applied
DFobsz <- st_set_geometry(DFobs2, NULL)
################################################################
################################################################


## transform these variables and put the new values into new columns in DFobs
for(var in emplog_vars){
  ### create new tranformed data column to preserve the original
  new_nm <- paste0(var,"_emplog")
  dat_vec_obs <- DFobsz[,var]
  # dat_vec_preds <- DFpreds[,var]
  
  #converting to 0-1 range
  dat_vec_obs <- dat_vec_obs/100
  # dat_vec_preds <- dat_vec_preds/100
  
  # dat_vec[dat_vec == 1] <- .9999
  # dat_vec[dat_vec == 0] <- .0001
  
  if(any(dat_vec_obs > 1 | dat_vec_obs < 0)){
    cat("ERROR: percentage variables outside logical bounds")
  }
  
  small_dat_vec_obs <- dat_vec_obs[dat_vec_obs <1 & dat_vec_obs >0]
  # small_dat_vec_preds <- dat_vec_preds[dat_vec_preds <1 & dat_vec_preds >0]
  op1_obs <- small_dat_vec_obs
  op2_obs<- 1-small_dat_vec_obs
  # op1_preds <- small_dat_vec_preds
  # op2_preds <- 1-small_dat_vec_preds
  
  ## minimum of op1 op2
  delt_obs <- min(c(op1_obs,op2_obs))
  # delt_preds <- min(c(op1_preds,op2_preds))
  
  ## getting set of frequencies
  freqs_obs <- NULL
  for(i in 1:length(dat_vec_obs)){
    if(dat_vec_obs[i] <= delt_obs){
      freqs_obs[i] <- delt_obs/2
    }else if(dat_vec_obs[i] >= 1- delt_obs){
      freqs_obs[i] <- 1-(delt_obs/2)
    }else{
      freqs_obs[i] <-dat_vec_obs[i]
    }
  }
  
 # freqs_preds <- NULL
 # for(i in 1:length(dat_vec_preds)){
 #   if(dat_vec_preds[i] <= delt_preds){
 #     freqs_preds[i] <- delt_preds/2
 #   }else if(dat_vec_preds[i] >= 1- delt_preds){
 #     freqs_preds[i] <- 1-(delt_preds/2)
 #   }else{
 #     freqs_preds[i] <-dat_vec_preds[i]
  #  }
 # }
  
  ##getting logits
  logits_obs <- log(freqs_obs/(1-freqs_obs))
  DFobsz[,new_nm] <- logits_obs
  
  # logits_preds <- log(freqs_preds/(1-freqs_preds))
  # DFpreds[,new_nm] <- logits_preds
}

```

## 4.20 Logit of Impervious Cover

On 04/30/2024 summaries of empirical logit of impervious cover show that minimum for impervious cover of catchment is greater than minima for watershed, watershed-riparian, and catchment-riparian.
I think the minimum non-zero proportions of 0.01 versus 0.02 is causing the different transformation values.
But, I can't get "hand" calculator versions to work.
Also see that impervious cover values of 0 and the second smallest value under each of the 4 geographies is assigned the empirical logit transformation.
Does result below have something to do with second smallest value for pct_imp_w, pct_imp_rp, and pct_imp_rp_c is 0.01, but for pct_imp_c it is 0.02?

Look at 2024 PowerPoint in presentations folder, and see VSCI and Impervious Cover: 4 Geographies, slide 11, and note that minimum values of impervious cover for Catchment differs from other 3 geographies.

```{r emplogit_of_impevious}
# names(DFobsz)
# summary(DFobsz$pct_imp_w_emplog)
# summary(DFobsz$pct_imp_rp_w_emplog)
# summary(DFobsz$pct_imp_c_emplog)
# summary(DFobsz$pct_imp_rp_c_emplog)
# 
# str(DFobsz$pct_imp_w_emplog)
# str(DFobsz$pct_imp_rp_w_emplog)
# str(DFobsz$pct_imp_c_emplog)
# str(DFobsz$pct_imp_rp_c_emplog)
# 
# zlogit <- dplyr::select(DFobsz, c(st_id_tren, pct_imp_w, pct_imp_w_emplog,pct_imp_rp_w, pct_imp_rp_w_emplog, pct_imp_c, pct_imp_c_emplog, pct_imp_rp_c, pct_imp_rp_c_emplog))

```

## 4.30 Logit variables into SSN

```{r logit_into_ssn}
names(DFobsz)

DFobsz2 <- dplyr::select(DFobsz, c(station_id_2, pct_for_w_emplog:pct_shrb_cr_emplog))
# put transformed covariates in an SF object
DFobs2a <- full_join(DFobs2, DFobsz2, by = join_by(station_id_2))

# put SF object into SSN
j_ssn3 <-  SSN2::ssn_put_data(DFobs2a,j_ssn3)

## produce eda.shape plots for each of the transformed variables
## I got error with code below
# for(var in emplog_vars_transformed){
#   ## eda 4 panel plot for all covariates
#   covariate_eda <- eda.shape(x = DFobsz[,var] , var_name = var, conf.bands=F)
#   print(covariate_eda)
# }

# commented out EDA code below
# unconvinced of impervious surface and wetland interaction
# ggplot(DFobs2a, aes(x=pct_imp_rp_c_emplog, y = vsci, colour = pct_wet_cr_f)) + geom_point() + geom_smooth()

# interaction plots of vsci with 3 covariates having high correlations to vsci crossed with wetland factor
# ggplot(DFobs2a, aes(x=pct_imp_rp_c_emplog, y = vsci, colour = pct_wet_cr_f)) + geom_point() + geom_smooth(method = lm, se = TRUE)

# ggplot(DFobs2a, aes(x=pct_for_cr_emplog, y = vsci, colour = pct_wet_cr_f)) + geom_point() + geom_smooth(method = lm, se = TRUE)
# 
# ggplot(DFobs2a, aes(x=tothab, y = vsci, colour = pct_wet_cr_f)) + geom_point() + geom_smooth(method = lm, se = TRUE)

# perhaps interaction between shrubs and ecoregions, but recall southeastern plains has n =7, but this line is similar to Piedmont line
# ggplot(DFobs2a, aes(x=pct_shrb_cr_emplog, y = vsci, colour = eco_region)) + geom_point() +geom_smooth(method =lm, se = TRUE)
# 
# ggplot(DFobs2a, aes(x=pct_shrb_cr_emplog, y = vsci, colour = eco_region)) + geom_point() +geom_smooth(method =lm, se = TRUE)

# Golden et al. 2016 JAWRA mentioned interaction of urban areas and TP. Slopes and intercepts overlap in plots with impervious surface area and forest.
# ggplot(DFobs2a, aes(x = pct_imp_c_emplog, y = vsci, colour = vahusb)) + geom_point() + geom_smooth(method = "lm", se = FALSE) + labs(x="Catchment Impervious (empirical logit)", y = "VSCI") + expand_limits(x=-10) + scale_colour_discrete(name="Subbasin")
# 
# ggplot(DFobs2a, aes(x = pct_for_w_emplog, y = vsci, colour = tp_f)) + geom_point() + geom_smooth(method = "lm", se = TRUE)

# this plot suggested by result from WsRpWq analysis
# ggplot(DFobs2a, aes(x = pct_for_wr_emplog, y = vsci, colour = eco_region)) + geom_point() + geom_smooth(method = "lm", se = TRUE)
# 
# ggplot(DFobs2a, aes(x = pct_for_wr_emplog, y = vsci, colour = vahusb)) + geom_point() + geom_smooth(method = "lm", se = TRUE)
# 
# ggplot(DFobs2a, aes(x = pct_imp_w_emplog, y = vsci, colour = vahusb)) + geom_point() + geom_smooth(method = "lm", se = TRUE)
# 
# ggplot(DFobs2a, aes(x = pct_imp_w_emplog, y = vscivcpmi)) + geom_point() + geom_smooth(method = "lm") + facet_wrap(vars(vahusb))
# 
# ggplot(DFobs2a, aes(x = do, y = vsci, colour = vahusb)) + geom_point() + geom_smooth(method = "lm", se = TRUE)
# 
# ggplot(DFobs2a, aes(x = vahusb, y = pct_imp_w_emplog)) + geom_boxplot()

# check suggested by Mike Dumelle if y-interecept at x = 0 looks different from y-intercept when x > 0
# ggplot(DFobsz, aes(x=pct_grs_c_emplog, y = vsci)) + geom_point() + geom_smooth(method = "lm")

```

## 5.50 SFS StreamCat Plots

```{r sfs_streamcat}
# scat_cat2 <- ggplot(DFobs2a, aes(x = pct_imp_c_emplog, y = vsci, colour = vahusb)) + geom_point() + geom_smooth(method = "lm", se = FALSE) + labs(x="Catchment Impervious (empirical logit)", y = "VSCI") + expand_limits(x=-10) + scale_colour_discrete(name="Subbasin")
# 
# scat_cat2
# png(file="figures_sfs/scat_cat2.png",width=4,height=3,units="in",res=150)
#   scat_cat2
# dev.off()

```

#6.0 ESDA Across Different Polygons Using vahusb as those 5 polygons are comprised of huc8s.
Commented out EDA code below
```{r esda_polygons}

#by regions
# ggplot(DFobs, aes(x=factor(vahusb), y=vscivcpmi)) + geom_boxplot(varwidth = TRUE)
# 
# ggplot(DFobs, aes(x=factor(eco_region), y=vscivcpmi)) + geom_boxplot(varwidth = TRUE)
# 
# ggplot(DFobs, aes(x=factor(bio_region), y=vscivcpmi)) + geom_boxplot(varwidth = TRUE)
# 
# ggplot(DFobsz, aes(x=pct_imp_rp_w_emplog, y = vscivcpmi, colour = eco_region)) + geom_point() + geom_smooth(se = TRUE, method = lm)
# 
# ggplot(DFobsz, aes(x = pct_imp_rp_w_emplog, y = l_tn)) + geom_point() + geom_smooth()
# 
# ggplot(DFobsz, aes(x = pct_imp_w_emplog, y = l_tn)) + geom_point() + geom_smooth()
# 
# ggplot(DFobsz, aes(x=tothab, y = vscivcpmi, colour = eco_region)) + geom_point() + geom_smooth(se = TRUE, method = lm)
# 
# ggplot(DFobsz, aes(x = l_tn, y = vscivcpmi, colour = eco_region)) + geom_point() + geom_smooth()
# 
# ggplot(DFobsz, aes(x = pct_imp_rp_w_emplog, y = l_tn)) + geom_point()

# ggplot(DFobsz, aes(x = PctImp_Ws_emplog, y = l_tn)) + geom_point() + geom_smooth(se = TRUE, method = lm)

# ggplot(DFobsz, aes(x=vahusb, y=vscivcpmi)) + geom_boxplot(varwidth = TRUE)
# 
# mapview(DFobs2)
```

## 7.0 Instream/Onsite Correlations

Commented out EDA code below.

```{r instream_onsite_corr}
# DFobs3 <- dplyr::select(DFobsz, st_id_tren, vsci, do:tss, tothab, l_no3, l_tn, l_tp, l_tds)

# instream
# ggpairs(DFobs3, columns = c(2:5,8,13:15))
# ggpairs(DFobs3, columns = c(2:5,8,15,18:19))
# nutrients
# ggpairs(DFobs3, columns = c(2,6:12,16:18))

# WQ and WS variables
# DFobs4 <- dplyr::select(DFobsz, st_id_tren, vsci, do,p_h, tothab, l_no3:l_spc, pct_imp_w_emplog, pct_grs_w_emplog, pct_hay_w_emplog, pct_shrb_w_emplog, pct_for_w_emplog)
# names(DFobs4)
# 
# ggpairs(DFobs4, columns = c(2:15))
# 
# # WQ and WS-Riparian
# DFobs5 <- dplyr::select(DFobsz, st_id_tren, vsci, do,p_h, tothab, l_no3:l_spc, pct_imp_rp_w_emplog, pct_grs_wr_emplog, pct_hay_wr_emplog, pct_shrb_wr_emplog, pct_for_wr_emplog)
# names(DFobs5)
# ggpairs(DFobs5, columns = c(2:15))

# WQ and Catchment
# DFobs6 <- dplyr::select(DFobsz, st_id_tren, vsci, do,p_h, tothab, l_no3:l_spc, pct_imp_c_emplog, pct_grs_c_emplog, pct_hay_c_emplog, pct_shrb_c_emplog, pct_for_c_emplog)
# 
# names(DFobs6)
# 
# ggpairs(DFobs6, columns = c(2:15))

# WQ and Catchment-Riparian
# DFobs7 <- dplyr::select(DFobsz, st_id_tren, vsci, do,p_h, tothab, l_no3:l_spc, pct_imp_rp_c_emplog, pct_grs_cr_emplog, pct_hay_cr_emplog, pct_shrb_cr_emplog, pct_for_cr_emplog)
# names(DFobs7)
# 
# ggpairs(DFobs7, columns = c(2:15))
# ggpairs(DFobs7, columns = c(2,11:15))

# Mix of Watershed, Catchment, Riparian
# DFobs8 <- dplyr::select(DFobsz, st_id_tren, vsci, pct_imp_w_emplog, pct_imp_c_emplog, pct_imp_rp_c_emplog, pct_for_wr_emplog, elev_ws, elev_cat)
# names(DFobs8)
# 
# ggpairs(DFobs8, columns = c(2:8))
```

# 8.0 Script 03 SSN Output
The ssn object with new covariates and transformed covariates is saved so it can be load into Script 04.
```{r script03_ssn_output}
save(j_ssn3, file = "outputs/j_ssn3.RData")
```


# ZZZ Subset Test

Cannot get subset working, return this later.

```{r subsetSSN_test}

# ssn_2018 <- ssn_subset(j_ssn3, path= "/ssn_objects/subset_test" ,subset = year == 2018)


```