Trend Analysis

Retrieve Discharge data

cat("\014") 

rm(list=ls())
library(dataRetrieval)
library(sf)
library(leaflet)
library(dplyr)
library(tidyr)
library(purrr)
library(tibble)
library(ggplot2)
library(base64enc)
library(lubridate)  
library(patchwork)
library(ggplot2)
library(grid)
library(stringr)
library(wql)
library(ggplot2)
library(grid)
library(lubridate)
library(measurements)
library(dplyr)
library(tidyverse)
library(openair)
library(heatwaveR)
library(data.table)
library(ggplot2)
library(grid)
library(lubridate)
library(measurements)
library(dplyr)
library(tidyverse)
library(openair)
library(heatwaveR)
library(data.table)
require(methods)
library(investr)
require(gridExtra)
library(Metrics)
library(scales)
library(htmlwidgets)
library(Kendall)
library(nhdplusTools)
library(tigris)
# change file_Path_Variable 
file_Path_Variable_IT <- "C:/Users/a905h226/OneDrive - University of Kansas/Desktop/TooBigForGithub_Steps_Workflow_Sept17/InputFiles"

file_Path_Variable_O<- "C:/Users/a905h226/OneDrive - University of Kansas/Desktop/Steps_Workflow_Sept17/Output"
file_Path_Variable_I<- "C:/Users/a905h226/OneDrive - University of Kansas/Desktop/Steps_Workflow_Sept17/InputFiles"


TargetDomain_shapefile_path <- file.path("C:/Users/a905h226/OneDrive - University of Kansas/Desktop/TooBigForGithub_Steps_Workflow_Sept17/InputFiles/EKSRBBoundary_shapefile/watershed_bndry.shp")
TargetDomain_shp <- sf::st_read(TargetDomain_shapefile_path, quiet = TRUE)
TargetDomain_shp <- sf::st_transform(TargetDomain_shp, crs = "EPSG:4326") 
TargetDomain_shp <- sf::st_as_sf(TargetDomain_shp)
combined_shp<-TargetDomain_shp
bounding_box <- st_bbox(combined_shp)

xmin <- sprintf("%.7f", bounding_box["xmin"])
ymin <- sprintf("%.7f", bounding_box["ymin"])
xmax <- sprintf("%.7f", bounding_box["xmax"])
ymax <- sprintf("%.7f", bounding_box["ymax"])

stations <- whatNWISsites(
  bBox = c(xmin, ymin, xmax, ymax)
)

stations_sf <- st_as_sf(stations, coords = c("dec_long_va", "dec_lat_va"), crs = st_crs(combined_shp))

stations_within_combined <- st_intersection(stations_sf, combined_shp)

site_numbers <- stations_within_combined$site_no
data_avail_df <- whatNWISdata(siteNumber = site_numbers)

stations_with_streamflow <- data_avail_df %>%
  filter(parm_cd == "00060") %>%
  distinct(site_no)

stations_with_streamflow_sf <- stations_within_combined %>%
  filter(site_no %in% stations_with_streamflow$site_no)

coords <- st_coordinates(stations_with_streamflow_sf)

streamflow_data_list <- map(stations_with_streamflow_sf$site_no, function(site) {
  data <- readNWISdv(siteNumber = site, parameterCd = "00060")
  if (nrow(data) > 0) {
    data <- data %>%
      select(Date, X_00060_00003) %>%
      rename(mean_streamflow = X_00060_00003)
    data <- tibble::tibble(Date = data$Date, mean_streamflow = data$mean_streamflow)
  } else {
    data <- tibble::tibble(Date = as.Date(character()), mean_streamflow = numeric())
  }
  return(data)
})


metrics_list <- map(streamflow_data_list, function(data) {
  if (nrow(data) > 0) {
    start_year <- min(year(data$Date), na.rm = TRUE)
    end_year <- max(year(data$Date), na.rm = TRUE)
    start_day <- min((data$Date), na.rm = TRUE)
    end_day <- max((data$Date), na.rm = TRUE)
    
    num_years <- end_year - start_year + 1
    total_days <- num_years * 365
    missing_days <- sum(is.na(data$mean_streamflow))
    percent_missing <- (missing_days / total_days) * 100
  } else {
    start_year <- NA
    end_year <- NA
    num_years <- 0
    percent_missing <- 100
    start_day<- NA
    end_day<- NA
  }
  tibble(start_year, end_year,start_day,end_day, num_years, percent_missing)
})
metrics_tibble <- bind_rows(metrics_list)

streamflow_tibble <- tibble(
  station_name = stations_with_streamflow_sf$station_nm,
  site_no = stations_with_streamflow_sf$site_no,
  station_lat = coords[, "Y"],
  station_lon = coords[, "X"],
  streamflow_data = streamflow_data_list,
  start_year = metrics_tibble$start_year,
  end_year = metrics_tibble$end_year,
  start_day = metrics_tibble$start_day,
  end_day = metrics_tibble$end_day,
  num_years = metrics_tibble$num_years,
  #percent_missing = round(metrics_tibble$percent_missing, 2)  # Round to 2 decimal places
)



streamflow_tibbles<- streamflow_tibble %>% 
  filter(map_lgl(streamflow_data, ~nrow(.x)>0))


streamflow_tibbles<- streamflow_tibbles %>% 
  mutate(streamflowUnit='cfs')

streamflow_tibbles<- streamflow_tibbles %>% select(station_name,site_no,station_lat,station_lon,streamflow_data,streamflowUnit)
StreamflowData <- streamflow_tibbles

extract_years<- function(df){
  year(as.Date(df$Date))
}

StreamflowData<- StreamflowData %>%  mutate(SamplingYears=map(streamflow_data,extract_years))

fill_missing_dates <- function(data) {
  if (nrow(data) > 0) {
    full_dates <- tibble(Date = seq(min(data$Date, na.rm = TRUE), max(data$Date, na.rm = TRUE), by = "day"))
    
    full_data <- full_dates %>%
      left_join(data, by = "Date")
    
    return(full_data)
  } else {
    return(tibble(Date = as.Date(character()), mean_streamflow = numeric()))
  }
}

StreamflowData2 <- StreamflowData %>%
  mutate(full_streamflow_data = map(streamflow_data, fill_missing_dates))
StreamflowData<- StreamflowData2


calculate_metrics <- function(data) {
  if (nrow(data) > 0) {
    start_year <- min(year(data$Date), na.rm = TRUE)
    end_year <- max(year(data$Date), na.rm = TRUE)
    start_day <- min(data$Date, na.rm = TRUE)
    end_day <- max(data$Date, na.rm = TRUE)
    
    num_years <- n_distinct(year(data$Date), na.rm = TRUE)
    total_days <- as.numeric(difftime(end_day, start_day, units = "days")) + 1
    missing_days <- sum(is.na(data$mean_streamflow)) # all dates are counted as we are taking the full_streamflow_data 
    percent_missing <- (missing_days / total_days) * 100
  } else {
    start_year <- NA
    end_year <- NA
    start_day <- NA
    end_day <- NA
    num_years <- NA
    total_days <- NA
    missing_days <- NA
    percent_missing <- NA
  }
  
  tibble(start_year, end_year, start_day, end_day, num_years, total_days, missing_days, percent_missing)
}

metrics <- StreamflowData %>%
  mutate(metrics = map(full_streamflow_data, calculate_metrics)) %>%
  unnest(metrics)

StreamflowData<- metrics

filter the data where percent missing is less than 10

Each year 95% data atleast

# filter the data where percent missing is less than 10



filtered_tibble <- StreamflowData %>%
  filter(percent_missing < 10)

filter_years_by_data_availability <- function(streamflow_data) {
  streamflow_data %>%
    mutate(year = year(Date)) %>%
    group_by(year) %>%
    summarise(
      days_in_year = n(),
      total_days_in_year = ifelse(leap_year(year), 366, 365),
      percent_data = (days_in_year / total_days_in_year) * 100
    ) %>%
    filter(percent_data >= 95) %>%
    select(year)
}

filtered_tibbleWP <- filtered_tibble %>%
  mutate(
    streamflow_data = map(streamflow_data, ~ {
      available_years <- filter_years_by_data_availability(.x)
      .x %>%
        mutate(year = year(Date)) %>%
        filter(year %in% available_years$year) %>%
        select(-year)  
    })
  )


filtered_tibbleWP

# A tibble: 59 × 16
   station_name   site_no station_lat station_lon streamflow_data streamflowUnit
   <chr>          <chr>         <dbl>       <dbl> <list>          <chr>         
 1 SMOKY HILL R … 068790…        39.0       -96.9 <tibble>        cfs           
 2 REPUBLICAN R … 068571…        39.0       -96.8 <tibble>        cfs           
 3 KANSAS R AT F… 068791…        39.1       -96.8 <tibble>        cfs           
 4 CLARK C NR JU… 068792…        39.0       -96.7 <tibble>        cfs           
 5 KANSAS R AT O… 068795…        39.1       -96.7 <tibble>        cfs           
 6 WILDCAT C AT … 068798…        39.2       -96.6 <tibble>        cfs           
 7 KINGS C NR MA… 068796…        39.1       -96.6 <tibble>        cfs           
 8 BIG BLUE R NR… 068870…        39.2       -96.6 <tibble>        cfs           
 9 KANSAS R AT W… 068875…        39.2       -96.3 <tibble>        cfs           
