---
title: "Temporal Mismatch Simulation"
output: html_document
date: "2026-04-06"
author: "Olivia Cronin-Golomb"
contact: "olivia.cronin-golomb@whoi.edu"
---

```{r setup, include=FALSE} 
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

#Set Up Environment
## Load libraries
```{r}
rm(list=ls()) #clear environment

library(tidyverse) 
library(RColorBrewer)
library(Metrics)
library(data.table)
library(cowplot)
library(gstat)
library(rlist)
library(scales)
library(ggforce)
library(extrafont)
library(dplyr)
library(ggplot2)
library(lubridate)
library('ggpubr')
library(grid)
```

## Set input and output directories
```{r}
#Set working directory
wd <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage1/Data/Processed/"

#Output directory
Paper_figures_filepath <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage1/Outputs/Simulation/"

Extra_figures_filepath <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage1/Outputs/Extra_Figures/"

#Table_filepath <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage/Outputs/Tables/Chlorophyll/"

```

#Select only days with >10 30min intervals
```{r}
#Read in 30min chl data
Chl <- read.csv(paste0(wd, "Chlorophyll/chlData30min.csv"))

#Make character date into a datetime object. Extract the date, time, year, and day. Add an index column.
Chl$Datetime <- parse_date_time(as.character(Chl$Datetime), order = "%y-%m-%d %H:%M:%S", tz = 'EST') 
Chl$Date <- format(Chl$Datetime,"%Y-%m-%d")
Chl$Time <- format(Chl$Datetime,"%H:%M")

Chl$year <- as.numeric(format(Chl$Datetime,"%Y"))
Chl$day <- as.numeric(format(Chl$Datetime,"%d"))
Chl$Index <- seq(1,nrow(Chl),1)

#Only select dates with days of mostly full sampling
allChl <- Chl %>%
  group_by(Date) %>%
  summarise(count = n())
allChl <- allChl[allChl$count <= max(allChl$count) & allChl$count >= max(allChl$count) - 1 ,]
allChl <- merge(Chl, allChl, by = "Date")
allChl <- allChl[order(allChl$Index),]
```

#Extract the noon sample, as that will be the field sample proxy that the values from the increasing time intervals will be compared against.
```{r}
Chl_noon <- allChl[allChl$Time == "12:00", ]
Chl_noon <- Chl_noon[,c(1,3)]
allChl <- merge(allChl, Chl_noon, by = "Date")
names(allChl)[names(allChl) == 'median_chl_mgL.x'] <- 'median_chl_mgL'
names(allChl)[names(allChl) == 'median_chl_mgL.y'] <- 'median_chl_noon'
```

# One day simulation

## Median chl with max and min for each half hour interval

### Function to calculate the stats in each half hour time interval
```{r}
timeInt <- as.data.frame(unique(allChl$Time))
names(timeInt) <- c("Time")
timeInt$numTime <- seq(9,15,0.5)

Daily_chl <- function(loc_name, interval){
  Hour_range <- c(12 - interval/2, 12 + interval/2)
  targetTime <- timeInt[timeInt$numTime >= Hour_range[1] & timeInt$numTime <= Hour_range[2],]
  targetChl <- merge(allChl, targetTime, by = "Time")
  targetChl <- targetChl[order(targetChl$Index),]
  
  locChl <- targetChl[targetChl$loc == loc_name,]
  Chl_range <- locChl %>%
    group_by(Date) %>%
    summarise(max = max(median_chl_mgL, na.rm=TRUE), 
              min = min(median_chl_mgL, na.rm=TRUE),
              range = max(median_chl_mgL, na.rm=TRUE) - min(median_chl_mgL, na.rm=TRUE),
              mean = mean(median_chl_mgL, na.rm = TRUE),
              median = median(median_chl_mgL, na.rm= TRUE),
              sd = sd(median_chl_mgL, na.rm = TRUE),
              quant_10 = quantile(median_chl_mgL, probs = c(.1)),
              quant_90 = quantile(median_chl_mgL, probs = c(.9)),
              count = mean(count),
              MAE = mae(median_chl_mgL, median_chl_noon),
              Chl_noon = mean(median_chl_noon))
  Chl_range$Loc <- loc_name
  Chl_range$Date <- as.Date(Chl_range$Date, format = "%Y-%m-%d")
  return(Chl_range)
}
```

### Plot for each half hour time interval
```{r}
Target_hours <- seq(1, 6, 1)
intDF <- data.frame()
my_legend <- NULL
fig_panels <- list("a","x","x","x","x","b")

# Use identical y label text for all plots, make it transparent after the first
ylab_text_const <- "Simulated in situ sample chl-a \nconcentration range (µg L⁻¹)"

