search

HARPgroup / HARParchive

This repo houses HARP code development items, resources, and intermediate work products.
1 stars 0 forks source link

FTABLEs for streams #331

Closed rburghol closed 1 year ago

rburghol commented 2 years ago

Objectives

I have updated our objectives, and split them out into bite-size sub-issues - @megpritch @durelles @jdkleiner please let me know your thoughts.

Overview

FTABLE Files

Data Model

test it! works like a champ

jsonlite::unserializeJSON(ftable_prop$proptext)

nicoledarling commented 2 years ago

Update:

Here are some questions Megan and I have after reviewing all the prep task documents:

Here are some key points we found from the readings:

rburghol commented 2 years ago

Just saw this, thanks for the email nudge:

Another note on the variables from VA Hydro: while you'll only need drainage area, province, and channel length to do the analysis, the other attributes that come from the regression model such as bank full depth, slope etc. Will be part of the validation that you perform to see if your procedure has the same outputs as that which we've been previously using in VAHydro. Also, these are to the best of my recollection, but it's been a long time since I did this analysis myself so if you see something that suggests that we need other data please let me know.

durelles commented 2 years ago

Hi all, When you do the deep dive into the hspf docs for F-tables, please send a link/summary to review. There may possibly be other props of interest. S

Sent from my iPhone

On Jul 25, 2022, at 8:24 AM, rburghol @.***> wrote:

 Just saw this, thanks for the email nudge:

Check out the local_channel object in the Sugar Hollow model for existing VA Hydro stream properties. Of all the variables that are in Local_channel currently you will need the drainage area, physiographic province, and channel length. The drainage area and physiographic province are needed for the Usgs regression model, but the channel length will be needed in order to compute the volume of the segment at different flow depth. Oh, this is where we need you all to dive deep into the HSPF manual and educate all of us on what those extra Ftable columns do. In my own experience, I've never used them. Another note on the variables from VA Hydro: while you'll only need drainage area, province, and channel length to do the analysis, the other attributes that come from the regression model such as bank full depth, slope etc. Will be part of the validation that you perform to see if your procedure has the same outputs as that which we've been previously using in VAHydro. Also, these are to the best of my recollection, but it's been a long time since I did this analysis myself so if you see something that suggests that we need other data please let me know.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.

megpritch commented 2 years ago

HSPF Manual Summary:

  1. Where in VAHydro do you want the FTABLE to go under? Just within the local_channel object?
  2. Are you able to clarify what the variables in the php code stand for? (hc, he, bfc, bfe, bc, be, z)

  3. And why are these values constant rather than by user-input?

    • We know from the USGS Publication that these equations correspond to lines 1112-1115 in the php code: image image
    • Therefore we know what h, bf, and b stand for, but still not hc, he, bfc, bfe, bc, be, z
rburghol commented 2 years ago

I used the equations in table 7 which had a different one for each physiographic province. Therefore at the beginning of our PHP code (lines 1054-1110) we set coefficients for the equation that were specific to the province. Those equations were in the form:

$C * DA^E$

So in the first lines if the PHP code I set the Coefficient and Exponent corresponding to the selected region, so:

Similarly, he, bfe, and be were the exponents for those respective equations.

This is not to suggest I used the right method, or the method we should use to match the current version of the CBP model, just that's the one I chose due to my objectives for the old model. That old model used a Muskingum method stream flow modeling approach that used iterative solving. Since we will want to validate our work against the FTABLEs in the CBP suite, we definitely need to determine which method they used. And if it used the single equation instead of the region specific, I definitely would like to know what the rationale was (assuming you all can find it in the documentation).

nicoledarling commented 2 years ago

Here are some questions Megan and I have after working on the script. If you are available tomorrow, can we set up a meeting to make sure we are on the right path.

Questions:

  1. Will the values we pull from Vahydro (drainage area, province, length) always be under “local_channel”? We found a river segment where it was under “river_channel” and are wondering if they will be under different names.
  2. What are the units of the variables he, bfe, and be? We are also wondering how to look at the units in Vahydro to make sure we can convert any units if needed.
  3. Megan and I think we need to pull the slope from Vahydro in order to calculate the discharge. Do you know if the slope in Sugar Hollow is a longitudinal or water surface slope?
  4. We made a rough draft of an Ftable in our script and wondering where in Vahydro the already made Ftables are that can confirm if we are doing our calculations right.
  5. Just to clarify, is the z variable the side slope of the channel?