10 MILL C NR PAX… 068885…        39.1       -96.2 <tibble>        cfs           
# ℹ 49 more rows
# ℹ 10 more variables: SamplingYears <list>, full_streamflow_data <list>,
#   start_year <int>, end_year <int>, start_day <date>, end_day <date>,
#   num_years <int>, total_days <dbl>, missing_days <int>,
#   percent_missing <dbl>

Cfs to MMD

filtered_tibbleWP_DataOmit<- filtered_tibbleWP %>% filter(map_lgl(streamflow_data,~!is.null(.x)&&nrow(.x)>0))


##which statons are missing #check data to see if it really does not have data
StationNotFound<- anti_join(filtered_tibbleWP,filtered_tibbleWP_DataOmit,by='station_name')
StationNotFound

# A tibble: 3 × 16
  station_name    site_no station_lat station_lon streamflow_data streamflowUnit
  <chr>           <chr>         <dbl>       <dbl> <list>          <chr>         
1 DEER C AT SE 6… 068896…        39.0       -95.6 <tibble>        cfs           
2 ROCK C NORTHEA… 068905…        39.2       -95.6 <tibble>        cfs           
3 L MILL C AT TO… 068925…        39.0       -94.8 <tibble>        cfs           
# ℹ 10 more variables: SamplingYears <list>, full_streamflow_data <list>,
#   start_year <int>, end_year <int>, start_day <date>, end_day <date>,
#   num_years <int>, total_days <dbl>, missing_days <int>,
#   percent_missing <dbl>

filtered_tibble_selected <- filtered_tibble %>%
  filter(station_name %in% StationNotFound$station_name)
filtered_tibble_selected

# A tibble: 3 × 16
  station_name    site_no station_lat station_lon streamflow_data streamflowUnit
  <chr>           <chr>         <dbl>       <dbl> <list>          <chr>         
1 DEER C AT SE 6… 068896…        39.0       -95.6 <tibble>        cfs           
2 ROCK C NORTHEA… 068905…        39.2       -95.6 <tibble>        cfs           
3 L MILL C AT TO… 068925…        39.0       -94.8 <tibble>        cfs           
# ℹ 10 more variables: SamplingYears <list>, full_streamflow_data <list>,
#   start_year <int>, end_year <int>, start_day <date>, end_day <date>,
#   num_years <int>, total_days <dbl>, missing_days <int>,
#   percent_missing <dbl>

filtered_tibble<-filtered_tibbleWP_DataOmit

hist(filtered_tibble$num_years)

AllYear_StreamflowData <- filtered_tibble


######Calculate mm per day
cfs_to_m3s <- 0.0283168  
seconds_per_day <- 86400 

retrieve_watershed_area <- function(site_no) {
  site_info <- readNWISsite(site_no)
  if (!is.null(site_info$drain_area_va)) {
    return(site_info$drain_area_va)  # Drainage area in square miles
  } else {
    return(NA)  
  }
}

AllYear_StreamflowData <- AllYear_StreamflowData %>%
  mutate(watershed_area_sqmi = map_dbl(site_no, retrieve_watershed_area)) %>%
  mutate(watershed_area_m2 = measurements::conv_unit(watershed_area_sqmi, from = "mi2", to = "m2"))
# 549 square miles 06863500
NoData<-AllYear_StreamflowData %>% filter( is.na(AllYear_StreamflowData$watershed_area_m2))
NoData$site_no

character(0)

convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * cfs_to_m3s * seconds_per_day
  streamflow_mm_per_day <- (streamflow_m3_per_day / watershed_area_m2) * 1000
  return(streamflow_mm_per_day)
}

AllYear_StreamflowData$streamflow_data <- AllYear_StreamflowData$streamflow_data %>%
  map2(AllYear_StreamflowData$watershed_area_m2, ~ .x %>%
         mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, .y)))
print(AllYear_StreamflowData)

# A tibble: 56 × 18
   station_name   site_no station_lat station_lon streamflow_data streamflowUnit
   <chr>          <chr>         <dbl>       <dbl> <list>          <chr>         
 1 SMOKY HILL R … 068790…        39.0       -96.9 <tibble>        cfs           
 2 REPUBLICAN R … 068571…        39.0       -96.8 <tibble>        cfs           
 3 KANSAS R AT F… 068791…        39.1       -96.8 <tibble>        cfs           
 4 CLARK C NR JU… 068792…        39.0       -96.7 <tibble>        cfs           
 5 KANSAS R AT O… 068795…        39.1       -96.7 <tibble>        cfs           
 6 WILDCAT C AT … 068798…        39.2       -96.6 <tibble>        cfs           
 7 KINGS C NR MA… 068796…        39.1       -96.6 <tibble>        cfs           
 8 BIG BLUE R NR… 068870…        39.2       -96.6 <tibble>        cfs           
 9 KANSAS R AT W… 068875…        39.2       -96.3 <tibble>        cfs           
10 MILL C NR PAX… 068885…        39.1       -96.2 <tibble>        cfs           
# ℹ 46 more rows
# ℹ 12 more variables: SamplingYears <list>, full_streamflow_data <list>,
#   start_year <int>, end_year <int>, start_day <date>, end_day <date>,
#   num_years <int>, total_days <dbl>, missing_days <int>,
#   percent_missing <dbl>, watershed_area_sqmi <dbl>, watershed_area_m2 <dbl>

filter the canals

canal_stations<- AllYear_StreamflowData[grepl('canal',AllYear_StreamflowData$station_name,ignore.case = TRUE),]
AllYear_StreamflowData_NoCanal<- AllYear_StreamflowData[!grepl('canal',AllYear_StreamflowData$station_name,ignore.case = TRUE),]
AllYear_StreamflowData_NoCanalF<-na.omit(AllYear_StreamflowData_NoCanal)
AllYear_StreamflowData_NoCanalF<- AllYear_StreamflowData_NoCanalF %>% rename(streamflow_data_CFS_MMD=streamflow_data)
AllYear_StreamflowData_NoCanalFk<- AllYear_StreamflowData_NoCanalF %>% 
  mutate(StreamflowMMD=map(streamflow_data_CFS_MMD,~select(.x,Date,mean_streamflow_mm_per_day)))
AllYear_StreamflowData_NoCanalFk<- AllYear_StreamflowData_NoCanalFk %>% mutate(streamflow_data = StreamflowMMD)
AllYear_StreamflowData_NoCanalFk

# A tibble: 56 × 20
   station_name           site_no station_lat station_lon streamflow_data_CFS_…¹
   <chr>                  <chr>         <dbl>       <dbl> <list>                
 1 SMOKY HILL R AT JUNCT… 068790…        39.0       -96.9 <tibble [1,827 × 3]>  
 2 REPUBLICAN R AT JUNCT… 068571…        39.0       -96.8 <tibble [22,281 × 3]> 
 3 KANSAS R AT FORT RILE… 068791…        39.1       -96.8 <tibble [22,280 × 3]> 
 4 CLARK C NR JUNCTION C… 068792…        39.0       -96.7 <tibble [2,557 × 3]>  
 5 KANSAS R AT OGDEN, KS  068795…        39.1       -96.7 <tibble [11,323 × 3]> 
 6 WILDCAT C AT SCENIC D… 068798…        39.2       -96.6 <tibble [4,738 × 3]>  
 7 KINGS C NR MANHATTAN,… 068796…        39.1       -96.6 <tibble [16,437 × 3]> 
 8 BIG BLUE R NR MANHATT… 068870…        39.2       -96.6 <tibble [27,029 × 3]> 
 9 KANSAS R AT WAMEGO, KS 068875…        39.2       -96.3 <tibble [38,351 × 3]> 
10 MILL C NR PAXICO, KS   068885…        39.1       -96.2 <tibble [25,933 × 3]> 
# ℹ 46 more rows
# ℹ abbreviated name: ¹​streamflow_data_CFS_MMD
# ℹ 15 more variables: streamflowUnit <chr>, SamplingYears <list>,
#   full_streamflow_data <list>, start_year <int>, end_year <int>,
#   start_day <date>, end_day <date>, num_years <int>, total_days <dbl>,
#   missing_days <int>, percent_missing <dbl>, watershed_area_sqmi <dbl>,
#   watershed_area_m2 <dbl>, StreamflowMMD <list>, streamflow_data <list>

MMD<- 'MMD'
AllYear_StreamflowData_NoCanalFk$streamflowUnit<- MMD                                                                         

Annual Runoff

USGS daily streamflow data in cubic feet per second (cfs).

Convert it to cubic meters per second (m³/s) → then to m³/day.

Divide by watershed area in m² and multiply by 1000 → this gives runoff depth in mm/day.

Aggregate daily values over the year (or season) to compute annual/seasonal runoff in mm

#Add run off to the column 
AllYear_StreamflowData_NoCanalFk <- AllYear_StreamflowData_NoCanalFk %>%
  mutate(
    annual_runoff = map(streamflow_data, ~ .x %>%
                          mutate(Year = year(Date)) %>%
                          group_by(Year) %>%
                          summarise(
                            runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
                            .groups = "drop"
                          ))
  )






annual_runoff_list <- AllYear_StreamflowData_NoCanalFk %>%
  mutate(
    annual_runoff = map(streamflow_data, ~ .x %>% ##it was in mmd 
                          mutate(Year = year(Date)) %>%
                          group_by(Year) %>%
                          summarise(runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
                                    .groups = "drop"))
  )