for (i in seq_along(Target_hours)) {

  # Get data per interval
  Chl_NH  <- Daily_chl("jlnewhope",  Target_hours[i])
  Chl_WO  <- Daily_chl("jlwhiteoak", Target_hours[i])
  Chl_Haw <- Daily_chl("jlhaw",      Target_hours[i])

  # Combine and assign categories
  combineChl <- as.data.frame(rbind(Chl_NH, Chl_WO, Chl_Haw))
  combineChl$timeInt <- Target_hours[i] / 2

  combineChl$Category[combineChl$Chl_noon >= 30] <- "high"
  combineChl$Category[combineChl$Chl_noon < 30 & combineChl$median > 7] <- "medium"
  combineChl$Category[combineChl$Chl_noon <= 7] <- "low"

  # Stats (median vs Chl_noon)
  Bias <- formatC(mean(combineChl$mean - combineChl$Chl_noon, na.rm = TRUE),
                  format = "f", digits = 2)
  CV   <- formatC(100 * sd(combineChl$mean, na.rm = TRUE) /
                     mean(combineChl$mean, na.rm = TRUE),
                  format = "f", digits = 2)
  nz   <- is.finite(combineChl$Chl_noon) & combineChl$Chl_noon > 0
  MAPE <- formatC(100 * mean(abs(combineChl$mean[nz] - combineChl$Chl_noon[nz]) /
                             combineChl$Chl_noon[nz], na.rm = TRUE),
                  format = "f", digits = 2)

  intDF <- rbind(intDF, combineChl)

  # Build plot
  int_plot <- ggplot(combineChl, aes(y = median, x = Chl_noon, color = factor(Category))) +
    geom_point(alpha = 0.8, size = 2) +
    geom_errorbar(aes(ymin = min, ymax = max), width = 1, alpha = 0.5, show.legend = FALSE) +
    geom_abline(linetype = "longdash") +
    geom_vline(xintercept = 7,  linetype = 2) +
    geom_vline(xintercept = 30, linetype = 2) +
    scale_color_manual(
      name = NULL, 
      breaks = c("low", "medium", "high"),
      labels = c("Oligotrophic (≤7 µg L⁻¹)", "Eutrophic (7–30 µg L⁻¹)", "Hypereutrophic (≥30 µg L⁻¹)"),
      values = c("low" = "#fc8961", "medium" = "#b73779", "high" = "#51127c")
    ) +
    scale_y_continuous(breaks = seq(0, 100, 20), limits = c(0, 110)) +
    scale_x_continuous(breaks = seq(0, 100, 20), limits = c(0, 110)) +
    labs(
      x = "Simulated satellite sample chl-a \nconcentration at time zero (µg L⁻¹)",
      y = ylab_text_const
    ) +
    coord_fixed(ratio = 1, xlim = c(0, 110), ylim = c(0, 110)) +
    theme_classic(base_family = "Times New Roman") +
    theme(
      axis.text.x = element_text(size = 18),
      axis.text.y = element_text(size = 18),
      axis.title.x = element_text(size = 20),
      axis.title.y = element_text(
        size = 20,
        margin = margin(r = 15),
        color = if (i == 1) "black" else "transparent"
      ),
      legend.text  = element_text(size = 20),
      legend.title = element_blank(),
      legend.spacing.y = unit(0.25, 'cm'),
      legend.position = c(0.5, 0.5),
      legend.background = element_rect(fill = 'transparent'),
      panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      text = element_text(family = "Times New Roman"),
      plot.margin  = margin(t = 10, r = 10, b = 10, l = 10)
    )

  # Grab legend once before removing it from the per-plot output
  if (is.null(my_legend)) {
    my_legend <- get_legend(int_plot)
  }
  int_plot <- int_plot + theme(legend.position = "none")

  int_plot <- int_plot +
    annotate("text", x = 62, y = 1, hjust = 0, vjust = 0,
             label = paste0(
               "Bias = ", Bias, " µg L⁻¹\n",
               # "CV = ", CV, "%\n",          
               "MAPE = ", MAPE, "%"
             ),
             size = 6, family = "Times New Roman")

  # Save
  ggsave(path = Paper_figures_filepath,
         filename = paste0("5",as.character(fig_panels[i]),"_chl_Interval", as.character(Target_hours[i] / 2), "hrs.png"),
         device = "png", plot = int_plot,
         width = 6, height = 6, dpi = 600, units = "in")
  
  #Optional if you want pdf for easy editing in adobe illustrator later
  #ggsave(filename = file.path(Paper_figures_filepath,  paste0("chl_Interval", as.character(Target_hours[i] / 2), "hrs.pdf")),
         #plot = int_plot,
         #device = cairo_pdf,
         #width = 6, height = 6, units = "in")
}

# find files that start with '5x'
files_to_move <- list.files(
  path = Paper_figures_filepath,
  pattern = "^5x",
  full.names = TRUE
)

# move them
file.rename(
  from = files_to_move,
  to = file.path(Extra_figures_filepath, basename(files_to_move))
)

```

### Create legend
```{r}
df <- data.frame(
  x = 1:3,
  y = 1:3,
  cat = factor(c("oligo", "eutro", "hyper"),
               levels = c("oligo", "eutro", "hyper")) 
)

p <- ggplot(df, aes(x, y, color = cat)) +
  geom_point(size = 6) +
scale_color_manual(
    breaks = c("oligo", "eutro", "hyper"),
    labels = c(
      expression(Oligotrophic~"(" * "\u2264"~7~mu*"g L"^{-1} * ")"),
      expression(Eutrophic~"(7-30"~mu*"g L"^{-1} * ")"),
      expression(Hypereutrophic~"(" * "\u2265"~30~mu*"g L"^{-1} * ")")
    ),
    values = c(
      "oligo" = "#fc8961",
      "eutro" = "#b73779",
      "hyper" = "#51127c"
    ),
    name = NULL
  ) +

  guides(color = guide_legend(override.aes = list(size = 2))) +
  theme_minimal(base_family = "Times New Roman") +
  theme(
    legend.text       = element_text(size = 12),
    legend.position   = "right",
    legend.spacing.y  = unit(0.4, "cm"),
    legend.background = element_rect(fill = "white", color = "white")
  )

my_legend <- get_legend(p)

png(
  filename = file.path(Paper_figures_filepath, "5_legend.png"),
  width = 2.7, height = 1.3, units = "in", res = 600,
  type = "cairo"   # important
)
grid.newpage()
grid.draw(my_legend)
dev.off()

```

## Plot MAD for each half hour time interval

### Regular
```{r}
residStats <- data.frame()

