Central Valley Enhanced

Acoustic Tagging Project

logo





6-Year Study San Joaquin River Steelhead - March Releases

2020-2021 Season (PROVISIONAL DATA)


1. Project Status

Study is complete, all tags are no longer active. All times in Pacific Standard Time.


Telemetry Study Template for this study can be found here

try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))

library(knitr)
library(kableExtra)
library(lubridate)
library(data.table)
library(ggplot2)
library(RMark)
library(scales)
library(viridis)
library(forcats)
library(reshape2)
library(png)
library(dataRetrieval)
library(rerddap)

##################################################################################################################
#### ASSIGN STUDY ID IN THIS NEXT LINE OF CODE ####
study <- "SJ_Steelhead_2021"
##################################################################################################################

detects_study <- fread("study_detections_archive.csv", stringsAsFactors = F, colClasses = c(DateTime_PST = "character", RelDT = "character"))
detects_study <- as.data.frame(detects_study[detects_study$Study_ID == study,])
detects_study$DateTime_PST <- as.POSIXct(detects_study$DateTime_PST, format = "%Y-%m-%d %H:%M:%S", tz="Etc/GMT+8")
detects_study$release_time <- as.POSIXct(detects_study$RelDT, format = "%Y-%m-%d %H:%M:%S", tz="Etc/GMT+8")

colnames(detects_study)[which(colnames(detects_study) == "Weight")] <- "weight"
colnames(detects_study)[which(colnames(detects_study) == "Length")] <- "length"
colnames(detects_study)[which(colnames(detects_study) == "Rel_rkm")] <- "release_rkm"
colnames(detects_study)[which(colnames(detects_study) == "Rel_loc")] <- "release_location"
colnames(detects_study)[which(colnames(detects_study) == "rkm")] <- "river_km"

latest <- read.csv("latest_download.csv", stringsAsFactors = F)$x

##################################################################################################################
#### TO RUN THE FOLLOWING CODE CHUNKS FROM HERE ON DOWN USING R ERDDAP, UN-COMMENT THESE NEXT 9 LINES OF CODE ####
##################################################################################################################
# cache_delete_all()
# query=paste('&',"Study_ID",'="',study,'"',sep = '')
# datafile=URLencode(paste("https://oceanview.pfeg.noaa.gov/erddap/tabledap/","FEDcalFishTrack",".csv?",query,sep = ''))
# options(url.method = "libcurl", download.file.method = "libcurl", timeout = 180)
# detects_study <- data.frame(read.csv(datafile,row.names = NULL, stringsAsFactors = F))
# detects_study <- detects_study[-1,]
# detects_study$DateTime_PST <- as.POSIXct(detects_study$local_time, format = "%Y-%m-%d %H:%M:%S", "Etc/GMT+8")
# detects_study$release_time <- as.POSIXct(detects_study$release_time, format = "%Y-%m-%d %H:%M:%S", "Etc/GMT+8")
# detects_study$river_km <- as.numeric(detects_study$river_km)
##################################################################################################################


if (nrow(detects_study) == 0){
  cat("Study has not yet begun")
}else{
  
  if (min(detects_study$release_time) > Sys.time()){
    cat("Study has not yet begun, below data is a placeholder:")
  }
  if (min(detects_study$release_time) < Sys.time()){
    cat(paste("Study began on ", min(detects_study$release_time), ", see tagging details below:", sep = ""))
  }
  
  ########################################################################
  #### ASSIGN RELEASE GROUPS HERE ####
  #######################################################################
    detects_study$Release <- detects_study$release_location
    #detects_study[detects_study$release_time > as.POSIXct("2021-03-09") & detects_study$release_time < as.POSIXct("2021-03-17"), "Release"] <- "Release 2"
    #detects_study[detects_study$release_time > as.POSIXct("2021-03-17"), "Release"] <- "Release 3"
  #######################################################################
  
  study_tagcodes <- unique(detects_study[,c("TagCode", "release_time", "weight", "length", "release_rkm", "release_location", "Release")])

  
  release_stats <- aggregate(list(First_release_time = study_tagcodes$release_time),
                             by= list(Release = study_tagcodes$Release),
                             FUN = min)
  release_stats <- merge(release_stats,
                         aggregate(list(Last_release_time = study_tagcodes$release_time),
                                   by= list(Release = study_tagcodes$Release),
                                   FUN = max),
                         by = c("Release"))
  
  
  release_stats <- merge(release_stats, aggregate(list(Number_fish_released =
                                                         study_tagcodes$TagCode),
                                                  by= list(Release = study_tagcodes$Release),
                                                  FUN = function(x) {length(unique(x))}),
                         by = c("Release"))
  
  release_stats <- merge(release_stats,
                         aggregate(list(Release_location = study_tagcodes$release_location),
                                   by= list(Release = study_tagcodes$Release),
                                   FUN = function(x) {head(x,1)}),
                         by = c("Release"))
  release_stats <- merge(release_stats,
                         aggregate(list(Release_rkm = study_tagcodes$release_rkm),
                                   by= list(Release = study_tagcodes$Release),
                                   FUN = function(x) {head(x,1)}),
                         by = c("Release"))
  release_stats <- merge(release_stats,
                         aggregate(list(Mean_length = as.numeric(study_tagcodes$length)),
                                   by= list(Release = study_tagcodes$Release),
                                   FUN = mean, na.rm = T),
                         by = c("Release"))
  release_stats <- merge(release_stats,
                         aggregate(list(Mean_weight = as.numeric(study_tagcodes$weight)),
                                   by= list(Release = study_tagcodes$Release),
                                   FUN = mean, na.rm = T),
                         by = c("Release"))
  

  release_stats[,c("Mean_length", "Mean_weight")] <- round(release_stats[,c("Mean_length", "Mean_weight")],1)
  
  release_stats$First_release_time <- format(release_stats$First_release_time, tz = "Etc/GMT+8")
  
  release_stats$Last_release_time <- format(release_stats$Last_release_time, tz = "Etc/GMT+8")
  
  release_stats <- release_stats[order(release_stats$First_release_time),]
  
  kable(release_stats, format = "html", row.names = F) %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left")

}
Study began on 2021-03-23 10:00:00, see tagging details below:
Release First_release_time Last_release_time Number_fish_released Release_location Release_rkm Mean_length Mean_weight
Durham_Ferry 2021-03-23 10:00:00 2021-03-26 11:30:00 200 Durham_Ferry 180.0 210.2 92.1
Stockton 2021-03-23 14:00:00 2021-03-26 17:15:00 100 Stockton 135.5 212.3 92.7
Head_of_Old_River 2021-03-23 16:00:00 2021-03-26 13:15:00 100 Head_of_Old_River 156.0 211.2 90.6


