Estimating Varying AGWRCs with Outflow/AGWS

[ilonah22]

On 2/10 we discussed using the formulas from these regression lines to estimate outflows with varying input AGWRCs.

• Here are the coefficients and intercepts from the plots above
• Here are the landseg weight for each gage
• Here are the model parameters
• Here are the area weighted model AGWRCs

The scripts that I am using for calculations are:

• flow_calcs.R
• landseg_weight_calcs.R

Steps/Options for Simulating varying AGWRC: @ilonah22

1. Use the line that we plot through our chart as a proposed baseflow function
   • 1a Compare Model regression line to HSPF AGWRC coefficient (weighted param value)
   • If our method of identifying events and capturing AGWRCs is accurate, then the values should look similar (numerically)
   • Link to CSV with model parameters used in 1a calcs:
   • On the server these params are in the file: /media/model/p6/input/param/pas/P620171001WQf/PWATER.csv
2. Quantile regression of baseflow events for better AGWRC estimate
   • Do this if we conclude than using the median line is producing a biased result (too high or too low AGWRC)
   • First check values quantitatively to determine what is "close enough"
3. Comparison of AGWRC vs USGS Flow to provide model parameters
   • After we have determined what the relationship between flow and AGWRC should be, translate into relationship between storage and AGWRC
4. R-based model using $AGWRC=f(AGWS)$ and HSPF equations #51
   • Eventually be able to use model to project into the future given current flow
5. Modify HSPF (special actions)#### Extracting formula for $C_{AGWRC} = f(AGWS)$ - Note: ultimately we will be simulating $C_{AGWR} = f(S_{AGW})$, and since we understand the relationship between AGWRC and Qout, this should be trivial -- but let us not forget this.

The first thing I wanted to do was set up the correct formula using these coefficients. In this example, from Cootes Store gage data.

lm(formula = event_AGWRC ~ log(median_flow))

$eventAGWRC = -0.0003047*log(Q) +0.9445$

$\frac{eventAGWRC - 0.9445}{-0.0003047} = log(Q)$

$Q = 10^{\frac{eventAGWRC - 0.9445}{-0.0003047}}$

[ilonah22]

Initial Results

Using the formula $Q = 10^{\frac{eventAGWRC - intercept}{slope}}$, I calculated some flows for Cootes Store and got some interesting results:

Obviously, some very large numbers for flow calculated. I went back and checked my math from above, and everything seems to be in order, but I'd appreciate if someone could check it too. @willprokopik, I guess it is also possible that I am misinterpreting the results of the regression, so maybe you could give it a try and let me know if you get a different formula.

[rburghol]

Hey @ilonah22 -- I have not checked your math here, though I'm assuming that it works out. But my question revolves around what question you're trying to answer here -- I may totally be missing something so correct me if so..

We're supposed to be moving towards a model where AGWRC varies at different AGWS levels, which we are estimating in this relationship as varying according to baseflow (note this is important, as storm-flows lead to spurious results), and I think you've solved for flow as a function of AGWRC?

Let's step back abd answer:

• Are you working on bullet 1 "Use the line that we plot through our chart as a proposed baseflow function -- Compare to HSPF AGWRC coefficient" from above.
• How does this approach fit into the goal of equipping our model with a varying AGWRC?
• We discussed doing an R simulation (bullet 4 from above), but we need to do the intermediate steps first I think. Also note that for that, we will need to predict flow, but we will do so from AGWS not Q.

[ilonah22]

@rburghol, I think I misunderstood what our goal was. In my notes from last weeks meeting, I had written down that we would be using our event calculated AGWRCs to estimate flow as a means of seeing whether our understanding of the model worked when using our newly calculated and varying AGWRCs. Reading back on the description now I see that you had intended working in the other direction and using flow to get AGWRC.

I guess I don't really understand how using flow to calculate a varying AGWRC is any different than what we have already done with the lm(eventAGWRC ~ log(medianFlow) or our event identification.

Also, just to answer your specific questions:

• I was working on bullet 1, using the coefficients and intercepts from Will's app to plug into the equation above.
• I was originally thinking that out plan was to input varying AGWRCs into the model or an R version of the model, and try and use that to calculate some of our key variables (outflow, storage, ET?). But, I now realize this isn't the plan, so I think you may have to explain to me again what our goal was with this step.

[rburghol]

Ultimately, we do wish to simulate outflow. But the intermediate steps are crucial, even though they appear simple on their surface.

So status:

• you completed the first bullet, but you should verify that you're getting reasonable values by plugging in the range of flows and plotting it.
• note: there is a sub bullet on the first task that has not yet been complete completed.
• For the sub bullet, you need to harvest the AGWRC values and do a comparison. I am pretty sure that we shared the CSV file with all of the AGWRC values for all landsegs early on in our project, so start by looking back through those earlier issues. If it leaves out out of me I will pop it back up here.

These bullets then found the foundation for the more sophisticated prototyping of the model in our subsequent bullets.

[ilonah22]

Forgot to put these up yesterday, but here are some plots to show the newly predicted AGWRCs calculated using outflow:
Cootes Store:

Mount Jackson:

Strasburg:

It sort of seems that the slopes get more extreme (both positive and negative) as you move through the watershed. I think today it would be good to decide whether I should move onto quantile analysis or trying to get area weighted values using multiple landsegs.

Use table with number of cells for area weighted values: HARPgroup/HARParchive#1427 (comment)

[rburghol]

@ilonah22 Question to guide you to answer your questions:

1. "whether I should move onto quantile analysis or trying to get area weighted values using multiple landseg" What are the benefits of skipping step 1a (compare to model AGWRC -- aka the weighted param value) in the issue task list?
2. How could the answer to 1a inform whether we need to do the quantile analysis?
3. Also check this comment from last week Estimating Varying AGWRCs with Outflow/AGWS #47 (comment)

[rburghol]

@ilonah22 - note, the goal IS to estimate "Outflow with Varying AWRCs", that 's what the final bullet in the task list will accomplish :) . The first step is indeed "Estimating Varying AGWRCs with Outflow/AGWS"

[rburghol]

@willprokopik @ilonah22 @benhdye -- if you are interested in learning a powerful system for doing spatial analysis and or just have something to double-check your R-GIS analysis, here is the area-weighted components (Table 1) as derived by an SQL query on our database (Code 1), and then finally an example of how to use the same query from Rstudio (Code 2). Notes:

• By my eye, the totals for this spatial query exceed 1.0 by a little bit -- and the sum total of overlap percents should always equal 1.0. This is most likely due to overlapping shapes, or some irregularity in the shapes that causes a double counting.
• So, to manage this would require normalizing each of these % overlaps by the sum of the % overlaps. See Code 3 for my attempt on that, which also introduces a bit more complicated SQL structures WITH, which is very useful IMO.
• The normalized pct_overlaps are presented in Table 2.

Code 1: Find the overlapping area percentage for each land segment over the full drainage shape of a USGS gage (the area should correspond well to the modeled river segments that we use in Shenandoah).

select a.hydrocode, b.hydrocode, 
  st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom))
    /st_area2d(a.dh_geofield_geom) 
from dh_feature_fielded as a 
left outer join dh_feature_fielded as b 
on (
  a.dh_geofield_geom && b.dh_geofield_geom 
  and b.ftype='cbp6_landseg'
) 
where a.hydrocode = 'usgs_ws_01634000' 
  and a.ftype='usgs_full_drainage' 
  and st_area2d(b.dh_geofield_geom) > 0 
  and st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom)) > 0;

Table 1: Results from Query 1.

hydrocode	hydrocode	% overlaps
usgs_ws_01634000	H51165	0.05728021177690966
usgs_ws_01634000	N51139	3.895692612701101e-05
usgs_ws_01634000	N51165	0.48421801110997803
usgs_ws_01634000	N51171	0.4816838283013124
usgs_ws_01634000	N51187	0.0001204666330892836
usgs_ws_01634000	N51660	0.0015151872697934045
usgs_ws_01634000	N54031	0.021819394748784433
usgs_ws_01634000	N54071	7.713128555783912e-05

Code 2: How to run this same query from Rstudio (must be connected to VT VPN or be on campus).

library("hydrotools")
library("sqldf")
basepath="/var/www/R"
source("/var/www/R/config.R")
weighted_landsegs <- sqldf(
  "select a.hydrocode, b.hydrocode, 
     st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom))
     / st_area2d(a.dh_geofield_geom) 
   from dh_feature_fielded as a 
   left outer join dh_feature_fielded as b 
   on (
     a.dh_geofield_geom && b.dh_geofield_geom 
     and b.ftype='cbp6_landseg'
   ) 
   where a.hydrocode = 'usgs_ws_01634000' 
   and a.ftype='usgs_full_drainage' 
   and st_area2d(b.dh_geofield_geom) > 0 
   and st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom)) > 0
  ",
  connection=ds$connection
)

Code 3: Find the overlapping area percentage for each land segment over the full drainage shape of a USGS gage, but normalize on the sum of overlap fractions in order to avoid double counting from shape weirdness.