annual_runoff_all <- annual_runoff_list %>%
  select(site_no, watershed_area_m2, annual_runoff) %>%
  unnest(annual_runoff)

Should i use area average or desoto?

# Area-weighted runoff by year across EKSrb
annual_runoff_eksrb <- annual_runoff_all %>%
  mutate(weight = watershed_area_m2 / sum(watershed_area_m2, na.rm = TRUE),
         weighted_runoff = runoff_mm * weight) %>%
  group_by(Year) %>%
  summarise(mean_runoff_mm = sum(weighted_runoff, na.rm = TRUE), .groups = "drop")



annual_runoff_eksrb_plot<- ggplot(data = annual_runoff_eksrb, aes(x = Year, y = mean_runoff_mm)) +
  geom_line(color = "blue", size = 1) +
  geom_point(color = "darkblue", size = 2) +
  labs(
    title = "Annual Area-Weighted Runoff for EKSrb",
    x = "Year",
    y = "Runoff (mm)"
  ) +
  theme_minimal(base_size = 14)

annual_runoff_eksrb_plot

##instead averaging, use desoto as De Soto sits near the very bottom of the EKSrb, so its record captures virtually all the flow generated throughout the entire basin upstream of that point.

de_soto_runoff <- AllYear_StreamflowData_NoCanalFk %>%
  filter(station_name == "KANSAS R AT DESOTO, KS")

annual_runoff_desoto <- de_soto_runoff %>%
  mutate(
    annual_runoff = map(streamflow_data, ~ .x %>%
                          mutate(Year = year(Date)) %>%
                          group_by(Year) %>%
                          summarise(runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
                                    .groups = "drop"))
  ) %>%
  select(site_no, watershed_area_m2, annual_runoff) %>%
  unnest(annual_runoff)

annual_runoff_eksrb_desoto <- annual_runoff_desoto %>%
  rename(mean_runoff_mm = runoff_mm) %>%
  select(Year, mean_runoff_mm)

annual_runoff_eksrb_desoto_plot<- ggplot(data = annual_runoff_eksrb_desoto, aes(x = Year, y = mean_runoff_mm)) +
  geom_line(color = "blue", size = 1) +
  geom_point(color = "darkblue", size = 2) +
  labs(
    title = "Annual Runoff for De Soto Station (EKSrb Outlet)",
    x = "Year",
    y = "Runoff (mm)"
  ) +
  theme_minimal(base_size = 14)


annual_runoff_eksrb_desoto_plot

annual_runoff_eksrb_desoto_plot/annual_runoff_eksrb_plot

Add precip data

TargetDomain_shp <- sf::st_read(TargetDomain_shapefile_path, quiet = TRUE)
TargetDomain_shp <- sf::st_transform(TargetDomain_shp, crs = "EPSG:32614") 
TargetDomain_shp <- sf::st_as_sf(TargetDomain_shp)

Precipitation_Data_Gridded_Filtered <- readRDS(file.path(file_Path_Variable_O, "step14_Precipitation_Data_Gridded_Filtered_extremes_Upuntil2023_feb28.rds"))
Precipitation_Data_Gridded_Filtered_TargetDomain <- Precipitation_Data_Gridded_Filtered
Precipitation_Data_Gridded_Filtered_TargetDomain <- sf::st_as_sf(Precipitation_Data_Gridded_Filtered_TargetDomain, coords = c("x", "y"), crs = 32614)
Precipitation_Data_Gridded_Filtered_TargetDomain <- sf::st_intersection(Precipitation_Data_Gridded_Filtered_TargetDomain, TargetDomain_shp)




annual_precip_trend <- Precipitation_Data_Gridded_Filtered_TargetDomain %>%
  unnest(annual_precip_data) %>%
  group_by(year) %>%
  summarise(avg_precip_mm = mean(precip_mm, na.rm = TRUE), .groups = "drop")

min_year <- min(annual_precip_trend$year, na.rm = TRUE)
max_year <- max(annual_precip_trend$year, na.rm = TRUE)

#Reference ET = How much water could be lost if the crop has unlimited water. #P − Runoff ET = “How much water was actually lost based on what came in (P) and what ran off (R). This method do not #include GW storage #Open ET considers gw

Add two type of ET, eq derived and open et

# ET values represent the total sum of evapotranspiration in millimeters for each month per pixel.
eksrb_waterbalance <- annual_runoff_eksrb %>%
  left_join(annual_precip_trend, by = c("Year" = "year")) %>%
  mutate(ET_mm_PminusR = avg_precip_mm - mean_runoff_mm)

eksrb_waterbalance <- eksrb_waterbalance %>%
                      select(Year,avg_precip_mm,mean_runoff_mm,ET_mm_PminusR)


#Get open ET data
ET_Precip_Monthly <- readRDS(file.path(file_Path_Variable_O, "Step19_precip_et_monthly_nearest_filtered_combined.rds"))
ET_Precip_Monthly_TargetDomain <- ET_Precip_Monthly
ET_Precip_Monthly_TargetDomain<- ET_Precip_Monthly_TargetDomain %>% rename(x=et_x,y=et_y)
ET_Precip_Monthly_TargetDomain <- sf::st_as_sf(ET_Precip_Monthly_TargetDomain, coords = c("x", "y"), crs = 32614)
ET_Precip_Monthly_TargetDomain <- sf::st_intersection(ET_Precip_Monthly_TargetDomain, TargetDomain_shp)

ET_Precip_Monthly_TargetDomain <- ET_Precip_Monthly_TargetDomain %>%
  mutate(grid_id = row_number())

openet_annual_et <- ET_Precip_Monthly_TargetDomain %>%
  select(grid_id, ET_ET_Data) %>%
  unnest(ET_ET_Data) %>%
  group_by(grid_id, Year) %>%
  summarise(annual_et = sum(et_ensemble_sam, na.rm = TRUE), .groups = "drop") %>%
  group_by(Year) %>%
  summarise(ET_OpenET_mm = mean(annual_et, na.rm = TRUE), .groups = "drop")

eksrb_waterbalance <- eksrb_waterbalance %>%
  left_join(openet_annual_et, by = "Year")


eksrb_waterbalance_Et<- eksrb_waterbalance %>% select(Year,ET_OpenET_mm,ET_mm_PminusR)
eksrb_waterbalance_Et<- na.omit(eksrb_waterbalance_Et)

Diff between two type of et

eksrb_waterbalance_Et_long <- eksrb_waterbalance_Et %>%
  pivot_longer(cols = c(ET_OpenET_mm, ET_mm_PminusR),
               names_to = "ET_source",
               values_to = "ET_mm")

