waterTempThresholds.Rmd

---
title: "Data sketches of water temperatures and thresholds"
knit: (function(input_file, encoding) {
  out_dir <- 'docs';
  rmarkdown::render(input_file,
 encoding=encoding,
 output_file=file.path(dirname(input_file), out_dir, 'index.html'))})

output: html_document
---

```{r setup, include=FALSE}
here::i_am("waterTempThresholds.Rmd")

knitr::opts_chunk$set(
	echo = FALSE,
	message = FALSE,
	warning = FALSE
)
library(plotly)
library(tidyverse)
library(dplyr)
library(magrittr)

library(chron)
library(stringr)
library(htmltools)

library(lubridate)
library(dataRetrieval)

# plotting tools
library(ggplot2)
library(ggthemes)
library(ggExtra)
library(ggtext)
library(gghighlight)
library(scales)
library(hrbrthemes)
# Spatial packages

library(sf)
library(maptools)
library(RColorBrewer)
library(mapview)
library(knitr)
library(rmarkdown)

# library(weatherData)
library(lunar)
# Tidal data
library(VulnToolkit) 
library(Tides)

# timeseries 
library(tsibble)
library(fable)
library(feasts)
library(xts)
library(here)
library(units)

source(here::here("../LI_SaltMarshMonitoring/R/truBlue.R")) # lots of code went into trublue utilities from another project. consider making a package.
library(lunar)
library(suncalc)
library(extrafont)

```

### Water temperature stations

Exploring temperature data in the Peconics.
2 USGS guage stations exist with continuous data collection. 


```{r sites}
h5("Data presented within this document have been made available through the",
      a(href = 'http://waterdata.usgs.gov/nwis', "USGS NWIS Web Interface:")
    )
# estuarySites <- whatNWISsites(stateCd = "NY", 
#               # water temp C
#               parameterCd = "00010") %>% 
#   filter(site_tp_cd == "ES") %>% 
#   pull(site_no)

peconicSiteIDs <- c("01304200", "01304562")

# pulling continuous data down from USGS takes some time so this will run 
# only when passing update as TRUE  
update = FALSE

if (update) {
  # get unit values for sites to explore
  continuousTemps <- readNWISuv(
    siteNumbers = peconicSiteIDs,
    startDate = "2013-01-01",
    #first full year of data collection 2013.
    endDate = Sys.Date(),
    # statCd = c("00001", "00002", "00003"),
    parameterCd = c("00010", "62620")
  ) %>%
    dataRetrieval::renameNWISColumns()
 
  # filter out coded -9999s and 9999s
  continuousTemps <-
    continuousTemps %>% filter(Wtemp_Inst < 100 &
                                 Wtemp_Inst > -10) 
  # save data to save time and processing each time.
  write_rds(continuousTemps, here("data/continuousTemps.rds"))
} else{
  continuousTemps <- read_rds(here("data/continuousTemps.rds"))
}


estuarySiteLocations <- readNWISsite(peconicSiteIDs) %>% 
  st_as_sf(coords = c("dec_long_va", "dec_lat_va"), crs = 4269) %>% 
  mutate(stationName = case_when(
    station_nm == "ORIENT HARBOR AT ORIENT NY" ~ "Orient Harbor",
    station_nm == "PECONIC RIVER AT COUNTY HWY 105 AT RIVERHEAD NY" ~ "Peconic River Rte 105"
  ))

# popup 

	popupLink <- function(siteName, siteNo){
		url <- paste('<a href="https://waterdata.usgs.gov/nwis/inventory?agency_code=USGS&site_no=', 
			     siteNo,'"','> Access data and site information from NWIS </a>', sep = "")
		
		popupFinal <- paste(sep = "<br/>", 
				    siteName, 
				    url)
		popupFinal
	}

stationMap <- mapview(estuarySiteLocations, popup = popupLink(siteName = estuarySiteLocations$stationName, siteNo = estuarySiteLocations$site_no))
cntrcords <- c(mean(st_coordinates(estuarySiteLocations)[,1]),
               mean(st_coordinates(estuarySiteLocations)[,2]))
stationMap@map %>% leaflet::setView(cntrcords[1], cntrcords[2], zoom = 10)

```