megpritch commented 2 years ago

Here are the ftable values we got so far and are looking to check:

image

rburghol commented 2 years ago

Reponse to questions:

nicoledarling commented 2 years ago

Update:

Megan and I tried many different things to correct the output of the Ftable:

megpritch commented 2 years ago

Using the depths from the uci, this is the ftable we calculate:

image

rburghol commented 2 years ago

OK - @megpritch @nicoledarling solving for DA is very smart. You got 106 sqmi, I just checked the VAHydro model for OR1_7700_7980 and the area in there is 97.75 sqmi, so, this is excellent work on your part! But, it looks like VAHydro has the correct area, so, where did the 2,796.95 value come from? Sounds like we have a data retrieval error somewhere -- maybe check our REST chain of data retrievals?

VAHydro River Channel object: http://deq1.bse.vt.edu:81/d.dh/om-model-info/4805858

megpritch commented 2 years ago

@rburghol @nicoledarling

I Have Uncovered a Sequence of Errors

In Summary:

image image

The River Channel Object I see: http://deq1.bse.vt.edu:81/d.dh/om-model-info/4805894

My Train of Thought:

[1] My search results for OR1_7700_7980: image

[2] Initial R code for pulling the model: image

[3] R-Generated Environment from pulling river segment feature: image

[4] The Different Spellings on VAHydro: image image

Note: if "Cub" is the typo and "Club" is correct, then I noticed "Lower Cub Creek" also has this typo in its feature name.

rburghol commented 2 years ago

Thanks for the details, couple of things:

megpritch commented 2 years ago

Thank you for the clarification about VAHydro. Exploring the website had been what helped me to understand the code better and know what to put for the different attributes of RomFeature$new.

image

rburghol commented 2 years ago

Thanks for the detail. My observations and questions fwiw:

megpritch commented 2 years ago

image

rburghol commented 2 years ago

Definitely look further -- this is exactly the kind of nuance we are here to understand and sort out.

Great work -- we should probably spend a few minutes going over these graphs tomorrow, when we get in let's look at our schedules and pop-up a zoom.

megpritch commented 2 years ago

On the left [specific n= 0.095 & generic n~0.05] and on the right [specific n= 0.04 and generic n = 0.037] image

nicoledarling commented 2 years ago

Notes from meeting with Dr. Scott

We had a meeting with Dr. Scott today discussing the floodplain geometry. He gave us many resources to start researching the different methods and equations.

3 methods to calculate floodplain geometry:

Megan and I started to research the HAND method to see how it worked and what GIS mapping was involved. Source: https://doi.org/10.1111/1752-1688.12661 Here is what we found:

Other sources we are looking at: https://nhess.copernicus.org/articles/19/2405/2019/ https://cfim.ornl.gov/data/ https://water.noaa.gov

nicoledarling commented 2 years ago

Changing floodplain side slope

Here is some documentation for changing the floodplain side slope based on visual inspection. I spent some time today trying different factors to multiply the slope by to match the uci.

For the Rockfish segment, the floodplain slope is about 6.10-6.11 times as big as the channel slope. Here is the Ftable and uci to compare.

Our Ftable: image Uci Ftable: image

The numbers are very close for area and only differ up to 5 acres. Volume they are a little more further apart and discharge is very far off. Here are the graphs: image

I tried to see if any of the other variables in our code could somehow produce the 6.1 factor by multiplying and dividing different numbers but had no luck. I am currently going back into our many resources and rereading about floodplain slope to see if there's any information helpful to this.

megpritch commented 2 years ago

Documenting Recent Work

R Estimation Method:

HAND with GIS Method:

2022-08-17 (2)

rburghol commented 2 years ago

Thanks for the update on this. I think we definitely want to avoid doing the whole DEM analysis ourselves, if we can possibly avoid it. IT will just be hard to automate, and since the bay program has a team of 1-2 people working on this full time, I think our efforts can be better spent on seeing if we can find high-efficiency methods, like finding them in NHD as you mentioned. To that:

megpritch commented 2 years ago

Yes, Joey showed us how he pulled values from NHD into R using the nhdplusTools package, and talked a little bit about using the same getNHDplus function to potentially search for values within a certain HUC. He showed us an example where he pulled an NHD feature from a coordinate, and there were a lot of different variables attached to that point. We can definitely look into what's there and what the values stand for. If we could find a floodplain side slope, that would be ideal.

rburghol commented 2 years ago

Awesome - thanks for the detail. I am sure that we can find some value to use for an estimate, and the question becomes is it good or good enough. I fully anticipate the floodplains FTABLEs to be replaced by something moving forward so good enough will probably be the bar we set :). In general:

jdkleiner commented 2 years ago

Just tuning into this thread - Agreed with Rob, you "should" be able to take the tools available in the nhdplusTools package to replace your manual ArcGIS workflow with a more automated method. I think there should be enough useful data there to get us moving forward.

You all can look further into National Water Model tools in R, a quick search yields this resource (though it doesn't look to be actively maintained) https://github.com/mikejohnson51/NWM

nicoledarling commented 2 years ago

Using get_elev_along_path function in R

Joey recommended we try using this function to get different elevation values along the Rockfish River. I included 10 points in the code (the example from nhdplus has 3) to make the graphs a little more accurate.

Graph 1: image

Graph 2: image

Here's an example of the code with the first 3 points:

points <-sf::st_as_sf(c(point1, point2, point3, point4, point5, point6, point7, point8))

(Table1 <-get_elev_along_path(points, 100))

if(!is.null(Table1)) { bbox <-sf::st_bbox(Table1) + c(-0.005, -0.005, 0.005, 0.005)

nhdplusTools::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_geometry(Table1), 3857), pch = ".", add = TRUE, col = "red") plot(sf::st_transform(sf::st_sfc(point1, crs = 4326), 3857), add = TRUE) plot(sf::st_transform(sf::st_sfc(point2, crs = 4326), 3857), add = TRUE) plot(sf::st_transform(sf::st_sfc(point3, crs = 4326), 3857), add = TRUE) }

jdkleiner commented 2 years ago

@nicoledarling this is very interesting indeed! One avenue for exploration, rather than looking at a line down a stretch of river as you've done above, you could try looking at a few river transects or cross sections. I think this will help get us closer to describing that relation between channel width/depth/area/slope.

get_elev example

megpritch commented 2 years ago

Documentation of Assumptions for FTABLE Floodplain Calculations:

Visual representation of channel cross-section [adapted from BASINS Technical Note 1]: image Variables:

Main Assumptions & Clarifications:

The Code for Reference:

#----Calculating FTABLE----
#depth
cdepth <- seq(0,h,length=10) #channel
fdepth <- seq(h+1, h*4 ,length=9) #floodplain

Abf <- ((sw+b)/2)*h #channel cross-sect area @ bankfull
Pbf <- b + 2*h*sqrt(z**2 +1) #channel wetted perimeter @ bankfull

#--Discharge Calculation Info:
# Q = V * A ; where A = cross sectional area (aka flow area)
# Manning's Eqn: V = (1.49/n) * R^(2/3) * S^(1/2) 
# (1.49/n) is English ; (1/n) is metric
# Hydraulic Radius = A/P , cross-sect.area/wetted perim.

fn_make_trap_ftable <- function(depth, clength, cslope, b, z, n, h, bf, fp, Abf, Pbf) { 
  sw <- b + 2*z*depth #surface width
  sw[depth == 0] <- 0
  A <- ((sw+b)/2)*depth #cross-sect. area (trapezoid)
  P <- b + 2*depth*sqrt(z**2 +1) #wetted perimeter (trapezoid)

  if (fp == TRUE) {
    A <- A + Abf #combine fp and channel@bankfull trapezoidal areas
    P <- P + Pbf - bf #fp + wetted perim.@bankfull - bf width of channel
    depth <- depth + h
  }

  area <- (sw * clength)/43560 #surface area
  vol <- (clength * A)/43560 # /43560 to convert to ft-acres
  disch <- (1.49/n) * (A/P)**(2/3) * cslope**0.5 * A
  ftable <- data.frame(depth, area, vol, disch)
  return(ftable)
}

