Purpose: join redd count and environmental covariate data and explore correlations between key environmental variables.
Load YCT redd count data and plot time series
redds <- read_csv("Data/ReddCounts_WGFD_1971-2021_cleaned.csv") %>%
filter(!stream %in% c("Cody", "Christiansen", "Dave", "Laker")) # drop Cody (too few data points) and Salt River tribs
redds %>% ggplot() + geom_point(aes(x = year, y = reddspermile)) + geom_line(aes(x = year, y = reddspermile)) + facet_wrap(~ stream, scales = "free_y") + theme_bw() + theme(panel.grid = element_blank())
Map sites:
reddlocs <- redds %>% group_by(stream) %>% summarize(lat = unique(lat), long = unique(long)) # summarize location data
reddlocs$stream_short <- c("THCH", "BLKT", "BLCR", "COCA", "FISH", "FLAT", "LTBC", "LOBC", "NOWL", "PRCE", "SRSC", "SPRG", "UPBC")
reddlocs_sp <- st_as_sf(reddlocs, coords = c("long", "lat"), crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84")
#st_write(reddlocs_sp, "Figures/Map/map_data/ReddCounts_SpatialLocations.shp", append = FALSE)
mapview(reddlocs_sp) # mapExplore autocorrelation, which may be due to stock-recruit dynamics, autocorrelated environmental drivers, or both
sites <- unique(redds$stream)
par(mfrow = c(4,4), mgp = c(2.5,1,0), mar = c(3.5,3.5,1,1))
for (i in 1:length(sites)) {
d <- filter(redds, stream == sites[i])
acf(d$reddsperkm, lag.max = 6, plot = TRUE, na.action = na.pass)
legend("topright", legend = sites[i], bty = "n")
}
For each site, add years in which data was not collected, then calculate lagged redds
lagg <- 4
rlist <- list()
sites <- unique(redds$stream)
for (i in 1:length(sites)) {
d <- filter(redds, stream == sites[i])
df1 <- tibble(year = seq(from = min(d$year, na.rm = TRUE), to = max(d$year, na.rm = T), by = 1))
rlist[[i]] <- df1 %>% left_join(d) %>% mutate(stream = sites[i],
miles = unique(d$miles),
kms = unique(d$kms),
lat = unique(d$lat),
long = unique(d$long),
reddsperkm_lag4 = lag(reddsperkm, n = lagg))
}
redds2 <- bind_rows(rlist)Explore the effect of spawners on productivity to motivate use of Ricker model, sensu Jones et al. (2020).
# recruits ~ spawners, combined
# redds2 %>% ggplot(aes(x = reddsperkm_lag4, y = reddsperkm)) +
# geom_point(aes(colour = stream)) +
# geom_smooth(method = "loess")
# recruits ~ spawners, by stream
# redds2 %>% ggplot(aes(x = reddsperkm_lag4, y = reddsperkm)) +
# geom_point(aes(colour = year)) +
# geom_smooth() +
# facet_wrap(~stream, scales = "free")
# # productivity ~ time, by stream
# redds2 %>% ggplot(aes(x = year, y = log(reddsperkm/reddsperkm_lag4))) +
# geom_point(aes(colour = year)) +
# geom_smooth() +
# facet_wrap(~stream, scales = "free_y")
# productivity ~ spawning stock, combined
# redds2 %>% ggplot(aes(x = reddsperkm_lag4, y = log(reddsperkm/reddsperkm_lag4))) +
# geom_smooth(method = "lm") +
# geom_point(aes(colour = stream))
# productivity ~ spawning stock, by stream
redds2 %>% ggplot(aes(x = reddsperkm_lag4, y = log(reddsperkm/reddsperkm_lag4), colour = year)) +
geom_smooth(method = "lm") +
geom_point(aes(colour = year)) +
facet_wrap(~stream, scales = "free")
Get range of years by site
redds2 %>% group_by(stream) %>% summarize(start = min(year, na.rm = TRUE), stop = max(year, na.rm = TRUE))# A tibble: 13 × 3
stream start stop
<chr> <dbl> <dbl>
1 3 Channel 1973 2014
2 Blacktail 1973 2016
3 Blue Crane 1985 2015
4 Cowboy Cabin 1980 2016
5 Fish 1984 2016
6 Flat 1976 2015
7 Little Bar BC 1978 2016
8 Lower Bar BC 1971 2021
9 Nowlin 1976 2015
10 Price 1975 2006
11 Snake River Side Channel 1986 2002
12 Spring 1988 2007
13 Upper Bar BC 1974 2016minby <- min(redds2$year) - 4 # minimum brood yearPer Murdoch et al (2023) CJFAS, only include variables for which Pearson’s r is < 0.6.
Notes: * drop jld_winmean and natq_winmean from same model (pearson r = 0.639 since 1960, 0.79 since 1989)
Load derived data files
jldramp <- read_csv("Data/Derived/JLD_RampDown_Summary_CalendarYear_1960-2022.csv") #%>% mutate(broodyr = year) %>% select(-year)
jldflow <- read_csv("Data/Derived/JLD_ManagedFlow_Covariates_BroodYear_1960-2022.csv")
natflow <- read_csv("Data/Derived/SnakeTribs_NaturalFlow_Covariates_BroodYear_1960-2022.csv") %>% rename(flowstn = site) #%>% filter(site %in% c("SnakeBeFlat", "Salt"))
airtemp <- read_csv("Data/Derived/AirTemperature_Covariates_BroodYear_1960-2022.csv") %>% filter(site %in% c("moose", "afton")) %>% rename(tempstn = site)
unique(natflow$flowstn)[1] "Buffalo" "Fish" "Flat" "Greys" "GrosVentre"
[6] "Salt" "SnakeBeFlat" "SnakeMoran" "SnakeNat" Use pairs plots to check correlation among key variables, within datasets:
ggpairs(jldramp %>% select(jld_rampdur, jld_rampratemindoy, jld_rampratemin_log), progress = FALSE ) # entire time series
ggpairs(jldramp %>% filter(broodyr >= 1989) %>% select(jld_rampdur, jld_rampratemindoy, jld_rampratemin_log) , progress = FALSE) # post-1989
ggpairs(jldflow %>% select(jld_winmean, jld_winmin, jld_winvar, jld_summean, jld_peakmag, jld_peaktime), progress = FALSE)
ggpairs(jldflow %>% filter(broodyr >= 1989) %>% select(jld_winmean, jld_winvar, jld_summean, jld_peakmag, jld_peaktime), progress = FALSE)
ggpairs(natflow %>% filter(flowstn == "SnakeNat") %>% select(natq_annmin, natq_winmean, natq_winvar, natq_peakmag, natq_peaktime), progress = FALSE)
ggpairs(natflow %>% filter(broodyr >= 1989, flowstn == "SnakeNat") %>% select(natq_winmean, natq_winvar, natq_peakmag, natq_peaktime), progress = FALSE)
ggpairs(natflow %>% filter(flowstn == "SnakeBeFlat") %>% select(natq_winmean, natq_winvar, natq_peakmag, natq_peaktime), progress = FALSE)
ggpairs(natflow %>% filter(broodyr >= 1989, flowstn == "SnakeBeFlat") %>% select(natq_winmean, natq_winvar, natq_peakmag, natq_peaktime), progress = FALSE)
ggpairs(airtemp %>% filter(tempstn == "moose") %>% select(temp_falmean, temp_winmean, temp_sprmean, temp_summean), progress = FALSE)
ggpairs(airtemp %>% filter(broodyr >= 1989, tempstn == "moose") %>% select(temp_falmean, temp_winmean, temp_sprmean, temp_summean), progress = FALSE)
Check correlation among variables across datasets:
p <- jldramp %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log) %>%
left_join(jldflow %>% select(broodyr, jld_winmin, jld_summean, jld_peakmag, jld_peaktime)) %>%
left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peakmag, natq_peaktime)) %>%
#left_join(natflow %>% filter(flowstn == "SnakeBeFlat") %>% select(broodyr, natq_winmean)) %>%
left_join(airtemp %>% filter(tempstn == "moose") %>% select(broodyr, temp_falmean, temp_winmean, temp_sprmean, temp_summean))
range(cor(p, use = "pairwise.complete.obs")[upper.tri(cor(p, use = "pairwise.complete.obs"))])[1] -0.5061091 0.5862950ggpairs(p, progress = FALSE)
Post 1989
p <- jldramp %>% filter(broodyr >= 1989) %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log) %>%
left_join(jldflow %>% filter(broodyr >= 1989) %>% select(broodyr, jld_winmin, jld_summean, jld_peakmag, jld_peaktime)) %>%
left_join(natflow %>% filter(broodyr >= 1989) %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peakmag, natq_peaktime)) %>%
#left_join(natflow %>% filter(broodyr >= 1989) %>% filter(flowstn == "SnakeBeFlat") %>% select(broodyr, natq_winmean)) %>%
left_join(airtemp %>% filter(broodyr >= 1989) %>% filter(tempstn == "moose") %>% select(broodyr, temp_falmean, temp_winmean, temp_sprmean, temp_summean))
range(cor(p, use = "pairwise.complete.obs")[upper.tri(cor(p, use = "pairwise.complete.obs"))])[1] -0.5481930 0.5453191ggpairs(p, progress = FALSE)
Experienced flow
p <- jldramp %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log) %>%
left_join(jldflow %>% select(broodyr, jld_winvar, jld_summean)) %>%
#left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peakmag, natq_peaktime)) %>%
left_join(natflow %>% filter(flowstn == "SnakeBeFlat") %>% select(broodyr, natq_winmean, natq_peakmag, natq_peaktime)) %>%
left_join(airtemp %>% filter(tempstn == "moose") %>% select(broodyr, temp_falmean, temp_winmean, temp_sprmean, temp_summean)) %>% select(-broodyr)
range(cor(p, use = "complete.obs")[upper.tri(cor(p, use = "complete.obs"))])[1] -0.5976980 0.6003041ggpairs(p, progress = FALSE)
JLD and natural winter flow highly correlated, especially after 1989
ddd <- jldflow %>% select(broodyr, jld_winmean) %>% left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_winmean))
cor.test(ddd$jld_winmean, ddd$natq_winmean)
Pearson's product-moment correlation
data: ddd$jld_winmean and ddd$natq_winmean
t = 5.6337, df = 46, p-value = 1.022e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4335332 0.7812567
sample estimates:
cor
0.6389597 ddd1 <- ddd %>% filter(broodyr >= 1989)
cor.test(ddd1$jld_winmean, ddd1$natq_winmean)
Pearson's product-moment correlation
data: ddd1$jld_winmean and ddd1$natq_winmean
t = 7.227, df = 32, p-value = 3.29e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6123541 0.8889125
sample estimates:
cor
0.7874557 ddd1 <- ddd %>% filter(broodyr < 1989)
cor.test(ddd1$jld_winmean, ddd1$natq_winmean)
Pearson's product-moment correlation
data: ddd1$jld_winmean and ddd1$natq_winmean
t = 1.563, df = 12, p-value = 0.144
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1526044 0.7731439
sample estimates:
cor
0.4112756 While the timing of peak natural vs. managed flows are not correlated during the full or either partial time period, this may be driven by years of missing data for JLD peak flows…where there was not a JLD peak. Time series plots are suspicious
# peak timing...ok...
ddd <- jldflow %>% select(broodyr, jld_peaktime) %>% left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peaktime))
cor.test(ddd$jld_peaktime, ddd$natq_peaktime)
Pearson's product-moment correlation
data: ddd$jld_peaktime and ddd$natq_peaktime
t = -0.43671, df = 36, p-value = 0.6649
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3833791 0.2529610
sample estimates:
cor
-0.07259273 plot(ddd$jld_peaktime, ddd$natq_peaktime, ylab = "brood year DOY", xlab = "year")
plot(jld_peaktime ~ broodyr, ddd, type = "b", col = "blue")
lines(natq_peaktime ~ broodyr, ddd, type = "b", col = "red")
legend("bottomright", legend = c("JLD", "natural"), lty = 1, col = c("blue", "red"))
ddd1 <- ddd %>% filter(broodyr >= 1989)
cor.test(ddd1$jld_peaktime, ddd1$natq_peaktime)
Pearson's product-moment correlation
data: ddd1$jld_peaktime and ddd1$natq_peaktime
t = -0.62341, df = 23, p-value = 0.5392
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4986378 0.2805140
sample estimates:
cor
-0.1289046 #plot(ddd1$jld_peaktime, ddd1$natq_peaktime)
#abline(a = 0, b = 1)
ddd1 <- ddd %>% filter(broodyr < 1989)
cor.test(ddd1$jld_peaktime, ddd1$natq_peaktime)
Pearson's product-moment correlation
data: ddd1$jld_peaktime and ddd1$natq_peaktime
t = 0.013513, df = 11, p-value = 0.9895
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5481414 0.5538164
sample estimates:
cor
0.004074438 #plot(ddd1$jld_peaktime, ddd1$natq_peaktime)No multicollinearity issues between the magnitude of peak natural vs. managed flows, although these variables are somewhat correlated prior to 1989.
# peak magnitude...ok...
ddd <- jldflow %>% select(broodyr, jld_peakmag) %>% left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peakmag))
cor.test(ddd$jld_peakmag, ddd$natq_peakmag)
Pearson's product-moment correlation
data: ddd$jld_peakmag and ddd$natq_peakmag
t = 4.9086, df = 46, p-value = 1.194e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3625576 0.7461337
sample estimates:
cor
0.586295 plot(jld_peakmag ~ broodyr, ddd, type = "b", col = "blue", ylim = c(3000,20000))
lines(natq_peakmag ~ broodyr, ddd, type = "b", col = "red")
legend("topleft", legend = c("JLD", "natural"), lty = 1, col = c("blue", "red"))
ddd1 <- ddd %>% filter(broodyr >= 1989)
cor.test(ddd1$jld_peakmag, ddd1$natq_peakmag)
Pearson's product-moment correlation
data: ddd1$jld_peakmag and ddd1$natq_peakmag
t = 3.6311, df = 32, p-value = 0.0009749
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2471686 0.7426824
sample estimates:
cor
0.5401836 ddd1 <- ddd %>% filter(broodyr < 1989)
cor.test(ddd1$jld_peakmag, ddd1$natq_peakmag)
Pearson's product-moment correlation
data: ddd1$jld_peakmag and ddd1$natq_peakmag
t = 3.2601, df = 12, p-value = 0.006828
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2431875 0.8916793
sample estimates:
cor
0.6853377 JLD peak flow magnitude and combined winter mean flow are correlated (>0.6)…perhaps an artifact of “good water years”
ddd <- jldflow %>% select(broodyr, jld_peakmag, jld_winmin)
cor.test(ddd$jld_peakmag, ddd$jld_winmin)
Pearson's product-moment correlation
data: ddd$jld_peakmag and ddd$jld_winmin
t = 2.2111, df = 61, p-value = 0.03079
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.02640982 0.48727185
sample estimates:
cor
0.2723925 plot(jld_peakmag ~ broodyr, ddd, type = "b", col = "blue")
par(new = TRUE)
plot(jld_winmin ~ broodyr, ddd, type = "b", col = "red")
legend("topleft", legend = c("JLD peak mag", "Combined winter mean"), lty = 1, col = c("blue", "red"))
ddd1 <- ddd %>% filter(broodyr >= 1989)
cor.test(ddd1$jld_peakmag, ddd1$jld_winmin)
Pearson's product-moment correlation
data: ddd1$jld_peakmag and ddd1$jld_winmin
t = 2.4703, df = 32, p-value = 0.01902
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.07173556 0.65034644
sample estimates:
cor
0.4001932 ddd1 <- ddd %>% filter(broodyr < 1989)
cor.test(ddd1$jld_peakmag, ddd1$jld_winmin)
Pearson's product-moment correlation
data: ddd1$jld_peakmag and ddd1$jld_winmin
t = 1.9257, df = 27, p-value = 0.06473
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.02176591 0.63335086
sample estimates:
cor
0.347512 Check correlations between variables in year t and t+1 (to separate age-0 and age-1 effects).
Winter variation is the only variable that is moderately correlated between time t and t+1 (0.65)
cor.test(jldramp$jld_rampdur, lead(jldramp$jld_rampdur))
Pearson's product-moment correlation
data: jldramp$jld_rampdur and lead(jldramp$jld_rampdur)
t = 0.96669, df = 60, p-value = 0.3376
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1299488 0.3623978
sample estimates:
cor
0.1238389 cor.test(jldramp$jld_rampratemindoy, lead(jldramp$jld_rampratemindoy))
Pearson's product-moment correlation
data: jldramp$jld_rampratemindoy and lead(jldramp$jld_rampratemindoy)
t = 4.2439, df = 60, p-value = 7.735e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2621889 0.6520108
sample estimates:
cor
0.4804914 cor.test(jldramp$jld_rampratemin_log, lead(jldramp$jld_rampratemin_log))
Pearson's product-moment correlation
data: jldramp$jld_rampratemin_log and lead(jldramp$jld_rampratemin_log)
t = 0.70935, df = 60, p-value = 0.4809
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1622692 0.3333699
sample estimates:
cor
0.09119508 cor.test(jldflow$jld_winmin, lead(jldflow$jld_winmin)) # pearsons r = 0.65
Pearson's product-moment correlation
data: jldflow$jld_winmin and lead(jldflow$jld_winmin)
t = 4.3865, df = 60, p-value = 4.729e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.2771030 0.6611604
sample estimates:
cor
0.4927661 cor.test(jldflow$jld_summean, lead(jldflow$jld_summean))
Pearson's product-moment correlation
data: jldflow$jld_summean and lead(jldflow$jld_summean)
t = 0.27338, df = 60, p-value = 0.7855
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2164039 0.2825498
sample estimates:
cor
0.03527073 cor.test(jldflow$jld_peakmag, lead(jldflow$jld_peakmag)) # pearsons r = 0.65
Pearson's product-moment correlation
data: jldflow$jld_peakmag and lead(jldflow$jld_peakmag)
t = 2.5116, df = 60, p-value = 0.01473
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.06357118 0.51828219
sample estimates:
cor
0.3084419 cor.test(jldflow$jld_peaktime, lead(jldflow$jld_peaktime))
Pearson's product-moment correlation
data: jldflow$jld_peaktime and lead(jldflow$jld_peaktime)
t = 1.0467, df = 44, p-value = 0.3009
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1407981 0.4268542
sample estimates:
cor
0.1558709 d <- natflow %>% filter(flowstn == "SnakeNat")
cor.test(d$natq_peakmag, lead(d$natq_peakmag))
Pearson's product-moment correlation
data: d$natq_peakmag and lead(d$natq_peakmag)
t = 0.17576, df = 45, p-value = 0.8613
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2629533 0.3110188
sample estimates:
cor
0.02619124 cor.test(d$natq_peaktime, lead(d$natq_peaktime))
Pearson's product-moment correlation
data: d$natq_peaktime and lead(d$natq_peaktime)
t = 0.85953, df = 45, p-value = 0.3946
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1661387 0.3996718
sample estimates:
cor
0.1270916 d <- airtemp %>% filter(tempstn == "moose")
cor.test(d$temp_falmean, lead(d$temp_falmean))
Pearson's product-moment correlation
data: d$temp_falmean and lead(d$temp_falmean)
t = 2.5647, df = 56, p-value = 0.01303
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.07193934 0.53749544
sample estimates:
cor
0.3242108 cor.test(d$temp_winmean, lead(d$temp_winmean))
Pearson's product-moment correlation
data: d$temp_winmean and lead(d$temp_winmean)
t = 0.82595, df = 58, p-value = 0.4122
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1502169 0.3521051
sample estimates:
cor
0.1078205 cor.test(d$temp_sprmean, lead(d$temp_sprmean))
Pearson's product-moment correlation
data: d$temp_sprmean and lead(d$temp_sprmean)
t = 0.96896, df = 52, p-value = 0.337
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1395630 0.3871302
sample estimates:
cor
0.1331742 cor.test(d$temp_summean, lead(d$temp_summean))
Pearson's product-moment correlation
data: d$temp_summean and lead(d$temp_summean)
t = 1.3246, df = 53, p-value = 0.191
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.09059664 0.42415727
sample estimates:
cor
0.1790029 # correlations among natural, managed, and combined winter flow
ddd <- jldflow %>% select(broodyr, jld_annmin) %>%
left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_annmin)) %>%
left_join(natflow %>% filter(flowstn == "SnakeBeFlat") %>% select(broodyr, natq_annmin) %>% rename(expq_annmin = natq_annmin))
ggpairs(ddd)
ggpairs(ddd %>% filter(broodyr >= 1989))
ddd <- jldflow %>% select(broodyr, jld_peakmag) %>%
left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peakmag)) %>%
left_join(natflow %>% filter(flowstn == "SnakeBeFlat") %>% select(broodyr, natq_peakmag) %>% rename(expq_peakmag = natq_peakmag))
ggpairs(ddd)
ggpairs(ddd %>% filter(broodyr >= 1989))
ddd <- jldflow %>% select(broodyr, jld_peaktime) %>%
left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_peaktime)) %>%
left_join(natflow %>% filter(flowstn == "SnakeBeFlat") %>% select(broodyr, natq_peaktime) %>% rename(expq_peaktime = natq_peaktime))
ggpairs(ddd)
ggpairs(ddd %>% filter(broodyr >= 1989))
# winter flow time series plots
ddd <- jldflow %>% select(broodyr, jld_annmin) %>%
left_join(natflow %>% filter(flowstn == "SnakeNat") %>% select(broodyr, natq_annmin)) %>%
mutate(tot_annmin = jld_annmin + natq_annmin,
prop_annmin = jld_annmin/tot_annmin) %>% filter(!is.na(tot_annmin))
ddd$cor <- NA
ddd$cor[5:44] <- rollapply(ddd, width = 9, function(x) abs(cor(x[,2], x[,3])), by.column = FALSE)
par(mfrow = c(3,1), mar = c(3,4,0.5,0.5), mgp = c(2.5,1,0))
plot(tot_annmin ~ broodyr, ddd, type = "n", ylim = c(0,2500), xlab = "", ylab = "Winter flow (cfs)", bty = "l")
polygon(x = c(ddd$broodyr, rev(ddd$broodyr)), y = c(rep(0, times = dim(ddd)[1]), rev(ddd$tot_annmin)), border = NA, col = "grey")
polygon(x = c(ddd$broodyr, rev(ddd$broodyr)), y = c(rep(0, times = dim(ddd)[1]), rev(ddd$natq_annmin)), border = NA, col = "white")
polygon(x = c(ddd$broodyr, rev(ddd$broodyr)), y = c(rep(0, times = dim(ddd)[1]), rev(ddd$jld_annmin)), border = NA, col = "white")
polygon(x = c(ddd$broodyr, rev(ddd$broodyr)), y = c(rep(0, times = dim(ddd)[1]), rev(ddd$natq_annmin)), border = NA, col = scales::alpha("darkorchid2", 1))
polygon(x = c(ddd$broodyr, rev(ddd$broodyr)), y = c(rep(0, times = dim(ddd)[1]), rev(ddd$jld_annmin)), border = NA, col = scales::alpha("goldenrod1", 1))
legend("topright", legend = c("Total", "Natural", "Managed"), fill = c("grey", "darkorchid2", "goldenrod1"), bty = "n")
plot(prop_annmin ~ broodyr, ddd, type = "l", bty = "l", xlab = "", ylab = "JLD contribution to winter flow")
points(prop_annmin ~ broodyr, ddd, pch = 16)
plot(cor ~ broodyr, ddd, type = "l", bty = "l", xlab = "", ylab = "Correlation (rolling 9-year window)", ylim = c(0,1))
points(cor ~ broodyr, ddd, pch = 16)
abline(h = 0.6, lty = 2, col = "red")Pull out relevant variables, trim to temporal extent of redd count data (by brood year. I.e., redd counts start at 1971 but environmental data starts at 1967), and standardize (mean centered and scaled).
# Managed flow variables
covsrc_jldq <- jldramp %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log) %>%
filter(broodyr >= minby) %>%
left_join(jldflow %>% select(broodyr, jld_winmean, jld_winmin, jld_winvar, jld_summean, jld_peakmag, jld_peaktime)) %>%
#mutate(jld_winmin_log = ifelse(jld_winmin_log < 4, NA, jld_winmin_log)) %>%
mutate(z_jld_rampdur = scale(jld_rampdur, center = TRUE, scale = TRUE)[,1],
z_jld_rampratemindoy = scale(jld_rampratemindoy, center = TRUE, scale = TRUE)[,1],
z_jld_rampratemin_log = scale(jld_rampratemin_log, center = TRUE, scale = TRUE)[,1],
z_jld_winmean = scale(jld_winmean, center = TRUE, scale = TRUE)[,1],
z_jld_winmin = scale(jld_winmin, center = TRUE, scale = TRUE)[,1],
z_jld_winvar = scale(jld_winvar, center = TRUE, scale = TRUE)[,1],
z_jld_summean = scale(jld_summean, center = TRUE, scale = TRUE)[,1],
z_jld_peakmag = scale(jld_peakmag, center = TRUE, scale = TRUE)[,1],
z_jld_peaktime = scale(jld_peaktime, center = TRUE, scale = TRUE)[,1])
#ggpairs(covsrc_jldq %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log, jld_winmean, jld_winmin, jld_winvar, jld_summean, jld_peakmag, jld_peaktime))
# Managed flow variables - ramp-down only
covsrc_jldq_rd <- jldramp %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log) %>%
filter(broodyr >= minby) %>%
#mutate(jld_rampratemin = ifelse(jld_rampratemin < -2000, NA, jld_rampratemin)) %>%
mutate(z_jld_rampdur = scale(jld_rampdur, center = TRUE, scale = TRUE)[,1],
z_jld_rampratemindoy = scale(jld_rampratemindoy, center = TRUE, scale = TRUE)[,1],
z_jld_rampratemin_log = scale(jld_rampratemin_log, center = TRUE, scale = TRUE)[,1])
#ggpairs(covsrc_jldq_rd %>% select(broodyr, jld_rampdur, jld_rampratemindoy, jld_rampratemin_log))
# Natural/wild flow variables
covsrc_natq <- natflow %>% filter(flowstn %in% c("SnakeNat")) %>% # make sure flow station is correct!
select(flowstn, broodyr, natq_falmean_log, natq_winmean_log, natq_sprmean_log, natq_summean_log, natq_winvar, natq_peakmag, natq_peaktime) %>% #rename(flowstn = site) %>%
filter(broodyr >= minby) %>%
mutate(z_natq_winvar = scale(natq_winvar, center = TRUE, scale = TRUE)[,1]) %>%
mutate(z_natq_falmean_log = scale(natq_falmean_log, center = TRUE, scale = TRUE)[,1],
z_natq_winmean_log = scale(natq_winmean_log, center = TRUE, scale = TRUE)[,1],
z_natq_sprmean_log = scale(natq_sprmean_log, center = TRUE, scale = TRUE)[,1],
z_natq_summean_log = scale(natq_summean_log, center = TRUE, scale = TRUE)[,1],
z_natq_peakmag = scale(natq_peakmag, center = TRUE, scale = TRUE)[,1],
z_natq_peaktime = scale(natq_peaktime, center = TRUE, scale = TRUE)[,1]) %>% ungroup()
#ggpairs(covsrc_natq %>% select(broodyr, natq_falmean_log, natq_winmean_log, natq_sprmean_log, natq_summean_log, natq_winvar, natq_peakmag, natq_peaktime))
# Experienced flow variables
covsrc_expq <- natflow %>% filter(flowstn %in% c("SnakeBeFlat")) %>%
select(flowstn, broodyr, natq_winmean, natq_winvar, natq_summean, natq_peakmag, natq_peaktime) %>%
filter(broodyr >= minby) %>%
mutate(z_natq_winmean = scale(natq_winmean, center = TRUE, scale = TRUE)[,1],
z_natq_winvar = scale(natq_winvar, center = TRUE, scale = TRUE)[,1],
z_natq_summean = scale(natq_summean, center = TRUE, scale = TRUE)[,1],
z_natq_peakmag = scale(natq_peakmag, center = TRUE, scale = TRUE)[,1],
z_natq_peaktime = scale(natq_peaktime, center = TRUE, scale = TRUE)[,1])
#ggpairs(covsrc_expq %>% select(broodyr, natq_winmean, natq_winvar, natq_summean, natq_peakmag, natq_peaktime))
# Temperature variables
covsrc_temp <- airtemp %>% filter(tempstn %in% c("moose")) %>%
select(tempstn, broodyr, temp_falmean, temp_winmean, temp_sprmean, temp_summean) %>% #rename(tempstn = site) %>% #group_by(tempstn) %>%
filter(broodyr >= minby) %>%
mutate(z_temp_falmean = scale(temp_falmean, center = TRUE, scale = TRUE)[,1],
z_temp_winmean = scale(temp_winmean, center = TRUE, scale = TRUE)[,1],
z_temp_sprmean = scale(temp_sprmean, center = TRUE, scale = TRUE)[,1],
z_temp_summean = scale(temp_summean, center = TRUE, scale = TRUE)[,1]) %>% ungroup()
#ggpairs(covsrc_temp %>% select(broodyr, temp_falmean, temp_winmean, temp_sprmean, temp_summean))Create and write out summary tables (means and standard deviations)
# JLD Flow
covsrc_jldq_summary <- tibble(cov = c("z_jld_rampdur", "z_jld_rampratemindoy", "z_jld_rampratemin_log", "z_jld_winmean", "z_jld_winmin", "z_jld_winvar", "z_jld_summean", "z_jld_peakmag", "z_jld_peaktime"),
mean = c(attr(scale(covsrc_jldq$jld_rampdur, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_rampratemindoy, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_rampratemin_log, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_winmean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_winmin, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_winvar, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_summean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_peakmag, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_jldq$jld_peaktime, center = TRUE, scale = TRUE), "scaled:center")),
sd = c(attr(scale(covsrc_jldq$jld_rampdur, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_rampratemindoy, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_rampratemin_log, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_winmean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_winmin, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_winvar, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_summean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_peakmag, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_jldq$jld_peaktime, center = TRUE, scale = TRUE), "scaled:scale")),
min_raw = c(min(covsrc_jldq$jld_rampdur, na.rm = TRUE),
min(covsrc_jldq$jld_rampratemindoy, na.rm = TRUE),
min(covsrc_jldq$jld_rampratemin_log, na.rm = TRUE),
min(covsrc_jldq$jld_winmean, na.rm = TRUE),
min(covsrc_jldq$jld_winmin, na.rm = TRUE),
min(covsrc_jldq$jld_winvar, na.rm = TRUE),
min(covsrc_jldq$jld_summean, na.rm = TRUE),
min(covsrc_jldq$jld_peakmag, na.rm = TRUE),
min(covsrc_jldq$jld_peaktime, na.rm = TRUE)),
max_raw = c(max(covsrc_jldq$jld_rampdur, na.rm = TRUE),
max(covsrc_jldq$jld_rampratemindoy, na.rm = TRUE),
max(covsrc_jldq$jld_rampratemin_log, na.rm = TRUE),
max(covsrc_jldq$jld_winmean, na.rm = TRUE),
max(covsrc_jldq$jld_winmin, na.rm = TRUE),
max(covsrc_jldq$jld_winvar, na.rm = TRUE),
max(covsrc_jldq$jld_summean, na.rm = TRUE),
max(covsrc_jldq$jld_peakmag, na.rm = TRUE),
max(covsrc_jldq$jld_peaktime, na.rm = TRUE)),
min_zscore = c(min(covsrc_jldq$z_jld_rampdur, na.rm = TRUE),
min(covsrc_jldq$z_jld_rampratemindoy, na.rm = TRUE),
min(covsrc_jldq$z_jld_rampratemin_log, na.rm = TRUE),
min(covsrc_jldq$z_jld_winmean, na.rm = TRUE),
min(covsrc_jldq$z_jld_winmin, na.rm = TRUE),
min(covsrc_jldq$z_jld_winvar, na.rm = TRUE),
min(covsrc_jldq$z_jld_summean, na.rm = TRUE),
min(covsrc_jldq$z_jld_peakmag, na.rm = TRUE),
min(covsrc_jldq$z_jld_peaktime, na.rm = TRUE)),
max_zscore = c(max(covsrc_jldq$z_jld_rampdur, na.rm = TRUE),
max(covsrc_jldq$z_jld_rampratemindoy, na.rm = TRUE),
max(covsrc_jldq$z_jld_rampratemin_log, na.rm = TRUE),
max(covsrc_jldq$z_jld_winmean, na.rm = TRUE),
max(covsrc_jldq$z_jld_winmin, na.rm = TRUE),
max(covsrc_jldq$z_jld_winvar, na.rm = TRUE),
max(covsrc_jldq$z_jld_summean, na.rm = TRUE),
max(covsrc_jldq$z_jld_peakmag, na.rm = TRUE),
max(covsrc_jldq$z_jld_peaktime, na.rm = TRUE)))
write_csv(covsrc_jldq_summary, "Data/Derived/ManagedFlow_SummaryMeanSD_1967-2022.csv")
# Natural Flow
covsrc_natq_summary <- tibble(cov = c("z_natq_falmean_log", "z_natq_winmean_log", "z_natq_sprmean_log", "z_natq_summean_log", "z_natq_winvar", "z_natq_peakmag", "z_natq_peaktime"),
mean = c(attr(scale(covsrc_natq$natq_falmean_log, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_natq$natq_winmean_log, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_natq$natq_sprmean_log, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_natq$natq_summean_log, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_natq$natq_winvar, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_natq$natq_peakmag, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_natq$natq_peaktime, center = TRUE, scale = TRUE), "scaled:center")),
sd = c(attr(scale(covsrc_natq$natq_falmean_log, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_natq$natq_winmean_log, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_natq$natq_sprmean_log, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_natq$natq_summean_log, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_natq$natq_winvar, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_natq$natq_peakmag, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_natq$natq_peaktime, center = TRUE, scale = TRUE), "scaled:scale")),
min_raw = c(min(covsrc_natq$natq_falmean_log, na.rm = TRUE),
min(covsrc_natq$natq_winmean_log, na.rm = TRUE),
min(covsrc_natq$natq_sprmean_log, na.rm = TRUE),
min(covsrc_natq$natq_summean_log, na.rm = TRUE),
min(covsrc_natq$natq_winvar, na.rm = TRUE),
min(covsrc_natq$natq_peakmag, na.rm = TRUE),
min(covsrc_natq$natq_peaktime, na.rm = TRUE)),
max_raw = c(max(covsrc_natq$natq_falmean_log, na.rm = TRUE),
max(covsrc_natq$natq_winmean_log, na.rm = TRUE),
max(covsrc_natq$natq_sprmean_log, na.rm = TRUE),
max(covsrc_natq$natq_summean_log, na.rm = TRUE),
max(covsrc_natq$natq_winvar, na.rm = TRUE),
max(covsrc_natq$natq_peakmag, na.rm = TRUE),
max(covsrc_natq$natq_peaktime, na.rm = TRUE)),
min_zscore = c(min(covsrc_natq$z_natq_falmean_log, na.rm = TRUE),
min(covsrc_natq$z_natq_winmean_log, na.rm = TRUE),
min(covsrc_natq$z_natq_sprmean_log, na.rm = TRUE),
min(covsrc_natq$z_natq_summean_log, na.rm = TRUE),
min(covsrc_natq$z_natq_winvar, na.rm = TRUE),
min(covsrc_natq$z_natq_peakmag, na.rm = TRUE),
min(covsrc_natq$z_natq_peaktime, na.rm = TRUE)),
max_zscore = c(max(covsrc_natq$z_natq_falmean_log, na.rm = TRUE),
max(covsrc_natq$z_natq_winmean_log, na.rm = TRUE),
max(covsrc_natq$z_natq_sprmean_log, na.rm = TRUE),
max(covsrc_natq$z_natq_summean_log, na.rm = TRUE),
max(covsrc_natq$z_natq_winvar, na.rm = TRUE),
max(covsrc_natq$z_natq_peakmag, na.rm = TRUE),
max(covsrc_natq$z_natq_peaktime, na.rm = TRUE)))
write_csv(covsrc_natq_summary, "Data/Derived/NaturalFlow_SummaryMeanSD_1975-2022.csv")
# Experienced Flow
covsrc_expq_summary <- tibble(cov = c("z_expq_winmean", "z_expq_winvar", "z_expq_summean", "z_expq_peakmag", "z_expq_peaktime"),
mean = c(attr(scale(covsrc_expq$natq_winmean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_expq$natq_winvar, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_expq$natq_summean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_expq$natq_peakmag, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_expq$natq_peaktime, center = TRUE, scale = TRUE), "scaled:center")),
sd = c(attr(scale(covsrc_expq$natq_winmean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_expq$natq_winvar, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_expq$natq_summean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_expq$natq_peakmag, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_expq$natq_peaktime, center = TRUE, scale = TRUE), "scaled:scale")),
min_raw = c(min(covsrc_expq$natq_winmean, na.rm = TRUE),
min(covsrc_expq$natq_winvar, na.rm = TRUE),
min(covsrc_expq$natq_summean, na.rm = TRUE),
min(covsrc_expq$natq_peakmag, na.rm = TRUE),
min(covsrc_expq$natq_peaktime, na.rm = TRUE)),
max_raw = c(max(covsrc_expq$natq_winmean, na.rm = TRUE),
max(covsrc_expq$natq_winvar, na.rm = TRUE),
max(covsrc_expq$natq_summean, na.rm = TRUE),
max(covsrc_expq$natq_peakmag, na.rm = TRUE),
max(covsrc_expq$natq_peaktime, na.rm = TRUE)),
min_zscore = c(min(covsrc_expq$z_natq_winmean, na.rm = TRUE),
min(covsrc_expq$z_natq_winvar, na.rm = TRUE),
min(covsrc_expq$z_natq_summean, na.rm = TRUE),
min(covsrc_expq$z_natq_peakmag, na.rm = TRUE),
min(covsrc_expq$z_natq_peaktime, na.rm = TRUE)),
max_zscore = c(max(covsrc_expq$z_natq_winmean, na.rm = TRUE),
max(covsrc_expq$z_natq_winvar, na.rm = TRUE),
max(covsrc_expq$z_natq_summean, na.rm = TRUE),
max(covsrc_expq$z_natq_peakmag, na.rm = TRUE),
max(covsrc_expq$z_natq_peaktime, na.rm = TRUE)))
write_csv(covsrc_expq_summary, "Data/Derived/ExperiencedFlow_SummaryMeanSD_1975-2022.csv")
# Temperature
covsrc_temp_summary <- tibble(cov = c("z_temp_falmean", "z_temp_winmean", "z_temp_sprmean", "z_temp_summean"),
mean = c(attr(scale(covsrc_temp$temp_falmean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_temp$temp_winmean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_temp$temp_sprmean, center = TRUE, scale = TRUE), "scaled:center"),
attr(scale(covsrc_temp$temp_summean, center = TRUE, scale = TRUE), "scaled:center")),
sd = c(attr(scale(covsrc_temp$temp_falmean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_temp$temp_winmean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_temp$temp_sprmean, center = TRUE, scale = TRUE), "scaled:scale"),
attr(scale(covsrc_temp$temp_summean, center = TRUE, scale = TRUE), "scaled:scale")),
min_raw = c(min(covsrc_temp$temp_falmean, na.rm = TRUE),
min(covsrc_temp$temp_winmean, na.rm = TRUE),
min(covsrc_temp$temp_sprmean, na.rm = TRUE),
min(covsrc_temp$temp_summean, na.rm = TRUE)),
max_raw = c(max(covsrc_temp$temp_falmean, na.rm = TRUE),
max(covsrc_temp$temp_winmean, na.rm = TRUE),
max(covsrc_temp$temp_sprmean, na.rm = TRUE),
max(covsrc_temp$temp_summean, na.rm = TRUE)),
min_zscore = c(min(covsrc_temp$z_temp_falmean, na.rm = TRUE),
min(covsrc_temp$z_temp_winmean, na.rm = TRUE),
min(covsrc_temp$z_temp_sprmean, na.rm = TRUE),
min(covsrc_temp$z_temp_summean, na.rm = TRUE)),
max_zscore = c(max(covsrc_temp$z_temp_falmean, na.rm = TRUE),
max(covsrc_temp$z_temp_winmean, na.rm = TRUE),
max(covsrc_temp$z_temp_sprmean, na.rm = TRUE),
max(covsrc_temp$z_temp_summean, na.rm = TRUE)))
write_csv(covsrc_temp_summary, "Data/Derived/Temperature_SummaryMeanSD_1967-2022.csv")Age-0 (4-year lag)
redds3a <- redds2
redds3a$broodyr <- redds3a$year - 4 # create broodyr variable for joining
redds3a$tempstn <- ifelse(redds3a$stream %in% c("Christiansen", "Dave", "Laker"), "afton", "moose") # basin ID to match to NWS temp station
redds3a$flowstn <- ifelse(redds3a$stream %in% c("Christiansen", "Dave", "Laker"), "Salt", "SnakeNat") # basin ID to match to USGS flow station
redds3a <- redds3a %>%
left_join(covsrc_jldq, by = "broodyr") %>%
left_join(covsrc_natq, by = c("broodyr", "flowstn")) %>%
left_join(covsrc_temp, by = c("broodyr", "tempstn"))Age-1 (3-year lag)
redds3b <- redds2
redds3b$broodyr <- redds3b$year - 3 # create broodyr variable for joining
redds3b$tempstn <- ifelse(redds3b$stream %in% c("Christiansen", "Dave", "Laker"), "afton", "moose") # basin ID to match to NWS temp station
redds3b$flowstn <- ifelse(redds3b$stream %in% c("Christiansen", "Dave", "Laker"), "Salt", "SnakeNat") # basin ID to match to USGS flow station
redds3b <- redds3b %>%
left_join(covsrc_jldq, by = "broodyr") %>%
left_join(covsrc_natq, by = c("broodyr", "flowstn")) %>%
left_join(covsrc_temp, by = c("broodyr", "tempstn"))Age-0 (4-year lag)
redds3c <- redds2
redds3c$broodyr <- redds3c$year - 4 # create broodyr variable for joining
redds3c$tempstn <- ifelse(redds3c$stream %in% c("Christiansen", "Dave", "Laker"), "afton", "moose") # basin ID to match to NWS temp station
redds3c$flowstn <- ifelse(redds3c$stream %in% c("Christiansen", "Dave", "Laker"), "Salt", "SnakeBeFlat") # basin ID to match to USGS flow station
redds3c <- redds3c %>%
left_join(covsrc_jldq_rd, by = "broodyr") %>%
left_join(covsrc_expq, by = c("broodyr", "flowstn")) %>%
left_join(covsrc_temp, by = c("broodyr", "tempstn"))Age-1 (3-year lag)
redds3d <- redds2
redds3d$broodyr <- redds3d$year - 3 # create broodyr variable for joining
redds3d$tempstn <- ifelse(redds3d$stream %in% c("Christiansen", "Dave", "Laker"), "afton", "moose") # basin ID to match to NWS temp station
redds3d$flowstn <- ifelse(redds3d$stream %in% c("Christiansen", "Dave", "Laker"), "Salt", "SnakeBeFlat") # basin ID to match to USGS flow station
redds3d <- redds3d %>%
left_join(covsrc_jldq_rd, by = "broodyr") %>%
left_join(covsrc_expq, by = c("broodyr", "flowstn")) %>%
left_join(covsrc_temp, by = c("broodyr", "tempstn"))write_csv(redds3a, "Data/Derived/ReddCounts_WGFD_1971-2021_cleaned_withFlowTempCovariates_SeparateFlowComponents_age0.csv")
write_csv(redds3b, "Data/Derived/ReddCounts_WGFD_1971-2021_cleaned_withFlowTempCovariates_SeparateFlowComponents_age1.csv")
write_csv(redds3c, "Data/Derived/ReddCounts_WGFD_1971-2021_cleaned_withFlowTempCovariates_ExperiencedFlow_age0.csv")
write_csv(redds3d, "Data/Derived/ReddCounts_WGFD_1971-2021_cleaned_withFlowTempCovariates_ExperiencedFlow_age1.csv")plot(natq_winmean ~ broodyr, covsrc_expq, type = "b")
plot(natq_winvar ~ broodyr, covsrc_expq, type = "b")
plot(natq_peakmag ~ broodyr, covsrc_expq, type = "b")
plot(natq_peaktime ~ broodyr, covsrc_expq, type = "b")