2. Real-time Fish Detections

try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))

## THIS CODE CHUNK WILL NOT WORK IF USING ONLY ERDDAP DATA, REQUIRES ACCESS TO LOCAL FILES
if (nrow(detects_study[is.na(detects_study$DateTime_PST)==F,]) == 0){
  "No detections yet"
} else {
  arrivals <- aggregate(list(DateTime_PST = detects_study$DateTime_PST), by = list(general_location = detects_study$general_location, TagCode = detects_study$TagCode), FUN = min)
    
  beacon_by_day <- fread("beacon_by_day.csv", stringsAsFactors = F)
  beacon_by_day$day <- as.Date(beacon_by_day$day)
  
  gen_locs <- read.csv("realtime_locs.csv", stringsAsFactors = F)
  
  arrivals$day <- as.Date(format(arrivals$DateTime_PST, "%Y-%m-%d", tz = "Etc/GMT+8"))
  
  arrivals_per_day <- aggregate(list(New_arrivals = arrivals$TagCode), by = list(day = arrivals$day, general_location = arrivals$general_location), length)
  arrivals_per_day$day <- as.Date(arrivals_per_day$day)

  ## Now subset to only look at data for the correct beacon for that day
  beacon_by_day <- as.data.frame(beacon_by_day[which(beacon_by_day$TagCode == beacon_by_day$beacon),])
  
  endtime <- min(as.Date(format(Sys.time(), "%Y-%m-%d")), max(as.Date(detects_study$release_time)+(as.numeric(detects_study$tag_life)*1.5)))
  ## Now only keep beacon by day for days since fish were released
  beacon_by_day <- beacon_by_day[beacon_by_day$day >= as.Date(min(study_tagcodes$release_time)) & beacon_by_day$day <= endtime,]  
  
  beacon_by_day <- merge(beacon_by_day, gen_locs[,c("location", "general_location","rkm")], by = "location", all.x = T)

  arrivals_per_day <- merge(unique(beacon_by_day[,c("general_location", "day", "rkm")]), arrivals_per_day, all.x = T, by = c("general_location", "day"))
  
  arrivals_per_day$day <- factor(arrivals_per_day$day)
  
  ## Remove bench test and other NA locations
  arrivals_per_day <- arrivals_per_day[!arrivals_per_day$general_location == "Bench_test",]
  arrivals_per_day <- arrivals_per_day[is.na(arrivals_per_day$general_location) == F,]

  ## Remove sites that were not operation the whole time
  gen_locs_days_in_oper <- aggregate(list(days_in_oper = arrivals_per_day$day), by = list(general_location = arrivals_per_day$general_location), FUN = length)
  gen_locs_days_in_oper <- gen_locs_days_in_oper[gen_locs_days_in_oper$days_in_oper == max(gen_locs_days_in_oper$days_in_oper),]
  
  arrivals_per_day_in_oper <- arrivals_per_day[arrivals_per_day$general_location %in% gen_locs_days_in_oper$general_location,]
  fish_per_site <- aggregate(list(fish_count = arrivals_per_day_in_oper$New_arrivals), by = list(general_location = arrivals_per_day_in_oper$general_location), FUN = sum, na.rm = T)
  
  active_gen_locs <- gen_locs[is.na(gen_locs$stop),]
  active_gen_locs <- active_gen_locs[active_gen_locs$general_location %in% fish_per_site$general_location,]
  ## estimate mean lat and lons for each genloc
  gen_locs_mean_coords <- aggregate(list(latitude = active_gen_locs$latitude), by = list(general_location = active_gen_locs$general_location), FUN = mean)
  gen_locs_mean_coords <- merge(gen_locs_mean_coords, aggregate(list(longitude = active_gen_locs$longitude), by = list(general_location = active_gen_locs$general_location), FUN = mean))
  
  fish_per_site <- merge(fish_per_site, gen_locs_mean_coords)
  
  library(leaflet)
  library(maps)
  library(htmlwidgets)
  library(leaflet.extras)

  icons <- awesomeIcons(iconColor = "lightblue",
                      #library = "ion",
                      text = fish_per_site$fish_count)
  
  leaflet(data = fish_per_site) %>%
      # setView(-72.14600, 43.82977, zoom = 8) %>% 
      addProviderTiles("Esri.WorldStreetMap", group = "Map") %>%
      addProviderTiles("Esri.WorldImagery", group = "Satellite") %>% 
      addProviderTiles("Esri.WorldShadedRelief", group = "Relief") %>%
      # Marker data are from the sites data frame. We need the ~ symbols
      # to indicate the columns of the data frame.
      addMarkers(~longitude, ~latitude, label = ~fish_count, group = "Receiver Sites", popup = ~general_location, labelOptions = labelOptions(noHide = T, textsize = "15px")) %>% 
      #addAwesomeMarkers(~longitude, ~latitude, icon = icons, labelOptions(textsize = "15px")) %>%
      addScaleBar(position = "bottomleft") %>%
      addLayersControl(
          baseGroups = c("Street Map", "Satellite", "Relief"),
          options = layersControlOptions(collapsed = FALSE))
}

