Oroville_temperature.Rmd

```{r}
library(dplyr)
library(imputeTS)
library(tidyr)
library(ggplot2)
library(cder) #list of sensors here https://cdec.water.ca.gov/reportapp/javareports?name=SensList
```
```{r}
Oroville_temp <- cdec_query("ORI", 25, "D", "2003-01-01", Sys.Date())
Oroville_temp
```


```{r}
#Get Oroville elevation from CDEC, sensor 6 is for elevetation, D is for daily data
Oroville_elev <- cdec_query("ORO", 6, "D", "2003-01-01", Sys.Date())
#making sure the date format is correct
Oroville_elev$DateTime <- as.Date(Oroville_elev$DateTime, format = "%Y/%m/%d")
Oroville_elev



# Define the custom renaming function
#rename_columns <- function(name) {
#  paste0("Temperature (\u00B0F) at ", name, " feet")
#}

# Read the Excel file and rename columns
Oroville <- readxl::read_excel("C:/Users/gdourado/Documents/Oroville_Temperature.xlsx", sheet = "Sheet1")%>%
  rename_with(~ sub(".*([0-9]{3})$", "\\1", .)) %>%    # Extract the last 3 digits
  mutate(across(where(is.numeric), ~ na_interpolation(.)),
         TIME = as.Date(TIME, "%Y-%m-%d"))#%>% # fill gaps with linear interpolation %>%
#rename_with(~ rename_columns(.), .cols = where(is.numeric))   # Apply the custom renaming function
  
Oroville


storage_curve <- readxl::read_excel("C:/Users/gdourado/Documents/Oroville_Temperature.xlsx", sheet = "Sheet2")
storage_curve

Oroville_elev %>% filter(DataFlag == "N")
```

```{r}
#sanity check
ggplot(Oroville%>%  filter(lubridate::year(TIME) == 2011)) + geom_line(aes(x=TIME, y = `630`)) 
#

```

```{r}
#Oroville temperature dataset is missing few days in 2016, let's make the two datasets the same lenght
# Find the missing dates
missing_dates <- as.Date(setdiff(Oroville_elev$DateTime, Oroville$TIME), format = "%Y-%m-%d")

# Create a data frame with missing dates and NA for temperature columns
missing_rows <- data.frame(
  TIME = missing_dates
)

# Add NA columns for each temperature measurement
temperature_columns <- colnames(Oroville)[grep("^[0-9]{3}$", colnames(Oroville))]
for (col in temperature_columns) {
  missing_rows[[col]] <- NA
}

# Combine the missing rows with the original dataset
Oroville_combined <- bind_rows(Oroville, missing_rows)

# Sort the combined dataset by date
Oroville_combined <- Oroville_combined %>% arrange(TIME) %>%
  mutate(across(all_of(temperature_columns), ~ na_interpolation(.))) # Interpolate the temperature columns after adding NA dates

Oroville_combined

Oroville_combined$`Elevation (feet)` <- Oroville_elev$Value
```


```{r}
# Install and load required library
library(imputeTS)

# Function to interpolate extreme temperature values
interpolate_extreme_temps <- function(temps) {
  interpolated_temps <- na_interpolation(temps, option = "spline")
  return(interpolated_temps)
}

# Function to drop extreme values and replace them through interpolation
interpolate_extreme_and_drop <- function(data) {
  # Extract the temperature columns
  temp_cols <- grep("^[0-9]{3}$", colnames(data), value = TRUE)
  
  for (col in temp_cols) {
    # Drop values above 100°F and below 32°F
    data[[col]][data[[col]] > 100 | data[[col]] < 32] <- NA
    
    # Drop values that vary more than 15 degrees from one day to another
    temp_diff <- abs(diff(data[[col]]))
    data[[col]][c(FALSE, temp_diff > 10)] <- NA
    
    # Perform interpolation for temperature columns only
    data[[col]] <- interpolate_extreme_temps(data[[col]])
  }
  
  return(data)
}

# Make a copy of Oroville
Oroville2 <- Oroville_combined

# Apply the interpolation function to your dataset, excluding TIME and Surface columns
Oroville2[ , -which(names(Oroville2) %in% c("TIME", "Surface"))] <- interpolate_extreme_and_drop(Oroville2[ , -which(names(Oroville2) %in% c("TIME", "Surface"))])

# Print the updated dataset
print(Oroville2)

```

