Open rburghol opened 1 year ago
In order for the Analysts to help summarize the watershed model summary, we have started on a document that outlines our work from the summer and fall. This document will hopefully give background information about the switch from hspf to hsp2, and then more specific information on the steps that we have taken to adapt data from h5's, analyzing it, and how data is in VAHydro. In the appendix we will provide example figures/graphs on examples about the analysis we have done. Finally, a summary on key statistics of the analysts' work and next steps in the project for the spring.
Table of contents for now:
Does this sound like a useful piece of documentation for the final report from the harp project? Is there anything that we should add? We can start working on this, and then at some point in January provide a draft for review.
Yes - I almost capitalized and bolder:). Thanks, S
On Dec 16, 2022, at 4:32 PM, juliabruneau @.***> wrote:
In order for the Analysts to help summarize the watershed model summary, we have started on a document that outlines our work from the summer and fall. This document will hopefully give background information about the switch from hspf to hsp2, and then more specific information on the steps that we have taken to adapt data from h5's, analyzing it, and how data is in VAHydro. In the appendix we will provide example figures/graphs on examples about the analysis we have done. Finally, a summary on key statistics of the analysts' work and next steps in the project for the spring.
Table of contents for now:
Introduction Watershed Modeling by the Virginia Department of Environmental Quality The Hydrological Simulation Program (HSP) 3.1 HSP Model Run and VAHydro Data Processing Workflow 3.1.1. Data from HSP Model Runs 3.1.2. Organizing Data 3.1.3. Summarizing Data 3.1.4. Data in VAHydro (ftables??) Statistics on HARP Analysts’ Work Next Steps in the Project References Appendix A. Example Figures B. List of Scripts and Their Uses Does this sound like a useful piece of documentation for the final report from the harp project? Is there anything that we should add? We can start working on this, and then at some point in January provide a draft for review.
— Reply to this email directly, view it on GitHub https://github.com/HARPgroup/HARParchive/issues/659#issuecomment-1355639634, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC4K423HZNUTJITJOJ7W4OLWNTNVFANCNFSM6AAAAAASQZWWHE. You are receiving this because you are subscribed to this thread.
@juliabruneau Yes, I think this list looks great and will form some excellent documentation of your work on this project!
For the watershed model summaries we've started outlining here, I believe the idea is to come up with a way to produce a dynamic R Markdown, capable of receiving inputs for watershed of interest and model runs of interest, then the script will produce a tidy document characterizing model validation and error analysis comparing scenario runs of HSPF to HSP2, or HSP2 to HSP2, or HSP2 to usgs gage and so on. The R Markdown document would contain maps, figures, tables and metric comparisons.
See latest mapping script: https://github.com/HARPgroup/HARParchive/issues/749
The scenario_vs_scenario.R
file is now renamed in HARParchive/HARP-2022-Summer/AutomatedScripts as ws_model_summary.R
(as discussed in initial meeting).
I have included these stats in addition to the pre-existing l90, l30, l7 & l1 stats:
All of these stats are the ones we have previously calculated and pushed to VAHydro (that rely solely on Qout).
I used the pre-existing riverseg: PM2_4860_4670 and gageid: 01643700 from the rmd file, since I'm not 100% sure how to match gages with river segments.
Automated testing: https://github.com/HARPgroup/vahydro/issues/309
Ran the chunk:
rmarkdown::render('C:/Users/VT_SA/Documents/GitHub/HARParchive/HARP-2022-Summer/AutomatedScripts/ws_model_summary.Rmd',
+ output_file = 'C:/Users/VT_SA/Documents/HARP/MarkdownSummaryTest',
+ params = list( doc_title = ("Test HSP2 Model Summary"),
+ rseg.file.path = c( "/media/model/p6/out/river/hsp2_2022/hydr/JA4_7280_7340_hydrd_wy.csv", "path/to/flow/data" ),
+ rseg.hydrocode = c("JA4_7280_7340","vahydrosw_wshed_JA4_7280_7340"),
+ rseg.ftype = c("cbp60","vahydro"),
+ rseg.model.version = c("cbp-6.0","cbp-6.1"),
+ runid.list = c("hsp2_2022","subshed"),
+ rseg.metric.list = c("7q10", "l90_Qout", "l30_year")
+ )
+ )
Received this error:
Quitting from lines 55-115 (ws_model_summary.Rmd)
Error in fromJSON(rseg_model_info.r) :
json_str must be a character string
This is the source of the error (line 4):
if (exists("json_obj_url")) {
rseg_obj_url.r <- paste(json_obj_url, rseg.model.r$pid, sep="/")
rseg_model_info.r <- ds$auth_read(rseg_obj_url.r, "text/json", "")
rseg_model_info.r <- fromJSON(rseg_model_info.r)
rseg_obj_url.r
[1] "http://deq1.bse.vt.edu:81/d.dh/node/62/7042044"
> rseg_model_info.r <- ds$auth_read(rseg_obj_url.r, "text/json", "")
> rseg_model_info.r
[1] FALSE
These parameters are for running the rmd:
rseg.file.path = c( "/media/model/p6/out/river/hsp2_2022/hydr/JL1_6770_6850_hydrd_wy.csv", "/media/model/p6/out/river/vadeq_2021/hydr/JL1_6770_6850_hydrd_wy.csv" ),
rseg.hydrocode = c("vahydrosw_wshed_JL1_6770_6850","vahydrosw_wshed_JL1_6770_6850"),
rseg.ftype = c("vahydro","vahydro"),
rseg.model.version = c("cbp-6.0","cbp-6.0”),
runid.list = c("hsp2_2022”, ”vadeq_2021”),
Currently using grid.arrange
to display images side-by-side, problems:
ggraw()
downloads a screengrab of the png![ ](url)
generates high-quality images, but they are not grob objectsQuality with GRID.ARRANGE()
:
VS Quality with
![ ](url)
for example
Matching the time series with match()
:
Current loop gives us NA for the sdate
and edate
match()
works as a function / the output
for (r in 1:length(rseg.model.version)) {
rseg.data.r <- data.table::fread(rseg.file.path[r])
assign(paste0('rseg.data.',r),rseg.data.r)
rseg.cols.r <- names(rseg.data.r)
if (r == 1){ data.match <- rseg.data.r$date sdate <- as.Date(min(data.match)) edate <- as.Date(max(data.match)) } else { data.match <- match(data.match, rseg.data.r) sdate <- as.Date(min(data.match)) edate <- as.Date(max(data.match)) } }
for (r in 1:length(rseg.model.version)) { rseg.flow.r <- zoo(as.numeric(as.character(rseg.data.r$Qout)), order.by = rseg.data.r$index) rseg.flow.r <- window(rseg.flow.r, start = sdate, end = edate) mode(rseg.flow.r) <- 'numeric'
assign(paste0('rseg.flow.',r),rseg.flow.r) }
- [x] Investigate `match()` further / another solution?
This is excellent progress @juliabruneau -- I have one suggestion for the width
parameter:
col_width <- 1.0 / num_models
match
seems intriguing, will this just automatically create an overlapping timeseries?The code to use match()
for having identical date ranges for comparison seems to work:
> head(rseg.flow.1)
1984-10-01 11:30:00 1984-10-02 11:30:00 1984-10-03 11:30:00 1984-10-04 11:30:00
178.3537 235.1851 225.2305 199.9935
1984-10-05 11:30:00 1984-10-06 11:30:00
175.5761 156.1415
> head(rseg.flow.2)
1984-10-01 11:30:00 1984-10-02 11:30:00 1984-10-03 11:30:00 1984-10-04 11:30:00
178.3537 235.1851 225.2305 199.9935
1984-10-05 11:30:00 1984-10-06 11:30:00
175.5761 156.1415
> tail(rseg.flow.1)
2005-12-26 11:30:00 2005-12-27 11:30:00 2005-12-28 11:30:00 2005-12-29 11:30:00
1005.2903 835.9024 664.5346 643.3333
2005-12-30 11:30:00 2005-12-31 11:30:00
610.4067 546.9785
> tail(rseg.flow.2)
2005-12-26 11:30:00 2005-12-27 11:30:00 2005-12-28 11:30:00 2005-12-29 11:30:00
1005.2903 835.9024 664.5346 643.3333
2005-12-30 11:30:00 2005-12-31 11:30:00
610.4067 546.9785
I took a random 3rd dataset; it kept the same start date, and changed the end date of all the datasets to 1991:
> tail(rseg.flow.1)
1991-03-29 11:30:00 1991-03-30 11:30:00 1991-03-31 11:30:00 1991-04-01 11:30:00
803.8933 738.3517 470.9966 338.3706
1991-04-02 11:30:00 1991-04-03 11:30:00
270.7234 233.8215
> tail(rseg.flow.2)
1991-03-29 11:30:00 1991-03-30 11:30:00 1991-03-31 11:30:00 1991-04-01 11:30:00
803.8933 738.3517 470.9966 338.3706
1991-04-02 11:30:00 1991-04-03 11:30:00
270.7234 233.8215
> tail(rseg.flow.3)
1991-03-29 11:30:00 1991-03-30 11:30:00 1991-03-31 11:30:00 1991-04-01 11:30:00
803.8933 738.3517 470.9966 338.3706
1991-04-02 11:30:00 1991-04-03 11:30:00
270.7234
The 90-day low flow bar plot generated using the IHA group2()
function (as it is in the output document):
I was able to generate scatterplots for all different low flow scenarios. Right now, the output figure height should still be edited so that the figures aren't as squished. But this is what it looks like right now:
Here's a closer look of just one plot of each:
Update: Created a %-Diff Boxplot:
This loop is not generating unique datasets for some reason:
for (r in 1:length(rseg.model.version)) {
rseg.flow.r <- zoo(as.numeric(as.character(rseg.data.r$Qout)), order.by = rseg.data.r$index)
rseg.flow.r <- window(rseg.flow.r, start = sdate, end = edate)
mode(rseg.flow.r) <- 'numeric'
assign(paste0('rseg.flow.',r),rseg.flow.r)
}
Result - named uniquely but the data is identical from model2:
> head(rseg.flow.1)
1984-10-01 11:30:00 1984-10-02 11:30:00 1984-10-03 11:30:00 1984-10-04 11:30:00 1984-10-05 11:30:00 1984-10-06 11:30:00
178.3537 235.1851 225.2305 199.9935 175.5761 156.1415
> head(rseg.flow.2)
1984-10-01 11:30:00 1984-10-02 11:30:00 1984-10-03 11:30:00 1984-10-04 11:30:00 1984-10-05 11:30:00 1984-10-06 11:30:00
178.3537 235.1851 225.2305 199.9935 175.5761 156.1415
input data in the zoo function:
> head(rseg.data.1$Qout)
[1] 180.3891 235.8186 224.8300 199.4006 175.1088 155.8118
> head(rseg.data.2$Qout)
[1] 178.3537 235.1851 225.2305 199.9935 175.5761 156.1415
> head(rseg.data.1$index)
[1] "1984-10-01 11:30:00" "1984-10-02 11:30:00" "1984-10-03 11:30:00" "1984-10-04 11:30:00" "1984-10-05 11:30:00"
[6] "1984-10-06 11:30:00"
> head(rseg.data.2$index)
[1] "1984-10-01 11:30:00" "1984-10-02 11:30:00" "1984-10-03 11:30:00" "1984-10-04 11:30:00" "1984-10-05 11:30:00"
[6] "1984-10-06 11:30:00"
Scatterplots
length() > 2
, so that Model 1 and Model 2 are compared
plot
legth() == 3
, and thus the comparison between Model 1 and Model 3 are added onto the previously generated plots from the first loop
plot <- plot + ggline() ...
Future Work
dataset
,x
, y
Creating functions would significantly shorten the markdown, since we do have to have many lines within these ggplot figure generation codes (different labels and axis names & a legend if 3 models).
Example JUST to generate l90 figures with 3 Models:
First Loop:
l90.scatter <- ggplot(l90.data, aes(x = get(names(l90.data)[2]), y = get(names(l90.data)[3]),
color = "Model 2")) +
geom_point() +
geom_smooth(formula = y ~ x, method=lm) +
theme_light() +
theme(axis.title.x = element_text(size = 9),
axis.title.y = element_text(size = 9),
legend.position = "none") +
ylab(colnames(l90.data[2])) +
xlab(colnames(l90.data[3])) +
geom_abline(intercept = 0, slope = 1, col = 'lightgray')
l90_diff <- unlist(l90.data[2] - l90.data[3]) # has to unlist in order to plot with ggplot!
l90_diff_pct <- 100.0 * (l90_diff) / l90.data[3]
l90.diff <- ggplot(l90.data, aes(x = get(names(l90.data)[1]), y = l90_diff,
color = "Model 2")) +
geom_point() +
geom_smooth(formula = y ~ x, method=lm) +
theme_light() +
theme(axis.title.x = element_text(size = 9),
axis.title.y = element_text(size = 9),
legend.position = "none") +
ylab(paste("Difference: Model 1 - Model 2")) +
xlab(paste("Year"))
Second Loop:
legend_colors <- c("Model 2" = "blue", "Model 3" = "red")
l90.scatter <- l90.scatter +
geom_point(l90.data, mapping = aes(x = get(names(l90.data)[2]), y = get(names(l90.data)[4]),
color = "Model 3")) +
geom_smooth(formula = y ~ x, method=lm) +
ylab(paste("Model 2 and 3")) +
xlab(paste("Model 1")) +
labs(color = "Model 1 vs:") +
scale_color_manual(values = legend_colors) +
theme(legend.title = element_text(size=10),
legend.text = element_text(size=10),
legend.position = "top",
legend.key.size = unit(0.25,"cm"))
l90_diff.2 <- unlist(l90.data[2] - l90.data[4])
l90_diff_pct.2 <- 100.0 * (l90_diff) / l90.data[4]
l90.diff <- l90.diff +
geom_point(l90.data, mapping = aes(x =get(names(l90.data)[1]), y = l90_diff.2,
color = "Model 3")) +
geom_smooth(formula = y ~ x, method=lm) +
ylab(paste("Model Difference")) +
xlab(paste("Year")) +
labs(color = "Model 1 vs:") +
scale_color_manual(values = legend_colors) +
theme(legend.title = element_text(size=10),
legend.text = element_text(size=10),
legend.position = "top",
legend.key.size = unit(0.25,"cm"))
I've been getting this error from rendering and chunk-running the rmd:
Error in gregexec(pattern = what, x$txt, ...) :
could not find function "gregexec"
and it shows up from this line:
model.info <- model.comparison %>% as.data.frame() %>% add_rownames() %>% qflextable()
I believe originating from the qflextable()
function
From research it seems like the function is in base R, but maybe called by qflextable
?
Additional variables now need to be set within the render command which are used to locate and read in the length needed by @megpritch
these are ftable.submod
, ftable.id
, and seg
more to be added after testing tmrw
Use:
rmarkdown::render([path to rmd], output_file = [output path], params = list( doc_title = ("Test HSP2 Model Summary"), rseg.file.path = c([hydr filepath1] , [hydr filepath2] ), rseg.hydrocode = c([hydrocode1],[hydrocode ]), rseg.ftype = c([ftype1],[ftype2]), rseg.model.version = c([version1],[version2]), runid.list = c([id1],[id2]), rseg.metric.list = c([list of metrics]), ftable.submod = c([submod1],[submod2]), ftable.id = c([fid1],[fid2]), seg = c([seg1],[seg2])))
Examples: (glenn's)
rmarkdown::render('~/Desktop/GitHub/HARParchive/HARP-2022-Summer/AutomatedScripts/ws_model_summary.Rmd', output_file = '~/Desktop/HARPteam22/ws_summ_doc1', params = list( doc_title = ("Test HSP2 Model Summary"), rseg.file.path = c("/media/model/p6/out/river/hsp2_2022/hydr/JA4_7280_7340_hydrd_wy.csv","/media/model/p6/out/river/subsheds/hydr/JA4_7280_7340_hydrd_wy.csv" ), rseg.hydrocode = c("JA4_7280_7340","vahydrosw_wshed_JA4_7280_7340"), rseg.ftype = c("cbp60","vahydro"), rseg.model.version = c("cbp-6.0","cbp-6.1"), runid.list = c("hsp2_2022","subsheds"), rseg.metric.list = c("Qout","l90_Qout","l30_Qout","l07_Qout","l01_Qout","consumptive_use_frac","wd_cumulative_mgd","ps_cumulative_mgd","wd_mgd","ps_mgd"), ftable.submod = c("p6", "p6"), ftable.id = c("vahydro_2022", "vahydro_2022"), seg = c("JA4_7280_7340", "JA4_7280_7340")))
(megan's)
rmarkdown::render('C:/aa_HARP/aa_GitHub/HARParchive/HARP-2022-Summer/AutomatedScripts/ws_model_summary.Rmd',
output_file = 'C:/aa_HARP/MarkdownFTABLEtest', params = list( doc_title = ("Test HSP2 Model Summary"), rseg.file.path = c("/media/model/p6/out/river/hsp2_2022/hydr/JA4_7280_7340_hydrd_wy.csv","/media/model/p6/out/river/subsheds/hydr/JA4_7280_7340_hydrd_wy.csv" ), rseg.hydrocode = c("JA4_7280_7340","vahydrosw_wshed_JA4_7280_7340"), rseg.ftype = c("cbp60","vahydro"), rseg.model.version = c("cbp-6.0","cbp-6.1"), runid.list = c("hsp2_2022","subsheds"), rseg.metric.list = c("Qout","l90_Qout","l30_Qout","l07_Qout","l01_Qout","consumptive_use_frac","wd_cumulative_mgd","ps_cumulative_mgd","wd_mgd","ps_mgd"), ftable.submod = c("p6", "p6"), ftable.id = c("vahydro_2022", "vahydro_2022"), seg = c("JA4_7280_7340", "JA4_7280_7340")))
Still some undesirable overlap on the y-axis with the channel length, but the depths are now more readable. We also removed the extra FTABLE plots showing the channel from different angles.
@megpritch Nice work!
The % difference for low metrics between scenarios (2020, 2040, climate change) was used throughout the 2020 state plan and will likely be incorporated into summary documents for the 2023 project
% difference columns were added to the metrics table at the end of the rmd and compare metric values to the first model supplied in the params (probably 2020 demand scenario -- runid_11)
Examples for 2 and 3-model comparisons
This may be better put in a separate table, and the % difference can be removed entirely if not useful
Test CIA routine on Rappahannock. Looks tight:
rmarkdown::render('/usr/local/home/git/HARParchive/HARP-2022-Summer/AutomatedScripts/ws_model_summary.Rmd',
output_file = '/WorkSpace/modeling/projects/rappahannock/cia/cia_hsp2022_subsheds.docx',
params = list( doc_title = ("Test HSP2 Model Summary"),
rseg.file.path = c("http://deq1.bse.vt.edu:81/p6/out/river/hsp2_2022/hydr/RU5_6030_0001_hydrd_wy.csv", "http://deq1.bse.vt.edu:81/p6/out/river/subsheds/hydr/RU5_6030_0001_hydrd_wy.csv" ),
rseg.hydrocode = c("RU5_6030_0001","vahydrosw_wshed_RU5_6030_0001"),
rseg.ftype = c("cbp60","vahydro"),
rseg.model.version = c("cbp-6.0","cbp-6.1"),
runid.list = c("hsp2_2022","subsheds"),
rseg.metric.list = c("7q10", "l90_Qout", "l30_year")
)
)
Testd beaver creek (subwatershed):
rmarkdown::render('/usr/local/home/git/vahydro/R/OWS_summaries/ws_model_summary.Rmd',
output_file = '/WorkSpace/modeling/projects/james_river/rivanna/beaver_hsp2/cia_vahydro_subsheds.docx',
params = list( doc_title = ("Test HSP2 Model Summary"),
rseg.file.path = c("http://deq1.bse.vt.edu:81/data/proj3/out/runlog600.351963.log", "http://deq1.bse.vt.edu:81/p6/out/river/subsheds/hydr/JL1_6562_6560_hydrd_wy.csv" ),
rseg.hydrocode = c("vahydrosw_wshed_JL1_6562_6560","vahydrosw_wshed_JL1_6562_6560"),
rseg.ftype = c("vahydro","vahydro"),
rseg.model.version = c("vahydro-1.0","cbp-6.1"),
runid.list = c("runid_600","subsheds"),
rseg.metric.list = c("7q10", "l90_Qout", "l30_year")
)
)
Checklist for completing the Rmd.
_hydrd_wy.csv
[ ] Setting the model comparison scatter plots to a 1:1 ratio - while displaying the entire figure well (labels, readability)
I want to think of how we'd like to present model results to characterize the validation of the new hsp2 system as well as general model error analysis. Last WSP we did a basic error analysis in an appendix, based on error between models and usgs gages but got feedback from public comment period that it would benefit from more context, graphics, etc. We also did appendices with "minor basin" summaries, tho these were heavy on narrative and light on data (just summary flow stats).
So I see this as an opportunity to work on a flexible watershed summary that includes usgs error analysis. Something that is usable in both wsp and permit analysis, something that benefits the public understanding but also can benefit us modelers, and also be used for model validation (it should be also able to compare model to model).
This will ultimately form the "final report deliverable" for the harp project.
Retrieving JSON objects with all model results in them
obj_url <- paste(json_obj_url, model$pid, sep="/") model_info <- ds$auth_read(obj_url, "text/json", "") model_info <- fromJSON(model_info ) elid <- model_info [[1]]$om_element_connection$value
3 models Including OM, hspf and hsp2
Resulting Document Layout: