2 Explore Data

Purpose: Load EcoDrought data and explore visually

2.1 Site information

Code

siteinfo <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_SiteInformation.csv")
datatable(siteinfo)

Code

siteinfo_sp <- st_as_sf(siteinfo, coords = c("long", "lat"), crs = 4326)
mapview(siteinfo_sp, zcol = "designation")

2.2 Flow and temp data

Code

dat <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_FlowTempData_DailyWeekly.csv")
dat

# A tibble: 202,014 × 31
   station_no site_name      site_id basin   subbasin region   lat  long elev_ft
   <chr>      <chr>          <chr>   <chr>   <chr>    <chr>  <dbl> <dbl>   <dbl>
 1 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 2 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 3 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 4 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 5 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 6 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 7 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 8 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 9 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
10 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
# ℹ 202,004 more rows
# ℹ 22 more variables: area_sqmi <dbl>, designation <chr>, date <date>,
#   DischargeReliability <dbl>, TempReliability <dbl>, flow_mean <dbl>,
#   flow_min <dbl>, flow_max <dbl>, tempc_mean <dbl>, tempc_min <dbl>,
#   tempc_max <dbl>, flow_mean_filled <dbl>, flow_mean_cms <dbl>,
#   flow_mean_filled_cms <dbl>, area_sqkm <dbl>, Yield_mm <dbl>,
#   Yield_filled_mm <dbl>, flow_mean_7 <dbl>, flow_mean_filled_7 <dbl>, …

Unique basins

Code

unique(dat$basin)

[1] "Flathead"       "West Brook"     "Donner Blitzen" "Paine Run"     
[5] "Piney River"    "Staunton River" "Shields River"  "Snake River"

2.2.1 View daily temp

2.2.2 View daily flow

2.3 Daymet data

Trim siteinfo to sites with temp/flow data

Code

siteinfo_sub <- siteinfo %>% filter(site_name %in% unique(dat$site_name))

Download Daymet daily air temperature data for each site (point locations)

Code

climlist <- vector("list", length = dim(siteinfo_sub)[1])
for (i in 1:dim(siteinfo_sub)[1]) {
  clim <- download_daymet(site = siteinfo_sub$site_name[i], lat = siteinfo_sub$lat[i], lon = siteinfo_sub$long[i], start = 2010, end = 2024, internal = T)
  climlist[[i]] <- tibble(clim$data) %>% 
    mutate(air_temp_mean = (tmax..deg.c. + tmin..deg.c.)/2, 
           date = as.Date(paste(year, yday, sep = "-"), "%Y-%j"),
           site_name = siteinfo_sub$site_name[i]) %>%
    select(12,2,11,10,4,6,3,5) %>% rename(precip_mmday = 5, swe_kgm2 = 6, daylength_sec = 7, shortrad_wm2 = 8)
  print(i)
}
climdf <- do.call(rbind, climlist)
write_csv(climdf, "data/Daymet_daily.csv")

Bind climate data to EcoD flow-temp data

Code

climdf <- read_csv("data/Daymet_daily.csv")
dat <- dat %>% left_join(climdf)

Plot an example air (orange) and water (blue) temperature time series

Code

dat %>% filter(site_name == "McGeeCreekTrib") %>% #select(date, tempc_mean, air_temp_mean) %>% 
  ggplot() + 
  geom_line(aes(x = date, y = air_temp_mean), color = "darkorange") + 
  geom_line(aes(x = date, y = tempc_mean), color = "blue") +
  theme_bw() + ylab("temperature (deg C)")

2.4 Breakpoints

Define periods of air-water temperature synchronization per Letcher et al. (2016, PeerJ)

Show how this works for a single site: McGeeCreekTrib

Code

sub <- dat %>% 
  filter(site_name == "NF Spread Creek Upper") %>%
  mutate(index = ifelse(tempc_mean > 0, (tempc_mean - air_temp_mean) / tempc_mean, NA),
         year = year(date))