2.1 Map of unique fish detections at operational realtime detection locations 


try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))

detects_benicia <- detects_study[detects_study$general_location %in% c("Benicia_west", "Benicia_east"),]

if (nrow(detects_benicia)>0) {
  detects_benicia <- merge(detects_benicia,aggregate(list(first_detect = detects_benicia$DateTime_PST), by = list(TagCode= detects_benicia$TagCode), FUN = min))
  
  detects_benicia$Day <- as.Date(detects_benicia$first_detect, "Etc/GMT+8")
  
  starttime <- as.Date(min(detects_benicia$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_benicia$release_time)+(as.numeric(detects_benicia$tag_life))))
  #wlk_flow <- cdec_query("COL", "20", "H", starttime, endtime+1)
  #wlk_flow$datetime <- as.Date(wlk_flow$datetime)
  #wlk_flow_day <- aggregate(list(parameter_value = wlk_flow$parameter_value), by = list(Day = wlk_flow$datetime), FUN = mean, na.rm = T)
  
  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_benicia$Release))])
  
  tagcount <- aggregate(list(unique_tags = detects_benicia$TagCode), by = list(Day = detects_benicia$Day, Release = detects_benicia$Release ), FUN = function(x){length(unique(x))})
  tagcount1 <- reshape2::dcast(tagcount, Day ~ Release)
                    
  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 <- merge(daterange1, wlk_flow_day, by = "Day", all.x = T)
  daterange2 <- daterange1
  
  rownames(daterange2) <- daterange2$Day
  daterange2$Day <- NULL
  
  par(mar=c(6, 5, 2, 5) + 0.1)
  # barp <- barplot(t(daterange2[,1:ncol(daterange2)]), plot = FALSE, beside = T)
  # barplot(t(daterange2[,1:ncol(daterange2)]), beside = T, col=brewer.pal(n = rel_num, name = "Dark2"), 
  #         xlab = "", ylab = "Number of fish arrivals per day", 
  #         ylim = c(0,max(daterange2[,1:ncol(daterange2)], na.rm = T)*1.2), 
  #         las = 2, xlim=c(0,max(barp)+1), cex.lab = 1.5, yaxt = "n", xaxt = "n", border = NA)#, 
  #         #legend.text = colnames(daterange2[,1:ncol(daterange2)-1]),
  #         #args.legend = list(x ='topright', bty='n', inset=c(-0.2,0)), title = "Release Group")
  # legend(x ='topleft', legend = colnames(daterange2)[1:ncol(daterange2)], fill= brewer.pal(n = rel_num, name = "Set1"), horiz = T, title = "Release")
  # ybreaks <- if(max(daterange2[,1:ncol(daterange2)], na.rm = T) < 4) {max(daterange2[,1:ncol(daterange2)], na.rm = T)} else {5}
  # xbreaks <- if(ncol(barp) > 10) {seq(1, ncol(barp), 2)} else {1:ncol(barp)}
  # barpmeans <- colMeans(barp)
  # axis(1, at = barpmeans[xbreaks], labels = rownames(daterange2)[xbreaks], las = 2)
  # axis(2, at = pretty(0:max(daterange2[,1:ncol(daterange2)], na.rm = T), ybreaks))
  # box()
  daterange2$Date <- as.Date(row.names(daterange2))
  daterange3 <- melt(daterange2, id.vars = "Date", variable.name = ".", )
  
  # p <- ggplot(data = daterange3, aes(x = Date, y = value, color = ., fill = .)) +
  #   geom_bar(stat='identity') +
  #   ylab("Number of fish arrivals per day") +
  #   #xlim(range(daterange$Day)) +
  #   #geom_line(data= daterange2_flow, aes(x = Date, y = parameter_value/500), color = alpha("#947FFF", alpha = 0.5))+
  #   #scale_x_date(date_breaks = "5 days") +
  #   #scale_y_continuous(name = "Number of fish arrivals per day",
  #     # Add a second axis and specify its features
  #   #  sec.axis = sec_axis(~.*500, name="Second Axis")) +
  #   theme_bw() +
  #   theme(panel.border = element_rect(colour = "black", fill=NA))

  
  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 items to isolate)", xref="paper", yref="paper",
                     x=0.01, xanchor="left",
                     y=1.056, yanchor="top",    # Same y as legend below
                     legendtitle=TRUE, showarrow=FALSE ) %>%
    #add_lines(x=~daterange2_flow$Date, y=~daterange2_flow$parameter_value, line = list(color = alpha("#947FFF", alpha = 0.5)), yaxis="y2", showlegend=FALSE, inherit=FALSE) %>%
    layout(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.066),
          margin=list(l = 50,r = 100,b = 50,t = 50)
           )



}else{
  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)
}