ggplot(eksrb_waterbalance_Et_long, aes(x = Year, y = ET_mm, color = ET_source)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  scale_color_manual(values = c("ET_OpenET_mm" = "blue", "ET_mm_PminusR" = "darkgreen"),
                     labels = c("OpenET", "P - Runoff")) +
  labs(
    title = "Annual Evapotranspiration (ET) in EKSrb",
    subtitle = "Comparison of OpenET and Water Balance (P - Runoff) Estimates",
    x = "Year",
    y = "ET (mm)",
    color = "ET Source"
  ) +
  theme_minimal(base_size = 14)

ggplot(eksrb_waterbalance_Et, aes(x = ET_OpenET_mm, y = ET_mm_PminusR)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(
    title = "Comparison of ET Estimates: OpenET vs P - Runoff",
    x = "ET from OpenET (mm/year)",
    y = "ET from P - Runoff (mm/year)"
  ) +
  theme_minimal(base_size = 14)

1. Irrigation Is Captured by OpenET but Not by P − Runoff

#In dry years, vegetation may use stored groundwater or soil moisture to sustain evapotranspiration: # # OpenET “sees” the ET from vegetation via satellite. # # P − R underestimates ET because the water source isn’t from that year’s precipitation.

From grebecht paper

Average Annual Precipitation, Streamflow, and Evapotranspiration for Dry ~1961–1980! and Wet Period ~1981–2001! as a Result of

Decade-Long Climate Variations and Difference between Dry and Wet Period in Percent

eksrb_waterbalance_Et_WithOpenEt<- eksrb_waterbalance %>% select(Year,avg_precip_mm ,mean_runoff_mm,ET_OpenET_mm,ET_mm_PminusR)
eksrb_waterbalance_Et_WithOutOpenEt<- eksrb_waterbalance %>% select(Year,avg_precip_mm ,mean_runoff_mm,ET_mm_PminusR)
eksrb_waterbalance_Et_WithOpenEt<- na.omit(eksrb_waterbalance_Et_WithOpenEt)
eksrb_waterbalance_Et_WithOutOpenEt<- na.omit(eksrb_waterbalance_Et_WithOutOpenEt)




eksrb_waterbalance_desoto <- annual_runoff_eksrb_desoto %>%
  left_join(annual_precip_trend, by = c("Year" = "year")) %>%
  mutate(ET_mm_PminusR = avg_precip_mm - mean_runoff_mm)

summary_table_desoto <- eksrb_waterbalance_desoto %>%
  filter(Year >= 1961, Year <= 2001) %>%
  mutate(period = case_when(
    Year >= 1961 & Year <= 1980 ~ "Dry",
    Year >= 1981 & Year <= 2001 ~ "Wet"
  )) %>%
  group_by(period) %>%
  summarise(
    Precip_mm = mean(avg_precip_mm, na.rm = TRUE),
    Runoff_mm = mean(mean_runoff_mm, na.rm = TRUE),
    ET_mm = mean(ET_mm_PminusR, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  pivot_wider(names_from = period, values_from = c(Precip_mm, Runoff_mm, ET_mm)) %>%
  mutate(
    ΔPrecip = Precip_mm_Wet - Precip_mm_Dry,
    ΔRunoff = Runoff_mm_Wet - Runoff_mm_Dry,
    ΔET = ET_mm_Wet - ET_mm_Dry,
    `%ΔPrecip` = 100 * ΔPrecip / Precip_mm_Dry,
    `%ΔRunoff` = 100 * ΔRunoff / Runoff_mm_Dry,
    `%ΔET` = 100 * ΔET / ET_mm_Dry,
    Runoff_to_P = (Runoff_mm_Dry + Runoff_mm_Wet) / (Precip_mm_Dry + Precip_mm_Wet),
Streamflow_sensitivity = `%ΔRunoff` / `%ΔPrecip`,
    ET_sensitivity = `%ΔET` / `%ΔPrecip`
  ) %>%
  select(
    Precip_mm_Dry, Precip_mm_Wet, ΔPrecip, `%ΔPrecip`,
    Runoff_mm_Dry, Runoff_mm_Wet, ΔRunoff, `%ΔRunoff`,
    ET_mm_Dry, ET_mm_Wet, ΔET, `%ΔET`,
    Runoff_to_P, Streamflow_sensitivity, ET_sensitivity
  )

print(summary_table_desoto)

# A tibble: 1 × 15
  Precip_mm_Dry Precip_mm_Wet ΔPrecip `%ΔPrecip` Runoff_mm_Dry Runoff_mm_Wet
          <dbl>         <dbl>   <dbl>      <dbl>         <dbl>         <dbl>
1          897.          928.    31.4       3.50          44.3          54.1
# ℹ 9 more variables: ΔRunoff <dbl>, `%ΔRunoff` <dbl>, ET_mm_Dry <dbl>,
#   ET_mm_Wet <dbl>, ΔET <dbl>, `%ΔET` <dbl>, Runoff_to_P <dbl>,
#   Streamflow_sensitivity <dbl>, ET_sensitivity <dbl>

summary_table_desoto <- summary_table_desoto %>%
  pivot_longer(cols = everything(),
               names_to = "Metric",
               values_to = "Value")
summary_table_desoto #looks much more reasonable to me 

# A tibble: 15 × 2
   Metric                    Value
   <chr>                     <dbl>
 1 Precip_mm_Dry          897.    
 2 Precip_mm_Wet          928.    
 3 ΔPrecip                 31.4   
 4 %ΔPrecip                 3.50  
 5 Runoff_mm_Dry           44.3   
 6 Runoff_mm_Wet           54.1   
 7 ΔRunoff                  9.85  
 8 %ΔRunoff                22.2   
 9 ET_mm_Dry              852.    
10 ET_mm_Wet              874.    
11 ΔET                     21.5   
12 %ΔET                     2.52  
13 Runoff_to_P              0.0540
14 Streamflow_sensitivity   6.35  
15 ET_sensitivity           0.722 

Add two type of ET, eq derived and open et

library(gt)

table_desoto_side_by_side <- summary_table_desoto %>%
  pivot_wider(names_from = Metric, values_from = Value)

gt_table <- gt(table_desoto_side_by_side) %>%
  tab_header(
    title = "Table 2. Average Annual Precipitation, Streamflow, and Evapotranspiration for Dry (1961–1980) and Wet (1981–2001) Periods",
    subtitle = "As a result of decade-long climate variations in EKSrb"
  ) %>%
  cols_label(
    Precip_mm_Dry = "Precipitation (Dry)",
    Precip_mm_Wet = "Precipitation (Wet)",
    Runoff_mm_Dry = "Streamflow (Dry)",
    Runoff_mm_Wet = "Streamflow (Wet)",
    ET_mm_Dry = "Evapotranspiration (Dry)",
    ET_mm_Wet = "Evapotranspiration (Wet)",
    Runoff_to_P = "Runoff/Precip Ratio",
    Streamflow_sensitivity = "Streamflow Sensitivity",
    ET_sensitivity = "ET Sensitivity"
  ) %>%
  fmt_number(
    columns = everything(),
    decimals = 4
  ) %>%
  cols_align(align = "center", columns = everything())

gt_table

Table 2. Average Annual Precipitation, Streamflow, and Evapotranspiration for Dry (1961–1980) and Wet (1981–2001) Periods
As a result of decade-long climate variations in EKSrb
Precipitation (Dry)	Precipitation (Wet)	ΔPrecip	%ΔPrecip	Streamflow (Dry)	Streamflow (Wet)	ΔRunoff	%ΔRunoff	Evapotranspiration (Dry)	Evapotranspiration (Wet)	ΔET	%ΔET	Runoff/Precip Ratio	Streamflow Sensitivity	ET Sensitivity
896.6167	927.9732	31.3565	3.4972	44.3023	54.1479	9.8456	22.2236	852.3144	873.8253	21.5110	2.5238	0.0540	6.3547	0.7217

format a func

Hydrologic Summary Metrics (Table 1)

• Precip(Dry): Average annual precipitation (mm) during the selected dry period.

• Precip(Wet): Average annual precipitation (mm) during the selected wet period.

• ΔPrecip: Absolute difference in precipitation between wet and dry periods:

  ΔPrecip=Precip(Wet)−Precip(Dry) = - ΔPrecip=Precip(Wet)−Precip(Dry)

• %ΔPrecip: Percent change in precipitation from dry to wet period:

  %ΔPrecip=ΔPrecipPrecip(Dry)×100% = %ΔPrecip=Precip(Dry)ΔPrecip​×100

• Q(Dry): Average annual runoff or streamflow (mm) during the dry period.

• Q(Wet): Average annual runoff during the wet period.

• ΔRunoff: Difference in runoff between wet and dry periods:

  ΔRunoff=Q(Wet)−Q(Dry) = - ΔRunoff=Q(Wet)−Q(Dry)

• %ΔRunoff: Percent change in runoff from dry to wet period.

• ET(Dry): Evapotranspiration (mm) during the dry period, derived from the water balance:

  ET=Precipitation−Runoff = - ET=Precipitation−Runoff

• ET(Wet): Evapotranspiration during the wet period.

• ΔET: Change in ET between wet and dry periods.

• %ΔET: Percent change in ET from dry to wet period.

Derived Metrics (Table 2)

• Runoff/Precip: The ratio of runoff to precipitation during the dry and wet period (averaged). A measure of runoff efficiency.

• QSensitivity: Streamflow sensitivity — change in runoff per unit change in precipitation:

  QSensitivity=ΔRunoffΔPrecip = QSensitivity=ΔPrecipΔRunoff​

  Values close to 1 suggest a strong streamflow response to precipitation changes; values near 0 indicate minimal hydrologic response.

• ETSensitivity: ET sensitivity — change in evapotranspiration per unit change in precipitation:

  ETSensitivity=ΔETΔPrecip = ETSensitivity=ΔPrecipΔET​

  Higher values suggest that increases in precipitation are largely consumed by evapotranspiration rather than contributing to runoff.

make_decadal_summary_table <- function(data, dry_years, wet_years, title_text) {
  
  summary_table <- data %>%
    filter(Year %in% c(dry_years, wet_years)) %>%
    mutate(period = case_when(
      Year %in% dry_years ~ "Dry",
      Year %in% wet_years ~ "Wet",
      TRUE ~ NA_character_
    )) %>%
    group_by(period) %>%
    summarise(
      Precip_mm = mean(avg_precip_mm, na.rm = TRUE),
      Runoff_mm = mean(mean_runoff_mm, na.rm = TRUE),
      ET_mm = mean(ET_mm_PminusR, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_wider(names_from = period, values_from = c(Precip_mm, Runoff_mm, ET_mm)) %>%
    mutate(
      ΔPrecip = Precip_mm_Wet - Precip_mm_Dry,
      ΔRunoff = Runoff_mm_Wet - Runoff_mm_Dry,
      ΔET = ET_mm_Wet - ET_mm_Dry,
      `%ΔPrecip` = 100 * ΔPrecip / Precip_mm_Dry,
      `%ΔRunoff` = 100 * ΔRunoff / Runoff_mm_Dry,
      `%ΔET` = 100 * ΔET / ET_mm_Dry,
      #Runoff_to_P = Runoff_mm_Dry / Precip_mm_Dry,
      Runoff_to_P = (Runoff_mm_Dry + Runoff_mm_Wet) / (Precip_mm_Dry + Precip_mm_Wet), #average 2 is cancelled 
      Retention_ratio = 1 - Runoff_to_P ,
   Streamflow_sensitivity = `%ΔRunoff` / `%ΔPrecip`,
    ET_sensitivity = `%ΔET` / `%ΔPrecip`     
    )
  
  table1 <- summary_table %>%
    select(
      Precip_mm_Dry, Precip_mm_Wet, ΔPrecip, `%ΔPrecip`,
      Runoff_mm_Dry, Runoff_mm_Wet, ΔRunoff, `%ΔRunoff`,
      ET_mm_Dry, ET_mm_Wet, ΔET, `%ΔET`
    ) %>%
    gt() %>%
    tab_header(
      title = title_text,
      subtitle = "Hydrologic metrics based on dry and wet years"
    ) %>%
    cols_label(
      Precip_mm_Dry = "Precip(Dry)",
      Precip_mm_Wet = "Precip(Wet)",
      ΔPrecip = "ΔPrecip", `%ΔPrecip` = "%ΔPrecip",
      Runoff_mm_Dry = "Q(Dry)",
      Runoff_mm_Wet = "Q(Wet)",
      ΔRunoff = "ΔRunoff", `%ΔRunoff` = "%ΔRunoff",
      ET_mm_Dry = "ET(Dry)",
      ET_mm_Wet = "ET(Wet)",
      ΔET = "ΔET", `%ΔET` = "%ΔET"
    ) %>%
    fmt_number(columns = everything(), decimals = 2) %>%
    cols_align("center", columns = everything())
  
  table2 <- summary_table %>%
    select(Runoff_to_P,Retention_ratio, Streamflow_sensitivity, ET_sensitivity) %>%
    gt() %>%
    tab_header(
      title = title_text,
      subtitle = "Derived Metrics"
    ) %>%
    cols_label(
      Runoff_to_P = "Runoff/Precip",
      Retention_ratio = "Retention Ratio" ,
      Streamflow_sensitivity = "QSensitivity",
      ET_sensitivity = "ETSensitivity"
    ) %>%
    fmt_number(columns = everything(), decimals = 2) %>%
    cols_align("center", columns = everything())
  
  return(list(Hydrologic_Table = table1, Derived_Metrics_Table = table2))
}

1. Dust Bowl vs Post-Dust Bowl

Periods 1930–1950 (Dust Bowl) vs 1951–1970 (Recovery/Post-Dust Bowl)

Why Compare extreme drought era to post-recovery; high relevance for Great Plains.

What to Compare Precipitation, runoff, ET, runoff/precip, sensitivities.

Later—- Match this with persistent dry events—should show strong contrast.

2. Pre-1970s vs Modern Climate Shift

Why Mid-century vs recent hydroclimatic regime shift (many changes in land use, irrigation, ET tech).

Main idea is to See if streamflow sensitivity changed with increasing ET demand and altered land cover.

3. Post-1980 vs Post-2000

Periods 1981–2000 vs 2001–2020

increased climate variability

Dust Bowl (1930–1950) vs Post-Dust Bowl (1951–1970)
Hydrologic metrics based on dry and wet years
Precip(Dry)	Precip(Wet)	ΔPrecip	%ΔPrecip	Q(Dry)	Q(Wet)	ΔRunoff	%ΔRunoff	ET(Dry)	ET(Wet)	ΔET	%ΔET
843.72	867.14	23.41	2.78	38.82	42.14	3.32	8.54	804.90	825.00	20.10	2.50

Mid-Century (1950–1969) vs Modern (1990–2009)
Hydrologic metrics based on dry and wet years
Precip(Dry)	Precip(Wet)	ΔPrecip	%ΔPrecip	Q(Dry)	Q(Wet)	ΔRunoff	%ΔRunoff	ET(Dry)	ET(Wet)	ΔET	%ΔET
865.12	936.01	70.89	8.19	43.80	45.80	2.00	4.57	821.32	890.21	68.89	8.39

Post-1980 (1981–2000) vs Post-2000 (2001–2020)
Hydrologic metrics based on dry and wet years
Precip(Dry)	Precip(Wet)	ΔPrecip	%ΔPrecip	Q(Dry)	Q(Wet)	ΔRunoff	%ΔRunoff	ET(Dry)	ET(Wet)	ΔET	%ΔET
922.47	950.92	28.45	3.08	54.46	39.58	−14.88	−27.33	868.01	911.34	43.33	4.99

Dust Bowl (1930–1950) vs Post-Dust Bowl (1951–1970)
Derived Metrics
Runoff/Precip	Retention Ratio	QSensitivity	ETSensitivity
0.05	0.95	3.08	0.90

Mid-Century (1950–1969) vs Modern (1990–2009)
Derived Metrics
Runoff/Precip	Retention Ratio	QSensitivity	ETSensitivity
0.05	0.95	0.56	1.02

Post-1980 (1981–2000) vs Post-2000 (2001–2020)
Derived Metrics
Runoff/Precip	Retention Ratio	QSensitivity	ETSensitivity
0.05	0.95	−8.86	1.62

Interpretation of Decadal Hydroclimatic Shifts in EKSrb The hydrologic summaries comparing dry and wet periods across three distinct eras reveal significant changes in water balance behavior in the Eastern Kansas River Basin (EKSrb), shaped by long-term climate variability and changing land–water interactions.

1. Dust Bowl (1930–1950) vs Post-Dust Bowl (1951–1970) This period marks the transition from one of the most extreme drought events in U.S. history to a phase of moderate climatic recovery. Precipitation increased modestly (~2.8%), with a proportionally higher increase in streamflow (~8.5%). Evapotranspiration (ET) also increased (~2.5%), reflecting improved vegetation water use. Streamflow sensitivity to precipitation was low but positive (0.14), indicating slight recovery in runoff generation. ET sensitivity remained moderate (0.86), suggesting that a significant portion of the added water supported atmospheric moisture loss.

2. Mid-Century (1950–1969) vs Modern (1990–2009) Despite a substantial rise in precipitation (~8.2%), streamflow only increased marginally (~2%), resulting in very low streamflow sensitivity (0.03). In contrast, ET increased markedly (~8.4%), implying that much of the additional water was consumed through evapotranspiration rather than contributing to runoff. These trends are consistent with intensified land use, expanded irrigation, or rising atmospheric demand. The basin appears to have become more ET-dominated during this period.

3. Post-1980 (1981–2000) vs Post-2000 (2001–2020) The most striking result emerges in the recent decades. Although precipitation rose (~3.1%), streamflow declined substantially (−27%), resulting in a negative streamflow sensitivity (−0.52). Meanwhile, ET increased by ~5%, surpassing the precipitation gain. This suggests that vegetation and land surfaces increasingly rely on stored moisture or groundwater to sustain ET during variable climate conditions. The elevated ET sensitivity (1.52) reinforces this pattern. These results reflect increasing decoupling between rainfall and runoff, possibly driven by enhanced evaporative demand, changes in land cover, or climate warming.

Overall Implications Across the 20th and early 21st centuries, the EKSrb has transitioned toward a less runoff-efficient, more evapotranspiration-driven hydrologic regime. These shifts raise important questions about future water availability and watershed resilience under climate change, especially in regions where even modest precipitation increases no longer translate to greater streamflow.

Seasonal component

assign_season_and_year <- function(date) {
  month <- as.numeric(format(date, "%m"))
  year_val <- year(date)
  
  if (month == 12) {
    return(list(season = "Winter", adjusted_year = year_val + 1))
  } else if (month %in% c(1, 2)) {
    return(list(season = "Winter", adjusted_year = year_val))
  } else if (month %in% c(3, 4, 5)) {
    return(list(season = "Spring", adjusted_year = year_val))
  } else if (month %in% c(6, 7, 8)) {
    return(list(season = "Summer", adjusted_year = year_val))
  } else {
    return(list(season = "Fall", adjusted_year = year_val))
  }
}




de_soto_runoff <- AllYear_StreamflowData_NoCanalFk %>%
  filter(station_name == "KANSAS R AT DESOTO, KS")

annual_runoff_desoto <- de_soto_runoff %>%
  mutate(
    annual_runoff = map(streamflow_data, ~ .x %>%
                          mutate(Year = year(Date)) %>%
                          group_by(Year) %>%
                          summarise(runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
                                    .groups = "drop"))
  ) %>%
  select(site_no, watershed_area_m2, annual_runoff) %>%
  unnest(annual_runoff)

annual_runoff_eksrb_desoto <- annual_runoff_desoto %>%
  rename(mean_runoff_mm = runoff_mm) %>%
  select(Year, mean_runoff_mm)











de_soto_runoff <- AllYear_StreamflowData_NoCanalFk %>%
  filter(station_name == "KANSAS R AT DESOTO, KS")

de_soto_seasonal <- de_soto_runoff %>%
  mutate(
    seasonal_streamflow = map(streamflow_data, ~ .x %>%
                                mutate(
                                  season_data = lapply(Date, assign_season_and_year),
                                  Season = sapply(season_data, `[[`, "season"),
                                  adjusted_year = sapply(season_data, `[[`, "adjusted_year")
                                ) %>%
                                select(-season_data)
    )
  )

seasonal_runoff_desoto <- de_soto_seasonal %>%
  select(site_no, watershed_area_m2, seasonal_streamflow) %>%
  unnest(seasonal_streamflow) %>%
  group_by(site_no, adjusted_year, Season) %>%
  summarise(
    seasonal_runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 91.3125,
    .groups = "drop"
  )
seasonal_runoff_desoto

# A tibble: 429 × 4
   site_no  adjusted_year Season seasonal_runoff_mm
   <chr>            <dbl> <chr>               <dbl>
 1 06892350          1918 Fall                 3.61
 2 06892350          1918 Spring               5.00
 3 06892350          1918 Summer               4.65
 4 06892350          1918 Winter               2.48
 5 06892350          1919 Fall                 6.02
 6 06892350          1919 Spring              19.0 
 7 06892350          1919 Summer              13.2 
 8 06892350          1919 Winter               3.91
 9 06892350          1920 Fall                 4.51
10 06892350          1920 Spring               7.34
# ℹ 419 more rows

seasonal_precip_area_avg <- Precipitation_Data_Gridded_Filtered_TargetDomain %>%
  unnest(data) %>%
  mutate(
    season_info = map(YearMonth, assign_season_and_year),
    Season = map_chr(season_info, "season"),
    adjusted_year = map_dbl(season_info, "adjusted_year")
  ) %>%
  group_by(adjusted_year, Season, geometry) %>%
  summarise(seasonal_precip_mm = sum(total_rainfall_monthly, na.rm = TRUE), .groups = "drop") %>%
  group_by(adjusted_year, Season) %>%
  summarise(avg_precip_mm = mean(seasonal_precip_mm, na.rm = TRUE), .groups = "drop")


ET_Precip_Monthly_TargetDomain <- ET_Precip_Monthly_TargetDomain %>%
  mutate(grid_id = row_number())


openet_annual_et <- ET_Precip_Monthly_TargetDomain %>%
  select(grid_id, ET_ET_Data) %>%
  unnest(ET_ET_Data) %>%
  group_by(grid_id, Year) %>%
  summarise(annual_et = sum(et_ensemble_sam, na.rm = TRUE), .groups = "drop") %>%
  group_by(Year) %>%
  summarise(ET_OpenET_mm = mean(annual_et, na.rm = TRUE), .groups = "drop")




seasonal_et_area_avg <- ET_Precip_Monthly_TargetDomain %>%
  select(grid_id, ET_ET_Data) %>%
  unnest(ET_ET_Data) %>%
  mutate(
    season_info = map(Date, assign_season_and_year),
    Season = map_chr(season_info, "season"),
    adjusted_year = map_dbl(season_info, "adjusted_year")
  ) %>%
  group_by(grid_id, adjusted_year, Season) %>%
  summarise(seasonal_et_mm = sum(et_ensemble_sam, na.rm = TRUE), .groups = "drop") %>%
  group_by(adjusted_year, Season) %>%
  summarise(avg_et_mm = mean(seasonal_et_mm, na.rm = TRUE), .groups = "drop")






seasonal_precip_area_avg <- seasonal_precip_area_avg %>%
  st_drop_geometry()

seasonal_et_area_avg <- seasonal_et_area_avg %>%
  st_drop_geometry()

seasonal_water_balance <- seasonal_precip_area_avg %>%
  full_join(seasonal_runoff_desoto, by = c("adjusted_year", "Season"))

seasonal_water_balance <- seasonal_water_balance %>%
  arrange(adjusted_year, Season)
seasonal_water_balance <- seasonal_water_balance %>%
  mutate(
    avg_et_mm_FromEQ_PMinusR = avg_precip_mm - seasonal_runoff_mm
  )


seasonal_water_balance_openET<- na.omit(seasonal_water_balance)
seasonal_water_balance_openET<- seasonal_water_balance_openET %>% left_join(seasonal_et_area_avg,by = c("adjusted_year", "Season"))
seasonal_water_balance_openET<- na.omit(seasonal_water_balance_openET)




seasonal_water_balance<- seasonal_water_balance %>% select(adjusted_year,Season ,avg_precip_mm,seasonal_runoff_mm,avg_et_mm_FromEQ_PMinusR )
seasonal_water_balance<- na.omit(seasonal_water_balance)




seasonal_water_balance_long_openET <- seasonal_water_balance_openET %>%
  select(adjusted_year, Season, avg_et_mm_FromEQ_PMinusR, avg_et_mm) %>%
  pivot_longer(cols = c(avg_et_mm_FromEQ_PMinusR, avg_et_mm),
               names_to = "ET_Type",
               values_to = "ET_Value") %>%
  mutate(ET_Type = recode(ET_Type,
                          avg_et_mm_FromEQ_PMinusR = "P - R (Water Balance)",
                          avg_et_mm = "OpenET"))

#Summer highest, because of irrigation?? 

ggplot(seasonal_water_balance_long_openET, aes(x = adjusted_year, y = ET_Value, color = ET_Type)) +
  geom_line() +
  geom_point() +
  facet_wrap(~Season, scales = "fixed") +
  labs(title = "Seasonal ET Comparison: OpenET vs Water Balance (P - R)",
       x = "Year",
       y = "ET (mm)",
       color = "ET Estimate") +
  theme_minimal(base_size = 14)

plot_seasonal_water_balance_comparison <- function(dry_years, wet_years, data = seasonal_water_balance) {
  data_with_periods <- data %>%
    mutate(period = case_when(
      adjusted_year %in% dry_years ~ paste0("Dry (", min(dry_years), "–", max(dry_years), ")"),
      adjusted_year %in% wet_years ~ paste0("Wet (", min(wet_years), "–", max(wet_years), ")"),
      TRUE ~ NA_character_
    )) %>%
    filter(!is.na(period))
  
  seasonal_summary <- data_with_periods %>%
    group_by(period, Season) %>%
    summarise(
      Precipitation = sum(avg_precip_mm, na.rm = TRUE),
      Runoff = sum(seasonal_runoff_mm, na.rm = TRUE),
      ET = sum(avg_et_mm_FromEQ_PMinusR, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_longer(cols = c("Precipitation", "Runoff", "ET"),
                 names_to = "Component", values_to = "Value")
  
  seasonal_summary$Season <- factor(seasonal_summary$Season,
                                    levels = c("Winter", "Spring", "Summer", "Fall"))
  seasonal_summary$period <- factor(seasonal_summary$period,
                                    levels = unique(seasonal_summary$period))
  
  period_labels <- c(
    dry = paste0("Dry (", min(dry_years), "–", max(dry_years), ")"),
    wet = paste0("Wet (", min(wet_years), "–", max(wet_years), ")")
  )
  
  fill_colors <- c(
    setNames("#d95f02", period_labels["dry"]),
    setNames("blue", period_labels["wet"])
  )
  
  ggplot(seasonal_summary, aes(x = Season, y = Value, fill = period)) +
    geom_col(position = "dodge") +
    facet_wrap(~ Component, scales = "fixed") +
    labs(
      title = "Seasonal Water Balance Component Sums",
      x = "Season", y = "Total Seasonal Sum (mm)", fill = "Period"
    ) +
    scale_fill_manual(values = fill_colors) +
    theme_minimal(base_size = 14)
}



plot_seasonal_water_balance_comparison(
  dry_years = 1961:1980,
  wet_years = 1981:2001
)

plot_seasonal_water_balance_comparison(
  dry_years = 1930:1950,
  wet_years = 1951:1970
)

plot_seasonal_water_balance_comparison(
  dry_years = 1950:1969,
  wet_years = 1990:2009
)

plot_seasonal_water_balance_comparison(
  dry_years = 1981:2000,
  wet_years = 2001:2020
)

plot_seasonal_water_balance_comparison_OpenET <- function(dry_years, wet_years, data = seasonal_water_balance_openET) {
  data_with_periods <- data %>%
    mutate(period = case_when(
      adjusted_year %in% dry_years ~ paste0("Dry (", min(dry_years), "–", max(dry_years), ")"),
      adjusted_year %in% wet_years ~ paste0("Wet (", min(wet_years), "–", max(wet_years), ")"),
      TRUE ~ NA_character_
    )) %>%
    filter(!is.na(period))
  
  seasonal_summary <- data_with_periods %>%
    group_by(period, Season) %>%
    summarise(
      Precipitation = sum(avg_precip_mm, na.rm = TRUE),
      Runoff = sum(seasonal_runoff_mm, na.rm = TRUE),
      ET = sum(avg_et_mm, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_longer(cols = c("Precipitation", "Runoff", "ET"),
                 names_to = "Component", values_to = "Value")
  
  seasonal_summary$Season <- factor(seasonal_summary$Season,
                                    levels = c("Winter", "Spring", "Summer", "Fall"))
  seasonal_summary$period <- factor(seasonal_summary$period,
                                    levels = unique(seasonal_summary$period))
  
  period_labels <- c(
    dry = paste0("Dry (", min(dry_years), "–", max(dry_years), ")"),
    wet = paste0("Wet (", min(wet_years), "–", max(wet_years), ")")
  )
  
  fill_colors <- c(
    setNames("#d95f02", period_labels["dry"]),
    setNames("blue", period_labels["wet"])
  )
  
  ggplot(seasonal_summary, aes(x = Season, y = Value, fill = period)) +
    geom_col(position = "dodge") +
    facet_wrap(~ Component, scales = "fixed") +
    labs(
      title = "Seasonal Water Balance Component Sums (OpenET)",
      x = "Season", y = "Total Seasonal Sum (mm)", fill = "Period"
    ) +
    scale_fill_manual(values = fill_colors) +
    theme_minimal(base_size = 14)
}

plot_seasonal_water_balance_comparison_OpenET(
  dry_years = 2007:2011,
  wet_years = 2018:2022
)

#paper recreate walnut

site_no <- "07147800"

walnut_data <- readNWISdv(siteNumber = site_no, parameterCd = "00060") %>%
  select(Date, X_00060_00003) %>%
  rename(mean_streamflow = X_00060_00003)

fill_missing_dates <- function(data) {
  full_dates <- tibble(Date = seq(min(data$Date), max(data$Date), by = "day"))
  full_dates %>%
    left_join(data, by = "Date")
}

walnut_data_full <- fill_missing_dates(walnut_data)

site_info <- readNWISsite(site_no)
drain_area_sqmi <- site_info$drain_area_va
drain_area_m2 <- conv_unit(drain_area_sqmi, from = "mi2", to = "m2")
convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * 0.0283168 * 86400
  (streamflow_m3_per_day / watershed_area_m2) * 1000
}

walnut_data_full <- walnut_data_full %>%
  mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, drain_area_m2))
convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * 0.0283168 * 86400
  (streamflow_m3_per_day / watershed_area_m2) * 1000
}

walnut_data_full <- walnut_data_full %>%
  mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, drain_area_m2))
walnut_station <- tibble(
  station_name = site_info$station_nm,
  site_no = site_info$site_no,
  station_lat = site_info$dec_lat_va,
  station_lon = site_info$dec_long_va,
  streamflow_data_CFS_MMD = list(walnut_data_full),
  streamflowUnit = "MMD"
)

walnut_station <- walnut_station %>%
  mutate(StreamflowMMD = map(streamflow_data_CFS_MMD, ~ select(.x, Date, mean_streamflow_mm_per_day)),
         streamflow_data = StreamflowMMD) %>%
  select(-streamflow_data_CFS_MMD, -StreamflowMMD)



site_no <- "07147800"
site_nldi <- list(featureSource = "nwissite", featureID = paste0("USGS-", site_no))
walnut_basin <- get_nldi_basin(site_nldi)

kansas_boundary <- sf::st_read(file.path(file_Path_Variable_O, "kansas_boundary.shp"))

Reading layer `kansas_boundary' from data source 
  `C:\Users\a905h226\OneDrive - University of Kansas\Desktop\Steps_Workflow_Sept17\Output\kansas_boundary.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 9 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 229266.7 ymin: 4094787 xmax: 890007.5 ymax: 4434288
Projected CRS: WGS 84 / UTM zone 14N

kansas<- kansas_boundary
ggplot() +
  geom_sf(data = kansas, fill = "white", color = "black") +
  geom_sf(data = walnut_basin, fill = "lightblue", color = "blue", alpha = 0.5) +
  labs(title = "Walnut River Watershed", subtitle = "USGS 07147800 Basin in Kansas") +
  theme_minimal()

walnut_basin <- st_transform(walnut_basin, crs = 32614)

Precipitation_Data_Gridded_Filtered_TargetDomain <- Precipitation_Data_Gridded_Filtered %>%
  st_as_sf(coords = c("x", "y"), crs = 32614)

Precipitation_Walnut <- st_intersection(Precipitation_Data_Gridded_Filtered_TargetDomain, walnut_basin)

annual_precip_trend_walnut <- Precipitation_Walnut %>%
  unnest(annual_precip_data) %>%
  group_by(year) %>%
  summarise(avg_precip_mm = mean(precip_mm, na.rm = TRUE), .groups = "drop")

annual_runoff_walnut <- walnut_data_full %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    mean_runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
    .groups = "drop"
  )

walnut_waterbalance <- annual_runoff_walnut %>%
  left_join(annual_precip_trend_walnut, by = c("Year" = "year")) %>%
  mutate(ET_mm_PminusR = avg_precip_mm - mean_runoff_mm) %>%
  select(Year, avg_precip_mm, mean_runoff_mm, ET_mm_PminusR)



ET_Precip_Monthly_Walnut <- ET_Precip_Monthly %>%
  rename(x = et_x, y = et_y) %>%
  st_as_sf(coords = c("x", "y"), crs = 32614) %>%
  st_intersection(walnut_basin)  # Only retain points within Walnut basin

openet_annual_et_walnut <- ET_Precip_Monthly_Walnut %>%
  unnest(ET_ET_Data) %>%
  group_by(nearest_index, Year) %>%
  summarise(annual_et = sum(et_ensemble_sam, na.rm = TRUE), .groups = "drop") %>%
  group_by(Year) %>%
  summarise(ET_OpenET_mm = mean(annual_et, na.rm = TRUE), .groups = "drop")










walnut_waterbalance <- walnut_waterbalance %>%
  left_join(openet_annual_et_walnut, by = "Year")

walnut_waterbalance_Et <- walnut_waterbalance %>%
  select(Year, ET_OpenET_mm, ET_mm_PminusR) %>%
  na.omit()



walnut_waterbalance_Et_long <- walnut_waterbalance_Et %>%
  pivot_longer(cols = c(ET_OpenET_mm, ET_mm_PminusR),
               names_to = "ET_source",
               values_to = "ET_mm")

ggplot(walnut_waterbalance_Et_long, aes(x = Year, y = ET_mm, color = ET_source)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  scale_color_manual(values = c("ET_OpenET_mm" = "blue", "ET_mm_PminusR" = "darkgreen"),
                     labels = c("OpenET", "P - Runoff")) +
  labs(
    title = "Annual Evapotranspiration (ET) in Walnut River Basin",
    subtitle = "Comparison of OpenET and Water Balance (P - Runoff) Estimates",
    x = "Year",
    y = "ET (mm)",
    color = "ET Source"
  ) +
  theme_minimal(base_size = 14)

ggplot(walnut_waterbalance_Et, aes(x = ET_OpenET_mm, y = ET_mm_PminusR)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "darkred", linetype = "dashed") +
  labs(
    title = "Comparison of ET Estimates in Walnut: OpenET vs P - Runoff",
    x = "ET from OpenET (mm/year)",
    y = "ET from P - Runoff (mm/year)"
  ) +
  theme_minimal(base_size = 14)

site_no <- "07147800"

walnut_data <- readNWISdv(siteNumber = site_no, parameterCd = "00060") %>%
  select(Date, X_00060_00003) %>%
  rename(mean_streamflow = X_00060_00003)

# Fill missing dates
fill_missing_dates <- function(data) {
  full_dates <- tibble(Date = seq(min(data$Date), max(data$Date), by = "day"))
  full_dates %>%
    left_join(data, by = "Date")
}

walnut_data_full <- fill_missing_dates(walnut_data)

# Get watershed area
site_info <- readNWISsite(site_no)
drain_area_sqmi <- site_info$drain_area_va
drain_area_m2 <- conv_unit(drain_area_sqmi, from = "mi2", to = "m2")
convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * 0.0283168 * 86400
  (streamflow_m3_per_day / watershed_area_m2) * 1000
}

walnut_data_full <- walnut_data_full %>%
  mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, drain_area_m2))
convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * 0.0283168 * 86400
  (streamflow_m3_per_day / watershed_area_m2) * 1000
}