```{r}
#sanity check
ggplot(Oroville2%>%
  filter(lubridate::year(TIME) == 2011)) + geom_line(aes(x=TIME, y = `630`))+ coord_cartesian(ylim = c(42,46.5))

#sanity check
ggplot(Oroville%>%
  filter(lubridate::year(TIME) == 2011)) + geom_line(aes(x=TIME, y = `630`)) + coord_cartesian(ylim = c(42,46.5))


ggplot(Oroville) + geom_line(aes(x=TIME, y = `630`)) + coord_cartesian(ylim = c(10,85))


ggplot(Oroville2) + geom_line(aes(x=TIME, y = `630`)) + coord_cartesian(ylim = c(10,85))

```
```{r}
# Define the columns and adjustments
columns <- c(612, 630, 649, 668, 686, 705, 723, 742, 751, 761, 770, 779, 789, 798, 807, 816, 826, 835, 844, 854, 863, 872, 882, 891)
adjustments <- c(-0.3, 0.6, -0.3, 0.6, -0.3, 0.5, -0.3, 0.1, -0.7, 0.2, -0.8, 0.1, -0.5, 0.6, 0.3, 0.9, -0.8, 0, 0, 0, 0, 0, 0, 0)

# Extract the elevation values
elevation_values <- Oroville2[["Elevation (feet)"]]

# Loop over each row
for (i in 1:nrow(Oroville2)) {
  # Get the elevation value for this row
  elevation <- elevation_values[i]
  
  # Loop over each temperature column and its adjustment
  for (j in seq_along(columns)) {
    col <- as.character(columns[j])
    adjustment <- adjustments[j]
    
    # Check if elevation is not missing
    if (!is.na(elevation)) {
      # Check if elevation is greater than the column name
      if (elevation > as.numeric(col)) {
        # Add adjustment to the current cell value
        Oroville2[i, col] <- Oroville2[i, col] + adjustment
      } else {
        # If elevation is smaller, set the cell value to NA
        Oroville2[i, col] <- NA
      }
    }
  }
}

# Print the updated dataset
print(Oroville2)

```


```{r}
Oroville3 <- Oroville2 %>%
  mutate(Year = lubridate::year(TIME)) %>%
  reshape2::melt(., id = c("Elevation (feet)", "TIME", "Year")) %>%
  rename(`Height of Measurement (feet)` = variable, 
         "Temperature (\u00B0F)" = value) %>%
  mutate(`Height of Measurement (feet)` = as.double(as.character(`Height of Measurement (feet)`)))
Oroville3

# Use approx() to find the closest match between elevation values
matched_storage <- approx(storage_curve$`Elevation (feet)`, storage_curve$`Storage (acre-feet)`, xout = Oroville3$`Height of Measurement (feet)`)

# Create a new column "Storage (acre-feet)" in Oroville3 and assign matched storage values
Oroville3$`Storage (acre-feet)` <- matched_storage$y
# Reorder Oroville3 dataframe based on TIME column
Oroville3 <- Oroville3[order(Oroville3$TIME), ]
Oroville3
# Print the updated Oroville3 dataframe
tail(Oroville3, 100)
#calculate the probability of exceedance with the 51-52 thresholds
```



```{r}
library(hrbrthemes)
library(RColorBrewer)
library(viridis)

#my_palette <- brewer.pal(n = 100, name = "RdYlBu")


  ggplot(Oroville3, aes(x = TIME, y = as.factor(`Height of Measurement (feet)`), fill= `Temperature (°F)`)) + 
    geom_tile() +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y", expand = c(0,0))+
 theme_bw()  +  
  #scale_fill_manual(values = my_palette) +
  scale_fill_viridis(option = "inferno", breaks = seq(32, 100, by = 8)) + 
    labs(x = NULL)+
    png("Heat map of water temperature at Oroville.png", width = 4000, height = 2000, res = 300)


```


# Load necessary libraries
library(dplyr)
library(tidyr)

# Drop rows with missing values in the Temperature column
Oroville3_clean <- Oroville3 %>%
  filter(!is.na(`Temperature (°F)`))

# Define the target temperature
target_temp <- 52