for (i in 1:length(Target_hours)){

  Chl_NH <- Daily_chl("jlnewhope", Target_hours[i])
  Chl_WO <- Daily_chl("jlwhiteoak", Target_hours[i])
  Chl_Haw <- Daily_chl("jlhaw", Target_hours[i])
  
  #Combine
  combineChl <- as.data.frame(rbind(Chl_NH, Chl_WO, Chl_Haw))
  combineChl$timeInt <- Target_hours[i]/2
  
  # >95 residuals
  combineChl100up <- combineChl[combineChl$Chl_noon >= 30,]
  lm_100up <- lm(combineChl100up$median - combineChl100up$Chl_noon ~ 0)
  resid100up_stats <- as.data.frame(as.list(summary(resid(lm_100up))))
  resid100up_stats$count <- nrow(combineChl100up)
  resid100up_stats$MAD <- sum(abs(combineChl100up$median - combineChl100up$Chl_noon))/nrow(combineChl100up)
  
  # 10 - 95 residuals
  combineChl10_100 <- combineChl[combineChl$Chl_noon < 30 & combineChl$Chl_noon > 7,]
  lm_10_100 <- lm(combineChl10_100$mean - combineChl10_100$Chl_noon ~ 0)
  resid10_100_stats <- as.data.frame(as.list(summary(resid(lm_10_100))))
  resid10_100_stats$count <- nrow(combineChl10_100)
  resid10_100_stats$MAD <- sum(abs(combineChl10_100$median - combineChl10_100$Chl_noon))/nrow(combineChl10_100)
  
  # < 10 residuals
  combineChl10down <- combineChl[combineChl$Chl_noon <= 7,]
  lm_10down <- lm(combineChl10down$mean - combineChl10down$Chl_noon ~ 0)
  resid10down_stats <- as.data.frame(as.list(summary(resid(lm_10down))))
  resid10down_stats$count <- nrow(combineChl10down)
  resid10down_stats$MAD <- sum(abs(combineChl10down$median - combineChl10down$Chl_noon))/nrow(combineChl10down)
  
  #Combine into df
  residDF <- as.data.frame(rbind(resid100up_stats, resid10_100_stats, resid10down_stats))
  residDF$category <- c("lm100", "lm10_100", "lm10")
  residDF$timeInt <- Target_hours[i]
  
  residStats <- rbind(residStats, residDF)
}

lm10_start <- c(0,0,0,0,0,0,1,0, "lm10", 0)
lm10_100_start <- c(0,0,0,0,0,0,1,0, "lm10_100", 0)
lm100_start <- c(0,0,0,0,0,0,1,0, "lm100", 0)

start_df <- as.data.frame(rbind(lm10_start,lm10_100_start,lm100_start))
names(start_df) <- names(residStats)
col.num <- c("Min.", "X1st.Qu.", "Median", "Mean", "X3rd.Qu.", "Max.", "count", "MAD", "timeInt")
start_df[col.num] <- sapply(start_df[col.num],as.numeric)

residStats <- as.data.frame(rbind(residStats, start_df))