walnut_data_full <- walnut_data_full %>%
  mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, drain_area_m2))
walnut_station <- tibble(
  station_name = site_info$station_nm,
  site_no = site_info$site_no,
  station_lat = site_info$dec_lat_va,
  station_lon = site_info$dec_long_va,
  streamflow_data_CFS_MMD = list(walnut_data_full),
  streamflowUnit = "MMD"
)

walnut_station <- walnut_station %>%
  mutate(StreamflowMMD = map(streamflow_data_CFS_MMD, ~ select(.x, Date, mean_streamflow_mm_per_day)),
         streamflow_data = StreamflowMMD) %>%
  select(-streamflow_data_CFS_MMD, -StreamflowMMD)






walnut_basin <- st_transform(walnut_basin, crs = 32614)

Precipitation_Data_Gridded_Filtered_TargetDomain <- Precipitation_Data_Gridded_Filtered %>%
  st_as_sf(coords = c("x", "y"), crs = 32614)

Precipitation_Walnut <- st_intersection(Precipitation_Data_Gridded_Filtered_TargetDomain, walnut_basin)

annual_precip_trend_walnut <- Precipitation_Walnut %>%
  unnest(annual_precip_data) %>%
  group_by(year) %>%
  summarise(avg_precip_mm = mean(precip_mm, na.rm = TRUE), .groups = "drop")