# Function to perform interpolation and handle NA cases
interpolate_temperature <- function(df, temp_target) {
  results <- df %>%
    group_by(TIME) %>%
    summarise(
      Height_Interp = {
        if (n() < 2) {
          NA
        } else {
          temp_vals <- `Temperature (°F)`
          height_vals <- `Height of Measurement (feet)`
          
          if (all(temp_vals < temp_target) || all(temp_vals > temp_target)) {
            NA
          } else {
            approx(x = temp_vals, y = height_vals, xout = temp_target)$y
          }
        }
      },
      Storage_Interp = {
        if (n() < 2) {
          NA
        } else {
          temp_vals <- `Temperature (°F)`
          storage_vals <- `Storage (acre-feet)`
          
          if (all(temp_vals < temp_target) || all(temp_vals > temp_target)) {
            NA
          } else {
            approx(x = temp_vals, y = storage_vals, xout = temp_target)$y
          }
        }
      }
    )
  
  results
}

# Perform interpolation
interpolated_results <- interpolate_temperature(Oroville3_clean, target_temp)

# Print the resulting data frame
print(interpolated_results)

# Optionally, write the results to a CSV file
#write.csv(interpolated_results, "interpolated_temperature_52_values.csv", row.names = FALSE)
Oroville2



library(dplyr)
library(tidyr)

# Create a sequence of dates from 2003-01-01 to 2024-06-05
date_sequence <- seq(as.Date("2003-01-01"), Sys.Date(), by = "day")

# Create a dataframe with the sequence of dates
df <- data.frame(TIME = rep(date_sequence, each = 280))  # 891 - 612 + 1 = 280

# Add a column for Height of Measurement (feet) ranging from 612 to 891 everyday
df <- df %>%
  mutate(`Height of Measurement (feet)` = rep(612:891, length.out = nrow(df)),
         Year = as.double(2003))

# Display the first few rows of the dataframe
df%>% filter(TIME == "2003-01-01")




```{r}
library(dplyr)

library(dplyr)

# Function to perform interpolation for a single day
interpolate_for_day <- function(df_day) {
  # Filter out rows where Elevation is smaller than Height of Measurement
  df_day <- df_day %>%
    filter(`Elevation (feet)` >= `Height of Measurement (feet)`)
  
  # Remove rows with NA in Temperature
  df_day <- df_day %>%
    filter(!is.na(`Temperature (°F)`))
  
  # Ensure there are at least two data points for interpolation
  if (nrow(df_day) < 2) {
    return(data.frame(
      `TIME` = df_day$TIME[1], 
      `Height_48` = NA, 
      `Storage_48` = NA,
      `Height_50` = NA,
      `Storage_50` = NA,
      `Height_52` = NA,
      `Storage_52` = NA
    ))
  }
  
  results <- data.frame(
    `TIME` = df_day$TIME[1], 
    `Height_48` = NA, 
    `Storage_48` = NA,
    `Height_50` = NA,
    `Storage_50` = NA,
    `Height_52` = NA,
    `Storage_52` = NA
  )
  
  target_temps <- c(48, 50, 52)
  
  for (temp in target_temps) {
    if (max(df_day$`Temperature (°F)`) < temp) {
      # Select the highest Height and corresponding Storage
      highest_row <- df_day %>%
        filter(`Height of Measurement (feet)` == max(`Height of Measurement (feet)`, na.rm = TRUE))
      results[paste0("Height_", temp)] <- highest_row$`Height of Measurement (feet)`
      results[paste0("Storage_", temp)] <- highest_row$`Storage (acre-feet)`
    } else {
      # Perform interpolation for Height of Measurement
      height_temp <- approx(x = df_day$`Temperature (°F)`, y = df_day$`Height of Measurement (feet)`, xout = temp)$y
      
      # Perform interpolation for Storage
      storage_temp <- approx(x = df_day$`Temperature (°F)`, y = df_day$`Storage (acre-feet)`, xout = temp)$y
      
      results[paste0("Height_", temp)] <- height_temp
      results[paste0("Storage_", temp)] <- storage_temp
    }
  }
  
  return(results)
}

# Function to apply interpolation to the entire dataset
interpolate_temperatures <- function(df) {
  # Split the dataframe by TIME
  df_split <- split(df, df$TIME)
  
  # Apply interpolation for each day and bind the results
  interpolated_results <- do.call(rbind, lapply(df_split, function(day_df) {
    tryCatch(interpolate_for_day(day_df), error = function(e) {
      data.frame(
        `TIME` = unique(day_df$TIME), 
        `Height_48` = NA, 
        `Storage_48` = NA,
        `Height_50` = NA,
        `Storage_50` = NA,
        `Height_52` = NA,
        `Storage_52` = NA
      )
    })
  }))
  
  return(interpolated_results)
}

# Apply the function to the Oroville3 dataset
interpolated_data <- interpolate_temperatures(Oroville3)

# View the interpolated data
print(interpolated_data)

# Apply the function to the Oroville3 dataset
interpolated_data <- interpolate_temperatures(Oroville3)

interpolated_data
```

