Introduce sample data sets and gages

USGS 01634000 Strasburg VA

Daymet: http://deq1.bse.vt.edu:81/files/met/usgs_ws_01634000-daymet-2022-2023.csv
PRISM: http://deq1.bse.vt.edu:81/files/met/usgs_ws_01634000-prism-2022-2023.csv
NLDAS: tbd
Rscript to match precip and USGS: https://github.com/HARPgroup/HARParchive/blob/master/HARP-2024-2025/access-file.R
- This script does various summaries of daily precip and flow, change in flow, change in next day flow
- Has variations for weekly flow, and analyzing a single month.
- Introduces sqldf for data frame summaries
- Uses dataRetrieval module
- Introduces use of lm for simple linear regression, and plotting of model regressions.
- Question 1: "why do certain days repeat non-zero values?"
  - see 2022/03/22 and 2022/03/23 data for repeating daymet and prism values
  - http://deq1.bse.vt.edu:81/met/daymet/2022/081/ORIGINAL_daymet_precip2022-03-22.bil_CBP.gtiff
  - http://deq1.bse.vt.edu:81/met/PRISM/2022/082/PRISM_ppt_stable_4kmD2_20220323_bil.bil_CBP.gtiff
  - http://deq1.bse.vt.edu:81/met/daymet/2022/081/ORIGINAL_daymet_precip2022-03-22.bil_CBP.gtiff
  - http://deq1.bse.vt.edu:81/met/daymet/2022/082/ORIGINAL_daymet_precip2022-03-23.bil_CBP.gtiff
- Question 2: "can we do a weekly precip versus change in flow?"
- Question 3: "how about weekly change in flow versus previous week precip?"

The storm separation code we developed in 2018 is here. Glancing through it, I believe it may be a good first step to working through this. I think the logic here is sound, but it needs some more recent testing. I'm going to try spinning it up early next week and will report back. For now, I'm going to pull baseflow separation from grwrat, a package that I like the looks of as it has storm separation functions in the package itself (although they seem complicated). Doesn't have to be our final answer on this, but it gives me a starting point.

[like] Scott, Durelle reacted to your message:

From: C. Brogan @.> Sent: Friday, June 7, 2024 4:30:50 PM To: HARPgroup/HARParchive @.> Cc: Scott, Durelle @.>; Assign @.> Subject: Re: [HARPgroup/HARParchive] Introduce sample data sets and gages (Issue #1270)

The storm separation code we developed in 2018 is herehttps://github.com/HARPgroup/HARParchive/tree/cd4c0380b65e61a6c0f0fea42dfa48cff369d99c/HARP-2018-2019/Flow%20Stats. Glancing through it, I believe it may be a good first step to working through this. I think the logic here is sound, but it needs some more recent testing. I'm going to try spinning it up early next week and will report back. For now, I'm going to pull baseflow separation from grwrat, a package https://tsamsonov.github.io/grwat/articles/baseflow.html that I like the looks of as it has storm separation functions in the package itself (although they seem complicated). Doesn't have to be our final answer on this, but it gives me a starting point.

— Reply to this email directly, view it on GitHubhttps://github.com/HARPgroup/HARParchive/issues/1270#issuecomment-2155166503, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AC4K427A765OWLXW4BBWGA3ZGHN3VAVCNFSM6AAAAABI5NGESSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJVGE3DMNJQGM. You are receiving this because you were assigned.Message ID: @.***>

R Code to look at changes in Q, on days when Q is increasing over prior day (evidence of rain increase).

# *** PRISM - Flow Change tomorrow, ONLY for days where Flow Change > 0, month of February only
mod_prism_mon_nz_ndd <- lm(nextday_d_cfs ~ prism_p_cfs, data=comp_data[which((comp_data$mo == 2) & (comp_data$nextday_d_cfs > 0)),])
summary(mod_prism_mon_nz_ndd)
plot(mod_prism_mon_nz_ndd$model$nextday_d_cfs ~ mod_prism_mon_nz_ndd$model$prism_p_cfs)

Look at Every month and plot their R-sq by month

# do all months and assemble a barplot of R^2
ndd_stats <- data.frame(row.names=c('month', 'rsquared_a'))
for (i in 1:12) {
  mod_prism_mon_nz_ndd <- lm(nextday_d_cfs ~ prism_p_cfs, data=comp_data[which((comp_data$mo == i) & (comp_data$nextday_d_cfs > 0)),])
  dsum <- summary(mod_prism_mon_nz_ndd)
  plot(mod_prism_mon_nz_ndd$model$nextday_d_cfs ~ mod_prism_mon_nz_ndd$model$prism_p_cfs)
  ndd_stats <- rbind(ndd_stats, data.frame(i, dsum$adj.r.squared))
}
barplot(ndd_stats$dsum.adj.r.squared ~ ndd_stats$i)

Week data

# do all months and assemble a barplot of R^2
nwd_stats <- data.frame(row.names=c('month', 'rsquared_a'))
for (i in 1:12) {
  # Weekly d cfs vs P
  mod_weekmo_prism_cfs <- lm(usgs_cfs ~ prism_p_cfs, data=week_data[which((week_data$mo == i)),])
  dsum <- summary(mod_weekmo_prism_cfs)
  plot(mod_weekmo_prism_cfs$model$usgs_cfs ~ mod_weekmo_prism_cfs$model$prism_p_cfs)

  #mod_week_daymet_d_cfs <- lm(today_d_cfs ~ daymet_p_cfs, data=week_data)
  #summary(mod_week_daymet_d_cfs)
  #plot(mod_week_daymet_d_cfs$model$today_d_cfs ~ mod_week_daymet_d_cfs$model$daymet_p_cfs)

  nwd_stats <- rbind(nwd_stats, data.frame(i, dsum$adj.r.squared))
}
barplot(nwd_stats$dsum.adj.r.squared ~ nwd_stats$i, main=paste(gage_info$station_nm ))

I tried two different ways to look at weekly data, one of the next_day flow and the other just usgs_cfs but for the week. In the end the graphs and r^2 is still the same for all the methods of looking at the months, could be an example of too many extra variables.

# do all months and assemble a barplot of R^2, looking at weekly average
week_ndd_stats <- data.frame(row.names=c('month', 'rsquared_a'))
for (i in 1:12) {
  mod_week_prism_mon_nz_ndd <- lm(usgs_cfs ~ prism_p_cfs, data=week_data[which((week_data$mo == i) & (week_data$usgs_cfs > 0)),])
  week_dsum <- summary(mod_week_prism_mon_nz_ndd)
  plot(mod_week_prism_mon_nz_ndd$model$usgs_cfs ~ mod_week_prism_mon_nz_ndd$model$prism_p_cfs)
  week_ndd_stats <- rbind(week_ndd_stats, data.frame(i, week_dsum$adj.r.squared))
}
barplot(week_ndd_stats$week_dsum.adj.r.squared ~ week_ndd_stats$i)
summary(mod_week_prism_mon_nz_ndd)

# do all months and assemble a barplot of R^2, seperated by week, and looking at next day flow
nex_week_ndd_stats <- data.frame(row.names=c('month', 'rsquared_a'))
for (i in 1:12) {
  nex_mod_week_prism_mon_nz_ndd <- lm(nextday_d_cfs ~ prism_p_cfs, data=week_data[which((week_data$mo == i) & (week_data$nextday_d_cfs > 0)),])
  nex_week_dsum <- summary(nex_mod_week_prism_mon_nz_ndd)
  plot(nex_mod_week_prism_mon_nz_ndd$model$nextday_d_cfs ~ nex_mod_week_prism_mon_nz_ndd$model$prism_p_cfs)
  nex_week_ndd_stats <- rbind(nex_week_ndd_stats, data.frame(i, nex_week_dsum$adj.r.squared))
}
barplot(nex_week_ndd_stats$nex_week_dsum.adj.r.squared ~ nex_week_ndd_stats$i)
summary(mod_week_prism_mon_nz_ndd)
mod_week_prism <- lm(usgs_cfs ~ prism_p_cfs, data=week_data)

fb09d129-1a3b-4ef6-b05a-58b4fa5f600e c08f13f4-d4aa-45d3-b595-12955aabc597

I also created a look at last weeks precipitation compared to this weeks streamflow, a similar overall graph, but in certain parts of the year it seems to have a much stronger correlation

last_week_data <- sqldf(
  "select min(week_begin) as week_begin, yr, wk, mo, min(dataset_day_begin) as dataset_day_begin,
     Lag(daymet_p_in,1) OVER(ORDER BY dataset_day_begin) as lastweek_daymet_p_in, 
     Lag(daymet_p_cfs,1) OVER(ORDER BY dataset_day_begin) as lastweek_daymet_p_cfs,
     Lag(prism_p_in,1) OVER(ORDER BY dataset_day_begin) as lastweek_prism_p_in, 
     Lag(prism_p_cfs,1) OVER(ORDER BY dataset_day_begin) as lastweek_prism_p_cfs,
     usgs_cfs as usgs_cfs, today_d_cfs as today_d_cfs, 
     nextday_d_cfs as nextday_d_cfs
   from week_data
   group by yr, wk
   order by yr, wk
  "
)

last_week_data = last_week_data[-1,]

las_week_ndd_stats <- data.frame(row.names=c('month', 'rsquared_a'))
for (i in 1:12) {
  las_week_prism_mon_nz_ndd <- lm(usgs_cfs ~ lastweek_prism_p_cfs, data=last_week_data[which((last_week_data$mo == i) & (last_week_data$usgs_cfs > 0)),])
  las_week_dsum <- summary(las_week_prism_mon_nz_ndd)
  plot(las_week_prism_mon_nz_ndd$model$usgs_cfs ~ las_week_prism_mon_nz_ndd$model$lastweek_prism_p_cfs)
  las_week_ndd_stats <- rbind(las_week_ndd_stats, data.frame(i, las_week_dsum$adj.r.squared))
}
barplot(las_week_ndd_stats$las_week_dsum.adj.r.squared ~ las_week_ndd_stats$i)
summary(las_week_prism_mon_nz_ndd)

3c428660-f4d7-4583-b9b8-ed6ebd25e10e

HARPgroup / HARParchive