gg_mad <- ggplot(residStats, aes(x = timeInt / 2, y = MAD, color = category)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(se = FALSE, method = "loess", formula = "y ~ x", linewidth = 1) +
  labs(
    x = "No. of hours since sampling \n",
    y = "\nMean absolute deviation (\u03bcg/L)"
  ) +
  scale_color_manual(
    name = NULL,
    labels = c("Oligotrophic (≤ 7 µg L⁻¹)", "Eutrophic (7–30 µg L⁻¹)", "Hypereutrophic ( ≥30 µg L⁻¹)"),
    values = c("lm10" = "#fc8961", "lm10_100" = "#b73779", "lm100" = "#51127c")
  ) +
  scale_x_continuous(breaks = seq(0, 3, 0.5), limits = c(0, 3.2)) +
  scale_y_continuous(breaks = seq(0, 4, 1),   limits = c(0, 4)) +
  coord_fixed(
    xlim  = c(0, 3.2),
    ylim  = c(0, 4),
    ratio = (3.2 - 0) / (4 - 0) 
  ) +
  theme_classic(base_family = "Times New Roman") +
  theme(
    axis.text.x  = element_text(size = 20),
    axis.text.y  = element_text(size = 20),
    axis.title.x = element_text(size = 22),
    axis.title.y = element_text(size = 24, margin = margin(r = 15), color = "black"),
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.margin  = margin(t = 10, r = 10, b = 10, l = 10)
  )

gg_mad

ggsave(
  path = Extra_figures_filepath,
  filename = "MAD_timeints.png",
  device = "png",
  plot = gg_mad,
  width = 7.5, height = 7.5, dpi = 600, units = "in"
)

```

### Log10(MAD) 
```{r}
residStats_log <- data.frame()

for (i in 1:length(Target_hours)){

  Chl_NH <- Daily_chl("jlnewhope", Target_hours[i])
  Chl_WO <- Daily_chl("jlwhiteoak", Target_hours[i])
  Chl_Haw <- Daily_chl("jlhaw", Target_hours[i])
  
  #Combine
  combineChl <- as.data.frame(rbind(Chl_NH, Chl_WO, Chl_Haw))
  combineChl$timeInt <- Target_hours[i]/2
  combineChl <- combineChl[combineChl$median > 0 & combineChl$Chl_noon > 0, ]
  
  # >95 residuals
  combineChl100up <- combineChl[combineChl$Chl_noon >= 30,]
  lm_100up <- lm(combineChl100up$median - combineChl100up$Chl_noon ~ 0)
  resid100up_stats <- as.data.frame(as.list(summary(resid(lm_100up))))
  resid100up_stats$count <- nrow(combineChl100up)
  resid100up_stats$MAD <- sum(abs(log10(combineChl100up$median) - log10(combineChl100up$Chl_noon)))/nrow(combineChl100up)
  
  # 10 - 95 residuals
  combineChl10_100 <- combineChl[combineChl$Chl_noon < 30 & combineChl$Chl_noon > 7,]
  lm_10_100 <- lm(combineChl10_100$median - combineChl10_100$Chl_noon ~ 0)
  resid10_100_stats <- as.data.frame(as.list(summary(resid(lm_10_100))))
  resid10_100_stats$count <- nrow(combineChl10_100)
  resid10_100_stats$MAD <- sum(abs(log10(combineChl10_100$median) - log10(combineChl10_100$Chl_noon)))/nrow(combineChl10_100)
  
  # < 10 residuals
  combineChl10down <- combineChl[combineChl$Chl_noon <= 7,]
  lm_10down <- lm(combineChl10down$median - combineChl10down$Chl_noon ~ 0)
  resid10down_stats <- as.data.frame(as.list(summary(resid(lm_10down))))
  resid10down_stats$count <- nrow(combineChl10down)
  resid10down_stats$MAD <- sum(abs(log10(combineChl10down$median) - log10(combineChl10down$Chl_noon)))/nrow(combineChl10down)
  
  #Combine into df
  residDF <- as.data.frame(rbind(resid100up_stats, resid10_100_stats, resid10down_stats))
  residDF$category <- c("lm100", "lm10_100", "lm10")
  residDF$timeInt <- Target_hours[i]
  
  residStats_log <- rbind(residStats_log, residDF)
}

lm10_start <- c(0,0,0,0,0,0,1,0, "lm10", 0)
lm10_100_start <- c(0,0,0,0,0,0,1,0, "lm10_100", 0)
lm100_start <- c(0,0,0,0,0,0,1,0, "lm100", 0)

start_df <- as.data.frame(rbind(lm10_start,lm10_100_start,lm100_start))
names(start_df) <- names(residStats_log)
col.num <- c("Min.", "X1st.Qu.", "Median", "Mean", "X3rd.Qu.", "Max.", "count", "MAD", "timeInt")
start_df[col.num] <- sapply(start_df[col.num],as.numeric)

residStats_log <- as.data.frame(rbind(residStats_log, start_df))

gg_mad_log <- ggplot(residStats_log, aes(x = timeInt / 2, y = MAD, color = category)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(se = FALSE, method = "loess", formula = "y ~ x", linewidth = 1) +
  labs(
    x = "No. of hours since sampling \n",
    y = "\nlog₁₀(mean absolute deviation) (µg L⁻¹)"
  ) +
  scale_color_manual(
    name = NULL,
    labels = c("Oligotrophic (≤ 7 µg L⁻¹)", "Eutrophic (7–30 µg L⁻¹)", "Hypereutrophic ( ≥30 µg L⁻¹)"),
    values = c("lm10" = "#fc8961", "lm10_100" = "#b73779", "lm100" = "#51127c")
  ) +
  scale_x_continuous(breaks = seq(0, 3, 0.5), limits = c(0, 3.2)) +
  scale_y_continuous(breaks = seq(0, 0.05, 0.01), limits = c(0, 0.05)) +
  coord_fixed(
    xlim  = c(0, 3.2),
    ylim  = c(0, 0.05),
    ratio = (3.2 - 0) / (0.05 - 0)   # 64 -> square panel in pixels
  ) +
  theme_classic(base_family = "Times New Roman") +
  theme(
    axis.text.x  = element_text(size = 18),
    axis.text.y  = element_text(size = 18),
    axis.title.x = element_text(size = 20),
    axis.title.y = element_text(size = 20, margin = margin(r = 15), color = "black"),
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.margin  = margin(t = 10, r = 10, b = 10, l = 10)
  )

gg_mad_log

ggsave(
  path = Paper_figures_filepath,
  filename = "5c_MAD_log_timeints.png",
  device = "png",
  plot = gg_mad_log,
  width = 6, height = 6, dpi = 600, units = "in"
)

```

## Get counts for each trophic category
```{r}
aggregate(count ~ category, data = residStats_log, sum)

```


# Multi-day simulation Set Up
## Calculate longest string of consecutive days
```{r}
Target_hours <- 6
Chl_NH <- Daily_chl("jlnewhope", Target_hours)
Chl_WO <- Daily_chl("jlwhiteoak", Target_hours)
Chl_Haw <- Daily_chl("jlhaw", Target_hours)

consecDays <- function(intData){
  intData %>% 
    group_by(cumsum(c(0, diff(Date) - 1))) %>%
    summarise(start_date = first(Date), 
              end_date = last(Date),
              ndays = n()) %>%
    select(start_date, end_date, ndays)
}

Chl_NH_days <- consecDays(Chl_NH); Chl_NH_days$Loc <- "jlnewhope"
Chl_WO_days <- consecDays(Chl_WO); Chl_WO_days$Loc <- "jlwhiteoak"
Chl_Haw_days <- consecDays(Chl_Haw); Chl_Haw_days$Loc <- "jlhaw"

combineDays <- rbind(Chl_NH_days, Chl_WO_days, Chl_Haw_days)
maxDays <- combineDays[combineDays$ndays == max(combineDays$ndays),]
maxDays$mid_date <- maxDays$start_date + floor((maxDays$end_date-maxDays$start_date)/2)
allChl$Date <- as.Date(format(allChl$Datetime,"%Y-%m-%d"))

targetDays <- allChl[allChl$Date >= maxDays$start_date & allChl$Date <= maxDays$end_date,]
row.names(targetDays) <- NULL

times <- c("09:00", "09:30", "10:00", "10:30", "11:00", "11:30", "12:00", "12:30", "13:00", "13:30", "14:00", "14:30", "15:00")

times_list <- rep(times, length(unique(targetDays$Date)))
dates_list <- as.data.frame(rep(unique(targetDays$Date), each = length(times)))
dates_list$time <- times_list
names(dates_list) <- c("Date", "Time")
dates_list$Date <- as.Date(format(dates_list$Date,"%Y-%m-%d"))

targetAll <- merge(dates_list, targetDays, by = c("Date", "Time"), all = T)
row.names(targetAll) <- NULL
targetAll$Datetime <- parse_date_time(paste(as.character(targetAll$Date), as.character(targetAll$Time)), order = "%y/%m/%d %H:%M", tz = 'EST') 
```

## Plot Consecutive Days Chlorophyll Time Series
```{r}
#Field samples
target_midDay <- targetAll[targetAll$Date == "2017-08-14",]
target_midDay$Index <- row.names(target_midDay)
row.names(target_midDay) <- NULL
middaySplit <- target_midDay[c(1,13),c(1:3)]
target_midpts <- target_midDay[target_midDay$Time == "10:00" | target_midDay$Time == "12:00" | target_midDay$Time == "14:00",]

midptIndex <- as.numeric(target_midpts$Index)

#Colors
TSSColors <- brewer.pal(7,"Greys")
TSSColors <- TSSColors[2:7]
names(TSSColors) <- c("Calender day","1 day", "2 days", "3 days", "1 week", "2 weeks")

#Plot time intervals
hourInts <- c(seq(480/2, 0, -1), seq(1,480/2, 1))
targetAll$timeSince <- hourInts
targetInts = subset(targetAll, subset = timeSince %in% c(0,6,13,26,39,91,169))
row.names(targetInts) <- NULL

targetInts$Group <- c("2 weeks","1 week","3 days","2 days","1 day","Calender day","Noon","Calender day","1 day","2 days","3 days","1 week","2 weeks")
targetInts$Color <- c(rlist::list.reverse(TSSColors),"black", TSSColors)
targetInts2 <- targetInts[-c(7),]

row.names(targetInts2) <- NULL

LineUp <- targetInts2$Datetime[1:6] - (60*60*3)
LineDown <- targetInts2$Datetime[7:12] + (60*60*3)
Lines <- c(LineUp, LineDown)
targetInts2$lines <- Lines

targetAll2 <- targetAll[-midptIndex,]

gg_chl_simple <- ggplot(na.omit(targetAll), aes(x = Datetime, y = median_chl_mgL)) +
  annotate("rect",
         xmin = as.POSIXct("2017-08-14 06:30:00", tz = "EST"),
         xmax = as.POSIXct("2017-08-14 19:00:00", tz = "EST"),
         ymin = -Inf,
         ymax = Inf,
         fill = "grey50",
         alpha = 0.5) + 
  geom_vline(data = targetInts2, aes(xintercept = lines, group = as.factor(Group), color = as.factor(Group)),
             show.legend = TRUE, linetype = "solid", key_glyph = "path") +
  scale_color_manual(
    values = TSSColors,
    breaks = c("Calender day", "1 day", "2 days", "3 days", "1 week", "2 weeks"),
    labels = c("Calender day", "1 day", "2 days", "3 days", "1 week", "2 weeks"),
    guide = guide_legend(order = 2, ncol = 2, title = NULL)
  ) +
  ggnewscale::new_scale_color() +
  geom_point(size = 2) +
  labs(
    x = NULL,
    y = "Chl-a concentration (µg L⁻¹)",
    color = ""
  ) +
  geom_point(data = target_midpts,
           aes(x = Datetime, y = median_chl_mgL, fill = Time),
           color = "black", shape = 24, size = 4, show.legend = FALSE) +
scale_fill_manual(values = c('#35b779', '#31688e', '#440154'), guide = "none") +
  scale_y_continuous(breaks = seq(0, 50, 10), limits = c(10, 43), expand = c(0,0)) +
  scale_x_datetime(
    breaks = c(as.POSIXct("2017-08-01 10:00:00 EST"),
               as.POSIXct("2017-08-07 10:00:00 EST"),
               as.POSIXct("2017-08-14 12:00:00 EST"),
               as.POSIXct("2017-08-21 16:00:00 EST"),
               as.POSIXct("2017-08-27 16:00:00 EST")),
    labels = date_format("%b %d", tz = "EST"),
    expand = c(0, 0),
    limits = c(as.POSIXct("2017-07-26 00:00:00 EST"),
               as.POSIXct("2017-09-02 23:59:00 EST"))) +
  annotate("text", label = "a.", x = as.POSIXct("2017-07-26 13:00:00 EST"), y = 42.5,
           hjust = 0, vjust = 1, size = 8, family = "Times") +
  #ggtitle(expression("Chlorophyll Time Series"^phantom("th") * "")) +
  theme_classic() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    axis.title.y = element_text(size = 20),
    plot.title = element_text(hjust = 0.5, size = 20),
    legend.text = element_text(size = 18),
    #legend.title = element_text(size = 17),
    #legend.background = element_rect(fill = 'transparent'),
    legend.position = c(0.2,0.1),
    text = element_text(family = "Times New Roman")
  )


gg_chl_simple

ggsave(path = Paper_figures_filepath, filename = paste0("4a_Interval_Chl_TimeSeries.png"), device = "png", plot = gg_chl_simple, width = 9.5, height = 6.5, dpi = 600, units = "in")

#Optional if you want pdf for easy editing in adobe illustrator later
#ggsave(
#  filename = file.path(Paper_figures_filepath, paste0("Interval_Chl_TimeSeries.pdf")),
#  plot = gg_chl_simple,
#  device = cairo_pdf,
#  width = 9.5, height = 6.5, units = "in"
#)
```
## Plot Midpoint Day Aug 14
```{r}
target_midDay2 <- target_midDay[ ! target_midDay$Index %in% as.character(midptIndex), ]

# add hour_offset column relative to that row’s Date
add_hour_offset <- function(df, date_col = "Date", dt_col = "Datetime", tz_fallback = "EST") {
  tz_use <- tryCatch(attr(df[[dt_col]], "tzone"), error = function(e) NULL)
  if (is.null(tz_use) || tz_use == "") tz_use <- tz_fallback
  noon_vec <- as.POSIXct(paste(df[[date_col]], "12:00:00"), tz = tz_use)
  df$hour_offset <- as.numeric(difftime(df[[dt_col]], noon_vec, units = "hours"))
  df
}

target_midDay2 <- add_hour_offset(target_midDay2)
target_midpts  <- add_hour_offset(target_midpts)
middaySplit    <- add_hour_offset(middaySplit)

gg_midday <- ggplot(na.omit(target_midDay2), aes(x = hour_offset, y = median_chl_mgL)) +
  # Shaded rectangle over satellite window 
  annotate("rect",
           xmin = min(middaySplit$hour_offset, na.rm = TRUE),
           xmax = max(middaySplit$hour_offset, na.rm = TRUE),
           ymin = -Inf, ymax = Inf, alpha = 0.3, fill = "grey50") +
  geom_vline(data = middaySplit, aes(xintercept = hour_offset),
             color = TSSColors[1], linetype = "solid", show.legend = FALSE) +
  # Error bars (±20% around median) and midpoints
  geom_errorbar(data = target_midpts,
                aes(x = hour_offset,
                    ymin = median_chl_mgL * 0.8, ymax = median_chl_mgL * 1.2, color = Time),
                width = 0.05, position = position_dodge(0.05), show.legend = FALSE) +
  geom_point(data = target_midpts,
             aes(x = hour_offset, fill = Time),
             color = "black", shape = 24, size = 4, show.legend = TRUE) +
  scale_fill_manual(
    labels = c("28.26 ± 5.65", "34.81 ± 6.96", "27.33 ± 5.47"),
    values = c('#35b779', '#31688e', '#440154'),
    guide = guide_legend(title = NULL)
  ) +
  scale_color_manual(values = c('#35b779', '#31688e', '#440154'), guide = "none") +
  labs(x = "No. of hours from time zero (noon)", y = NULL, fill = NULL) +
  scale_x_continuous(breaks = seq(-3, 3, 1), limits = c(-3.2, 3.2), expand = c(0, 0)) +
  scale_y_continuous(breaks = seq(10, 40, 10), limits = c(10, 43), expand = c(0, 0)) +
  annotate("text", label = "b.", x = -3.1, y = 42.5,
           hjust = 0, vjust = 1, size = 8, family = "Times") +
  theme_classic() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    axis.text.x = element_text(size = 18),
    axis.text.y = element_text(size = 18),
    axis.title  = element_text(size = 20),
    legend.text = element_text(size = 18),
    legend.title = element_blank(),
    legend.position = c(0.2, 0.1),
    legend.background = element_rect(fill = 'transparent'),
    plot.title = element_text(hjust = 0.5, size = 20),
    text = element_text(family = "Times New Roman")
  ) +
  guides(fill = guide_legend(title.position = "top", title.hjust = 0.5))