annual_runoff_walnut <- walnut_data_full %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    mean_runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
    .groups = "drop"
  )

walnut_waterbalance <- annual_runoff_walnut %>%
  left_join(annual_precip_trend_walnut, by = c("Year" = "year")) %>%
  mutate(ET_mm_PminusR = avg_precip_mm - mean_runoff_mm) %>%
  select(Year, avg_precip_mm, mean_runoff_mm, ET_mm_PminusR)



ET_Precip_Monthly_Walnut <- ET_Precip_Monthly %>%
  rename(x = et_x, y = et_y) %>%
  st_as_sf(coords = c("x", "y"), crs = 32614) %>%
  st_intersection(walnut_basin)  # Only retain points within Walnut basin

openet_annual_et_walnut <- ET_Precip_Monthly_Walnut %>%
  unnest(ET_ET_Data) %>%
  group_by(nearest_index, Year) %>%
  summarise(annual_et = sum(et_ensemble_sam, na.rm = TRUE), .groups = "drop") %>%
  group_by(Year) %>%
  summarise(ET_OpenET_mm = mean(annual_et, na.rm = TRUE), .groups = "drop")










walnut_waterbalance <- walnut_waterbalance %>%
  left_join(openet_annual_et_walnut, by = "Year")

walnut_waterbalance_Et <- walnut_waterbalance %>%
  select(Year, ET_OpenET_mm, ET_mm_PminusR) %>%
  na.omit()



