air_water_temps

Set custom plotting themes

EPSCoR SCTC Aquatic Ecology

Analysis in progress

R Session Info

sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.2
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plotrix_3.7-4      lme4_1.1-19        Matrix_1.2-15     
##  [4] ggridges_0.5.1     anytime_0.3.3      knitr_1.21        
##  [7] nnet_7.3-12        cowplot_0.9.4      FSA_0.8.22        
## [10] stringi_1.2.4      janitor_1.1.1      ggpmisc_0.3.0     
## [13] ggrepel_0.8.0      broom_0.5.1        DT_0.5            
## [16] lubridate_1.7.4    forcats_0.3.0      stringr_1.3.1     
## [19] dplyr_0.7.8        purrr_0.2.5        readr_1.3.1       
## [22] tidyr_0.8.2        tibble_2.0.0       ggplot2_3.1.0     
## [25] tidyverse_1.2.1    googlesheets_0.3.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.0         lattice_0.20-38    assertthat_0.2.0  
##  [4] digest_0.6.18      R6_2.3.0           cellranger_1.1.0  
##  [7] plyr_1.8.4         backports_1.1.3    RApiDatetime_0.0.4
## [10] evaluate_0.12      httr_1.4.0         pillar_1.3.1      
## [13] rlang_0.3.1        lazyeval_0.2.1     readxl_1.2.0      
## [16] rstudioapi_0.9.0   minqa_1.2.4        nloptr_1.2.1      
## [19] rmarkdown_1.11     splines_3.5.2      htmlwidgets_1.3   
## [22] munsell_0.5.0      compiler_3.5.2     modelr_0.1.2      
## [25] xfun_0.4           pkgconfig_2.0.2    htmltools_0.3.6   
## [28] tidyselect_0.2.5   crayon_1.3.4       withr_2.1.2       
## [31] MASS_7.3-51.1      grid_3.5.2         nlme_3.1-137      
## [34] jsonlite_1.6       gtable_0.2.0       magrittr_1.5      
## [37] scales_1.0.0       cli_1.0.1          bindrcpp_0.2.2    
## [40] xml2_1.2.0         generics_0.0.2     tools_3.5.2       
## [43] glue_1.3.0         hms_0.4.2          yaml_2.2.0        
## [46] colorspace_1.3-2   rvest_0.3.2        bindr_0.1.1       
## [49] haven_2.0.0

This document explores relationships between air and water temperature in the Kenai EPSCoR Aquatic Ecology Project.

• a.) Import water temperature data.
• b.) Import air temperature data.
• c.) Determine air-water sensitivity values to be used in bioenergetics modeling.

Water temperatures

Import Water Temperature Data Note: this code differs somewhat from that in “season_temp_stats.Rmd”

# links?

### Identify source sheets
# Lower reach sites
# Water temperatures for lower sites have been curated so as to provide maximum constinuous temporal extent
lower_sites_2015 <- gs_title("lower_sites_2015")
lower_sites_2016 <- gs_title("lower_sites_2016")

## Special procedure for Ptarmigan Creek 2015 & 2016
# a.) Ptarmigan Creek 2015: We want to use the file generated/discussed from "2015 Water Temperature Analaysis.Rmd" for both (lower and middle) reaches.  Lower and Middle sites are similar enough (and close enough) that we want to use the same dataset for both sites.

# b.) Ptarmigan Creek 2016: we want to use the MPC-2016 site for both middle and lower sites; they are within tolerance margin (see "2016 Water Temperature Analaysis.Rmd").  MPC-2016 site has most continuous dataset.

# Methods narration:  We replaced missing water temperature data with data from the closest available logger site, if the two sites averaged an annual difference of < 0.2 C .  

# see "2015" and "2016" Water Temperature Data .Rmd files for methods.


# 2015
## lower sites
LBC15 <- lower_sites_2015 %>%
  gs_read(ws = "LBC15", colnames = TRUE) %>%
  transform(DateTime = mdy_hm(DateTime))
LRR15 <- lower_sites_2015 %>%
  gs_read(ws = "LRR15", colnames = TRUE) %>%
  transform(DateTime = mdy_hm(DateTime))
LPC15 <- lower_sites_2015 %>%
  gs_read(ws = "LPC15", colnames = TRUE) %>%
  transform(DateTime = mdy_hm(DateTime))
LKR15 <- lower_sites_2015 %>%
  gs_read(ws = "LKR15", colnames = TRUE) %>%
  transform(DateTime = mdy_hm(DateTime))


# 2016
## lower sites
LBC16 <- lower_sites_2016 %>%
  gs_read(ws = "LBC16", colnames = TRUE) %>%
  transform(Date = as.character(Date),
            DateTime = as.POSIXct(DateTime))
LRR16 <- lower_sites_2016 %>%
  gs_read(ws = "LRR16", colnames = TRUE) %>%
  transform(Date = as.character(Date),
            DateTime = as.POSIXct(DateTime))

LPC16 <- lower_sites_2016 %>%
  gs_read(ws = "MPC16", colnames = TRUE) %>%
  transform(Date = as.character(Date),
            DateTime = as.POSIXct(DateTime)) %>%
 select(-Hydrologic_Segment,-Site) %>%
  mutate(Hydrologic_Segment = "Lower",
         Site = paste(Hydrologic_Segment, Stream_Name))

LKR16 <- lower_sites_2016 %>%
  gs_read(ws = "LKR16", colnames = TRUE) %>%
  transform(Date = as.character(Date),
            DateTime = as.POSIXct(DateTime))

## middle kenai river
MKR15 <- gs_title("2015_MKR_H2OTemp_Final.csv") %>%
  gs_read(ws = "2015_MKR_H2OTemp_Final.csv") %>%
  transform(DateTime = mdy_hm(DateTime))
MKR16 <- gs_title("2016_MKR_H2OTemp_Final.csv") %>%
  gs_read(ws = "2016_MKR_H2OTemp_Final.csv") %>%
  transform(Date = as.character(Date),
            DateTime = as.POSIXct(DateTime))

## middle ptarmigan creek 2015
MPC15 <- lower_sites_2015 %>%
  gs_read(ws = "LPC15", colnames = TRUE) %>%
  select(-Hydrologic_Segment,-Site) %>%
  mutate(Hydrologic_Segment = "Middle",
         Site = paste(Hydrologic_Segment, Stream_Name)) %>%
  transform(Date = as.character(Date),
            DateTime = mdy_hm(DateTime))

## middle ptarmigan creek 2016
MPC16 <- lower_sites_2016 %>%
  gs_read(ws = "MPC16", colnames = TRUE) %>%
  transform(Date = as.character(Date),
            DateTime = as.POSIXct(DateTime))

# bind together lower site dataframes (+ middle kenai river)
lower_sites <- bind_rows(
                         # lower sites
                         LBC15,LRR15,LPC15,LKR15,
                         LBC16,LRR16,LPC16,LKR16,
                         # middle kenai sites
                         MKR15,MKR16,
                         # those special ones
                         MPC15,MPC16) %>%
  transform(DateTime = as.character(DateTime)) %>%
  select(-X1)

Import middle and upper site water temp data

# links to google sheet data sources here?

# other sites (middle and upper reaches)
BC15 <- gs_title("Beaver_Creek_2015_Water_Temperatures")
BC16 <- gs_title("Beaver_Creek_2016_Water_Temperatures")
RR15 <- gs_title("Russian_River_2015_Water_Temperatures")
RR16 <- gs_title("Russian_River_2016_Water_Temperatures")
PC15 <- gs_title("Ptarmigan_Creek_2015_Water_Temperatures")
PC16 <- gs_title("Ptarmigan_Creek_2016_Water_Temperatures")

Sys.sleep(6)

## Beaver Creek
#2015
MBC2015 <- BC15 %>% 
  gs_read(ws = "MBC_WL_10324669",col_names=TRUE)
UBC2015 <- BC15 %>% 
  gs_read(ws = "UBC_WL_10537880",col_names=TRUE)
#2016
MBC2016 <- BC16 %>% 
  gs_read(ws = "MBC-WL 2004069",col_names=TRUE)