gg_midday

ggsave(path = Paper_figures_filepath,
       filename = "4b_Aug14_TimeSeries.png",
       device = "png", plot = gg_midday,
       width = 7, height = 6.8, dpi = 600, units = "in")


```
## Longest consectutive day stats
```{r}
df <- na.omit(targetAll)

# keep only positive chl values
df <- df[df$median_chl_mgL > 0, ]

mean_chl <- round(mean(df$median_chl_mgL), 2)
max_chl  <- round(max(df$median_chl_mgL), 2)
min_chl  <- round(min(df$median_chl_mgL), 2)
cv       <- sd(df$median_chl_mgL) / mean(df$median_chl_mgL) * 100

cat("Mean Chlorophyll-a:", mean_chl, "mg/L\n")
cat("Max Chlorophyll-a:", max_chl, "mg/L\n")
cat("Min Chlorophyll-a:", min_chl, "mg/L\n")
cat("36 day Coef of Var:", round(cv, 2), "\n")

# Calculate percent change
initial_value <- df$median_chl_mgL[1]      # First positive value
final_value   <- df$median_chl_mgL[nrow(df)]  # Last positive value
percent_change <- ((final_value - initial_value) / initial_value) * 100

cat("The percent change is", round(percent_change, 2), "%")
```

## Aug 14th stats
```{r}
df_sub <- targetAll %>%
  filter(Date == "2017-08-17",
         Time >= "09:00" & Time <= "15:00",
         median_chl_mgL > 0)