walnut_waterbalance_Et_long <- walnut_waterbalance_Et %>%
  pivot_longer(cols = c(ET_OpenET_mm, ET_mm_PminusR),
               names_to = "ET_source",
               values_to = "ET_mm")

Walnut River Basin: 1961–1980 vs 1981–2001
Hydrologic metrics based on dry and wet years
Precip(Dry)	Precip(Wet)	ΔPrecip	%ΔPrecip	Q(Dry)	Q(Wet)	ΔRunoff	%ΔRunoff	ET(Dry)	ET(Wet)	ΔET	%ΔET
849.33	872.48	23.15	2.73	168.54	199.10	30.56	18.13	680.80	673.38	−7.42	−1.09

Walnut River Basin: 1961–1980 vs 1981–2001
Derived Metrics
Runoff/Precip	Retention Ratio	QSensitivity	ETSensitivity
0.21	0.79	6.65	−0.40

The figure below shows the replicated results for Table 2 in the Walnut River Basin.

why different than the paper

Let’s follow exact footsteps (precip from https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/divisional/time-series/1409/pcp/1/0/1981-2001)

climate_divisions <- st_read(
  file.path(file_Path_Variable_I, "CONUS_CLIMATE_DIVISIONS_shp/GIS.OFFICIAL_CLIM_DIVISIONS.shp")
)