# sub_indexci <- sub %>% 
#   filter(yday >= 125, yday <= 275) %>% 
#   group_by(site_name) %>% 
#   summarize(index_low = quantile(index, probs = 0.001, na.rm = TRUE),
#             index_high = quantile(index, probs = 0.999, na.rm = TRUE)) %>%
#   ungroup()

sub_indexci <- sub %>% 
  filter(yday >= 125, yday <= 275) %>% 
  group_by(site_name, year) %>% 
  summarize(index_low = t.test(index, conf.level = 0.999)$conf.int[1],
            index_high = t.test(index, conf.level = 0.999)$conf.int[2]) %>%
  ungroup()

# 
# sub %>% 
#   ggplot() + 
#   geom_point(aes(x = date, y = index)) + 
#   theme_bw() + ylab("temperature index") + 
#   ylim(-20,20)
#   # geom_hline(yintercept = c(-1,1), color = "red")

sub %>% 
  ggplot() + 
  geom_point(aes(x = yday, y = index)) + 
  theme_bw() + ylab("temperature index") + 
  ylim(-20,20) +
  # geom_hline(yintercept = c(-1,1), color = "red") +
  geom_vline(xintercept = c(125,275), color = "blue") +
  geom_hline(data = sub_indexci, aes(yintercept = index_low), color = "red") +
  geom_hline(data = sub_indexci, aes(yintercept = index_high), color = "red") +
  facet_wrap(~year)

Code

dat

# A tibble: 202,014 × 37
   station_no site_name      site_id basin   subbasin region   lat  long elev_ft
   <chr>      <chr>          <chr>   <chr>   <chr>    <chr>  <dbl> <dbl>   <dbl>
 1 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 2 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 3 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 4 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 5 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 6 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 7 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 8 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
 9 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
10 12355347   Big Creek NWIS BIG     Flathe… Big Cre… Flat    48.6 -114.   3528.
# ℹ 202,004 more rows
# ℹ 28 more variables: area_sqmi <dbl>, designation <chr>, date <date>,
#   DischargeReliability <dbl>, TempReliability <dbl>, flow_mean <dbl>,
#   flow_min <dbl>, flow_max <dbl>, tempc_mean <dbl>, tempc_min <dbl>,
#   tempc_max <dbl>, flow_mean_filled <dbl>, flow_mean_cms <dbl>,
#   flow_mean_filled_cms <dbl>, area_sqkm <dbl>, Yield_mm <dbl>,
#   Yield_filled_mm <dbl>, flow_mean_7 <dbl>, flow_mean_filled_7 <dbl>, …