mean_chl <- round(mean(df_sub$median_chl_mgL, na.rm = TRUE), 2)
max_chl  <- round(max(df_sub$median_chl_mgL, na.rm = TRUE), 2)
min_chl  <- round(min(df_sub$median_chl_mgL, na.rm = TRUE), 2)
cv       <- sd(df_sub$median_chl_mgL, na.rm = TRUE) / mean(df_sub$median_chl_mgL, na.rm = TRUE) * 100

cat("Mean Chlorophyll-a:", mean_chl, "mg/L\n")
cat("Max Chlorophyll-a:", max_chl, "mg/L\n")
cat("Min Chlorophyll-a:", min_chl, "mg/L\n")
cat("Coef of Var:", round(cv, 2), "%\n")

# Percent change
initial_value <- df_sub$median_chl_mgL[1]
final_value   <- df_sub$median_chl_mgL[nrow(df_sub)]
percent_change <- ((final_value - initial_value) / initial_value) * 100

cat("The percent change is", round(percent_change, 2), "%")
```

# Simulation

## Plot 6hr Variogram
```{r}
df_all <- targetAll %>%
  mutate(Datetime = as.POSIXct(Datetime, tz = "EST"),
         value = as.numeric(median_chl_mgL)) %>%
  filter(is.finite(Datetime), is.finite(value))

# Adaptive variogram: bin all pairwise diffs by lag
make_adaptive_variogram <- function(df, bin_width_hours = 0.5, max_lag_hours = 72, min_pairs = 3) {
  cmb <- utils::combn(nrow(df), 2)
  i <- cmb[1, ]; j <- cmb[2, ]
  
  lag_hr  <- abs(as.numeric(difftime(df$Datetime[j], df$Datetime[i], units = "hours")))
  semivar <- 0.5 * (df$value[j] - df$value[i])^2
  
  out <- data.frame(lag_hr = lag_hr, semivar = semivar)
  
  out %>%
    filter(lag_hr <= max_lag_hours) %>%              # keep within max lag
    mutate(lag_bin = floor(lag_hr / bin_width_hours) * bin_width_hours) %>%
    group_by(lag_bin) %>%
    summarise(
      gamma = mean(semivar, na.rm = TRUE),
      n_pairs = n(),
      .groups = "drop"
    ) %>%
    filter(n_pairs >= min_pairs) %>%
    arrange(lag_bin)
}

make_adaptive_variogram <- function(df, bin_width_hours = 1, max_lag_hours = 72, min_pairs = 3) {
  if (nrow(df) < 2) return(tibble(lag_bin = numeric(), gamma = numeric(), n_pairs = numeric(), se = numeric()))
  
  cmb <- utils::combn(nrow(df), 2)
  i <- cmb[1, ]; j <- cmb[2, ]
  
  lag_hr  <- as.numeric(difftime(df$Datetime[j], df$Datetime[i], units = "hours"))
  semivar <- 0.5 * (df$value[j] - df$value[i])^2
  
  tibble(lag_hr = abs(lag_hr), semivar = semivar) %>%
    filter(lag_hr <= max_lag_hours) %>%
    mutate(lag_bin = round(lag_hr / bin_width_hours) * bin_width_hours) %>%
    group_by(lag_bin) %>%
    summarise(
      gamma = mean(semivar, na.rm = TRUE),
      n_pairs = n(),
      se = sd(semivar, na.rm = TRUE) / sqrt(n_pairs),
      .groups = "drop"
    ) %>%
    filter(n_pairs >= min_pairs) %>%
    arrange(lag_bin)
}


# Short-lag: every 0.5 hr up to 6 hr
vg_adapt_short <- make_adaptive_variogram(df_all, bin_width_hours = 0.5, max_lag_hours = 6)