WITH g AS (
  select a.hydrocode as gage, b.hydrocode as landseg, 
    st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom))
      /st_area2d(a.dh_geofield_geom) as pct_overlaps
  from dh_feature_fielded as a 
  left outer join dh_feature_fielded as b 
  on (
    a.dh_geofield_geom && b.dh_geofield_geom 
    and b.ftype='cbp6_landseg'
  ) 
  where a.hydrocode = 'usgs_ws_01634000' 
    and a.ftype='usgs_full_drainage' 
    and st_area2d(b.dh_geofield_geom) > 0 
    and st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom)) > 0
),
gsum as ( 
  select sum(g.pct_overlaps) as total 
  from g
) 
select g.gage, g.landseg, (g.pct_overlaps / gsum.total) as pct_overlaps
from g, gsum;

Table 2: Results from Code 2, with normalized percents overlapped to sum to 1.0. (Note: this markdown rendered table can be done as kableExtra::kable(weighted_landsegs,format="markdown")

gage	landseg	pct_overlaps
usgs_ws_01634000	H51165	0.05472179347600768
usgs_ws_01634000	N51139	3.7216916625328114e-05
usgs_ws_01634000	N51165	0.4625904335780543
usgs_ws_01634000	N51171	0.4601694399402104
usgs_ws_01634000	N51187	0.00011508599588172518
usgs_ws_01634000	N51660	0.0014475114927653629
usgs_ws_01634000	N54031	0.020844832380591555
usgs_ws_01634000	N54071	7.36862198637416e-05

[rburghol]

@ilonah22 - I shared these results with my colleague rom the bay program today and we were both super excited by the prospects, and also lookoing forward to the weighted AGWRC from the model files. Also, a reminder to:

• link the CSV that you uploaded that has all of those model parameters in the main body of this issue
• Link the R code that you are using to do your intermediate calculation in the main body of this issue

Thanks!

[ilonah22]

@rburghol - The code you sent above, for getting the percent overlaps was super helpful. I was able to get the area weighted AGWRCs using your code from above and the model parameters. I have linked those, the model parameters, and the scripts I am using in the main body of the issue as well.

These weighted AGWRCs are slightly lower than I would have expected. I know we have been talking about 0.97 but these are all <0.96.

This is an example of the sqldf that I used for each gage, just combined code from your comment and changed gage number in hydrocode, all included in landseg_weight_calcs.R:

CS_weight <- sqldf(
  "WITH g AS (
    select a.hydrocode as gage, b.hydrocode as landseg, 
      st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom))
        / st_area2d(a.dh_geofield_geom) as pct_overlaps
    from dh_feature_fielded as a 
    left outer join dh_feature_fielded as b 
    on (
      a.dh_geofield_geom && b.dh_geofield_geom 
      and b.ftype='cbp6_landseg'
    ) 
    where a.hydrocode = 'usgs_ws_01632000' 
      and a.ftype='usgs_full_drainage' 
      and st_area2d(b.dh_geofield_geom) > 0 
      and st_area2d(st_intersection(b.dh_geofield_geom, a.dh_geofield_geom)) > 0
  ),
  gsum as ( 
    select sum(g.pct_overlaps) as total 
    from g
  ) 
  select g.gage, g.landseg, (g.pct_overlaps / gsum.total) as pct_overlaps
  from g, gsum;
  ",
  connection = ds$connection
)

Since we have the area weighted AGWRCs now, I wanted to visually compare them to the AGWRCs that we calculated using the regression coefficients. Also, as you will be able to see from the new plots, it looks like before I used the same event flow data for all graphs, I think this was a problem with using csvs from the wrong branch, but @benhdye is currently pushing and merging the updated data to the main branch, so this should be cleared up for the future.

As a reminder, here are the coefficients gathered from Will's regressions:

Cootes Store:

Mount Jackson:

Strasburg:

I wasn't sure which observed event data you meant, but I made this boxplot to show the AGWRCs at different points in our methodology for Strasburg. The first image is at full scale, the second, zoomed in to (0.85,1.05).
Key for x axis:

• 'Calculated from observed Flow' - $Q_t/Q_{t-1}$
• 'Calculated using Model Coefs' - Using the regression slope and intercept from analyzing model AGWRCs against median event flow
• 'Calculated using Gage Coefs' - Using the regression slope and intercept from analyzing gage AGWRCs against median event flow
• 'Event Value after Trimming' - Event AGWRC determined after being run through Ben's trimming function