UBC2016 <- BC16 %>%
  gs_read(ws = "UBC-WL 10718566",col_names=TRUE)
UBC1_2016 <- BC16 %>%
  gs_read(ws = "UBC1-WL 10886777")
UBC2_2016 <- BC16 %>%
  gs_read(ws = "UBC2-WL 10886775")

Sys.sleep(6)

# Russian River
# 2015
MRR2015 <- RR15 %>% 
  gs_read(ws = "MRR_WL_10324639",col_names=TRUE)
URR2015 <- RR15 %>% 
  gs_read(ws = "URR_WL_10324638",col_names=TRUE)
#2016
MRR2016 <- RR16 %>% 
  gs_read(ws = "MRR-WL 10886776",col_names=TRUE)  # 9/11/18 confirmed no unintentional data dupication in import
MRR1_2016 <- RR16 %>%
  gs_read(ws = "MRR1-WL 10324644", col_names = TRUE)
MRR2_2016 <- RR16 %>%
  gs_read(ws = "MRR2-WL 10324638", col_names = TRUE)
#URR2016 <- RR16 %>% 
#  gs_read(ws = "URR-WL 10886778", col_names=TRUE)     # hashtagged out for now; will use MRR2 site instead; see below
URR1_2016 <- RR16 %>%
  gs_read(ws = "URR1-WL 10767518", col_names = TRUE)

Sys.sleep(6)

##Ptarmigan Creek
#2015
# 2015 MPC already imported above
UPC2015 <- PC15 %>% 
  gs_read(ws = "UPC_WL_10324642",col_names=TRUE)
#2016
# 2016 MPC already imported above
#MPC2016 <- PC16 %>% 
#  gs_read(ws = "MPC-WL 10324639",col_names=TRUE)
UPC2016 <- PC16 %>% 
  gs_read(ws = "UPC-WL 10324648",col_names=TRUE)


# Bind middle & upper sites data
middle_upper_sites <- bind_rows(
  #2015
  MBC2015,UBC2015,
  MRR2015,URR2015,
  UPC2015,
  #2016
  MBC2016,UBC2016,UBC1_2016,UBC2_2016, 
  MRR2016,MRR1_2016,MRR2_2016,URR1_2016,
  UPC2016) 

# format middle & upper site dataframe to match lower sites dataframe
middle_upper_sites <- middle_upper_sites %>%
  select(-one_of("Latitude","Longitude","Logger Type","Sample Interval (min)",
         "Logger Serial #","Site Description","Logger Model")) %>%
  rename(Stream_Name = `Stream Name`,
         Hydrologic_Segment = `Hydrologic Segment`,
         Site_ID = `Site ID`)  %>%
  transform(DateTime = mdy_hm(DateTime)) %>%
  mutate(Date = as.character(date(DateTime)),
         Week = week(DateTime),
         Site = paste(Hydrologic_Segment,Stream_Name)) %>%
  transform(DateTime = as.character(DateTime)) %>%
  select(DateTime,Temp_C,Stream_Name,Hydrologic_Segment,
         Site_ID,Agency,Site,Date,Week)

# bind lower site dataframe with middle & upper site dataframe
fnl_temps <- bind_rows(lower_sites,middle_upper_sites) %>%
  select(-Date) %>%
  mutate(Date = date(DateTime),
         year = year(DateTime)) %>%
  # remove extraneous site label info from some parts of "Site" column
  separate(Site,sep = " ", into = c("reach","river","waterbody")) %>%
  unite(reach,river,waterbody, sep = " ",col = "Site") 


# remove temporary objects
rm(lower_sites,middle_upper_sites,
  MBC2015,UBC2015,
  MRR2015,URR2015,
  MPC2015,UPC2015,
  MBC2016,UBC2016,UBC1_2016,UBC2_2016, 
  MRR2016,MRR1_2016,MRR2_2016,URR2016,URR1_2016,
  MPC2016,UPC2016,
  BC15,BC16,RR15,RR16,PC15,PC16,MKR15,MKR16)
  
# remove temporary objects
rm(lower_sites_2015,lower_sites_2016,LBC15,LRR15,LPC15,LKR15,LBC16,LRR16,MPC16,LKR16)

Eliminate loggers unwanted for growth analyses (e.g. water Temp loggers with no associated air temp dataset nearby)

# Issue: in 2016, there exist multiple temperature loggers in a given reach but only one fish sampling site.
# Resolved: for growth analyses, use only the one logger that was closest to sampling site.

# create function to exclude multiple criteria from same column, opposite of "%in%"
'%notin%' <- Negate('%in%')

# filter out water temperature loggers non-adjacent to fish sampling sites.

# Also: Upper Russian River 2016 (URR-2016) water temp is not available due to logger malfunction.  Replace with closest downstream logger; MRR2-2016:
loggers_excluded <- c("UBC1-WL-2016",
                      "UBC2-WL-2016",
                      "MRR1-WL-2016", # just above lower russian lake
                      #"MRR2-2016",   # mid-way between lower and upper russian lakes; keep this site; will use for upper russian river 2016 because that dataset too is short to use (2016 logger failure).
                      # leaving upper ptarmigan creek site in for now.  lentic site, may not be comparable to others.
                      "URR-WL-2016",   
                      "URR1-WL-2016"
                      #,
                      # "UPC-WL-2015", # upper ptarmigan creek site; lentic
                      # "UPC-WL-2016" # upper ptarmigan creek site; lentic
                      )

# exclude unwanted loggers
fnl_temps_sub <- fnl_temps %>%
  filter(Site_ID %notin% loggers_excluded) %>%
  filter(!is.na(Stream_Name),
         !is.na(Temp_C)) %>%
  transform(DateTime = as.POSIXct(DateTime))

## Special procedure for Upper Russian River 2016
# change MRR2-WL-2016 site from middle to upper reach.
# (upper russian river 2016 temp data largely missing due to logger malfunction)
z <- fnl_temps_sub %>%
  filter(Site_ID == "MRR2-WL-2016") %>%
  mutate_all(funs(str_replace_all(., "Middle", "Upper"))) %>%
  mutate(DateTime = as.POSIXct(DateTime),
         Temp_C = as.numeric(Temp_C),
         Week = as.numeric(Week),
         Date = as.Date(Date),
         year = as.numeric(year))

# reunite MRR2 data with overall dataframe
fnl_temps_sub <- bind_rows(fnl_temps_sub,z) %>%
  # clean up the "Site" column so that only name is present; not Site_ID, is present in this column
  separate(Site,sep = " ", into = c("reach","river","waterbody")) %>%
  unite(reach,river,waterbody, sep = " ",col = "Site")

# remove presence of rows with Site_ID = MRR2-WL-2016 and Site = 'Middle Russian River' to eliminate duplication (we are using MRR-WL-2015 site for MRR2016 data)

# create table of unwanted data and anti_join it to original table
anti_fnl_temps_sub <- fnl_temps_sub %>%
  filter(Site == "Middle Russian River" & Site_ID == "MRR2-WL-2016")

fnl_temps_sub <- anti_join(fnl_temps_sub,anti_fnl_temps_sub)

rm(LPC16,MPC15)

Plot water temperature data

Sys.sleep(6)

# plot raw data to confirm it is sensical...?
fnl_temps_x <- fnl_temps_sub %>%
  mutate(doy = yday(DateTime)) %>%
  separate(Site, sep = " ","reach", remove = F) %>%
  group_by(Stream_Name,reach,Site,Date,year,doy) %>%
  summarise(daily_mean_temp = mean(Temp_C))

ggplot(fnl_temps_x, aes(doy,daily_mean_temp, color = as.factor(year))) +
  geom_line() +
  facet_grid(Stream_Name ~ reach) +
  theme_bw() + 
  ggtitle("2015-2016 Water Temperature Data")

An annual water temperature dataset associated with each site is present in the “fnl_temps_sub” dataframe (subsetted final temperatures).

Experiment with alternative plot formats for water temp data to visualize contrast

fnl_temps_sub <- fnl_temps_sub %>%
  #rename(Temp_C = water_temp_C) %>%
  mutate(Week = week(DateTime),
         month = month(DateTime))

# plot A
fnl_temps_sub %>%
  mutate(doy = yday(Date)) %>%
  group_by(Stream_Name,year,Date,doy) %>%
  summarise(mean_temp_h20_stream = mean(Temp_C),
         h2o_se = std.error(Temp_C)) %>%
  ggplot(aes(doy,mean_temp_h20_stream,color = Stream_Name)) +
  geom_line() +
  facet_grid(.~year) 

# plot B
fnl_temps_sub %>%
  filter(month >= 5 & month <= 9) %>%
  group_by(Stream_Name,year,Week) %>%
  summarise(mean_temp_h20_stream = mean(Temp_C),
         h2o_se = std.error(Temp_C)) %>%
  ggplot(aes(Week,mean_temp_h20_stream,color = Stream_Name)) +
  geom_line(size = 2.1) +
  facet_grid(.~year) +
  bm_theme_1 +
  theme(strip.text = element_text(face="bold", size=18),
        axis.text.x = element_text(size = 16,face="bold"),
        axis.text.y = element_text(size = 16,face="bold"),
        axis.title.x = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold")) +
  xlab("Week") +
  ylab("Mean Weekly Water Temp (C)")

# plot B
fnl_temps_sub %>%
  mutate(month = month(Date)) %>%
  filter(month >= 5 & month <= 9) %>%
  group_by(Stream_Name,year,month) %>%
  summarise(mean_temp_h20_stream = mean(Temp_C),
         h2o_se = std.error(Temp_C)) %>%
  ggplot(aes(month,mean_temp_h20_stream,color = Stream_Name)) +
  geom_line(size = 2.1) +
  bm_theme_1 +
  theme(strip.text = element_text(face="bold", size=18),
        axis.text.x = element_text(size = 16,face="bold"),
        axis.text.y = element_text(size = 16, face="bold"),
        axis.title.x = element_text(face = "bold")) +
  xlab("Month") +
  ylab("Mean Monthly Water Temp (C)") +
  facet_grid(.~year) 

# can't get generalized stream name facet labeller to work here

# reorder legend
# consider which one would offer best contrast w/ future temps

Create a plot constrasting different sites within a tributary

# plot raw data to confirm it is sensical...?
fnl_temps_x <- fnl_temps_sub %>%
  mutate(doy = yday(DateTime)) %>%
  separate(Site, sep = " ","reach", remove = F) %>%
  group_by(Stream_Name,reach,Site,Date,year,doy) %>%
  summarise(daily_mean_temp = mean(Temp_C))

# plot

#2015
fnl_temps_x %>%
  filter(year == "2015") %>%
ggplot(aes(doy,daily_mean_temp, color = as.factor(reach))) +
  geom_line() +
  facet_grid(Stream_Name ~ .) +
  theme_bw() + 
  ggtitle("2015 Water Temperature Data")

# 2016
fnl_temps_x %>%
  filter(year == "2016") %>%
ggplot(aes(doy,daily_mean_temp, color = as.factor(reach))) +
  geom_line() +
  facet_grid(Stream_Name ~ .) +
  theme_bw() + 
  ggtitle("2016 Water Temperature Data")

Create a plot contrasting select water temperature metrics

# import water temperature metrics from each individual logger
tempmetrics <- gs_read(gs_title("water_temperature_metrics")) %>%
  select(year,river,reach,MWMT,MWAT)

# prep facet strip labels
metrics_labels = c("MWAT" = "Max Weekly\nAverage Temp",
                   "MWMT" = "Max Weekly\nMax Temp")

# plot
tempmetrics %>%
  gather(metric,value,MWMT,MWAT) %>%
  filter(!is.na(value)) %>%
  ggplot(aes(river,value, color = river)) +
  scale_x_discrete(limits=c("Beaver Creek",
                            "Russian River",
                            "Ptarmigan Creek",
                            "Kenai River"),
                   labels=c("Lowland","Montane","Glacial","Main Stem")) +
  geom_boxplot() +
  facet_grid(.~metric, labeller = labeller(metric = metrics_labels)) +
  xlab("") +
  ylab("Water Temperature (C)") +
  bm_theme_2 +
  theme(strip.text.x = element_text(size = 16, margin = margin(.2,0,.2,0, "cm")))

Air temperatures

Import air temperatures

Sys.sleep(6)

# 2015 Air temps
## Beaver Creek
LBC15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "LBC-AL 10324664") 
MBC15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "MBC-AL 2004069")
UBC15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "UBC-AL 2004060")

Sys.sleep(6)

## Russian River
LRR15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "LRR-AL 10537877")
MRR15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "MRR-AL 10324659")
URR15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "URR-AL 1065516")

Sys.sleep(6)

## Ptarmigan Creek
LPC15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "LPC-AL 10324665")
MPC15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "MPC-AL 10537882")
UPC15_air <- gs_read(gs_title("Air_Temperatures_2015_EPSCoR_SCTC.xlsx"),ws = "UPC-AL 10324647")


# 2016 Air Temps

## Beaver Creek
LBC16_air <- gs_read(gs_title("Beaver_Creek_2016_Air_Temperatures.xlsx"),ws = "LBC-AL 10324647")
MBC16_air <- gs_read(gs_title("Beaver_Creek_2016_Air_Temperatures.xlsx"),ws = "MBC-AL 10324645")
UBC16_air <- gs_read(gs_title("Beaver_Creek_2016_Air_Temperatures.xlsx"),ws = "UBC-AL 2004072")

Sys.sleep(6)

## Russian River
LRR16_air <- gs_read(gs_title("Russian_River_2016_Air_Temperatures.xlsx"),ws = "LRR-AL 10324665")
MRR16_air <- gs_read(gs_title("Russian_River_2016_Air_Temperatures.xlsx"),ws = "MRR-AL 10324640")
URR16_air <- gs_read(gs_title("Russian_River_2016_Air_Temperatures.xlsx"),ws = "URR-AL 10537882")

Sys.sleep(6)

## Ptarmigan Creek
LPC16_air <- gs_read(gs_title("Ptarmigan_Creek_2016_Air_Temperatures.xlsx"),ws = "LPC-AL 2004060")
MPC16_air <- gs_read(gs_title("Ptarmigan_Creek_2016_Air_Temperatures.xlsx"),ws = "MPC-AL 10537875")
UPC16_air <- gs_read(gs_title("Ptarmigan_Creek_2016_Air_Temperatures.xlsx"),ws = "UPC-AL 10537880")

Sys.sleep(6)

# Kenai River sites associated air temps

# Lower Kenai River
# downloaded csv to local directory from NWS LCD service on 12/17/18
LKR_air <- read.csv("/Users/bmeyer/Google Drive/Thesis/R Analysis/Temperature Data/NWS_download/1578886.csv") %>%
  select(DATE,HOURLYDRYBULBTEMPC,LATITUDE,LONGITUDE) %>%
  transform(DATE = round_date(mdy_hm(DATE),"15 minutes")) %>%  # round to nearest 0.25 hrs
  mutate(Agency = "NWS",
         `Stream Name` = "Kenai River",
         Hydrologic.Segment = "Lower",
         Site.ID = "LKR-AL-2015",
         Logger.Type = "Air",
         site = "Lower Kenai River",
         Date = date(DATE),
         year = year(DATE)) %>%
  rename(DateTime = DATE,
         Temp_C = HOURLYDRYBULBTEMPC,
         Latitude = LATITUDE,
         Longitude = LONGITUDE) %>%
  transform(DateTime = as.character(DateTime))

LKR15_air <- LKR_air %>%
  filter(year == 2015)

LKR16_air <- LKR_air %>%
  filter(year == 2016)

# Middle Kenai River
# USGS site does not record useable air temperature date, nearby LRR-AL-2015 site should be an appropriate subsitute -- near intersection of lower russian and mainstem kenai
# Examine LRR-AL-2015 site against Skilak Guard Station RAWS data later to check appropriateness of choice (https://raws.dri.edu/cgi-bin/rawMAIN.pl?akASKI)
MKR_air <- bind_rows(LRR15_air,LRR16_air) %>%
  # rename stream name from russian river to kenai river
  mutate(`Stream Name` = str_replace(`Stream Name`,"Russian River","Kenai River"),
         `Hydrologic Segment` = str_replace(`Hydrologic Segment`,"Lower","Middle"),
         year = year(mdy_hm(DateTime))) 

MKR15_air <- MKR_air %>%
  filter(year == 2015)

MKR16_air <- MKR_air %>%
  filter(year == 2016)

# bind 2015 air temps
air15 <- bind_rows(LBC15_air,MBC15_air,UBC15_air,
                   LRR15_air,MRR15_air,URR15_air,
                   LPC15_air,MPC15_air,UPC15_air,
                   MKR15_air) %>%
  select(-Notes)

# bind 2016 air temps
air16 <- bind_rows(LBC16_air,MBC16_air,UBC16_air,
                   LRR16_air,MRR16_air,URR16_air,
                   LPC16_air,MPC16_air,UPC16_air,
                   MKR16_air) %>%
  select(-X10,-X15,-X16)

# bind 2015 and 2016 air temps
airtemps <- bind_rows(air15,air16) %>%
  transform(DateTime = mdy_hm(DateTime)) %>%
  select(-Logger.Model,-Site.Description,-Sample.Interval..min.,-Logger.Serial..,-Contact) %>%
  mutate(site = as.factor(paste(Hydrologic.Segment,Stream.Name)),
         Date = date(DateTime),
         year = year(DateTime)) %>%
  filter(!is.na(DateTime)) %>%
  transform(DateTime = as.character(DateTime)) %>%
  # include lower mainstem kenai site
  bind_rows(LKR15_air,LKR16_air)

# subset to date range of interest
## May 1 - Nov 1, 2015 and 2016

z15 <- airtemps %>%
  filter(Date >= as.Date("2015-05-01") & Date <= as.Date("2015-11-01"))
z16 <- airtemps %>%
  filter(Date >= as.Date("2016-05-01") & Date <= as.Date("2016-11-01"))
airtemps <- bind_rows(z15,z16)

# remove temporary objects
#rm(LBC15_air,MBC15_air,UBC15_air,
#   LRR15_air,MRR15_air,URR15_air,
#   LPC15_air,MPC15_air,UPC15_air,
#   LBC16_air,MBC16_air,UBC16_air,
#   LRR16_air,MRR16_air,URR16_air,
#   LPC16_air,MPC16_air,UPC16_air,
#   LKR_air,LKR15_air,LKR16_air)

rm(z,z1,z15,z16)

Create plot showing extent of available data air temperature logger deployments

# calculate min and max extent of air temperature logger deployments
air_logger_times <- airtemps %>%
  group_by(site,year,Stream.Name,Agency) %>% 
  summarize(logger.deploy.start = min(Date),
            logger.deploy.end = max(Date),
            days = logger.deploy.end - logger.deploy.start) %>%
  mutate(logger.deploy.start.doy = yday(logger.deploy.start),
         logger.deploy.end.doy = yday(logger.deploy.end)) %>%
  transform(year = as.factor(year))


#  make plot that visualizes extent of each logger deployment faceted by year and river

# order facets
air_logger_times$site <- factor(air_logger_times$site, levels = 
                                    c("Lower Beaver Creek",
                                      "Middle Beaver Creek",
                                      "Upper Beaver Creek",
                                      "Lower Russian River",
                                      "Middle Russian River",
                                      "Upper Russian River",
                                      "Lower Ptarmigan Creek",
                                      "Middle Ptarmigan Creek",
                                      "Lower Kenai River",
                                      "Middle Kenai River"))


# color plot
 p <- air_logger_times %>%
   filter(!is.na(site)) %>%
   ggplot() +
  geom_segment(aes(x = logger.deploy.start.doy, xend= logger.deploy.end.doy,
                   y = year, yend = year, color = Agency),
               size = 3) + 
  facet_grid(site~.) +
  scale_x_continuous(breaks = c(152,182,213,244),
                     labels = c("June","July","August","Sept")) +
  theme_bw() +
  ggtitle("Air Temperature Logger\nDeployment Periods") +
  xlab("") +
  ylab("") +
  theme(strip.text.y = element_text(angle = 0, size = 14)) +
  theme(plot.title = element_text(size= 30, face = "bold", hjust = 0.3)) +
  theme(legend.title=element_text(size=24, face = "bold")) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        legend.text = element_text(size = 14)) +
  labs(colour="Agency")
 
 p

Plot raw air temperature data May 1 - Nov 1

# NOTE: it appears that there is once a day (~midnight) every day that has a missing air temp observation from the kenai airport air temp data set.  These must be filtered out in order to calculate daily means.

# plot raw data to confirm it is sensical...?
airtemps_means <- airtemps %>%
  filter(!is.na(Temp_C)) %>%
  mutate(doy = yday(DateTime)) %>%
  group_by(Stream.Name,Hydrologic.Segment,site,Date,year,doy) %>%
  summarise(daily_mean_airtemp = mean(Temp_C))

ggplot(airtemps_means, aes(doy,daily_mean_airtemp, color = as.factor(year))) +
  geom_line() +
  facet_grid(Stream.Name ~ Hydrologic.Segment) +
  theme_bw() + 
  ggtitle("2015-2016 Air Temperature Data")

NOTE: Air temperatures look visually anomalous for middle/upper Beaver Creek and Middle/Upper Russian river. BM verified on 1/17/19 that these data were recorded during logger deployment by manually checking deployment/retrieval dates. Consistency across multiple sites/rivers indicates this is likely “real”.

Join water and air temperature data in to same dataframe

# join water and air temps in same dataframe by DateTime and site

## prep water temps dataframe for join
fnl_temps_sub <- fnl_temps_sub %>%
  transform(DateTime = as.POSIXct(DateTime)) %>%
  rename(water_temp_C = Temp_C,
         Agency_water = Agency,
         Site_ID_water = Site_ID) %>%
  select(-Week)

## prep air temps dataframe for join
airtemps <- airtemps %>%
  transform(DateTime = as.POSIXct(DateTime)) %>%
  rename(air_temp_C = Temp_C,
         Stream_Name = Stream.Name,
         Hydrologic_Segment = Hydrologic.Segment,
         Site_ID_air = Site.ID,
         Site = site,
         Agency_air = Agency) %>%
  select(-Latitude,-Longitude,-Logger.Type)


# left_join of water and airtemps dfs 
temps <- left_join(fnl_temps_sub,airtemps, by = c("DateTime","year","Date","Hydrologic_Segment","Stream_Name","Site"))

Plot air-water data to visualize extent and overall shape(s)

# plot
temps %>%
  ggplot(aes(air_temp_C,water_temp_C, color = as.factor(year))) +
  geom_jitter(width = 2) +
  facet_grid(Stream_Name ~ Hydrologic_Segment) +
  theme_bw()

# summary table
temps %>%
  filter(!is.na(air_temp_C)) %>%
  group_by(Site) %>%
  summarise(n_obs = n())

Notes: Lower Kenai River site is less data-rich because air temperature data from the Kenai Airport is hourly rather than every 0.25 hrs.

a.) restrict to common time period for all sites (annual values for latest deployment and earliest retrieval.)

• For this exercise at least, it will be necessary to exclude lower russian river 2015 DateTime when calculating earliest/lastest deployment. LRR-2015 did not record water temp dat until 6/22/15.

# calculate latest deployment and earliest retrieval date for each year
logger_dates <- temps %>%
  filter(Site != "Lower Russian River") %>%
  group_by(Site,year) %>%
  summarise(deploy_date = min(DateTime),
            retrieve_date = max(DateTime)) %>%
  ungroup() %>%
  group_by(year) %>%
  summarise(max_deploy = max(deploy_date),
            min_retrieve = min(retrieve_date))

# restrict air/water temperature data to common time period calculated above
temps_mod  <- temps %>% 
    rowwise() %>% 
    mutate(present = any()) %>%
    mutate(present = any(DateTime >= logger_dates$max_deploy & DateTime <= logger_dates$min_retrieve)) %>% 
    filter(present) %>% 
    data.frame() %>%
    select(-present) 

Basic tables of ranges of temperatures throughout study area(s) and study period

# total range of instantaneous temps
temps_mod %>%
  summarise(maxtemp = max(water_temp_C),
            mintemp = min(water_temp_C))

# total range of daily mean temps
z <- temps_mod %>%
  group_by(Date,Site) %>%
  summarise(meantemp = mean(water_temp_C)) %>%
  ungroup() %>%
  summarise(max_daily = max(meantemp),
            min_daily = min(meantemp))
z

Sensitivity relationships

Calculate sensitivity and r2 at annual site level using time frame common to all sites

  
# formula structure
formula <- y ~ x

# formula and r2 labels
formula_labels <- stat_poly_eq(aes(label = paste("atop(",..eq.label.., ",", ..adj.rr.label.., ")", sep = "")),
               formula = formula, 
               rr.digits = 2, 
               coef.digits = 2, 
               parse = TRUE,
               size = 3,
               geom = "label", 
               alpha = 0.8)

# plot water vs air temp for each site/year
temps_mod %>%
  ggplot(aes(air_temp_C,water_temp_C)) +
  geom_jitter(width = 2) +
  facet_wrap(year~Site) +
 # formula_labels +
  theme_bw()

Notes: sensitivity varies greatly when raw data rather than daily or weekly means used (e.g same site can have strong or weak relationship from year to year when plotted at time scale of whole year). Not a meaningful scale of interpretation.

Plot overall relationship between daily mean water temperature and air temperature, using time frame common to all sites

# manual order of facets
temps_mod <- temps_mod %>%
  transform(Stream_Name = factor(Stream_Name, levels=c("Beaver Creek",
                                                       "Russian River",
                                                       "Ptarmigan Creek",
                                                       "Kenai River")))

# calc r2 values for overall watershed
temps_mod %>%
  filter(!is.na(water_temp_C),
         !is.na(air_temp_C)) %>%
  mutate(month = month(DateTime)) %>%
  group_by(Site,Date,month,year,Stream_Name,Hydrologic_Segment) %>%
  summarise(daily_air_C = mean(air_temp_C),
            daily_h2o_C = mean(water_temp_C)) %>%
  ggplot(aes(daily_air_C,daily_h2o_C)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(.~Stream_Name) +
  xlab("Daily Mean Air Temperature (C)") +
  ylab("Daily Mean Water Temperature (C)") +
  bm_theme_1 

#+
 # formula_labels

# calc r2 values by site
temps_mod %>%
  filter(!is.na(water_temp_C),
         !is.na(air_temp_C)) %>%
  mutate(month = month(DateTime)) %>%
  group_by(Site,Date,month,year,Stream_Name,Hydrologic_Segment) %>%
  summarise(daily_air_C = mean(air_temp_C),
            daily_h2o_C = mean(water_temp_C)) %>%
  ggplot(aes(daily_air_C,daily_h2o_C)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(Hydrologic_Segment~Stream_Name) +
  xlab("Daily Mean Air Temperature (C)") +
  ylab("Daily Mean Water Temperature (C)") +
  theme_bw() +
  bm_theme_1 

#+
 # formula_labels

Weekly mean values are more likely to be integrative of air-water temps. Calculate site-specific sensitivity values based on weekly means.

Calculate sensitivity (slope) and r2 for each site/week/year (Hastagged out 2/27/19)

# reduce data to tibble of site/year/week; use data within to calculate slope, r2, and n()

# air-water temp sensitivity formula
# air_h20_sens <- function(df) {
#  lm(water_temp_C ~ air_temp_C, data = df)}

# calculate slope of air-water relationship at site/year/week level
## note that here we are not using the dataset that is retricted to a common time frame for all sites
# sens_vals <- temps %>%
#  filter(!is.na(water_temp_C),
#         !is.na(air_temp_C))%>%
#  mutate(week = week(DateTime),
#         doy = yday(DateTime)) %>%
#  group_by(Site,Hydrologic_Segment,year,week) %>%
#  nest() %>%
  # calculate slope for weekly relationship of air-water temps
#  mutate(sens = map(data,air_h20_sens),
#         tidy_sens = map(sens,broom::tidy),
#         weekly_slope = map_dbl(tidy_sens,~.$estimate[2])) %>%
  # calculate r2 for weekly relationship of air-water temps
#  mutate(fit = map(data,air_h20_sens),
#         summary = map(fit,glance),
#         r_sq = map_dbl(summary,"r.squared")) %>%
#  select(Site,Hydrologic_Segment,year,week,weekly_slope,r_sq)

# table of values
# datatable(sens_vals, caption = "Air-Water Sensitivity Values by Site/Year/Week", 
#          filter = 'top', options = list(
#  pageLength = 5, autoWidth = TRUE)) %>%
#  formatSignif(c("weekly_slope","r_sq"), digits = 3) 

#Notes: The above approach (site/year/week sensitivity values) may employ a temporal scale too specific to be useful in relating to future projected air temperatures.  Consider other approaches, e.g.:  Mauger et al 2017 uses "3 years of weekly mean July temperatures," calculating sensitivity from a linear mixed-effects model with year as a random intercept.

#"We used a linear mixed effects model with a random intercept for year to account for the temporal correlation between weekly values within a year. We used Akaike’s information criterion (AIC) to compare results from a mixed effects model to one without a random intercept; the mixed effects model resulted in better model fit at 29 of the 48 sites."

#For my work, that would produce only 8 data points (2 years x 4 weeks per month) from which to acquire a slope (monthly site sensitivity value).  Thought: using weekly mean values is probably the minimum timeframe we want for average values.

#Other approaches: Dynamic Factor Analysis (DFA) approach used in Lisi 2015 and Winfree 2018 looks robust, but fairly daunting for me to figure out how to understand and execute.

#Let's try the identical approach as in Mauger 2017 to see how things look. (Do we get rational/useful results?):

2/22/19

Conversation w/ E. schoen, take-away: I have too few years of data to do site/month -specific sensitivity values (max 8 points per site), but site-specific values (inclusive of all months) could work (20-30 points per site; 2 yrs x ~14 weeks per site). If this approach is later though to be insufficiently robust, create 3 models and compare via AIC:

1.) r2 ~ site 2.) r2 ~ site + month 3.) r2 ~ site + season (early vs. late summer)

---

---

2/24/19 TO DO: convert work below here to reflect calculations at scale of site rather than site/month.

---

---

A.) Calculate model fit (r2) values for sensitivity relationships NOT using year as a random intercept.

# Create generic names for facet labels
stream_facets <- c(
  "Beaver Creek" = "Lowland", 
    "Russian River" = "Montane",
    "Ptarmigan Creek" = "Glacial", 
    "Kenai River" = "Main \nStem")

# calculate weekly mean water and air temp values (if at least 70% of data in that week captured)

# note: a full week of observations for records every 0.25 hrs is 672 observations, a full week of observations every 1.0 hrs is 168 (e.g. mainstem kenai river sites).  Create a column specifying 672 obs per week for all sites except mainstem kenai, 168 for mainstem kenai sites

# create intermediate dataframe to join to specify number of max observations per week
# using time frame common to all sites
weekly_obs <- temps_mod %>%
  distinct(Site,Stream_Name) %>%
  mutate(max_weeklyobs = ifelse(Stream_Name == "Kenai River",168,672))
  
# calculate mean weekly temperatures by site/year
meantemps <- temps %>% 
  filter(!is.na(air_temp_C)) %>%
  mutate(month = month(DateTime),
         week = week(DateTime)) %>%
  group_by(month,year,week,Site,Hydrologic_Segment) %>%
  mutate(n_week_obs = n()) %>%
  group_by(month,week,year,Site,Stream_Name,Hydrologic_Segment,n_week_obs) %>%
  summarise(mean_h20 = mean(water_temp_C),
            mean_air = mean(air_temp_C))
  
# it appears that at the lower kenai site (maybe the middle one too?) observation frequency was not always just hourly (some have > 168 obs/week).  That should be OK; specifify diff in methods though.