# Plot short-lag
p_short <- ggplot(vg_adapt_short, aes(x = lag_bin, y = gamma)) +
  geom_point(size = 2, color = "blue") +
  geom_line(color = "blue") +
  geom_text(aes(label = n_pairs), vjust = -1, size = 16/.pt) +  # match legend font size
  scale_x_continuous(breaks = seq(0, 6, 1)) +
  ylim(2,10) +
  labs(
    x = "Lag (hours)",
    y = expression(paste("Semivariance, ", gamma(h)))#,    title = "One-day Variogram (0.5–6 hr)"
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(size = 17, angle = 0, hjust = 0.5),
    axis.text.y = element_text(size = 17),
    axis.title.x = element_text(size = 18),
    axis.title.y = element_text(size = 18),
    plot.title = element_text(hjust = 0.5, size = 20),
    legend.text = element_text(size = 16),
    legend.title = element_text(size = 17),
    legend.spacing.y = unit(0.25, 'cm'),
    legend.position = "right",
    legend.background = element_rect(fill = "transparent"),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(family = "Times New Roman")
  )

ggsave(
  filename = paste0(Paper_figures_filepath, "6a_Chl_Variogram_6hr.png"),
  plot = p_short,
  device = "png",
  width = 6.5, height = 5.5, dpi = 300, units = "in"
)

p_short

```

## 10AM, Noon, and 2PM Sample Time Specific Intervals
```{r}
Sample_date <- "2017-08-14"
full_day <- c(0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6)

timeSubset <- function(dayRange, timeRange){
  days <- c(1:dayRange)
  day_starts <- (days*24) - 1
  days_reps <- rep(day_starts, each = length(full_day))
  days_add <- rep(full_day, dayRange)
  days_sum <- as.data.frame(cbind(days_reps, days_add))
  days_sum$days_sum <- days_sum$days_reps + days_sum$days_add
  
  day_ints <- days_sum$days_sum
  
  day_ints <- c(timeRange, day_ints)
  daysDF <- as.data.frame(day_ints)
  daysDF$Index <- seq(1, nrow(daysDF), 1)
  names(daysDF) <- c("Int", "Index")
  return(daysDF)
}

odds <- function(x) subset(x, x %% 2 == 1)
evens <- function(x) subset(x, x %% 2 == 0)

get_intervals <- function(Sample_time){
  hour <-as.numeric(sub("\\:.*", "", Sample_time))
  minute <- as.numeric(sub('.*:', '', Sample_time))
  
  midpoint <- targetAll[targetAll$Date == Sample_date & targetAll$Time == Sample_time,]
  
  morning <- seq(9,hour, 0.5) - 9
  evening <- seq(hour, 15, 0.5) - hour

  target_forward <- targetAll[which(targetAll$Date == Sample_date & targetAll$Time == Sample_time) : nrow(targetAll),]
  target_backwards <- targetAll[1:which(targetAll$Date == Sample_date & targetAll$Time == Sample_time)-1,]
  
  dayRangeF <- (length(unique(target_forward$Date))-1)
  dayRangeB <- (length(unique(target_backwards$Date))-1)
  
  #After sample time
  daysDFpost <- timeSubset(dayRangeF, evening)

  #Before sample time
  daysDFpre <- timeSubset(dayRangeB, morning)
  daysDF2 <- daysDFpre[-1,]
  daysDF2 <- daysDF2[order(daysDF2$Index, decreasing = TRUE), ]
  daysDF2$Index <- -1*daysDF2$Index
  
  #Combine
  daysDF <- as.data.frame(rbind(daysDF2, daysDFpost))
  row.names(daysDF) <- NULL

  targetAll$timePassed <- daysDF$Int
  targetAll$timePassedInt <- daysDF$Index

  #Remove duplicated rows and filter only for fully symmetrical time intervals 
  Duplicated <- as.numeric(abs(length(morning) - length(evening)) -1)
  
  if (length(morning) < length(evening)) {
    hour_intervals <- daysDF[-seq(nrow(daysDF),nrow(daysDF) - Duplicated, -1), ]
    hour_intervals <- sort(unique(abs(hour_intervals$Index)))
    hour_intervals <- odds(hour_intervals)
  } else if (length(morning) > length(evening)) {
    hour_intervals <- daysDF[-seq(1, Duplicated, 1), ]
    hour_intervals <- sort(unique(abs(hour_intervals$Index)))
    hour_intervals <- odds(hour_intervals)
  } else {
    hour_intervals <- sort(unique(abs(daysDF$Index)))
    hour_intervals <- odds(hour_intervals)
  }
  
  names(daysDF) <- c("timePassed", "timePassedInt")
  return_list <- list(hour_intervals, daysDF)
  return(return_list)
}

hour_intervals10 <- get_intervals(Sample_time = "10:00")
hour_intervals12 <- get_intervals(Sample_time = "12:00")
hour_intervals14 <- get_intervals(Sample_time = "14:00")

targetAll$Index <- seq(1,nrow(targetAll),1)
targetAll$loc <- "jlnewhope"
```

## Error Stats for Each Sample Time
```{r}
percent <- 0.2

Sample_time = "10:00"
timeStep <- hour_intervals10[[1]][4]