cftab <- fn_make_trap_ftable(cdepth, clength, cslope, b, z, n, h, bf, FALSE, Abf, Pbf) #in-channel
fptab <- fn_make_trap_ftable(fdepth-h, clength, cslope, 5*bf, z, nf, h, bf, TRUE, Abf, Pbf) #floodplain 
# ^base of floodplain = 5x bankfull width

ftable <- rbind(cftab, fptab)
durelles commented 2 years ago

Great documentation.The one thing I think I would add is both non-dimensional units (e.g., L2) and the units used in our specific model. The figure really helps illustrate the concept. You could also consider adding a statement related to how theres a drastic change in area once you get to line 10 in the table as expected.CheersDSSent from my iPadOn Oct 29, 2022, at 5:22 PM, megpritch @.***> wrote: Documentation of Assumptions for FTABLE Floodplain Calculations: Visual representation of channel cross-section [adapted from BASINS Technical Note 1]:

Variables:

depth = depth of the water sw = surface width of the water at certain depth channel: h=bankfull depth, b=base width, bf=bankfull width, z=side slope, cslope=longitudinal slope along channel length floodplain: 5*bf=base width, zfp=side slope=same as channel

Main Assumptions & Clarifications:

The channel and the floodplain were treated as two trapezoids. When choosing the depths for the 19 rows of the table: rows 1-10 account for the channel and go from zero to h in even increments; rows 11-19 are into the floodplain and go from (h+1) to (4*h) in even increments. The base of the floodplain trapezoid is 5x the channel bankfull width, and the side slope of the floodplain was assumed to be equal to that of the channel. Cross-sectional area of the water beyond bankfull is calculated as the sum of (1) the bankfull channel trapezoid area and (2) the area of water in the floodplain trapezoid at a specific depth past bankfull. Wetted perimeter into the floodplain is calculated as the sum of (1) bankfull wetted perimeter and (2) floodplain trapezoid wetted perimeter at a depth past bankfull, minus the bankfull width.

The Code for Reference:

----Calculating FTABLE----

depth

cdepth <- seq(0,h,length=10) #channel fdepth <- seq(h+1, h*4 ,length=9) #floodplain

Abf <- ((sw+b)/2)h #channel cross-sect area @ bankfull Pbf <- b + 2h*sqrt(z**2 +1) #channel wetted perimeter @ bankfull

--Discharge Calculation Info:

Q = V * A ; where A = cross sectional area (aka flow area)

Manning's Eqn: V = (1.49/n) R^(2/3) S^(1/2)

(1.49/n) is English ; (1/n) is metric

Hydraulic Radius = A/P , cross-sect.area/wetted perim.

fn_make_trap_ftable <- function(depth, clength, cslope, b, z, n, h, bf, fp, Abf, Pbf) { sw <- b + 2zdepth #surface width sw[depth == 0] <- 0 A <- ((sw+b)/2)depth #cross-sect. area (trapezoid) P <- b + 2depth*sqrt(z**2 +1) #wetted perimeter (trapezoid)

if (fp == TRUE) { A <- A + Abf #combine fp and @. trapezoidal areas P <- P + Pbf - bf #fp + wetted @. - bf width of channel depth <- depth + h }

area <- (sw clength)/43560 #surface area vol <- (clength A)/43560 # /43560 to convert to ft-acres disch <- (1.49/n) * (A/P)*(2/3) cslope*0.5 A ftable <- data.frame(depth, area, vol, disch) return(ftable) }

cftab <- fn_make_trap_ftable(cdepth, clength, cslope, b, z, n, h, bf, FALSE, Abf, Pbf) #in-channel fptab <- fn_make_trap_ftable(fdepth-h, clength, cslope, 5*bf, z, nf, h, bf, TRUE, Abf, Pbf) #floodplain

^base of floodplain = 5x bankfull width

ftable <- rbind(cftab, fptab)

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>