# join and calculate fraction of maximum possible observations present per week/site.  Retain weeks with >70% of observations present
meantemps <- left_join(meantemps,weekly_obs) %>%
  mutate(frac_week_obs = n_week_obs/max_weeklyobs) %>%
  # filter out weeks with less than 70% of observations present
  filter(frac_week_obs >= .7) %>%
  # exclude site/months with less than 4 week-values
  group_by(Site,month,Hydrologic_Segment) %>%
  mutate(n_weeks = n()) %>%
  filter(n_weeks > 3) %>%
# manually order appearance of sites
  transform(Stream_Name = factor(Stream_Name, levels=c("Beaver Creek",
                                         "Russian River",
                                         "Ptarmigan Creek",
                                         "Kenai River"))) 


## prep formula labels for each plot

# formula structure
formula <- y ~ x

# formula and r2 labels
formula_labels <- stat_poly_eq(aes(label = paste("atop(",..eq.label.., ",", ..adj.rr.label.., ")", sep = "")),
               formula = formula, 
               rr.digits = 2, 
               coef.digits = 2, 
               parse = TRUE,
               size = 3,
               geom = "label", 
               alpha = 0.8)


# plot seperated by site
meantemps %>%
  filter(month != 9,
         month != 10) %>%
  ggplot(aes(mean_air,mean_h20), color = year) + 
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(Hydrologic_Segment ~ Stream_Name, labeller = labeller(Stream_Name = stream_facets) ) +
  theme_bw() +
  theme(strip.text.y = element_text(angle = 360)) +
  ggtitle("Air-Water Temperature Sensitivity Relationships") +
 # formula_labels +
  stat_quadrant_counts(quadrants = 0L, labels.range.x = 11, labels.range.y = 18) 

# plot with sites within each drainage, segregated
z <- meantemps %>%
  filter(month != 9,
         month != 10) %>%
  ggplot(aes(mean_air,mean_h20, color = Hydrologic_Segment)) + 
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ Stream_Name, labeller = labeller(Stream_Name = stream_facets) ) +
  theme_bw() +
  theme(strip.text.y = element_text(angle = 360)) +
  #ggtitle("Air-Water Temperature Sensitivity") +
  stat_quadrant_counts(quadrants = 0L, labels.range.x = 11, labels.range.y = 18) +
  xlab("Weekly Mean Air Temperature (C)") +
  ylab("Weekly Mean Water Temperature (C)") +
  bm_theme_1 +
  theme(strip.text = element_text(face="bold", size=18)) +
  scale_y_continuous(breaks=seq(6, 16, 2), limits = c(8,16)) 
## with legend
z

## without legend
z + theme(legend.position = "none")

# plot with sites within each drainage, pooled
meantemps %>%
  filter(month != 9,
         month != 10) %>%
  ggplot(aes(mean_air,mean_h20)) + 
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(. ~ Stream_Name, labeller = labeller(Stream_Name = stream_facets) ) +
  theme_bw() +
  theme(strip.text.y = element_text(angle = 360)) +
 #ggtitle("Air-Water Temperature Sensitivity") +
 #formula_labels +
 #stat_quadrant_counts(quadrants = 0L, labels.range.x = 11, labels.range.y = 18) +
  xlab("Weekly Mean Air Temperature (C)") +
  ylab("Weekly Mean Water Temperature (C)") +
  bm_theme_1 +
  theme(strip.text = element_text(face="bold", size=18)) +
  scale_y_continuous(breaks=seq(8, 16, 2), limits = c(8,16)) +
  stat_quadrant_counts(quadrants = 0L, labels.range.x = 11, labels.range.y = 18)

Plot air-water senstivity value comparison

# air-water temp sensitivity formula
 air_h20_sens <- function(df) {
  lm(mean_h20 ~ mean_air, data = df)}

# calculate slope of air-water relationship at site level
## note that here we are using the dataset that is retricted to a common time frame for all sites
 sens_vals <- meantemps %>%
   filter(!is.na(mean_air),
          !is.na(mean_h20)) %>%
  group_by(Stream_Name,Hydrologic_Segment) %>%
  nest() %>%
  # calculate slopes for site-specific air-water temps
  mutate(sens = map(data,air_h20_sens),
         tidy_sens = map(sens,broom::tidy),
         site_slope = map_dbl(tidy_sens,~.$estimate[2])) %>%
  # calculate r2 for weekly relationship of air-water temps
  mutate(fit = map(data,air_h20_sens),
         summary = map(fit,glance),
         r_sq = map_dbl(summary,"r.squared")) %>%
  select(Stream_Name,Hydrologic_Segment,site_slope,r_sq)

# table of values
 datatable(sens_vals, caption = "Air-Water Sensitivity Values by Site", 
          filter = 'top', options = list(
  pageLength = 5, autoWidth = TRUE)) %>%
  formatSignif(c("site_slope","r_sq"), digits = 3)

visualize site-specific slope and r2 values

# create generic labels for streams
stream_labs <- scale_x_discrete(breaks=c("Beaver Creek", "Russian River", "Ptarmigan Creek", "Kenai River"),
                      labels=c("Lowland", "Montane", "Glacial", "Main Stem"))

## slopes
sens_vals %>%
  ggplot(aes(Stream_Name,site_slope, fill = Hydrologic_Segment)) +
  geom_col(position = position_dodge2(preserve = "single", padding = 0, width = 3)) +
  theme_bw() +
 # ggtitle("Air-Water Sensitivity Values") +
  xlab("") +
  ylab("Sensitivity Slope Value") +
  stream_labs +
  bm_theme_2

## r2 values
sens_vals %>%
  ggplot(aes(Stream_Name,r_sq, fill = Hydrologic_Segment)) +
  geom_col(position = position_dodge2(preserve = "single", padding = 0, width = 3)) +
  theme_bw() +
  #ggtitle("Air-Water Sensitivity Values") +
  xlab("") +
  ylab("Sensitivity r2 Value") +
  stream_labs +
  bm_theme_2

Calculate slope,r2, and n values for each site/month air water sensitivity relationship (year NOT a random intercept)

# calculate slope, r2, and n for each site from weekly mean air/water temp values
sens_model <- function(df) {
  lm(mean_h20 ~   mean_air, data = df)
}

# create values from nested tibble
sens <- meantemps %>%
  group_by(Stream_Name,Hydrologic_Segment,Site) %>%
  nest() %>%
  mutate(
    sens = map(data, sens_model),
    tidy_sens = map(sens, broom::tidy),
    sens_slope = map_dbl(tidy_sens, ~.$estimate[2]),
    sens_int = map_dbl(tidy_sens, ~.$estimate[[1]]),
    n_weeks = map_dbl(data, nrow),
    summary = map(sens,summary),
    r_sq = map_dbl(summary,"r.squared"))

# table of values
sens_vals <- sens %>%
  select(Site,r_sq,Stream_Name,Hydrologic_Segment,sens_slope,sens_int,n_weeks) %>%
  arrange(n_weeks)
sens_vals

Calculate slope,r2, and n values for each site air water sensitivity relationship (year AS a random intercept)

# calculate slope, r2, and n for each site/month from weekly mean air/water temp values
sens_model_lme <- function(df) {
  lm(mean_h20 ~ (1|year) + mean_air, data = df)
}

# create values from nested tibble
sens_lme <- meantemps %>%
  group_by(Stream_Name,Hydrologic_Segment) %>%
  nest() %>%
  mutate(
    sens = map(data, sens_model_lme),
    tidy_sens = map(sens, broom::tidy),
    sens_slope = map_dbl(tidy_sens, ~.$estimate[2]),
    sens_int = map_dbl(tidy_sens, ~.$estimate[[1]]),
    n_weeks = map_dbl(data, nrow),
    summary = map(sens,summary),
    r_sq = map_dbl(summary,"r.squared"))

# table of values
sens_vals_lme <- sens_lme %>%
  select(Stream_Name,r_sq,Hydrologic_Segment,sens_slope,sens_int,n_weeks) %>%
  arrange(n_weeks)
sens_vals_lme

NOTE: values for slope and r2 are IDENTICAL regardless of whether or not year is set as a random intercept. Not quite sure what the deal is….?!

We calculated thermal sensitivity for all sites using data from a seasonal time frame common to all locations and years, which we expressed as the slope of the stream-air temperature relationship based on weekly mean temperatures. We calculated weekly mean air/water temperatures from weeks with at least 70% of data present and took the slope of the resulting linear relationship as the site’s sensitivity.