2.4 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 = "")))
  
try(benicia <- read.csv("benicia_surv.csv", stringsAsFactors = F))

detects_benicia <- detects_study[detects_study$general_location %in% c("Benicia_west", "Benicia_east"),]
endtime <- min(as.Date(format(Sys.time(), "%Y-%m-%d")), max(as.Date(detects_study$release_time)+(as.numeric(detects_study$tag_life)*1.5)))

if (nrow(detects_benicia) == 0){
  if(as.numeric(difftime(Sys.time(), min(detects_study$RelDT), units = "days"))>30){
    WR.surv <- data.frame("Release"="ALL", "estimate"=0, "se"=NA, "lcl"=NA, "ucl"=NA, "Detection_efficiency"=NA)
  }else{
    WR.surv <- data.frame("Release"=NA, "estimate"="NO DETECTIONS YET", "se"=NA, "lcl"=NA, "ucl"=NA, "Detection_efficiency"=NA)
  }
  WR.surv1 <- WR.surv
  colnames(WR.surv1) <- c("Release Group", "Survival (%)", "SE", "95% lower C.I.", "95% upper C.I.", "Detection efficiency (%)")
  print(kable(WR.surv1, row.names = F, "html", caption = "3.5 Minimum survival to Benicia Bridge East Span (using CJS survival model)") %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))  
}else if (length(table(detects_benicia$general_location)) == 1){
  if(as.numeric(difftime(Sys.time(), min(detects_study$RelDT), units = "days"))>30){
    WR.surv <- data.frame("Release"="ALL", "estimate"=round(length(unique(detects_benicia$TagCode))/length(unique(detects_study$TagCode))*100,1), "se"=NA, "lcl"=NA, "ucl"=NA, "Detection_efficiency"=NA)
  }else{
    WR.surv <- data.frame("Release"=NA, "estimate"="NOT ENOUGH DETECTIONS", "se"=NA, "lcl"=NA, "ucl"=NA, "Detection_efficiency"=NA)
  }
    WR.surv1 <- WR.surv
    colnames(WR.surv1) <- c("Release Group", "Survival (%)", "SE", "95% lower C.I.", "95% upper C.I.", "Detection efficiency (%)")
    print(kable(WR.surv1, row.names = F, "html", caption = "3.5 Minimum survival to Benicia Bridge East Span (using CJS survival model)") %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))
} else {  
  ## Only do survival to Benicia here
  test3 <- detects_study[which(detects_study$river_km < 53),]
  
  ## Create inp for survival estimation
  
  inp <- as.data.frame(reshape2::dcast(test3, TagCode ~ river_km, fun.aggregate = length))
  
  ## Sort columns by river km in descending order
  # Count number of genlocs
  gen_loc_sites <- ncol(inp)-1
  
  inp <- inp[,c(1,order(names(inp[,2:(gen_loc_sites+1)]), decreasing = T)+1)]

  inp <- merge(study_tagcodes, inp, by = "TagCode", all.x = T)
  
  inp2 <- inp[,(ncol(inp)-gen_loc_sites+1):ncol(inp)]
  inp2[is.na(inp2)] <- 0
  inp2[inp2 > 0] <- 1
  
  inp <- cbind(inp, inp2)
  groups <- as.character(sort(unique(inp$Release)))
  groups_w_detects <- names(table(test3$Release))

  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 = inp$Release, stringsAsFactors = F)

    WR.process <- process.data(inp.df, 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)
    
    inp.df <- inp.df[inp.df$rel %in% groups_w_detects,]
    inp.df$rel <- factor(inp.df$rel, levels = groups_w_detects)
    if(length(groups_w_detects) > 1){
      WR.process <- process.data(inp.df, model="CJS", begin.time=1, groups = "rel") 
    
      WR.ddl <- make.design.data(WR.process)
      WR.mark.rel <- mark(WR.process, WR.ddl, model.parameters=list(Phi=list(formula=~time*rel),p=list(formula=~time)), silent = T, output = F)
    }else{  
      WR.process <- process.data(inp.df, model="CJS", begin.time=1) 
    
      WR.ddl <- make.design.data(WR.process)
      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 <- cbind(Release = "ALL",round(WR.mark.all$results$real[1,c("estimate", "se", "lcl", "ucl")] * 100,1))
    WR.surv.rel <- cbind(Release = groups_w_detects, round(WR.mark.rel$results$real[seq(from=1,to=length(groups_w_detects)*2,by = 2),c("estimate", "se", "lcl", "ucl")] * 100,1))
    WR.surv.rel <- merge(WR.surv.rel, data.frame(Release = groups), all.y = T)
    WR.surv.rel[is.na(WR.surv.rel$estimate),"estimate"] <- 0
    WR.surv <- rbind(WR.surv, WR.surv.rel)
    
  }else{
    inp.df <- data.frame(ch = as.character(inp$inp_final), freq = 1, stringsAsFactors = F)

    WR.process <- process.data(inp.df, 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.surv <- cbind(Release = c("ALL", groups),round(WR.mark.all$results$real[1,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.surv1 <- WR.surv
  colnames(WR.surv1) <- c("Release Group", "Survival (%)", "SE", "95% lower C.I.", "95% upper C.I.", "Detection efficiency (%)")

  print(kable(WR.surv1, row.names = F, "html", caption = "3.5 Minimum survival to Benicia Bridge East Span (using CJS survival model)") %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))    
  
}
3.5 Minimum survival to Benicia Bridge East Span (using CJS survival model)
Release Group Survival (%) SE 95% lower C.I. 95% upper C.I. Detection efficiency (%)
ALL 3 0.9 1.7 5.2 83.3
Durham_Ferry 3 1.2 1.4 6.5 NA
Head_of_Old_River 1 1.0 0.1 6.8 NA
Stockton 5 2.2 2.1 11.5 NA 
if(exists("benicia")==T & is.numeric(WR.surv1[1,2])){
  ## Find mean release time per release group, and ALL
  reltimes <- aggregate(list(RelDT = study_tagcodes$release_time), by = list(Release = study_tagcodes$Release), FUN = mean)
  reltimes <- rbind(reltimes, data.frame(Release = "ALL", RelDT = mean(study_tagcodes$release_time)))

  ## Assign whether the results are tentative or final
  quality <- "tentative"
  if(endtime < as.Date(format(Sys.time(), "%Y-%m-%d"))) { quality <- "final"}

  WR.surv <- merge(WR.surv, reltimes, by = "Release", all.x = T)
  
  WR.surv$RelDT <- as.POSIXct(WR.surv$RelDT, origin = '1970-01-01')
  benicia$RelDT <- as.POSIXct(benicia$RelDT)
  ## remove old benicia record for this studyID
  benicia <- benicia[!benicia$StudyID == unique(detects_study$Study_ID),]
  
  benicia <- rbind(benicia, data.frame(WR.surv, StudyID = unique(detects_study$Study_ID), data_quality = quality))
  
  write.csv(benicia, "benicia_surv.csv", row.names = F, quote = F) 
}


4. Detections statistics at all realtime receivers

try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))

if (nrow(detects_study[is.na(detects_study$DateTime_PST)==F,]) == 0){
  "No detections yet"
} else {

  arrivals <- aggregate(list(DateTime_PST = detects_study$DateTime_PST), by = list(general_location = detects_study$general_location, TagCode = detects_study$TagCode), FUN = min)
  
  tag_stats <- aggregate(list(First_arrival = arrivals$DateTime_PST), 
                         by= list(general_location = arrivals$general_location),
                         FUN = min)
  tag_stats <- merge(tag_stats, 
                     aggregate(list(Mean_arrival = arrivals$DateTime_PST), 
                         by= list(general_location = arrivals$general_location),
                         FUN = mean), 
                     by = c("general_location"))
  tag_stats <- merge(tag_stats, 
                     aggregate(list(Last_arrival = arrivals$DateTime_PST), 
                         by= list(general_location = arrivals$general_location),
                         FUN = max), 
                     by = c("general_location"))
  tag_stats <- merge(tag_stats, 
                     aggregate(list(Fish_count = arrivals$TagCode), 
                         by= list(general_location = arrivals$general_location), 
                         FUN = function(x) {length(unique(x))}), 
                     by = c("general_location"))
  tag_stats$Percent_arrived <- round(tag_stats$Fish_count/nrow(study_tagcodes) * 100,2)
      
  tag_stats <- merge(tag_stats, unique(detects_study[,c("general_location", "river_km")]))
  
  tag_stats <- tag_stats[order(tag_stats$river_km, decreasing = T),]
  
  tag_stats[,c("First_arrival", "Mean_arrival", "Last_arrival")] <- format(tag_stats[,c("First_arrival", "Mean_arrival", "Last_arrival")], tz = "Etc/GMT+8")
  
  tag_stats <- tag_stats[is.na(tag_stats$First_arrival)==F,]

  print(kable(tag_stats, row.names = F, 
              caption = "4.1 Detections for all releases combined",
              "html") %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))
  
  for (j in sort(unique(study_tagcodes$Release))) {
    
    if(nrow(detects_study[detects_study$Release == j,]) > 0 ) {
    
      temp <- detects_study[detects_study$Release == j,]
      
        arrivals1 <- aggregate(list(DateTime_PST = temp$DateTime_PST), by = list(general_location = temp$general_location, TagCode = temp$TagCode), FUN = min)
  
      rel_count <- nrow(study_tagcodes[study_tagcodes$Release == j,])
  
      tag_stats1 <- aggregate(list(First_arrival = arrivals1$DateTime_PST), 
                             by= list(general_location = arrivals1$general_location), 
                             FUN = min)
      tag_stats1 <- merge(tag_stats1, 
                         aggregate(list(Mean_arrival = arrivals1$DateTime_PST), 
                             by= list(general_location = arrivals1$general_location), 
                             FUN = mean), 
                         by = c("general_location"))
      tag_stats1 <- merge(tag_stats1, 
                   aggregate(list(Last_arrival = arrivals1$DateTime_PST), 
                       by= list(general_location = arrivals1$general_location), 
                       FUN = max), 
                   by = c("general_location"))
      tag_stats1 <- merge(tag_stats1, 
                         aggregate(list(Fish_count = arrivals1$TagCode), 
                                   by= list(general_location = arrivals1$general_location), 
                                   FUN = function(x) {length(unique(x))}), 
                         by = c("general_location"))
      
      tag_stats1$Percent_arrived <- round(tag_stats1$Fish_count/rel_count * 100,2)
    
      tag_stats1 <- merge(tag_stats1, unique(detects_study[,c("general_location", "river_km")]))
    
      tag_stats1 <- tag_stats1[order(tag_stats1$river_km, decreasing = T),]
      
      tag_stats1[,c("First_arrival", "Mean_arrival", "Last_arrival")] <- format(tag_stats1[,c("First_arrival", "Mean_arrival", "Last_arrival")], tz = "Etc/GMT+8")
      
      tag_stats1 <- tag_stats1[is.na(tag_stats1$First_arrival)==F,]
      
      final_stats <- kable(tag_stats1, row.names = F, 
            caption = paste("4.2 Detections for",j,"release groups", sep = " "),
            "html")
      
      print(kable_styling(final_stats, bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"), full_width = F, position = "left"))
      
    } else {
      cat("\n\n\\pagebreak\n")
      print(paste("No detections for",j,"release group yet", sep=" "), quote = F)
      cat("\n\n\\pagebreak\n")
    }
  }
}
4.1 Detections for all releases combined
general_location First_arrival Mean_arrival Last_arrival Fish_count Percent_arrived river_km
Old River 2021-03-27 22:52:38 2021-04-05 05:43:39 2021-04-26 19:31:10 10 2.50 153.001
MiddleRiver 2021-03-27 23:39:34 2021-04-07 14:25:29 2021-04-19 20:51:53 10 2.50 150.000
Holland_Cut_Quimby 2021-03-28 13:29:47 2021-04-08 04:13:19 2021-04-23 13:03:08 12 3.00 145.000
CVP_Trash_Rack_1 2021-05-30 06:12:44 2021-06-11 12:14:33 2021-06-18 09:08:23 3 0.75 144.500
SWP_intake 2021-04-07 09:12:03 2021-04-18 21:55:36 2021-04-28 08:05:12 11 2.75 142.721
SWP_radial_gates_DS 2021-04-20 05:41:39 2021-04-26 08:17:26 2021-05-02 04:27:53 4 1.00 142.721
SWP_radial_gates_US 2021-03-29 17:42:58 2021-04-10 21:13:08 2021-04-27 09:40:07 16 4.00 142.721
Old_River_Quimby 2021-03-30 20:59:49 2021-04-11 23:59:44 2021-04-24 02:59:40 2 0.50 141.000
Benicia_east 2021-04-03 09:01:07 2021-04-14 07:52:04 2021-04-27 08:32:41 10 2.50 52.240
Benicia_west 2021-04-03 09:03:35 2021-04-14 12:42:41 2021-04-27 08:35:40 12 3.00 52.040
4.2 Detections for Durham_Ferry release groups
general_location First_arrival Mean_arrival Last_arrival Fish_count Percent_arrived river_km
Old River 2021-03-30 23:15:53 2021-04-09 11:12:45 2021-04-26 19:31:10 6 3.0 153.001
MiddleRiver 2021-04-05 04:22:53 2021-04-12 12:39:21 2021-04-19 20:51:53 4 2.0 150.000
Holland_Cut_Quimby 2021-03-31 10:20:56 2021-04-12 04:07:15 2021-04-23 13:03:08 8 4.0 145.000
CVP_Trash_Rack_1 2021-05-30 06:12:44 2021-05-30 06:12:44 2021-05-30 06:12:44 1 0.5 144.500
SWP_intake 2021-04-07 09:12:03 2021-04-18 21:22:55 2021-04-28 08:05:12 8 4.0 142.721
SWP_radial_gates_DS 2021-04-20 05:41:39 2021-04-26 14:39:02 2021-05-02 04:27:53 3 1.5 142.721
SWP_radial_gates_US 2021-03-30 14:52:16 2021-04-12 05:40:12 2021-04-27 09:40:07 11 5.5 142.721
Old_River_Quimby 2021-04-24 02:59:40 2021-04-24 02:59:40 2021-04-24 02:59:40 1 0.5 141.000
Benicia_east 2021-04-12 09:02:05 2021-04-20 05:00:59 2021-04-27 08:32:41 4 2.0 52.240
Benicia_west 2021-04-11 17:20:15 2021-04-18 15:35:30 2021-04-27 08:35:40 6 3.0 52.040
4.2 Detections for Head_of_Old_River release groups
general_location First_arrival Mean_arrival Last_arrival Fish_count Percent_arrived river_km
Old River 2021-03-27 22:52:38 2021-03-28 09:04:50 2021-03-28 20:49:28 3 3 153.001
MiddleRiver 2021-03-27 23:39:34 2021-03-29 07:04:59 2021-03-30 14:30:25 2 2 150.000
Holland_Cut_Quimby 2021-03-28 13:29:47 2021-03-29 11:06:03 2021-03-30 10:55:40 3 3 145.000
CVP_Trash_Rack_1 2021-06-16 21:22:33 2021-06-16 21:22:33 2021-06-16 21:22:33 1 1 144.500
SWP_intake 2021-04-10 12:04:11 2021-04-18 10:22:53 2021-04-26 08:41:36 2 2 142.721
SWP_radial_gates_DS 2021-04-25 13:12:39 2021-04-25 13:12:39 2021-04-25 13:12:39 1 1 142.721
SWP_radial_gates_US 2021-03-29 17:42:58 2021-04-03 14:13:20 2021-04-09 07:22:11 3 3 142.721
Old_River_Quimby 2021-03-30 20:59:49 2021-03-30 20:59:49 2021-03-30 20:59:49 1 1 141.000
Benicia_east 2021-04-19 15:13:22 2021-04-19 15:13:22 2021-04-19 15:13:22 1 1 52.240
Benicia_west 2021-04-19 15:17:27 2021-04-19 15:17:27 2021-04-19 15:17:27 1 1 52.040
4.2 Detections for Stockton release groups
general_location First_arrival Mean_arrival Last_arrival Fish_count Percent_arrived river_km
Old River 2021-04-03 10:45:34 2021-04-03 10:45:34 2021-04-03 10:45:34 1 1 153.001
MiddleRiver 2021-04-03 09:50:31 2021-04-07 07:51:51 2021-04-14 20:29:05 4 4 150.000
Holland_Cut_Quimby 2021-04-05 08:23:42 2021-04-05 08:23:42 2021-04-05 08:23:42 1 1 145.000
CVP_Trash_Rack_1 2021-06-18 09:08:23 2021-06-18 09:08:23 2021-06-18 09:08:23 1 1 144.500
SWP_intake 2021-04-20 01:22:32 2021-04-20 01:22:32 2021-04-20 01:22:32 1 1 142.721
SWP_radial_gates_US 2021-04-12 06:29:27 2021-04-14 09:14:01 2021-04-16 11:58:35 2 2 142.721
Benicia_east 2021-04-03 09:01:07 2021-04-08 13:28:41 2021-04-19 12:16:09 5 5 52.240
Benicia_west 2021-04-03 09:03:35 2021-04-08 13:32:20 2021-04-19 12:19:41 5 5 52.040


4.3 Fish arrivals per day (“NA” means receivers were non-operational)

try(setwd(paste(file.path(Sys.getenv("USERPROFILE"),"Desktop",fsep="\\"), "\\Real-time data massaging\\products", sep = "")))

## THIS CODE CHUNK WILL NOT WORK IF USING ONLY ERDDAP DATA, REQUIRES ACCESS TO LOCAL FILES
if (nrow(detects_study[is.na(detects_study$DateTime_PST)==F,]) == 0){
  "No detections yet"
} else {
  arrivals <- aggregate(list(DateTime_PST = detects_study$DateTime_PST), by = list(general_location = detects_study$general_location, TagCode = detects_study$TagCode), FUN = min)
    
  beacon_by_day <- fread("beacon_by_day.csv", stringsAsFactors = F)
  beacon_by_day$day <- as.Date(beacon_by_day$day)
  
  gen_locs <- read.csv("realtime_locs.csv", stringsAsFactors = F)
  
  arrivals$day <- as.Date(format(arrivals$DateTime_PST, "%Y-%m-%d", tz = "Etc/GMT+8"))
  
  arrivals_per_day <- aggregate(list(New_arrivals = arrivals$TagCode), by = list(day = arrivals$day, general_location = arrivals$general_location), length)
  arrivals_per_day$day <- as.Date(arrivals_per_day$day)

  ## Now subset to only look at data for the correct beacon for that day
  beacon_by_day <- as.data.frame(beacon_by_day[which(beacon_by_day$TagCode == beacon_by_day$beacon),])
  
  endtime <- min(as.Date(format(Sys.time(), "%Y-%m-%d")), max(as.Date(detects_study$release_time)+(as.numeric(detects_study$tag_life)*1.5)))
  ## Now only keep beacon by day for days since fish were released
  beacon_by_day <- beacon_by_day[beacon_by_day$day >= as.Date(min(study_tagcodes$release_time)) & beacon_by_day$day <= endtime,]  
  
  beacon_by_day <- merge(beacon_by_day, gen_locs[,c("location", "general_location","rkm")], by = "location", all.x = T)

  arrivals_per_day <- merge(unique(beacon_by_day[,c("general_location", "day", "rkm")]), arrivals_per_day, all.x = T, by = c("general_location", "day"))
  
  arrivals_per_day$day <- factor(arrivals_per_day$day)
  
  ## Remove bench test and other NA locations
  arrivals_per_day <- arrivals_per_day[!arrivals_per_day$general_location == "Bench_test",]
  arrivals_per_day <- arrivals_per_day[is.na(arrivals_per_day$general_location) == F,]
  
  ## Change order of data to plot decreasing river_km
  arrivals_per_day <- arrivals_per_day[order(arrivals_per_day$rkm, decreasing = T),]
  arrivals_per_day$general_location <- factor(arrivals_per_day$general_location, unique(arrivals_per_day$general_location))
  
  # 
  # ggplot(data=arrivals_per_day, aes(x=general_location, y=fct_rev(as_factor(day)))) +
  # geom_tile(fill = "lightgray", color = "black") + 
  # geom_text(aes(label=New_arrivals)) +
  # labs(x="General Location", y = "Date") +
  # theme(panel.background = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1))
    crosstab <- xtabs(formula = arrivals_per_day$New_arrivals ~ arrivals_per_day$day + arrivals_per_day$general_location, addNA =T)
  crosstab[is.na(crosstab)] <- ""
  crosstab[crosstab==0] <- NA
  crosstab <- as.data.frame.matrix(crosstab)
  #colnames(crosstab) <- c("Butte Br", "Tower Br", "I8050 Br", "Old River", "Middle River", "CVP Tanks", "Georg Slough1", "Sac_Blw Georg1", "Georg Slough2", "Sac_Blw Georg2", "Benicia East", "Benicia West")

 kable(crosstab, align = "c") %>%
  kable_styling(c("striped", "condensed"), font_size = 11, full_width = F, position = "left") %>%
  #row_spec(0, angle = -45) %>%
  column_spec(column = 1:ncol(crosstab),width_min = "50px",border_left = T, border_right = T) %>%
  column_spec(1, bold = T, width_min = "75px")%>%
  scroll_box(height = "700px")
}
ButteBrRT TowerBridge I80-50_Br Old River MiddleRiver Holland_Cut_Quimby CVP_Tank CVP_Trash_Rack_1 SWP_intake SWP_radial_gates_DS SWP_radial_gates_US Old_River_Quimby Georgiana_Slough1 Sac_BlwGeorgiana Georgiana_Slough2 Sac_BlwGeorgiana2 Benicia_east Benicia_west
2021-03-23 NA NA
2021-03-24 NA NA
2021-03-25 NA NA
2021-03-26 NA NA
2021-03-27 1 1 NA NA
2021-03-28 2 1 NA NA
2021-03-29 1 NA NA 1
2021-03-30 1 1 1 NA NA 1 1
2021-03-31 1 1 NA NA 1
2021-04-01 NA NA
2021-04-02 NA NA 2
2021-04-03 1 2 NA NA 1 1
2021-04-04 1 1 NA NA
2021-04-05 1 2 NA NA 1 1
2021-04-06 NA NA 1 1 1
2021-04-07 1 NA 2 NA NA 1 1
2021-04-08 1 NA NA NA
2021-04-09 1 NA NA 2 NA NA
2021-04-10 NA 1 NA NA NA
2021-04-11 1 NA NA NA NA 1
2021-04-12 NA NA 1 NA NA 1 1
2021-04-13 NA NA 1 NA NA
2021-04-14 1 1 NA NA NA NA NA NA
2021-04-15 1 NA NA NA NA 1 NA NA
2021-04-16 1 1 NA NA NA NA 1 NA NA 1 1
2021-04-17 NA NA NA NA NA NA
2021-04-18 NA NA NA NA NA NA
2021-04-19 1 NA 1 2 NA NA 2 3
2021-04-20 1 NA 3 1 1 NA NA
2021-04-21 NA NA NA
2021-04-22 NA 1 NA NA
2021-04-23 1 NA NA NA
2021-04-24 NA 1 NA NA 1 1
2021-04-25 NA 1 1 NA NA
2021-04-26 1 NA 1 NA NA
2021-04-27 NA 1 1 NA NA 1 1
2021-04-28 NA 1 NA NA
2021-04-29 NA NA NA
2021-04-30 NA NA NA
2021-05-01 NA NA NA
2021-05-02 NA 1 NA NA
2021-05-03 NA NA NA
2021-05-04 NA NA NA
2021-05-05 NA NA NA
2021-05-06 NA NA NA
2021-05-07 NA NA NA
2021-05-08 NA NA NA
2021-05-09 NA NA NA
2021-05-10 NA NA NA
2021-05-11 NA NA NA
2021-05-12 NA NA NA
2021-05-13 NA NA NA
2021-05-14 NA NA
2021-05-15 NA NA
2021-05-16 NA NA
2021-05-17 NA NA
2021-05-18 NA NA
2021-05-19 NA NA
2021-05-20 NA NA
2021-05-21 NA NA
2021-05-22 NA NA
2021-05-23 NA NA
2021-05-24 NA NA NA
2021-05-25 NA NA
2021-05-26 NA NA NA
2021-05-27 NA NA NA
2021-05-28 NA NA NA
2021-05-29 NA NA NA
2021-05-30 1 NA NA NA
2021-05-31 NA NA NA
2021-06-01 NA NA NA
2021-06-02 NA NA NA
2021-06-03 NA NA NA
2021-06-04 NA NA NA
2021-06-05 NA NA NA
2021-06-06 NA NA NA
2021-06-07 NA NA NA
2021-06-08 NA NA NA
2021-06-09 NA NA NA
2021-06-10 NA NA NA
2021-06-11 NA NA NA
2021-06-12 NA NA NA
2021-06-13 NA NA NA
2021-06-14 NA NA NA
2021-06-15 NA NA NA
2021-06-16 1 NA NA NA
2021-06-17 NA NA NA
2021-06-18 1 NA NA NA
2021-06-19 NA NA NA
2021-06-20 NA NA NA
2021-06-21 NA NA NA
2021-06-22 NA NA NA
2021-06-23 NA NA NA
2021-06-24 NA NA NA
2021-06-25 NA NA NA
2021-06-26 NA NA NA
2021-06-27 NA NA NA
2021-06-28 NA NA NA
2021-06-29 NA NA NA
2021-06-30 NA NA NA
2021-07-01 NA NA NA
2021-07-02 NA NA NA
2021-07-03 NA NA NA
2021-07-04 NA NA NA
2021-07-05 NA NA NA
2021-07-06 NA NA NA
2021-07-07 NA NA NA NA
2021-07-08 NA NA NA NA
2021-07-09 NA NA NA NA
2021-07-10 NA NA NA NA
2021-07-11 NA NA NA NA
2021-07-12 NA NA NA NA
2021-07-13 NA NA NA NA
2021-07-14 NA NA NA NA
2021-07-15 NA NA NA NA
2021-07-16 NA NA NA NA
2021-07-17 NA NA NA NA
2021-07-18 NA NA NA NA
2021-07-19 NA NA NA NA
2021-07-20 NA NA NA NA NA
2021-07-21 NA NA NA NA NA
2021-07-22 NA NA NA NA NA
2021-07-23 NA NA NA NA NA 
rm(list = ls())
cleanup(ask = F)



For questions or comments, please contact