-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathHYOS_Hydropeaking.R
126 lines (113 loc) · 5.92 KB
/
HYOS_Hydropeaking.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
############################
# Hydropeaking Sensitivity HYOS
#############################
source("1spFunctions.R")
source("HYOS_1sp.R")
# read in LF temp and discharge data from 2007 to 2023
temp <- readNWISdv("09380000", "00010", "2007-10-01", "2023-05-01")
discharge <- readNWISdv("09380000", "00060", "2007-10-01", "2023-05-01")
# calculate average yearly flows
flow <- average.yearly.flows(flowdata = discharge, flow.column_name = "X_00060_00003", date.column_name = "Date" )
# calculate average yearly temperatures
temps <- average.yearly.temp(tempdata = temp, temp.column_name = "X_00010_00003", date.column_name = "Date")
# create a timeseries of average temperatures 100 years long
temps <- rep.avg.year(temps, n = 100, change.in.temp = 0, years.at.temp = 0)
# create a timeseries of average flows 100 years long
flows <- do.call("rbind", replicate(100, flow, simplify = FALSE))
# match dates
flows$dts <- as.Date(temps$dts)
# get discharge magnitude by dividing by bankfull discharge
flows$Discharge <- flows$Discharge/85000
# create sequence of hydropeaking intensities
hydropeak <- seq(0.0, 0.7, by = 0.05)
# makes some vectors for data to go into
means <- vector()
sd <- vector()
sizemeans <- vector()
sizesd <- vector()
S3Yrprod <- vector()
# cycle though hydropeaking scenarios
for (hyd in 1:length(hydropeak)){
set.seed(123) # make reproducible
# model sizes
sizes <- HYOSmodel(flow.data = flows$Discharge, temp.data = temps, baselineK = 10000, disturbanceK = 9000 , Qmin = 0.3, extinct = 50, iteration = 2, peaklist = hydropeak[hyd], peakeach = length(temps$Temperature), stage_output = "size")
# model abundances
out <- HYOSmodel(flow.data = flows$Discharge, temp.data = temps, baselineK = 10000, disturbanceK = 9000 , Qmin = 0.3, extinct = 50, iteration = 2, peaklist = hydropeak[hyd], peakeach = length(temps$Temperature))
# for each stage, calculate mean biomass
s1s <- mean(out[-c(1:260), 1, ]) * (0.0046 * (mean(sizes[-c(1:260)]))^2.926)
s2s <- mean(out[-c(1:260), 2, ]) * (0.0046 * (mean(sizes[-c(1:260)]+3))^2.926)
s3s <- mean(out[-c(1:260), 3, ]) * (0.0046 * (mean(sizes[-c(1:260)]+3))^2.926)
# sum the mean biomass of each stage to get mean timestep biomass
sizemeans[hyd] <- sum(c(s1s, s2s, s3s))
# produce standard devation
sizesd[hyd] <- sd(c(s1s, s2s, s3s))
# now instead of getting the mean, calculate biomass at every timestep
s1ss <- as.data.frame(cbind(rowMeans(out[-c(1:260), 1, ]) * (0.0046 * ((sizes[-c(1:260)]))^2.926), year(temps$dts[-c(1:259)])))
s2ss <- as.data.frame(cbind(rowMeans(out[-c(1:260), 2, ]) * (0.0046 * ((sizes[-c(1:260)]+3))^2.926), year(temps$dts[-c(1:259)])))
s3ss <- as.data.frame(cbind(rowMeans(out[-c(1:260), 3, ]) * (0.0046 * ((sizes[-c(1:260)]+3))^2.926), year(temps$dts[-c(1:259)])))
#
# find annual stage specific biomass based on year
s1sYr <- aggregate(V1 ~ V2, data = s1ss, FUN = sum, na.rm = TRUE)
s2sYr <- aggregate(V1 ~ V2, data = s2ss, FUN = sum, na.rm = TRUE)
s3sYr <- aggregate(V1 ~ V2, data = s3ss, FUN = sum, na.rm = TRUE)
# add all stages together to get average annual biomass (aka secondary production)
#Yrprod[hyd] <- sum(mean(c(s1sYr$V1, s2sYr$V1, s3sYr$V1), na.rm = T))
# average annual biomass of only stage 3 (emergent adults)
S3Yrprod[hyd] <- mean( s3sYr$V1, na.rm = T)
# calculate mean abundances at each timestep
means.list.HYOS <- mean.data.frame(out, burnin = 260, iteration = 2)
# calculate the average of mean abundances at each hydropeaking intensity
means[hyd] <- mean(means.list.HYOS$mean.abund)
# calculate the standard deviation of mean abundances at each hydropeaking intensity
sd[hyd] <- sd(means.list.HYOS$mean.abund, na.rm = T)
}
# compile abundance data
HYOS_hyd_means <- as.data.frame(cbind(hydropeak, means, sd, rep("HYOS", length(means))))
# compile timeseries biomass data
HYOS_hyd_size <- as.data.frame(cbind(hydropeak, sizemeans, sizesd, rep("HYOS", length(sizemeans))))
# compile timeseries biomass data
HYOS_hyd_yrprod <- as.data.frame(cbind(hydropeak, S3Yrprod, rep("HYOS", length(sizemeans))))
#
# ggplot(data = HYOS_hyd_means, aes(x = hydropeak, y = means, group = 1, color = "HYOS")) +
# geom_ribbon(aes(ymin = means - sd,
# ymax = means + sd),
# colour = 'transparent',
# alpha = .15,
# show.legend = T) +
# geom_point(show.legend = T, linewidth = 1, alpha = 0.8) +
# #geom_line(data = flow.magnitude, aes(x = as.Date(dts), y = X_00060_00003), color = "blue") +
# #geom_line(data = temps, aes(x = as.Date(dts), y = Temperature*1000), color = "green")+
# #coord_cartesian(ylim = c(0,6000)) +S1 and S2 (inds/m2)
# ylab('HYOS spp. Abund.') +
# xlab("")+
# labs(colour=" ")+
# theme_bw()+
# scale_color_manual(values = colors)+
# #scale_y_continuous(
# # sec.axis = sec_axis(~., name="HYOSidae Larvae (inds/m2)"
# # ))+
# theme(text = element_text(size = 13), axis.text.x = element_text(angle=45, hjust = 1, size = 12.5),
# axis.text.y = element_text(size = 13), )
#
#
# ggplot(data = HYOS_hyd_size, aes(x = hydropeak, y = sizemeans, group = 1, color = "HYOS")) +
# geom_ribbon(aes(ymin = sizemeans - sizesd,
# ymax = sizemeans + sizesd),
# colour = 'transparent',
# alpha = .15,
# show.legend = T) +
# geom_line(show.legend = T, linewidth = 1, alpha = 0.8) +
# #geom_line(data = flow.magnitude, aes(x = as.Date(dts), y = X_00060_00003), color = "blue") +
# #geom_line(data = temps, aes(x = as.Date(dts), y = Temperature*1000), color = "green")+
# #coord_cartesian(ylim = c(0,6000)) +S1 and S2 (inds/m2)
# ylab('HYOS spp. Biomass (mg)') +
# xlab("")+
# labs(colour=" ")+
# theme_bw()+
# scale_color_manual(values = colors)+
# #scale_y_continuous(
# # sec.axis = sec_axis(~., name="HYOSidae Larvae (inds/m2)"
# # ))+
# theme(text = element_text(size = 13), axis.text.x = element_text(angle=45, hjust = 1, size = 12.5),
# axis.text.y = element_text(size = 13), )
#