library(tidyr)
try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))
detects_6 <- detects_study %>% filter(general_location == "Benicia_west" | general_location == "Benicia_east")
if(nrow(detects_6) == 0){
plot(1:2, type = "n", xlab = "",xaxt = "n", yaxt = "n", ylab = "Number of fish arrivals per day")
text(1.5,1.5, labels = "NO DETECTIONS YET", cex = 2)
} else {
detects_6 <- detects_6 %>%
dplyr::left_join(., detects_6 %>%
group_by(TagCode) %>%
summarise(first_detect = min(DateTime_PST))) %>%
mutate(Day = as.Date(as.Date(first_detect, "Etc/GMT+8")))
starttime <- as.Date(min(detects_6$release_time), "Etc/GMT+8")
# Endtime should be either now, or end of predicted tag life, whichever comes first
endtime <- min(as.Date(format(Sys.time(), "%Y-%m-%d")),
max(as.Date(detects_study$release_time)+(as.numeric(detects_study$tag_life))))
daterange <- data.frame(Day = seq.Date(from = starttime, to = endtime, by = "day"))
rels <- unique(study_tagcodes$Release)
rel_num <- length(rels)
rels_no_detects <- as.character(rels[!(rels %in% unique(detects_6$Release))])
tagcount1 <- detects_6 %>%
group_by(Day, Release) %>%
summarise(unique_tags = length(unique(TagCode))) %>%
spread(Release, unique_tags)
daterange1 <- merge(daterange, tagcount1, all.x=T)
daterange1[is.na(daterange1)] <- 0
if(length(rels_no_detects)>0){
for(i in rels_no_detects){
daterange1 <- cbind(daterange1, x=NA)
names(daterange1)[names(daterange1) == "x"] <- paste(i)
}
}
## reorder columns in alphabetical order so its coloring in barplots is consistent
daterange1 <- daterange1[,order(colnames(daterange1))]
daterange2 <- daterange1
rownames(daterange2) <- daterange2$Day
daterange2$Day <- NULL
par(mar=c(6, 5, 2, 5) + 0.1)
daterange2$Date <- as.Date(row.names(daterange2))
daterange3 <- melt(daterange2, id.vars = "Date", variable.name = ".", )
daterange3$. <- factor(daterange3$., levels = sort(unique(daterange3$.), decreasing = T))
plot_ly(daterange3, width = 900, height = 600, dynamicTicks = TRUE) %>%
add_bars(x = ~Date, y = ~value, color = ~.) %>%
add_annotations(text="Release (click on legend<br> items to isolate)", xref="paper", yref="paper",
x=0.01, xanchor="left",
y=1.02, yanchor="bottom", # Same y as legend below
legendtitle=TRUE, showarrow=FALSE ) %>%
layout(yaxis2 = ay,showlegend = T,
barmode = "stack",
xaxis = list(title = "Date", mirror=T,ticks="outside",showline=T),
yaxis = list(title = "Number of fish arrivals per day", mirror=T,ticks="outside",showline=T),
legend = list(orientation = "h",x = 0.34, y = 1.01, yanchor="bottom", xanchor="left"),
margin=list(l = 50,r = 100,b = 50,t = 50))
}
2.6 Detections at Benicia Bridge for duration of tag life
3. Survival and Routing Probability
try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))
detects_tower <- detects_study %>% filter(general_location == "TowerBridge")
if(nrow(detects_tower) == 0){
WR.surv <- data.frame("Release"=NA, "Survival (%)"="NO DETECTIONS YET", "SE"=NA, "95% lower C.I."=NA,
"95% upper C.I."=NA, "Detection efficiency (%)"=NA)
colnames(WR.surv) <- c("Release", "Survival (%)", "SE", "95% lower C.I.",
"95% upper C.I.", "Detection efficiency (%)")
print(kable(WR.surv, row.names = F, "html", caption = "3.1 Minimum survival to Tower Bridge (using CJS
survival model). If Yolo Bypass Weirs are overtopping during migration, fish may have taken
that route, and therefore this is a minimum estimate of survival") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"),
full_width = F, position = "left"))
} else {
study_count <- nrow(study_tagcodes)
# Only do survival to Sac for now
surv <- detects_study %>% filter(river_km > 168 & river_km < 175)
# calculate mean and SD travel time
travel <- aggregate(list(first_detect = surv$DateTime_PST), by = list(Release = surv$Release, TagCode = surv$TagCode, RelDT = surv$RelDT), min)
travel$days <- as.numeric(difftime(travel$first_detect, travel$RelDT, units = "days"))
travel_final <- aggregate(list(mean_travel_time = travel$days), by = list(Release = travel$Release), mean)
travel_final <- merge(travel_final, aggregate(list(sd_travel_time = travel$days), by = list(Release = travel$Release), sd))
travel_final <- merge(travel_final, aggregate(list(n = travel$days), by = list(Release = travel$Release), length))
travel_final <- rbind(travel_final, data.frame(Release = "ALL", mean_travel_time = mean(travel$days), sd_travel_time = sd(travel$days),n = nrow(travel)))
# Create inp for survival estimation
inp <- as.data.frame(reshape2::dcast(surv, TagCode ~ river_km, fun.aggregate = length))
# Sort columns by river km in descending order
gen_loc_sites <- ncol(inp)-1 # Count number of genlocs
if(gen_loc_sites < 2){
WR.surv <- data.frame("Release"=NA, "Survival (%)"="NOT ENOUGH DETECTIONS", "SE"=NA, "95% lower C.I."=NA,
"95% upper C.I."=NA, "Detection efficiency (%)"=NA)
colnames(WR.surv) <- c("Release", "Survival (%)", "SE", "95% lower C.I.", "95% upper C.I.",
"Detection efficiency (%)")
print(kable(WR.surv, row.names = F, "html", caption = "3.1 Minimum survival to Tower Bridge (using CJS
survival model). If Yolo Bypass Weirs are overtopping during migration, fish may
have taken that route, and therefore this is a minimum estimate of survival") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"),
full_width = F,position = "left"))
} else {
inp <- inp[,c(1,order(names(inp[,2:(gen_loc_sites+1)]), decreasing = T)+1)] %>%
dplyr::left_join(study_tagcodes, ., by = "TagCode")
inp2 <- inp[,(ncol(inp)-gen_loc_sites+1):ncol(inp)] %>%
replace(is.na(.), 0) %>%
replace(., . > 0, 1)
inp <- cbind(inp, inp2)
groups <- as.character(sort(unique(inp$Release)))
surv$Release <- factor(surv$Release, levels = groups)
inp[,groups] <- 0
for (i in groups) {
inp[as.character(inp$Release) == i, i] <- 1
}
inp$inp_final <- paste("1",apply(inp2, 1, paste, collapse=""),sep="")
if(length(groups) > 1){
# make sure factor levels have a release that has detections first. if first release in factor order
# has zero detectins, model goes haywire
inp.df <- data.frame(ch = as.character(inp$inp_final), freq = 1,
rel = factor(inp$Release, levels = names(sort(table(surv$Release),decreasing = T))),
stringsAsFactors = F)
WR.process <- process.data(inp.df, model="CJS", begin.time=1, groups = "rel")
WR.ddl <- make.design.data(WR.process)
WR.mark.all <- mark(WR.process, WR.ddl,
model.parameters=list(Phi=list(formula=~time),p=list(formula=~time)),
silent = T, output = F)
WR.mark.rel <- mark(WR.process, WR.ddl,
model.parameters=list(Phi=list(formula=~time*rel),p=list(formula=~time)),
silent = T, output = F)
WR.surv <- round(WR.mark.all$results$real[1,c("estimate", "se", "lcl", "ucl")] * 100,1)
WR.surv <- rbind(WR.surv, round(WR.mark.rel$results$real[seq(from=1,to=length(groups)*2,by = 2),
c("estimate", "se", "lcl", "ucl")] * 100,1))
WR.surv$Detection_efficiency <- NA
WR.surv[1,"Detection_efficiency"] <- round(WR.mark.all$results$real[gen_loc_sites+1,"estimate"] * 100,1)
WR.surv <- cbind(c("ALL", names(sort(table(surv$Release),decreasing = T))), WR.surv)
}
if(length(intersect(colnames(inp),groups)) < 2){
inp$inp_final <- paste("1",apply(inp2, 1, paste, collapse=""), " ", 1,sep = "")
write.table(inp$inp_final,"WRinp.inp",row.names = F, col.names = F, quote = F)
WRinp <- convert.inp("WRinp.inp")
WR.process <- process.data(WRinp, model="CJS", begin.time=1)
WR.ddl <- make.design.data(WR.process)
WR.mark.all <- mark(WR.process, WR.ddl,
model.parameters=list(Phi=list(formula=~time),p=list(formula=~time)),
silent = T, output = F)
WR.mark.rel <- mark(WR.process, WR.ddl,
model.parameters=list(Phi=list(formula=~time),p=list(formula=~time)),
silent = T, output = F)
WR.surv <- round(WR.mark.all$results$real[1,c("estimate", "se", "lcl", "ucl")] * 100,1)
WR.surv <- rbind(WR.surv, round(WR.mark.rel$results$real[seq(from=1,to=length(groups)*2,by = 2),
c("estimate", "se", "lcl", "ucl")] * 100,1))
WR.surv$Detection_efficiency <- NA
WR.surv[1,"Detection_efficiency"] <- round(WR.mark.all$results$real[gen_loc_sites+1,"estimate"] * 100,1)
WR.surv <- cbind(c("ALL", groups), WR.surv)
}
colnames(WR.surv)[1] <- "Release"
WR.surv <- merge(WR.surv, travel_final, by = "Release", all.x = T)
WR.surv$mean_travel_time <- round(WR.surv$mean_travel_time,1)
WR.surv$sd_travel_time <- round(WR.surv$sd_travel_time,1)
colnames(WR.surv) <- c("Release", "Survival (%)", "SE", "95% lower C.I.",
"95% upper C.I.", "Detection efficiency (%)", "Mean time to Tower (days)", "SD of time to Tower (days)","Count")
WR.surv <- WR.surv %>% arrange(., Release)
print(kable(WR.surv, row.names = F, "html", caption = "3.1 Minimum survival to Tower Bridge (using CJS
survival model), and travel time. If Yolo Bypass Weirs are overtopping during migration, fish may have taken
that route, and therefore this is a minimum estimate of survival") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"),
full_width = F, position = "left"))
}
}
3.1 Minimum survival to Tower Bridge (using CJS survival model). If Yolo Bypass Weirs are overtopping during migration, fish may have taken that route, and therefore this is a minimum estimate of survival
Release
Survival (%)
SE
95% lower C.I.
95% upper C.I.
Detection efficiency (%)
NA
NOT ENOUGH DETECTIONS
NA
NA
NA
NA
try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))
route_results_possible <- FALSE
if(nrow(detects_study[is.na(detects_study$DateTime_PST)==F,]) == 0){
results_short <- data.frame("Measure"=NA, "Estimate"="NO DETECTIONS YET", "SE"=NA, "95% lower C.I."=NA,
"95% upper C.I."=NA)
colnames(results_short) <- c("Measure", "Estimate", "SE", "95% lower C.I.", "95% upper C.I.")
print(kable(results_short, row.names = F, "html", caption = "3.2 Reach-specific survival and probability
of entering Georgiana Slough") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"),
full_width = F, position = "left"))
} else {
# Only do survival to Georg split for now
test2 <- detects_study %>%
filter(general_location %in% c("TowerBridge", "I80-50_Br","SacRiverWalnutGrove_2", "Sac_BlwGeorgiana",
"Sac_BlwGeorgiana2", "Georg_Sl_1", "Georgiana_Slough2"))
# We can only do a multi-state model if there is at least one detection in each route
if(nrow(test2[test2$general_location %in% c("Sac_BlwGeorgiana", "Sac_BlwGeorgiana2"),]) == 0 |
nrow(test2[test2$general_location %in% c("Georg_Sl_1", "Georgiana_Slough2"),]) == 0){
results_short <- data.frame("Measure"=NA, "Estimate"="NOT ENOUGH DETECTIONS", "SE"=NA, "95% lower C.I."=NA,
"95% upper C.I."=NA)
colnames(results_short) <- c("Measure", "Estimate", "SE", "95% lower C.I.", "95% upper C.I.")
print(kable(results_short, row.names = F, "html", caption = "3.2 Reach-specific survival and probability of
entering Georgiana Slough") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"),
full_width = F, position = "left"))
} else {
# Make tagcode character
study_tagcodes$TagCode <- as.character(study_tagcodes$TagCode)
# Make a crosstab query with frequencies for all tag/location combination
test2$general_location <- factor(test2$general_location,
levels = c("TowerBridge", "I80-50_Br","SacRiverWalnutGrove_2", "Sac_BlwGeorgiana",
"Sac_BlwGeorgiana2", "Georg_Sl_1", "Georgiana_Slough2"))
test2$TagCode <- factor(test2$TagCode, levels = study_tagcodes$TagCode)
mytable <- table(test2$TagCode, test2$general_location) # A will be rows, B will be columns
# Change all frequencies bigger than 1 to 1. Here you could change your minimum cutoff to 2 detections,
# and then make another command that changes all detections=1 to 0
mytable[mytable>0] <- "A"
# Order in order of rkm
mytable2 <- mytable[, c("TowerBridge", "I80-50_Br","SacRiverWalnutGrove_2", "Sac_BlwGeorgiana", "Sac_BlwGeorgiana2",
"Georg_Sl_1", "Georgiana_Slough2")]
# Now sort the crosstab rows alphabetically
mytable2 <- mytable2[order(row.names(mytable2)),]
mytable2[which(mytable2[, "Sac_BlwGeorgiana"]=="A"), "Sac_BlwGeorgiana"] <- "A"
mytable2[which(mytable2[, "Sac_BlwGeorgiana2"]=="A"), "Sac_BlwGeorgiana2"] <- "A"
mytable2[which(mytable2[, "Georg_Sl_1"]=="A"), "Georg_Sl_1"] <- "B"
mytable2[which(mytable2[, "Georgiana_Slough2"]=="A"), "Georgiana_Slough2"] <- "B"
# Now order the study_tagcodes table the same way
study_tagcodes <- study_tagcodes[order(study_tagcodes$TagCode),]
# Paste together (concatenate) the data from each column of the crosstab into one string per row, add to tagging_meta.
# For this step, make sure both are sorted by FishID
study_tagcodes$inp_part1 <- apply(mytable2[,1:3],1,paste,collapse="")
study_tagcodes$inp_partA <- apply(mytable2[,4:5],1,paste,collapse="")
study_tagcodes$inp_partB <- apply(mytable2[,6:7],1,paste,collapse="")
# find last detection at each genloc
departure <- aggregate(list(depart = test2$DateTime_PST), by = list(TagCode = test2$TagCode, last_location = test2$general_location), FUN = max)
# subset for just juncture locations
departure <- departure[departure$last_location %in% c("Sac_BlwGeorgiana", "Sac_BlwGeorgiana2", "Georg_Sl_1", "Georgiana_Slough2"),]
# Find genloc of last known detection per tag
last_depart <- aggregate(list(depart = departure$depart), by = list(TagCode = departure$TagCode), FUN = max)
last_depart1 <- merge(last_depart, departure)
study_tagcodes <- merge(study_tagcodes, last_depart1[,c("TagCode", "last_location")], by = "TagCode", all.x = T)
# Assume that the Sac is default pathway, and for fish that were detected in neither route, it would get a "00" in inp so does not matter anyway
study_tagcodes$inp_final <- paste("A",study_tagcodes$inp_part1, study_tagcodes$inp_partA," 1 ;", sep = "")
# now put in exceptions...fish that were seen in georgiana last
study_tagcodes[study_tagcodes$last_location %in% c("Georg_Sl_1", "Georgiana_Slough2"), "inp_final"] <-
paste("A",
study_tagcodes[study_tagcodes$last_location %in% c("Georg_Sl_1", "Georgiana_Slough2"), "inp_part1"],
study_tagcodes[study_tagcodes$last_location %in% c("Georg_Sl_1", "Georgiana_Slough2"), "inp_partB"],
" 1 ;",
sep = "")
# At this point, some fish might not have been deemed to ever take a route based on last visit analysis. If so, model cannot be run
if(any(grepl(pattern = "A", study_tagcodes$inp_final)==T) & any(grepl(pattern = "B", study_tagcodes$inp_final)==T)){
write.table(study_tagcodes$inp_final,"WRinp_multistate.inp",row.names = F, col.names = F, quote = F)
WRinp <- convert.inp("WRinp_multistate.inp")
dp <- process.data(WRinp, model="Multistrata")
ddl <- make.design.data(dp)
#### p ####
# Cannott be seen at 2B or 3B or 4B (tower or I80 or walnut grove) and currently second georg is missing so set to 0 as well and to 1 for first line
ddl$p$fix = NA
ddl$p$fix[ddl$p$stratum == "B" & ddl$p$time %in% c(2,3,4,6)] = 0
ddl$p$fix[ddl$p$stratum == "B" & ddl$p$time %in% c(5)] = 1
#### Psi ####
# Only 1 transition allowed:
# from A to B at time interval 4 to 5
ddl$Psi$fix = 0
# A to B can only happen for interval 3-4
ddl$Psi$fix[ddl$Psi$stratum == "A"&
ddl$Psi$tostratum == "B" &
ddl$Psi$time == 4] = NA
#### Phi a.k.a. S ####
ddl$S$fix = NA
# None in B for reaches 1,2,3,4 and fixing it to 1 for 5 (between two georg lines). All getting fixed to 1
ddl$S$fix[ddl$S$stratum == "B" & ddl$S$time %in% c(1,2,3,4,5)] = 1
# For route A, fixing it to 1 for 5 (between two blw_georg lines)
ddl$S$fix[ddl$S$stratum == "A" & ddl$S$time == 5] = 1
# We use -1 at beginning of formula to remove intercept. This is because different routes probably should not share the same intercept
p.timexstratum = list(formula=~-1+stratum:time)
Psi.stratumxtime = list(formula=~-1+stratum:time)
S.stratumxtime = list(formula=~-1+stratum:time)
# Run model a first time
S.timexstratum.p.timexstratum.Psi.timexstratum = mark(dp, ddl,
model.parameters = list(S = S.stratumxtime,p = p.timexstratum,Psi = Psi.stratumxtime),
realvcv = T, silent = T, output = F)
# Identify any parameter estimates at 1, which would likely have bad SE estimates.
profile.intervals <- which(S.timexstratum.p.timexstratum.Psi.timexstratum$results$real$estimate %in% c(0,1) &
!S.timexstratum.p.timexstratum.Psi.timexstratum$results$real$fixed == "Fixed")
# Rerun model using profile interval estimation for the tricky parameters
S.timexstratum.p.timexstratum.Psi.timexstratum = mark(dp, ddl,
model.parameters=list(S=S.stratumxtime,p= p.timexstratum,Psi=Psi.stratumxtime),
realvcv = T, profile.int = profile.intervals, silent = T, output = F)
results <- S.timexstratum.p.timexstratum.Psi.timexstratum$results$real
results_short <- results[rownames(results) %in% c("S sA g1 c1 a0 o1 t1",
"S sA g1 c1 a1 o2 t2",
"S sA g1 c1 a2 o3 t3",
"S sA g1 c1 a3 o4 t4",
"p sA g1 c1 a1 o1 t2",
"p sA g1 c1 a2 o2 t3",
"p sA g1 c1 a3 o3 t4",
"p sA g1 c1 a4 o4 t5",
"p sB g1 c1 a4 o4 t5",
"Psi sA toB g1 c1 a3 o4 t4"),]
results_short <- round(results_short[,c("estimate", "se", "lcl", "ucl")] * 100,1)
# Now find estimate and CIs for AtoA route at junction
Psilist = get.real(S.timexstratum.p.timexstratum.Psi.timexstratum,"Psi",vcv=TRUE)
Psivalues = Psilist$estimates
routes <- TransitionMatrix(Psivalues[Psivalues$time==4 & Psivalues$cohort==1,],vcv.real=Psilist$vcv.real)
results_short$Measure <- c("Survival from Release to TowerBridge (minimum estimate since fish may have taken Yolo
Bypass)", "Survival from TowerBridge to I80-50_Br", "% arrived from I80-50_Br to Walnut Grove (not survival because
fish may have taken Sutter/Steam)","Survival from Walnut Grove to Georgiana Slough confluence","Detection probability at TowerBridge",
"Detection probability at I80-50_Br", "Detection probability at Walnut Grove", "Detection probability at Blw_Georgiana",
#"Detection probability at Georgiana Slough (fixed at 1)",
"Routing probability into Georgiana Slough (Conditional on fish arriving to junction)")
results_short <- results_short[,c("Measure", "estimate", "se", "lcl", "ucl")]
colnames(results_short) <- c("Measure", "Estimate", "SE", "95% lower C.I.", "95% upper C.I.")
print(kable(results_short, row.names = F, "html", caption = "3.2 Reach-specific survival and probability of entering Georgiana Slough") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))
route_results_possible <- TRUE
} else {
results_short <- data.frame("Measure"=NA, "Estimate"="NOT ENOUGH DETECTIONS YET", "SE"=NA, "95% lower C.I."=NA, "95% upper C.I."=NA)
colnames(results_short) <- c("Measure", "Estimate", "SE", "95% lower C.I.", "95% upper C.I.")
print(kable(results_short, row.names = F, "html", caption = "3.2 Reach-specific survival and probability of entering Georgiana Slough") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))
}
}
}
3.2 Reach-specific survival and probability of entering Georgiana Slough
Measure
Estimate
SE
95% lower C.I.
95% upper C.I.
Survival from Release to TowerBridge (minimum estimate since fish may have taken Yolo Bypass)
100.0
4.6
0.0
100.0
Survival from TowerBridge to I80-50_Br
67.1
2409.8
0.0
100.0
% arrived from I80-50_Br to Walnut Grove (not survival because fish may have taken Sutter/Steam)
2.4
84.8
0.0
100.0
Survival from Walnut Grove to Georgiana Slough confluence
99.5
747.5
0.0
100.0
Detection probability at TowerBridge
0.4
0.4
0.1
2.8
Detection probability at I80-50_Br
0.0
0.0
0.0
0.0
Detection probability at Walnut Grove
50.0
35.4
5.9
94.1
Detection probability at Blw_Georgiana
33.6
338.9
0.0
100.0
Routing probability into Georgiana Slough (Conditional on fish arriving to junction)
25.1
190.6
0.0
100.0
# If you do not have access to local files, uncomment and run next lines of code
#download.file("https://raw.githubusercontent.com/CalFishTrack/real-time/master/data/georg.png",destfile = "georg.png", quiet = T, mode = "wb")
georg <- readPNG("georg.png")
par(mar = c(2,0,0,0))
# Set up the plot area
plot(1:2, type="n", xlab="", ylab="", xaxt = "n", yaxt = "n")
# Get the plot information so the image will fill the plot box, and draw it
lim <- par()
rasterImage(georg, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
if(nrow(detects_study[is.na(detects_study$DateTime_PST) == F,]) == 0){
legend(x = 1.55,y = 1.6,legend = "No detections yet",col = "white", box.col = "light gray", bg = "light gray")
legend(x = 1.55,y = 1.45,legend = "No detections yet",col = "white", box.col = "light gray", bg = "light gray")
} else if (route_results_possible == F){
legend(x = 1.55,y = 1.6,legend = "Too few detections",col = "white", box.col = "light gray", bg = "light gray")
legend(x = 1.55,y = 1.45,legend = "Too few detections",col = "white", box.col = "light gray", bg = "light gray")
} else {
legend(x = 1.55,y = 1.6,legend = paste(round(routes$TransitionMat["A","A"],3)*100,
"% (", round(routes$lcl.TransitionMat["A","A"],3)*100, "-",
round(routes$ucl.TransitionMat["A","A"],3)*100,")", sep =""),
col = "white", box.col = "light gray", bg = "light gray")
legend(1.55,1.45, legend = paste(round(routes$TransitionMat["A","B"],3)*100,
"% (", round(routes$lcl.TransitionMat["A","B"],3)*100, "-",
round(routes$ucl.TransitionMat["A","B"],3)*100,")", sep =""),
box.col = "light gray", bg = "light gray")
}
mtext(text = "3.3 Routing Probabilities at Georgiana Slough Junction (with 95% C.I.s)", cex = 1.3, side = 1, line = 0.2, adj = 0)
if (nrow(detects_study[is.na(detects_study$DateTime_PST)==F,]) == 0){
plot(1:2, type = "n",xaxt = "n", yaxt = "n",
xlab = "Range of days study fish were present at Georgiana Sl Junction",
ylab = "Routing probability into Georgiana Slough at the junction")
text(1.5,1.5, labels = "NO DETECTIONS YET", cex = 2)
} else if(route_results_possible == F){
plot(1:2, type = "n",xaxt = "n", yaxt = "n",
xlab = "Range of days study fish were present at Georgiana Sl Junction",
ylab = "Routing probability into Georgiana Slough at the junction")
text(1.5,1.5, labels = "TOO FEW DETECTIONS", cex = 2)
} else {
library(repmis)
trytest <- try(source_data("https://code.usgs.gov/crrl_qfes/Enhanced_Acoustic_Telemetry_Project/raw/master/EAT_data_2024.Rdata?raw=True"))
if (inherits(trytest, "try-error")){
plot(1:2, type = "n",xaxt = "n", yaxt = "n",
xlab = "Range of days study fish were present at Georgiana Sl Junction",
ylab = "Routing probability into Georgiana Slough at the junction")
text(1.5,1.5, labels = "ERROR DOWNLOADING STARS", cex = 2)
} else {
detects_5 <- detects_study %>% filter(general_location == "SacRiverWalnutGrove_2")
detects_5 <- detects_5 %>%
dplyr::left_join(., detects_5 %>%
group_by(TagCode) %>%
summarise(first_detect = min(DateTime_PST))) %>%
mutate(Day = as.Date(as.Date(first_detect, "Etc/GMT+8")))
tagcount5 <- detects_5 %>%
group_by(Day) %>%
summarise(unique_tags = length(unique(TagCode)))
tagcount5$density <- tagcount5$unique_tags / sum(tagcount5$unique_tags)
test2 <- as.data.frame(test2)
# first, find min and max arrivals at georg for a study
min_georg <- as.Date(format(min(test2[test2$general_location %in% c("Sac_BlwGeorgiana", "Sac_BlwGeorgiana2","Georg_Sl_1",
"Georgiana_Slough2"),"DateTime_PST"]), "%Y-%m-%d"))
max_georg <- as.Date(format(max(test2[test2$general_location %in% c("Sac_BlwGeorgiana", "Sac_BlwGeorgiana2","Georg_Sl_1",
"Georgiana_Slough2"),"DateTime_PST"]), "%Y-%m-%d"))
psi_study <- psi_GeoCond[psi_GeoCond$Date <= max_georg & psi_GeoCond$Date >=min_georg-1,]
ggplot()+
geom_bar(data = tagcount5, aes(x = Day, y = density), alpha = 0.4, stat = "identity", fill = "red") +
geom_ribbon(data = psi_study, aes(x=Date, ymin = psi_geo.10, ymax = psi_geo.90), alpha = 0.3) +
geom_line(data = psi_study, aes(x=Date, y=psi_geo.50), size = 1.2) +
geom_point(aes(x = median(rep(tagcount5$Day, tagcount5$unique_tags)), y = tail(results_short$Estimate,1)/100), size = 4)+
geom_errorbar(width = 0.2, size = 1.2, aes(ymin = tail(results_short$`95% lower C.I.`,1)/100, ymax = tail(results_short$`95% upper C.I.`,1)/100, x = median(rep(tagcount5$Day, tagcount5$unique_tags))))+
xlim(c(min_georg, max_georg))+
xlab("Range of days study fish were present at Georgiana Sl Junction") +
geom_hline(yintercept=0, linetype="dashed") +
geom_segment(aes(x=median(rep(tagcount5$Day, tagcount5$unique_tags)),y =-0.01,xend=median(rep(tagcount5$Day, tagcount5$unique_tags)),yend =0.01)) +
geom_text(aes(x=median(rep(tagcount5$Day, tagcount5$unique_tags)),y =-0.02, label = "median"), color = "red")+
scale_y_continuous(name = "Routing probability into Georgiana Slough at the junction",limits = c(-0.02,1), sec.axis = dup_axis(name="Histogram of arrivals to junction")) +
theme_bw()+
theme(axis.text.y.right = element_text(color = "red"), axis.title.y.right = element_text(color = "red"))
}
}
3.4 STARS prediction (ribbon) vs. empirical (point) estimate of Routing Probability at Georgiana Slough Junction