range 20 - 40 weeks total for number of weeks

Quick exploration of predictors of r_sq values

sens_lm <- lm(r_sq ~ Stream_Name + Hydrologic_Segment, data = sens_vals)
summary(sens_lm)
## 
## Call:
## lm(formula = r_sq ~ Stream_Name + Hydrologic_Segment, data = sens_vals)
## 
## Residuals:
##         1         2         3         4         5         6         7 
## -0.230024  0.032133  0.058144 -0.055376 -0.029366 -0.002767 -0.009808 
##         8         9        10        11 
## -0.022325  0.259390 -0.181689  0.181689 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                 0.86522    0.14163   6.109  0.00170 **
## Stream_NameRussian River   -0.40932    0.16131  -2.537  0.05205 . 
## Stream_NamePtarmigan Creek -0.72687    0.16131  -4.506  0.00636 **
## Stream_NameKenai River     -0.48639    0.18627  -2.611  0.04760 * 
## Hydrologic_SegmentMiddle   -0.02153    0.13970  -0.154  0.88357   
## Hydrologic_SegmentUpper    -0.07874    0.15619  -0.504  0.63558   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1976 on 5 degrees of freedom
## Multiple R-squared:  0.8075, Adjusted R-squared:  0.615 
## F-statistic: 4.195 on 5 and 5 DF,  p-value: 0.0708
layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page 
plot(sens_lm)

require(jtools)
plot_summs(sens_lm)

Import site-specific downscaled projected air temperatures extracted from SNAP files (generated from script “raster_extract_ml_v2.Rmd”)

# read csv from local directory
projected_air <- read.csv("/Users/bmeyer/Google Drive/Thesis/R Analysis/air_water_temp_analysis/SNAP_downloads/projected_air_temps.csv") %>%
  select(-X,-raster_file) %>%
  rename(proj_air_Temp_C = Temp_C)

We can use our air-water sensitivity relationships generated from weekly means to generate weekly mean projected water temperatures

3/30/19: For purposes of comparing present water temperatures to future water temperatures, we will use air temperatures from downscaled SNAP data (monthly means) rather than our empirical air temperature.

# rename columns in sensitivity values dataframe
sens_vals <- sens_vals %>%
  rename(River = Stream_Name,
         Reach = Hydrologic_Segment)


# match air-water sensitivity values to projected air temps data frame
projected_temps <- left_join(projected_air,sens_vals) %>%
  filter(!is.na(n_weeks)) %>%
  # linear relationship using site/month specific sensitivity parameters
  mutate(proj_h2o_Temp_C = sens_slope*proj_air_Temp_C + sens_int) %>%
  # choose projected decades of interest
 mutate(decade = case_when(year >= 2030 & year <= 2039 ~"2030_2039",
                            year >= 2060 & year <= 2069 ~"2060_2069",
                            year == 2015|2016 ~ "2015_2016"))
    # alter this line of code if we end up using different decadal selection

Plot 2015-2016 temperatures juxtaposed with projected air temps by decade

## prep dfs of present-day and projected temps for join

# modify column names in projected_temps df 
projected_temps_mod <- projected_temps %>%
  select(-proj_air_Temp_C,-sens_slope,-sens_int,-r_sq,-n_weeks) %>%
  rename(mean_month_h2o = proj_h2o_Temp_C,
         Time_Period = decade,
         Stream_Name = River,
         Hydrologic_Segment = Reach) %>%
  group_by(scenario,Site,Stream_Name,Hydrologic_Segment,Time_Period,month) %>%
  summarise(mean_month_h2o = mean(mean_month_h2o))

# modify column names in empirical temps df 
temps_monthly <- temps %>%
  mutate(month = month(DateTime),
         scenario = "present",
         Time_Period = "2015_2016") %>%
  transform(Time_Period = as.character(Time_Period)) %>%
  group_by(scenario,month,Site,Stream_Name,Hydrologic_Segment,Time_Period) %>%
  # calculate monthly means of empirical temps
  summarise(mean_month_h2o = mean(water_temp_C)) 

# bind 2015/2016 monthly mean temps with projected monthly means
temps_monthly_fnl <- bind_rows(temps_monthly,projected_temps_mod) %>%
  # reorder appearance of rivers
  transform(Stream_Name = factor(Stream_Name, levels=c("Beaver Creek",
                                                       "Russian River",
                                                       "Ptarmigan Creek",
                                                       "Kenai River")))

Plot 2015-2016 water temps adjacent to projected water temps

# Create generic names for facet labels
stream_facets <- c(
  "Beaver Creek" = "Lowland", 
    "Russian River" = "Montane",
    "Ptarmigan Creek" = "Glacial", 
    "Kenai River" = "Main\nStem")

# plot A1B scenario water temps
sn <- "sresa1b" 
temps_monthly_fnl %>%
  filter(scenario %in% sn) %>%
  ggplot(aes(month,mean_month_h2o, fill = Time_Period)) +
  geom_col(position = position_dodge2(preserve = "single", padding = 0, width = 3)) +
  facet_grid(Stream_Name ~ Hydrologic_Segment, labeller = labeller(Stream_Name = stream_facets)) +
  theme_bw() +
  xlab("Month") +
  ylab("Mean Mothly Water Temp (C)") +
  theme(strip.text.y = element_text(angle = 360)) +
  bm_theme_1 +
  ggtitle("Projected Water Temperatures\nA1B Scenario (Mid-Range)")

# plot A2 scenario water temps
sn <- "sresa2"
temps_monthly_fnl %>%
  filter(scenario %in% sn) %>%
  ggplot(aes(month,mean_month_h2o, fill = Time_Period)) +
  geom_col(position = position_dodge2(preserve = "single", padding = 0, width = 3)) +
  facet_grid(Stream_Name ~ Hydrologic_Segment, labeller = labeller(Stream_Name = stream_facets)) +
  theme_bw() +
  xlab("Month") +
  ylab("Mean Mothly Water Temp (C)") +
  theme(strip.text.y = element_text(angle = 360)) +
  bm_theme_1 +
  ggtitle("Projected Water Temperatures\nA2 Scenario (Rapid Increase)")

Alternative plot style(s) for 2015-2016 vs future water temp data

# summarise temp data to group by watershed
z <- temps_monthly_fnl %>%
 # filter(scenario == "sresa1b") %>%
  group_by(Stream_Name,Time_Period,month,scenario) %>%
  summarise(mean_month_h2o = mean(mean_month_h2o)) 

# create plot object
p <- z %>%
  ggplot(aes(month,mean_month_h2o,color = Time_Period)) +
  geom_line(size = 1.3) +
  facet_grid(.~Stream_Name, labeller = labeller(Stream_Name = stream_facets)) +
  xlab("Month") +
  ylab("Mean Monthly Water Temperature (C)") +
  bm_theme_1 +
  theme(axis.text.x = element_text(size =16),
        axis.text.y = element_text(size =16),
        strip.text.x = element_text(size = 18,margin = margin(.2,0,.2,0, "cm"))) +
  scale_y_continuous(breaks=seq(8, 16, 2), limits = c(8,16)) +
 scale_color_manual(values=c("red","blue","green3"))
  

# a1b overall plot
p %+% subset(z, scenario %in% "sresa1b") 

# projections
p %+% subset(z, scenario %in% "sresa1b" & Time_Period %in% c("2015_2016")) +
  ggtitle("Mid-Range (A1B) 2015-2016")

p %+% subset(z, scenario %in% "sresa1b" & Time_Period %in% c("2015_2016","2030_2039")) +
  ggtitle("Mid-Range (A1B) 2030-2039") 

p %+% subset(z, scenario %in% "sresa1b" & Time_Period %in% c("2015_2016","2030_2039","2060_2069")) +
  ggtitle("Mid-Range (A1B) 2060-2069") 

# a2 overall plot
p %+% subset(z, scenario %in% "sresa2") 

# projections
p %+% subset(z, scenario %in% "sresa2" & Time_Period %in% c("2015_2016")) +
  ggtitle("Rapid Increase (A2)") 