errorTime <- function(timeStep, Sample_time){
  targetSub <- as.data.frame(cbind(targetAll, data))
  midpt_value <- targetSub[targetSub$Date == Sample_date & targetSub$Time == Sample_time,]
  
  targetTime <- targetSub[targetSub$timePassedInt <= timeStep & targetSub$timePassedInt >= -1*timeStep,]

  intMax = round(max(targetTime$median_chl_mgL, na.rm=TRUE),2)
  intMin = round(min(targetTime$median_chl_mgL, na.rm=TRUE),2)
  intRange = round(max(targetTime$median_chl_mgL, na.rm=TRUE) - min(targetTime$median_chl_mgL, na.rm=TRUE),2)
  intMean = round(mean(targetTime$median_chl_mgL, na.rm = TRUE),2)
  intMedian = round(median(targetTime$median_chl_mgL, na.rm= TRUE),2)
  intSD = round(sd(targetTime$median_chl_mgL, na.rm = TRUE),2)
  intQuant10 = round(as.numeric(quantile(targetTime$median_chl_mgL, probs = c(.1), na.rm = TRUE),2))
  intQuant90 = round(as.numeric(quantile(targetTime$median_chl_mgL, probs = c(.9), na.rm = TRUE),2))
  intCount = as.numeric(nrow(targetTime))
  intNA = as.numeric(sum(is.na(targetTime$median_chl_mgL)))
  intMAE = round(sum(abs(na.omit(targetTime$median_chl_mgL) - na.omit(targetTime$mean_chl_noon))) / nrow(targetTime),2)
  
  TimesUp <- nrow(targetTime[ targetTime$median_chl_mgL > (midpt_value$median_chl_mgL* (1 + percent)), ])
  TimesDown <- nrow(targetTime[ targetTime$median_chl_mgL < (midpt_value$median_chl_mgL*(1 - percent)), ])
  Times = nrow(targetTime) - (TimesUp + TimesDown)
  
  timePassed <- max(targetTime$timePassed)*2
  probability <- (Times/nrow(targetTime)) * 100

  statsList <- c(intMax, intMin, intRange, intMean, intMedian, intSD, intQuant10, intQuant90, intCount, intNA, intMAE, Times, timePassed, probability)
  statsDF <- as.data.frame(t(as.data.frame(statsList)))
  names(statsDF) <- c("Max", "Min", "Range", "Mean", "Median", "SD", "Quant10", "Quant90", "Count", "Count_NA", "MAE", "nTimes", "timePassed", "Probability")
  row.names(statsDF) <- paste0("int", as.character(timeStep))
  return(statsDF)
}
```

## Apply Functions
```{r}
intError10 <- data.frame()
data <- hour_intervals10[[2]]
ints <- hour_intervals10[[1]]

for(i in 1:length(ints)){
  intErr <- errorTime(timeStep = ints[i], Sample_time = "10:00" )
  intErr[sapply(intErr, is.infinite)] <- 100
  intErr$sampleTime <- "10:00"
  intError10 <- rbind(intError10, intErr)
}

intError12 <- data.frame()
data <- hour_intervals12[[2]]
ints <- hour_intervals12[[1]]

for(i in 1:length(ints)){
  intErr <- errorTime(timeStep = ints[i], Sample_time = "12:00" )
  intErr[sapply(intErr, is.infinite)] <- 100
  intErr$sampleTime <- "12:00"
  intError12 <- rbind(intError12, intErr)
}

intError14 <- data.frame()
data <- hour_intervals14[[2]]
ints <- hour_intervals14[[1]]

for(i in 1:length(ints)){
  intErr <- errorTime(timeStep = ints[i], Sample_time = "14:00" )
  intErr[sapply(intErr, is.infinite)] <- 100
  intErr$sampleTime <- "14:00"
  intError14 <- rbind(intError14, intErr)
}

intError2 <- as.data.frame(rbind(intError10, intError12, intError14))
intError <- as.data.frame(rbind(intError10[c(6,11,13:15)], intError12[c(6,11,13:15)], intError14[c(6,11,13:15)]))
```

## Plot Sampling Probabilities
```{r}
daysSegs <- as.data.frame(cbind(c("Calender day", "1 day", "2 days", "3 days", "1 week", "2 weeks"), c(3,24,48,72,7*24,14*24)))
names(daysSegs) <- c("Group", "xint")

intError$Group <- NA
intError$Group[intError$sampleTime == "10:00"] <- "A"
intError$Group[intError$sampleTime == "12:00"] <- "B"
intError$Group[intError$sampleTime == "14:00"] <- "C"

gg_err <- ggplot(intError, aes(x = (timePassed/2), y = Probability, color = as.factor(Group))) + 
  geom_vline(data = daysSegs, aes(xintercept = as.numeric(xint), 
                                  group = factor(Group), 
                                  color = factor(Group)), 
             show.legend = FALSE, 
             linetype = "solid") +
  
  # First legend for time-based grouping
  scale_color_manual(name = "Time since sample",
                     breaks = c("Calender day", "1 day", "2 days", "3 days", "1 week", "2 weeks"),
                     values = TSSColors,
                     guide = guide_legend(order = 2)) +
  
  ggnewscale::new_scale_color() +
  
  # Main data points with second legend for chl-a groups
  geom_point(size = 3, aes(color = Group)) +
  
  scale_color_manual(name = "Chl-a concentration of in situ\nsample ±20% (µg L⁻¹)",
                     breaks = c("A", "B", "C"),
                     labels = c("28.26 ± 5.65 (10 AM)", "34.81 ± 6.96 (12 PM)", "27.33 ± 5.47 (2 PM)"),
                     values = c('#35b779', '#31688e', '#440154'),
                     guide = guide_legend(order = 1)) +
  
  scale_x_continuous(name = "Time since sampling",
                     breaks = c(3, 24, 48, 72, 168, 336),
                     labels = c("Calendar day", "1 day", "2 days", "3 days", "1 week", "2 weeks")) +
  
  scale_y_continuous(name = "Probability of in situ sample match (%)",
                     limits = c(0, 100),
                     breaks = seq(0, 100, 20)) +
  
  theme_classic() +
  theme(
    axis.text.x = element_text(size = 18, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 18),
    axis.title.x = element_text(size = 19),
    axis.title.y = element_text(size = 19),
    plot.title = element_text(hjust = 0.5, size = 20),
    legend.text = element_text(size = 18),
    legend.title = element_text(size = 19),
    legend.spacing.y = unit(0.25, 'cm'),
    legend.position = c(0.8, 0.83),
    legend.background = element_rect(fill = "white"),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    text = element_text(family = "Times New Roman")
  )

gg_err

ggsave(path = Paper_figures_filepath, filename = paste0("6b_Probability_TimeSeries_Simulation.png"), device = "png", plot = gg_err, width = 9.5, height = 6.5, dpi = 300, units = "in")

#Optional if you want pdf for easy editing in adobe illustrator later
#ggsave(
#  filename = file.path(Paper_figures_filepath, "6b_Probability_TimeSeries_Simulation.pdf"),
#  plot = gg_err,
#  device = cairo_pdf,
#  width = 9.5, height = 6.5, units = "in"
#)
```

## Sampling Probabilities Stats
```{r}
prob_10am <- intError[intError$sampleTime == "10:00",]
prob_12pm <- intError[intError$sampleTime == "12:00",]
prob_14pm <- intError[intError$sampleTime == "14:00",]
```