```{r}
library(dplyr)
library(lubridate)
library(reshape2)

# Ensure the TIME column is of Date type
interpolated_data$TIME <- as.Date(interpolated_data$TIME)

# Extract the day of the year from the TIME column
interpolated_data <- interpolated_data %>%
  mutate(day_of_year = yday(TIME))

# Function to calculate percentiles
calculate_percentiles <- function(data) {
  percentiles <- c(5, 10, 25, 50, 75, 90, 95)
  sapply(percentiles, function(p) quantile(data, probs = 1 - p / 100, na.rm = TRUE))
}

# Calculate percentiles for each day of the year and each storage column
percentile_data <- interpolated_data %>%
  group_by(day_of_year) %>%
  summarise(
    Storage_48_5 = calculate_percentiles(Storage_48)[1],
    Storage_48_10 = calculate_percentiles(Storage_48)[2],
    Storage_48_25 = calculate_percentiles(Storage_48)[3],
    Storage_48_50 = calculate_percentiles(Storage_48)[4],
    Storage_48_75 = calculate_percentiles(Storage_48)[5],
    Storage_48_90 = calculate_percentiles(Storage_48)[6],
    Storage_48_95 = calculate_percentiles(Storage_48)[7],
    Storage_50_5 = calculate_percentiles(Storage_50)[1],
    Storage_50_10 = calculate_percentiles(Storage_50)[2],
    Storage_50_25 = calculate_percentiles(Storage_50)[3],
    Storage_50_50 = calculate_percentiles(Storage_50)[4],
    Storage_50_75 = calculate_percentiles(Storage_50)[5],
    Storage_50_90 = calculate_percentiles(Storage_50)[6],
    Storage_50_95 = calculate_percentiles(Storage_50)[7],
    Storage_52_5 = calculate_percentiles(Storage_52)[1],
    Storage_52_10 = calculate_percentiles(Storage_52)[2],
    Storage_52_25 = calculate_percentiles(Storage_52)[3],
    Storage_52_50 = calculate_percentiles(Storage_52)[4],
    Storage_52_75 = calculate_percentiles(Storage_52)[5],
    Storage_52_90 = calculate_percentiles(Storage_52)[6],
    Storage_52_95 = calculate_percentiles(Storage_52)[7]
  )



```

```{r}

# Melt the data to rearrange it into the desired format
melted_data <- melt(percentile_data, id.vars = "day_of_year", 
                    variable.name = "percentile", value.name = "storage_value")

# Separate the temperature and percentile into different columns
melted_data <- melted_data %>%
  mutate(temperature = gsub("Storage_(\\d+)_.*", "\\1", percentile),
         percentile = gsub("Storage_\\d+_(.*)", "\\1", percentile))

# View the resulting data
print(melted_data)

melted_data$percentile <- factor(melted_data$percentile, levels = c("5", "10", "25", "50", "75", "90", "95"))

melted_data
```


```{r}
melted_data2 <- melted_data %>% arrange(day_of_year)
# Calculate 7-day moving average and transform storage_value
library(dplyr)

# Find the last 21 rows
last_21_rows <- melted_data2[7624:7686,]

# Find the first 21 rows
first_21_rows <- melted_data2[1:63,]

# Append the last 21 rows at the beginning of the dataset
melted_data3 <- bind_rows(last_21_rows, melted_data2)

# Append the first 21 rows at the end of the dataset
melted_data4 <- bind_rows(melted_data3, first_21_rows)

melted_data4 <- melted_data4[-c(7750:7770),]
melted_data4
```


```{r}

melted_data5 <- melted_data4 %>%
rename(Exceedances = percentile) %>%
group_by(temperature, Exceedances) %>%
#arrange(day_of_year) %>%
mutate(`Volume (TAF)` = zoo::rollmean(storage_value, 7, fill = NA, align = "center") / 1000) %>%
ungroup() %>%
filter(!is.na(`Volume (TAF)`)) %>%
mutate(date = as.Date(day_of_year - 1, origin = "2000-01-01")) # Convert day_of_year to actual dates
melted_data5

```