## Plots
 
Data has been collected since 2012, with 2013 being the first 'full' year of data collection.


```{r tempdata, echo=FALSE}
# go tsibble.
theme_set(theme_ipsum())

temps <- continuousTemps %>% 
  # join to get station names
  left_join(estuarySiteLocations %>% 
              st_drop_geometry() %>% 
              select(site_no, stationName), by = "site_no") %>% 
  as_tsibble(key = stationName, index = dateTime) %>% 
  rename(tidalHt = X_NAVD88_62620_Inst) %>% 
  # fill_gaps() %>% 
  mutate(month = month(dateTime),
         year = year(dateTime),
         day_of_year = yday(dateTime),
         commonDate = as.Date(paste0("2000-",format(dateTime, "%j")), "%Y-%j"),
         year2019 = ifelse(year == 2019, 1, 0)) %>% 
  # Calc ecdf first -
  group_by(stationName, year) %>% 
  mutate(calc_ecdf = ecdf(Wtemp_Inst)(Wtemp_Inst),
         inv_calc_ecdf = 1- calc_ecdf,
         totTime = length(Wtemp_Inst) * 6,
         yearmins = 60*24*365.25, # minutes in a year
         mins_yr = units::set_units(yearmins, minutes),
         hours_yr = units::set_units(mins_yr, hours),
         hours_above = hours_yr * inv_calc_ecdf,
         hours_num = drop_units(hours_above),
         days_num = drop_units(set_units(hours_above, days)))

```


```{r tempplot, echo=FALSE}

temps  %>% 
  ggplot(aes(x = commonDate, y = Wtemp_Inst, group = year)) +
  geom_line() + 
  gghighlight(year == 2022) +
  scale_color_gradient(low = "grey", high = "blue", aesthetics = "color") +
  # scale_color_viridis_c(option = "inferno") +
  # facet_grid(stationName ~ year) +
  facet_grid(stationName ~ .) +
  scale_x_date(labels = function(x) format(x, "%d-%b")) +
  theme(panel.spacing = unit(.5, "lines")) +
  # scale_x_continuous(labels = function(x) format(as.Date(as.character(x), "%j"), "%d-%b")) +
  theme(legend.position = "none", axis.text.x = element_text(angle = 90, hjust = 1)) + 
  labs(y = "Temp (c)", x = NULL, title = "Water temperature through the year", subtitle = "2019 shown in blue")


```


```{r}


temps  %>% 
  ggplot(aes(x = commonDate, y = Wtemp_Inst, group = year)) +
  geom_line() + 
  gghighlight(year == 2022) +
  # scale_color_gradient(low = "grey", high = "blue", aesthetics = "color") +
  # scale_color_viridis_d(option = "inferno") +
  # facet_grid(stationName ~ year) +
  facet_grid(stationName ~ .) +
  scale_x_date(labels = function(x) format(x, "%d-%b")) +
  theme(panel.spacing = unit(.5, "lines")) +
  scale_x_continuous(labels = function(x) format(as.Date(as.character(x), "%j"), "%d-%b")) +
  theme(legend.position = "none", axis.text.x = element_text(angle = 90, hjust = 1)) + 
  labs(y = "Temp (c)", x = NULL, title = "Water temperature through the year", subtitle = "2019 shown in blue")



```


### Cumulative Distribution Approach - 
More to do... And think about... 

Focus on the temps above threshold temps?

```{r message=FALSE, warning=FALSE, paged.print=FALSE}
plot_ecdf_temps <- temps %>% filter(year != 2022) %>% 
  ggplot(aes(x =  Wtemp_Inst, color = as.character(year))) + 
  stat_ecdf(pad = F, geom = "line") + 
  scale_y_reverse() +
  facet_wrap(~ stationName, nrow = 2) +
  labs(color = NULL) +
  # scale_y_continuous(trans = "time")
  # scale_x_continuous(limits = c(min = 25, max = 30)) + 
  NULL

plot_ecdf_temps

ggplotly(plot_ecdf_temps)



plot_time_adjustedcdf <- temps %>% group_by(stationName, year) %>% arrange(hours_num) %>%  
  filter(Wtemp_Inst > 20) %>% 
  ggplot(aes(x = Wtemp_Inst, y = days_num, group = as.character(year))) + 
  geom_line() + 
  geom_ribbon(data ={ . %>% filter(Wtemp_Inst > 25)}, aes(ymin = 0, ymax = days_num), fill = 'red') +
  ylim(c(0, 100)) + 
  geom_hline(aes(yintercept = 50)) + 
  facet_grid(stationName ~ year) +
  labs(color = NULL, fill = NULL) +
  scale_x_continuous(breaks = c(20, 22, 24, 26, 28)) + 
  theme(legend.position = "none") + 
  labs(y = "Number of days at or below temperature X",
       x = "Temperature")


plot_time_adjustedcdf

# ggplotly(plot_time_adjustedcdf)

```


```{r temp_raster}
# Not enough daily variation to make this effective.

# 
# 
water_temp_raster <- temps %>% 
  filter(month %in% c(4:10)) %>% 
  filter(year != 2022) %>% 
  as_tibble() %>% 
  filter(stationName ==  "Peconic River Rte 105") %>%
  separate(dateTime, into = c("Date", "measTime"), sep = " ") %>%
  	mutate(measTime = as.POSIXct(.$measTime, format = "%H:%M:%S", tz = 'GMT'),
  	       time = as.character(.$measTime, format = "%H:%M:%S", tz = 'GMT'),
  	       date = as.Date(Date)) %>% # timezone of raw data is GMT
  append_truBlu(date_col = "date") %>% 
  ggplot(aes(y = commonDate, 
             x = measTime, 
             fill = Wtemp_Inst, 
             interpolate = T)) + #, 
             # alpha = 1 - (max(Wtemp_Inst) - Wtemp_Inst)/max(Wtemp_Inst))) +
  geom_raster() +
	scale_fill_viridis_c(option = "inferno") +
  	# scale_fill_viridis_c(option = "viridis") +
  	# scale_fill_viridis_c(option = "plasma") +
	scale_x_datetime(breaks = NULL) + 
  scale_y_date(date_labels = "%B", date_breaks = "1 month") +
  theme_few() +
  theme(legend.position = "left", 
        panel.spacing.y = unit(5, "points"),
        panel.spacing.x = unit(15, "points")) +
  guides(alpha = "none") + 
  labs(y = NULL, x = NULL, fill = "Water temperature C")


# perhaps more intuitive longitudinal look with date across horizontal.
water_temp_raster + 
  facet_grid(year ~ stationName, switch = "y") + 
  coord_flip() 

water_temp_raster + 
  facet_grid(year ~ stationName, switch = "y") + 
  coord_flip() 


# # vertical days - time across horizontal 
# water_temp_raster +
#   facet_grid(. ~ year)
## Tidal height ----
tidal_ht_raster <- temps %>% filter(month %in% c(5:10)) %>% 
  as_tibble() %>% filter(stationName !=  "Peconic River Rte 105") %>% 
  separate(dateTime, into = c("Date", "Time"), sep = " ") %>%
  	mutate(Time = as.POSIXct(.$Time, format = "%H:%M:%S", tz = 'GMT')) %>% # timezone of raw data is GMT
  # ggplot(aes(y = commonDate, x = Time, fill = Wtemp_Inst)) +
  ggplot(aes(y = commonDate, x = Time, fill = tidalHt)) +
		theme_bw() +
  geom_raster(interpolate = T) +
	scale_fill_viridis_c(option = "cividis", aesthetics = c("color", "fill")) +
  # scale_fill_gradient(low = "blue", high = "red") +
		# scale_x_datetime(breaks = date_breaks('1 hour'), labels = date_format("%H")) +
    scale_x_datetime(labels = NULL) +
  scale_y_date(breaks = date_breaks('1 month'), labels = date_format("%b")) + 
  facet_grid(. ~ year) + 
  theme(legend.position = "bottom") +
  labs(y = NULL, x = NULL)

tidal_ht_raster


temps %>% filter(month %in% c(5:10)) %>% 
  as_tibble() %>% filter(stationName !=  "Peconic River Rte 105") %>% 
  separate(dateTime, into = c("Date", "Time"), sep = " ") %>%
  	mutate(Time = as.POSIXct(.$Time, format = "%H:%M:%S", tz = 'GMT')) %>% # timezone of raw data is GMT
  HL(level = .$tidalHt, time = .$dateTime) %>% HL.plot()


library(patchwork)
water_temp_raster + tidal_ht_raster

## Moon phasses----

temps %>% 
  filter(month %in% c(4:10)) %>% 
  filter(year != 2022) %>% 
  as_tibble() %>% 
  filter(stationName ==  "Peconic River Rte 105") %>%
  separate(dateTime, into = c("Date", "measTime"), sep = " ") %>%
  	mutate(measTime = as.POSIXct(.$measTime, format = "%H:%M:%S", tz = 'GMT'),
  	       time = as.character(.$measTime, format = "%H:%M:%S", tz = 'GMT'),
  	       date = as.Date(Date)) %>% # timezone of raw data is GMT
  filter(moon.font == "O") %>% 
  append_truBlu(date_col = "date") %>% 
  ggplot(aes(x = commonDate, y = time)) +
  geom_text(aes(label = moon.font)) + #TODO: Moon font not showing..

# library(rayshader)
# 
# a <- water_temp_raster + 
#   facet_grid(year ~ stationName, switch = "y") + 
#   coord_flip()
# 
# plot_gg(ggobj = a, width = 5, height = 5, multicore = TRUE, scale = 200 )
  
```