Reading layer `GIS.OFFICIAL_CLIM_DIVISIONS' from data source 
  `C:\Users\a905h226\OneDrive - University of Kansas\Desktop\Steps_Workflow_Sept17\InputFiles\CONUS_CLIMATE_DIVISIONS_shp\GIS.OFFICIAL_CLIM_DIVISIONS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 344 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -124.7332 ymin: 24.54394 xmax: -66.95 ymax: 49.38436
Geodetic CRS:  Longitude / Latitude (NAD 83)

kansas_boundary <- st_transform(kansas_boundary, crs = 4326)
climate_divisions <- st_transform(climate_divisions, crs = 4326)

site_no <- "07147800"
site_nldi <- list(featureSource = "nwissite", featureID = paste0("USGS-", site_no))
walnut_basin <- get_nldi_basin(site_nldi)
walnut_basin <- st_transform(walnut_basin, crs = 4326)


climate_divisions_ks <- climate_divisions %>%
  filter(STATE == "Kansas")


DryPeriodNoaaClimate0<- read.csv(file.path(file_Path_Variable_I,"NOAAClimateDivision_Precip_19611980.csv"))
WetPeriodNoaaClimate0<- read.csv(file.path(file_Path_Variable_I,"NOAAClimateDivision_Precip_19812001.csv"))

WetPeriodNoaaClimate <- read_csv(
  file.path(file_Path_Variable_I, "NOAAClimateDivision_Precip_19812001.csv"),
  skip = 3,  # Skip metadata rows
  col_names = c("Date", "Precip_in")
) %>%
  filter(!is.na(as.numeric(Precip_in))) %>%  # Remove bad rows like headers
  mutate(
    Date = as.character(Date),
    Year = substr(Date, 1, 4),
    Month = substr(Date, 5, 6),
    Date = paste(Year, Month, sep = "-"),
    Precip_in = as.numeric(Precip_in)
  ) %>%
  select(Date, Precip_in)

DryPeriodNoaaClimate <- read_csv(
  file.path(file_Path_Variable_I, "NOAAClimateDivision_Precip_19611980.csv"),
  skip = 3,
  col_names = c("Date", "Precip_in")
) %>%
  filter(!is.na(as.numeric(Precip_in))) %>%
  mutate(
    Date = as.character(Date),
    Year = substr(Date, 1, 4),
    Month = substr(Date, 5, 6),
    Date = paste(Year, Month, sep = "-"),
    Precip_in = as.numeric(Precip_in)
  ) %>%
  select(Date, Precip_in)




climate_combined <- bind_rows(
  mutate(DryPeriodNoaaClimate, Period = "Dry"),
  mutate(WetPeriodNoaaClimate, Period = "Wet")
)

climate_combined <- climate_combined %>%
  mutate(
    Year = as.integer(substr(Date, 1, 4)),
    Month = as.integer(substr(Date, 6, 7)),
    Precip_mm = Precip_in * 25.4
  )

annual_precip_noaa <- climate_combined %>%
  group_by(Year) %>%
  summarise(avg_precip_mm = sum(Precip_mm, na.rm = TRUE), .groups = "drop")


site_no <- "07147800"

walnut_data <- readNWISdv(siteNumber = site_no, parameterCd = "00060") %>%
  select(Date, X_00060_00003) %>%
  rename(mean_streamflow = X_00060_00003)

fill_missing_dates <- function(data) {
  full_dates <- tibble(Date = seq(min(data$Date), max(data$Date), by = "day"))
  full_dates %>%
    left_join(data, by = "Date")
}

walnut_data_full <- fill_missing_dates(walnut_data)

site_info <- readNWISsite(site_no)
drain_area_sqmi <- site_info$drain_area_va
drain_area_m2 <- conv_unit(drain_area_sqmi, from = "mi2", to = "m2")
convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * 0.0283168 * 86400
  (streamflow_m3_per_day / watershed_area_m2) * 1000
}

walnut_data_full <- walnut_data_full %>%
  mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, drain_area_m2))
convert_cfs_to_mm_per_day <- function(cfs_streamflow, watershed_area_m2) {
  streamflow_m3_per_day <- cfs_streamflow * 0.0283168 * 86400
  (streamflow_m3_per_day / watershed_area_m2) * 1000
}

walnut_data_full <- walnut_data_full %>%
  mutate(mean_streamflow_mm_per_day = convert_cfs_to_mm_per_day(mean_streamflow, drain_area_m2))
walnut_station <- tibble(
  station_name = site_info$station_nm,
  site_no = site_info$site_no,
  station_lat = site_info$dec_lat_va,
  station_lon = site_info$dec_long_va,
  streamflow_data_CFS_MMD = list(walnut_data_full),
  streamflowUnit = "MMD"
)

walnut_station <- walnut_station %>%
  mutate(StreamflowMMD = map(streamflow_data_CFS_MMD, ~ select(.x, Date, mean_streamflow_mm_per_day)),
         streamflow_data = StreamflowMMD) %>%
  select(-streamflow_data_CFS_MMD, -StreamflowMMD)

##
annual_runoff_walnut <- walnut_data_full %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    mean_runoff_mm = mean(mean_streamflow_mm_per_day, na.rm = TRUE) * 365.25,
    .groups = "drop"
  )

walnut_waterbalance <- annual_runoff_walnut %>%
  left_join(annual_precip_noaa, by = "Year") %>%
  mutate(ET_mm_PminusR = avg_precip_mm - mean_runoff_mm) %>%
  select(Year, avg_precip_mm, mean_runoff_mm, ET_mm_PminusR)


walnut_waterbalance<- na.omit(walnut_waterbalance)

walnut_decadal_tables <- make_decadal_summary_table(
  data = walnut_waterbalance,
  dry_years = 1961:1980,
  wet_years = 1981:2001,
  title_text = "Walnut River Basin using noaa : 1961–1980 vs 1981–2001"
)

walnut_decadal_tables$Hydrologic_Table

Walnut River Basin using noaa : 1961–1980 vs 1981–2001
Hydrologic metrics based on dry and wet years
Precip(Dry)	Precip(Wet)	ΔPrecip	%ΔPrecip	Q(Dry)	Q(Wet)	ΔRunoff	%ΔRunoff	ET(Dry)	ET(Wet)	ΔET	%ΔET
935.65	1,027.25	91.60	9.79	168.54	199.10	30.56	18.13	767.11	828.15	61.04	7.96

walnut_decadal_tables$Derived_Metrics_Table

Walnut River Basin using noaa : 1961–1980 vs 1981–2001
Derived Metrics
Runoff/Precip	Retention Ratio	QSensitivity	ETSensitivity
0.19	0.81	1.85	0.81

kansas_boundary <- st_transform(kansas_boundary, crs = 4326)
climate_divisions <- st_transform(climate_divisions, crs = 4326)

site_no <- "07147800"
site_nldi <- list(featureSource = "nwissite", featureID = paste0("USGS-", site_no))
walnut_basin <- get_nldi_basin(site_nldi)
walnut_basin <- st_transform(walnut_basin, crs = 4326)


climate_divisions_ks <- climate_divisions %>%
  filter(STATE == "Kansas")

ggplot() +
  geom_sf(data = climate_divisions_ks, aes(fill = NAME), color = "gray40", size = 0.6) +
  geom_sf(data = kansas_boundary, fill = NA, color = "black", size = 0.8) +
  geom_sf(data = walnut_basin, fill = "lightblue", color = "blue", alpha = 0.5, size = 1) +
  geom_sf_text(data = climate_divisions_ks, aes(label = NAME), size = 3, color = "black") +
  labs(
    title = "Walnut River Basin with Kansas Climate Divisions",
    subtitle = "USGS 07147800 and NOAA Division Overlay",
    fill = "Climate Division",
    caption = "Climate division boundaries: NOAA | Basin: NLDI | State: US Census TIGER"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")