```{r}
# Convert the date column to Date format
melted_data5$date <- as.Date(melted_data5$date)

melted_data5 <- melted_data5 %>% 
  mutate(temperature = paste0("Oroville Lake Cold Water Pool Volume ≤", temperature, "°F - Percent Exceedances (2003-2023)",sep=""))

melted_data5$temperature <- factor(melted_data5$temperature, levels = c("Oroville Lake Cold Water Pool Volume ≤48°F - Percent Exceedances (2003-2023)", "Oroville Lake Cold Water Pool Volume ≤50°F - Percent Exceedances (2003-2023)", "Oroville Lake Cold Water Pool Volume ≤52°F - Percent Exceedances (2003-2023)"))
melted_data5
```

```{r}
year2024 <- interpolated_data%>%
  mutate(year = year(TIME)) %>%
  filter(year == 2024) %>%
  select(day_of_year, Storage_48, Storage_50, Storage_52) %>%
  rename(`Oroville Lake Cold Water Pool Volume ≤48°F - Percent Exceedances (2003-2023)` = Storage_48,  
         `Oroville Lake Cold Water Pool Volume ≤50°F - Percent Exceedances (2003-2023)` = Storage_50,
         `Oroville Lake Cold Water Pool Volume ≤52°F - Percent Exceedances (2003-2023)` = Storage_52)%>%
  melt(., id = "day_of_year") %>%
  rename(`Volume (TAF)` = value, temperature = variable) %>%
  arrange(day_of_year)
  
year2024


# Find the last 21 rows
last_21_rows <- year2024[1:9,]

# Find the first 21 rows
first_21_rows <- year2024[358:366,]

# Append the last 21 rows at the beginning of the dataset
year2024 <- bind_rows(last_21_rows, year2024)

# Append the first 21 rows at the end of the dataset
year2024 <- bind_rows(year2024, first_21_rows)


# Calculate the moving average
year2024 <- year2024 %>%
  group_by(temperature) %>%
  mutate(`Volume (TAF)` = zoo::rollmean(`Volume (TAF)`, 7, fill = NA, align = "center")/1000) %>%
  drop_na()%>%
  mutate(date = as.Date(day_of_year - 1, origin = "2000-01-01"))
year2024

year2024$temperature <- factor(year2024$temperature, levels = c("Oroville Lake Cold Water Pool Volume ≤48°F - Percent Exceedances (2003-2023)", "Oroville Lake Cold Water Pool Volume ≤50°F - Percent Exceedances (2003-2023)", "Oroville Lake Cold Water Pool Volume ≤52°F - Percent Exceedances (2003-2023)"))
```


```{r}
# Plot the data
hline_data <- data.frame(Exceedances = "Bottom of\nlowest shutter", yintercept = 678.79952)

# Plot the data
ggplot() +
  facet_wrap(~temperature, nrow = 3) +   
  scale_y_continuous(limit = c(0, 3000), expand = c(0,0)) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b %d", expand = c(0, 0)) +
  theme_bw() +
    geom_ribbon(data = year2024, mapping = aes(ymax = `Volume (TAF)`, ymin = 0, x = date, fill = "Year 2024"), alpha = 0.45) +
    guides(fill = guide_legend(title = "")) +
  scale_fill_manual(values = "gray")+
  
  geom_line(data = melted_data5, mapping = aes(x = date, y = `Volume (TAF)`, color = Exceedances)) +
  scale_color_manual(values = c( "dodgerblue4", "deepskyblue2", "chartreuse3", "gold","darkorange","red","darkred", "black")) +
  geom_hline(data = hline_data, aes(yintercept = yintercept, color = Exceedances), linetype = "dashed") +
 # annotate("text", x = max(melted_data5$date), y = 3760, label = "Storage capacity", hjust = 1, vjust = 1, color = "black") +
  theme(
    plot.title = element_text(hjust = 0.5),
    strip.placement = "outside",
    strip.background = element_blank(),
    panel.spacing.y = unit(0, "lines"),
    axis.title.x = element_blank(),  # Remove x-axis label
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)  # Use Arial font
  ) +
 png("Oroville_temperature_excedprob3.png", height = 2500, width = 2000, res = 300)

```