## Differencing
```{r}



temps %>% filter(month %in% c(5:10)) %>% 
  as_tibble() %>%
  filter(stationName !=  "Peconic River Rte 105") %>% 
  separate(dateTime, into = c("Date", "Time"), sep = " ") %>%
  	mutate(Time = as.POSIXct(.$Time, format = "%H:%M:%S", tz = 'GMT'),
  	       diffed = difference(Wtemp_Inst, lag = 1)) %>% 
   ggplot(aes(y = commonDate, x = Time, fill = diffed)) +
		theme_bw() +
  geom_raster(interpolate = T) +
	scale_fill_viridis_c(option = "cividis", aesthetics = c("color", "fill")) +
  # scale_fill_gradient(low = "blue", high = "red") +
		# scale_x_datetime(breaks = date_breaks('1 hour'), labels = date_format("%H")) +
    scale_x_datetime(labels = NULL) +
  scale_y_date(breaks = date_breaks('1 month'), labels = date_format("%b")) + 
  facet_grid(. ~ year) + 
  theme(legend.position = "bottom") +
  labs(y = NULL, x = NULL)

  

```


## Tidal characteristics


```{r}
library(Tides)

temps %>% select(dateTime, tidalHt) %>% 
  filter_index("2020-05-01" ~ "2020-06-30") %>% 
  ggplot(aes(x = dateTime, y = tidalHt)) +
  geom_path()

temps %>% select(dateTime, tidalHt, Wtemp_Inst) %>% 
  filter_index("2020-05-01" ~ "2020-06-30") %>% autoplot(vars(tidalHt, difference(tidalHt)))

temps %>% 
  filter(stationName !=  "Peconic River Rte 105") %>% 
  filter_index("2020-05-01" ~ "2020-06-30") %>% 
  as_tibble() %>%
  ungroup() %>% 
  select(time = dateTime, h = tidalHt) %>% # must rename the columns to time and h. 
  tidyr::drop_na() %>%
  extrema(h0 = 0)

tideOrient <- temps %>% 
  filter(stationName !=  "Peconic River Rte 105") %>% 
  filter_index("2020-06-01" ~ "2020-06-15") %>% 
  as_tibble() %>%
  ungroup() %>% 
  select(time = dateTime, h = tidalHt) %>% # must rename the columns to time and h. 
  tidyr::drop_na() %>% extrema(h0 = 0) 

View(tideOrient$h)
View(tideOrient$HL)

tidal_c_orient <-  temps %>%
  filter(stationName !=  "Peconic River Rte 105") %>%
  # filter_index("2020-06-01" ~ "2020-06-15") %>%
  as_tibble() %>%
  ungroup() %>%
  select(time = dateTime, h = tidalHt) %>% # must rename the columns to time and h.
  tidyr::drop_na() %>%
  TidalCharacteristics(h0 = 0, hoffset = .5)

```