--- title: "Explore Data" format: html --- **Purpose:** Load EcoDrought data and explore visually ```{r include=FALSE} library(tidyverse) library(dygraphs) library(mapview) library(DT) library(sf) library(daymetr) ``` ## Site information ```{r} siteinfo <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_SiteInformation.csv") datatable(siteinfo) siteinfo_sp <- st_as_sf(siteinfo, coords = c("long", "lat"), crs = 4326) mapview(siteinfo_sp, zcol = "designation") ``` ## Flow and temp data ```{r} dat <- read_csv("C:/Users/jbaldock/OneDrive - DOI/Documents/USGS/EcoDrought/EcoDrought Working/Data/EcoDrought_FlowTempData_DailyWeekly.csv") dat ``` Unique basins ```{r} unique(dat$basin) ``` ### View daily temp ::: panel-tabset #### West Brook ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "West Brook") %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` #### Paine Run ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Paine Run", date >= as_date("2018-01-01")) %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` #### Staunton River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Staunton River", date >= as_date("2018-01-01")) %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` #### Flathead River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Flathead") %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` #### Yellowstone River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Shields River") %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` #### Snake River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Snake River") %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` #### Donner und Blitzen River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Donner Blitzen") %>% select(date, site_name, tempc_mean) %>% spread(key = site_name, value = tempc_mean) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "Daily mean temperature (deg. C)") ``` ::: ### View daily flow ::: panel-tabset #### West Brook ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "West Brook") %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` #### Paine Run ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Paine Run", date >= as_date("2018-01-01")) %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` #### Staunton River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Staunton River", date >= as_date("2018-01-01")) %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` #### Flathead River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Flathead") %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` #### Yellowstone River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Shields River") %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` #### Snake River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Snake River") %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` #### Donner und Blitzen River ```{r, echo=FALSE} dat %>% filter(designation %in% c("little", "medium"), basin == "Donner Blitzen") %>% mutate(logYield = log(Yield_filled_mm)) %>% select(date, site_name, logYield) %>% spread(key = site_name, value = logYield) %>% dygraph() %>% dyRangeSelector() %>% dyAxis("y", label = "ln(Specific discharge, mm)") ``` ::: ## Daymet data Trim siteinfo to sites with temp/flow data ```{r} siteinfo_sub <- siteinfo %>% filter(site_name %in% unique(dat$site_name)) ``` Download Daymet daily air temperature data for each site (point locations) ```{r eval=FALSE} climlist <- vector("list", length = dim(siteinfo_sub)[1]) for (i in 1:dim(siteinfo_sub)[1]) { clim <- download_daymet(site = siteinfo_sub$site_name[i], lat = siteinfo_sub$lat[i], lon = siteinfo_sub$long[i], start = 2010, end = 2024, internal = T) climlist[[i]] <- tibble(clim$data) %>% mutate(air_temp_mean = (tmax..deg.c. + tmin..deg.c.)/2, date = as.Date(paste(year, yday, sep = "-"), "%Y-%j"), site_name = siteinfo_sub$site_name[i]) %>% select(12,2,11,10,4,6,3,5) %>% rename(precip_mmday = 5, swe_kgm2 = 6, daylength_sec = 7, shortrad_wm2 = 8) print(i) } climdf <- do.call(rbind, climlist) write_csv(climdf, "data/Daymet_daily.csv") ``` Bind climate data to EcoD flow-temp data ```{r} climdf <- read_csv("data/Daymet_daily.csv") dat <- dat %>% left_join(climdf) ``` Plot an example air (orange) and water (blue) temperature time series ```{r} dat %>% filter(site_name == "McGeeCreekTrib") %>% #select(date, tempc_mean, air_temp_mean) %>% ggplot() + geom_line(aes(x = date, y = air_temp_mean), color = "darkorange") + geom_line(aes(x = date, y = tempc_mean), color = "blue") + theme_bw() + ylab("temperature (deg C)") ``` ## Breakpoints Define periods of air-water temperature synchronization per Letcher et al. (2016, PeerJ) Show how this works for a single site: *McGeeCreekTrib* ```{r} sub <- dat %>% filter(site_name == "NF Spread Creek Upper") %>% mutate(index = ifelse(tempc_mean > 0, (tempc_mean - air_temp_mean) / tempc_mean, NA), year = year(date)) # sub_indexci <- sub %>% # filter(yday >= 125, yday <= 275) %>% # group_by(site_name) %>% # summarize(index_low = quantile(index, probs = 0.001, na.rm = TRUE), # index_high = quantile(index, probs = 0.999, na.rm = TRUE)) %>% # ungroup() sub_indexci <- sub %>% filter(yday >= 125, yday <= 275) %>% group_by(site_name, year) %>% summarize(index_low = t.test(index, conf.level = 0.999)$conf.int[1], index_high = t.test(index, conf.level = 0.999)$conf.int[2]) %>% ungroup() # # sub %>% # ggplot() + # geom_point(aes(x = date, y = index)) + # theme_bw() + ylab("temperature index") + # ylim(-20,20) # # geom_hline(yintercept = c(-1,1), color = "red") sub %>% ggplot() + geom_point(aes(x = yday, y = index)) + theme_bw() + ylab("temperature index") + ylim(-20,20) + # geom_hline(yintercept = c(-1,1), color = "red") + geom_vline(xintercept = c(125,275), color = "blue") + geom_hline(data = sub_indexci, aes(yintercept = index_low), color = "red") + geom_hline(data = sub_indexci, aes(yintercept = index_high), color = "red") + facet_wrap(~year) ``` ```{r} dat ```