p %+% subset(z, scenario %in% "sresa2" & Time_Period %in% c("2015_2016","2030_2039")) +
  ggtitle("Rapid Increase (A2)")

p %+% subset(z, scenario %in% "sresa2" & Time_Period %in% c("2015_2016","2030_2039","2060_2069")) +
  ggtitle("Rapid Increase (A2)")

Calculate table of mean % change in water temp under each scenario from 2015-2016 monthly mean

# create dataframe where re-casted present-day temps are adjacent to projected temps

# by overall watershed
pct_chg_h2o <- temps_monthly_fnl %>%
  filter(scenario != "present") %>%
  group_by(scenario,Time_Period,Stream_Name) %>%
  summarise(mean_h2o = mean(mean_month_h2o)) %>%
  spread(Time_Period,mean_h2o) %>%
  mutate(`2030_2039_pct` = ((`2030_2039` - `2015_2016`)/`2015_2016`),
         `2060_2069_pct` = ((`2060_2069` - `2015_2016`)/`2015_2016`)) %>%
  select(scenario,Stream_Name,`2030_2039_pct`,`2060_2069_pct`) %>%
  gather(years,pct,`2030_2039_pct`,`2060_2069_pct`) %>%
  mutate(years = str_replace(years, "_pct", ""),
         years = str_replace(years, "2030_2039","2030-2039"),
         years = str_replace(years, "2060_2069","2060-2069")) 

pct_chg_h2o

library(scales)
#
p <- pct_chg_h2o %>%
  #filter(scenario %in% sn) %>%
  ggplot(aes(years,pct, fill = scenario)) +
  geom_col(position = position_dodge2(preserve = "single", padding = 0, width = 3)) +
  facet_grid(.~Stream_Name, labeller = labeller(Stream_Name = stream_facets)) +
  theme_bw() +
  xlab("") +
  ylab("") +
  bm_theme_2 +
  theme(strip.text.y = element_text(angle = 360,size = 18),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(face = "bold",size = 18),
        legend.text = element_text(size = 14)) +
  # change legend
  scale_fill_discrete(name="Scenario",
                         breaks=c("sresa1b","sresa2"),
                         labels=c("Mid-Range","Rapid Increase")) +
  scale_y_continuous(labels=scales::percent_format(accuracy = 2)) 

# with legend
p

# without legend
p + theme(legend.position = "none")

We have projected mean monthly water temperatures. We want to use these as inputs for bioenergetics simulations, which means we need them structured as daily mean water temperatures.

To do: create data frame that just uses mean monthly temperature value for daily mean (constant temp), for months 05-09, decades 2030-2039, 2060-2069.

Note: for these simulations, we are keeping everything the same except water temps…

Possible startegy: export the “projected_temps” file created here, do the creation of files for simulations in “create_biol_input_FB4_v2.Rmd”

# export projected water and air temperatures to google sheet
## very large file; thus using gs_upload() rather than gs_new()
setwd("/Users/bmeyer/Google Drive/Thesis/R Analysis/air_water_temp_analysis")
write.csv(projected_temps, file = "projected_temps.csv")

# upload as gs_sheet
## remove previous version of sheet
gs_delete(gs_title("projected_temps"))
gs_upload("projected_temps",file = "/Users/bmeyer/Google Drive/Thesis/R Analysis/air_water_temp_analysis/projected_temps.csv")

The above google sheet is used and imported near the beginning of the scripts “create_biol_input_FB4_v2_YEAR_SCENARIO.Rmd”.

Anomaly values EDA

What degree of temperature anomaly exists between water temperature loggers and fish trap sites?

# A.) excise water temp logger data from durations of site vists
# B.) match probe instantaneous data from site visits
# superimpose

# import instantaneous water temp values
## 2015
trap_deploy15 <- gs_read(gs_title("2015 EPSCoR SCTC.xlsx"), ws = "B Fishing Data",skip = 2) %>%
  select(Reach,River,deploy_dt,collect_dt,Deploy_Temp,Total_Catch,Sample_Event) %>%
  mutate(Site = paste(Reach,River),
         year = "2015")

## 2016
trap_deploy16 <- gs_read(gs_title("2016_EPSCoR_Aquatic_Ecology_Database"), ws = "B Fishing Data",skip = 2) %>%
  select(Reach,River,deploy_dt,collect_dt,Deploy_Temp,Total_Catch,Sample_Event) %>%
  mutate(Site = paste(Reach,River),
         year = "2016")

# combine years
trap_deploy <- bind_rows(trap_deploy15,trap_deploy16) %>%
  transform(deploy_dt = mdy_hm(deploy_dt),
            collect_dt = mdy_hm(collect_dt)) %>%
  filter(!is.na(Deploy_Temp)) %>%
  mutate(deploy_dt = round_date(deploy_dt,"15 minutes"),
        collect_dt = round_date(collect_dt,"15 minutes")) %>%
  separate(Sample_Event, sep = "-", into = c("a","b","Sample_Event_Num"),remove = F) %>%
  select(-a,-b)
 

# select water logger temp data from temporal periods of sample events 

## calculate latest deployment and earliest retrieval time from each sample event
sample_events <- trap_deploy %>%
  #filter(Site != "Lower Russian River") %>%
  group_by(Site,year,Sample_Event,Sample_Event_Num) %>%
  summarise(event_start = min(deploy_dt),
            event_end = max(collect_dt))

## restrict air/water temperature data to time periods calculated above, by site
logger_temps_subset  <- temps %>% 
    rowwise() %>% 
    mutate(present = any()) %>%
    mutate(present = any(DateTime >= sample_events$event_start & DateTime <= sample_events$event_end & Site == sample_events$Site)) %>% 
    filter(present) %>% 
    data.frame() %>%
    select(-present) 

# match instantaneous handheld probe water temp data to logger water temp data from nearest site
## prep data frame for join
trap_deploy <- trap_deploy %>%
  select(deploy_dt,Deploy_Temp,Site,Sample_Event_Num) %>%
  rename(DateTime = deploy_dt) 


# left join instantaneous temps with logger temps
anomaly_temps <- left_join(logger_temps_subset,trap_deploy,by = c("Site","DateTime")) %>%
  fill(Sample_Event_Num, .direction = "down") %>%
  fill(Sample_Event_Num, .direction = "up")

rm(trap_deploy15,trap_deploy16)

# example plots of anomaly values

# lower beaver creek 2016
anomaly_temps %>%
  mutate(doy = yday(DateTime)) %>%
  filter(Site == "Lower Beaver Creek",
        # Sample_Event_Num == c("2","3"),
         year == "2016"
         ) %>%
  ggplot(aes(DateTime,water_temp_C)) +
  geom_point() +
  geom_jitter(aes(DateTime,Deploy_Temp, color = "red"), width = 0.3) +
  facet_grid(.~Sample_Event_Num, scales = "free_x") +
  theme_bw() +
  bm_theme +
  ggtitle("2016 Lower Beaver Creek\nWater Temperatures")

  
  
# upper russian river
anomaly_temps %>%
  mutate(doy = yday(DateTime)) %>%
  filter(Site == "Upper Russian River",
        # Sample_Event_Num == c("2","3"),
         year == "2016"
         ) %>%
  ggplot(aes(DateTime,water_temp_C)) +
  geom_point() +
  geom_jitter(aes(DateTime,Deploy_Temp, color = "red"),width = 0.1) +
  facet_grid(.~Sample_Event_Num, scales = "free_x") +
  theme_bw() +
  bm_theme +
  ggtitle("2016 Upper Russian River\nWater Temperatures")

# middle kenai river
anomaly_temps %>%
  mutate(doy = yday(DateTime)) %>%
  filter(Site == "Lower Ptarmigan Creek",
        # Sample_Event_Num == c("2","3"),
         year == "2016"
         ) %>%
  ggplot(aes(DateTime,water_temp_C)) +
  geom_point() +
  geom_jitter(aes(DateTime,Deploy_Temp, color = "red"),width = 0.1) +
  facet_grid(.~Sample_Event_Num, scales = "free_x") +
  theme_bw() +
  bm_theme +
  ggtitle("2016 Lower Ptarmigan Creek\nWater Temperatures")