Gradient and Sinuosity Calculations

[philipbaileynar]

We need gradient and sinuosity calculations for reach typing.

The goal of this ticket is to produce standalone Python modules for calculating gradient and sinuosity. (Other calculations will follow soon...) Essentially this code takes a polyline feature class and adds the gradient and sinuosity attributes (dropping them first if they are already present) and populates them with the new values.

We already have gradient as part of the BRAT reach geometry module. So I recommend that this is simply moved up to rscommons and called from there. Then refactor to BRAT to call it from rscommons as well.

Sinuosity does not use an underlying raster, but it can essentially be written as a new standalone module that mimics some of the pattersn from the reach geometry code.

Important Notes

1. How did GNAT calculate sinuosity and gradient in the past? Is it different?
2. What scale does it make sense to calculate these attributes? We will get different values for very short reaches as opposed to longer lines.
3. The gradient code linked above samples the elevations as the end of reach polylines. It does use a buffer (30m I think?) but it samples the raw DEM. USU in the past has used an averaging filter on the DEM first to avoid anomalous values.
4. Sinuosity and gradient are the first of many calculations that we will need for reach typing. Keep in mind that we will be extending this code. Code these current features as standalone code but think about how we will combine methods into a coherent GNAT tool.

[KellyMWhitehead]

arcGNAT previously generated:

• Planform Sinuosity (Channel Segment length/distance of channel segment endpoints). This is easy and only needs the line network as an input.
• Channel Sinuosity (Channel Segment length/ valley bottom Centerline length). We will need vb centerline and a method to spatially join the two (i.e. probably use vb sausages for this).
• Valley Bottom Sinuosity (Valley Bottom Centerline Segment length/distance of Valley Bottom Centerline Segment endpoints)
• Gradient: dem was sampled directly at the endpoints, no averaging.

These tools can be run at any scale, but tend to work better for longer reaches (less 'noise' and fewer geometry special cases). We never quite figured out the idea segmentation to use during the arcGNAT era.

Shopping list for GNAT inputs (as we develop the GNAT rs project)

• Flowline Network (Segmented?)
• DEM (or use py3dep?)
• Valley Bottom (with sausages?)
• Valley Bottom Centerline

[joewheaton]

@lauren-herbine you should chime in on this one...

Gradient

We should differentiate

• Segment Gradient
• Valley Bottom Gradient

We should also look into what we did in BRAT a long time ago to smooth and get reasonable gradients. I vaguely recall:

• Simple Search Window in vicinity of top of reach elevation point and bottom of reach elevation point to find nearest minimum and another check to make sure upstream was higher than downstream (if not force no slope). I think we also looked at upstream and downstream slopes for some sort of assurances.

Sinuosity

A few things:

For ideas on symbolization...

And, way back in 2016 Kirstie and I drew this on my whiteboard when we were trying to clarify these things.

Primary Planform Measures

Out of it comes the idea that we need direct measures of:

• Valley StraightLine Distance V SL
• Valley Bottom CenterLine Distance V CL
• Valley Bottom StraightLine Distance VB SL
• Valley Bottom CenterLine Distance VB CL
• Channel StraightLine Distance Ch SL
• Channel CenterLine Distance Ch CL
• Channel Thalweg Distance - Ch Th

A big challenge is whether we calculate channel distances as total of all anabranches in multithreaded system or mainstems?

Primary Sinuosity Metrics

To calculate primary sinuosity metrics:

• Valley Sinuosity V CL/ VSL
• Valley Bottom Sinuosity - VB CL/ VB SL
• Planform Sinuosity - Ch CL/ Ch SL
• Channel Sinuosity - Ch CL/ VB CL
• Thalweg Sinuosity - Ch Th/ Ch CL

Note similarity to @KellyMWhitehead :

• Planform Sinuosity (Channel Segment length/distance of channel segment endpoints). This is easy and only needs the line network as an input.
• Channel Sinuosity (Channel Segment length/ valley bottom Centerline length). We will need vb centerline and a method to spatially join the two (i.e. probably use vb sausages for this).
• Valley Bottom Sinuosity (Valley Bottom Centerline Segment length/distance of Valley Bottom Centerline Segment endpoints)

Why do we need all this?

Well the various sinuosity measures are useful by themselves for discriminating different conditions. However, it is the intercomparing of some of these sinuosity values at the riverscape setting that helps us discriminate different settings really well.

For example:

• A high channel sinuosity and low valley bottom sinuosity suggest
• If we have a high planform sinuosity, is it explained or inherited by a high valley bottom sinuosity?
• If we have a high valley bottom sinuosity, and low valley sinuosity, that suggest presence of lots of fans pushing valley bottom path around.
• If valley bottom sinuosity and valley sinuosity are basically equal, it suggests lack of fans. If true, but both values are moderate to high, it suggests structural geologic controls.
• If channel sinuosity is moderate to high, but VB sinuosity is low, it suggests active meandering (i.e. planform controlled or laterally unconfined).
• If planform sinuosity is high, but channel sinuosity is low, it suggests inheritance of sinuosity from valley bottom sinuosity.
• If thalweg sinuosity is low, likely very planar morphology in channel units.

[philipbaileynar]

What length scale are wee going to measure reach sinuosity? All rivers are straight if measured at some silly minute scale.

We are already segmenting NHD+ HR at 300m, land ownership, roads and rail. There's also talk of segmenting at ecoregions (#236 ).

The length over which we measure sinuosity affects the result. How do we pick the right length if our reaches are so fragmented that they are nearly all straight?

[lauren-herbine]

Sinuosity

Inputs:

From @joewheaton :
Out of it comes the idea that we need direct measures of:

• Valley StraightLine Distance V SL
• Valley Bottom CenterLine Distance V CL
• Valley Bottom StraightLine Distance VB SL
• Valley Bottom CenterLine Distance VB CL
• Channel StraightLine Distance Ch SL
• Channel CenterLine Distance Ch CL
• Channel Thalweg Distance - Ch Th

Metrics:

from @joewheaton
Valley Sinuosity V CL/ VSL
Valley Bottom Sinuosity - VB CL/ VB SL
Planform Sinuosity - Ch CL/ Ch SL
Channel Sinuosity - Ch CL/ VB CL
Thalweg Sinuosity - Ch Th/ Ch CL

Scale

@philipbaileynar let's start segmenting at 500m. I think level IV ecoregions would be too long in some cases and could "pinch" a reach in others:

an example from CO

Multithreaded systems

... coming soon...

Gradient

Question for you all (@philipbaileynar @KellyMWhitehead @joewheaton)- the NHD+ has slope as an attribute. Could this be used?
Ex from something I've already put together:

[lauren-herbine]

an update for moving forward:

Let's calculate those things that are easy:
(from @KellyMWhitehead)

• Planform Sinuosity (Channel Segment length/distance of channel segment endpoints). This is easy and only needs the line network as an input.

• Channel Sinuosity (Channel Segment length/ valley bottom Centerline length). We will need vb centerline and a method to spatially join the two (i.e. probably use vb sausages for this).

• Valley Bottom Sinuosity (Valley Bottom Centerline Segment length/distance of Valley Bottom Centerline Segment endpoints)

• Segment gradient

• Valley bottom gradient

For this first run I am interested in seeing how well channel sinuosity discriminates between river styles. This is what we were using in the MHFD, and for that "small" area it worked well. For multi-threaded systems, we will use the mainstem (if it is anastomosing/the NHD shows multiple channels). For braided systems, we will use the length of the channel belt. See figure below:

I am still thinking segmenting at 500 m will be a good place to start.

[joewheaton]

Slope

Please see how we calculate and display symbology for reach slope in BRAT. The slope in NHD+HR could be good, but it will be based on its segmentation, which may suffice.

We also need to think about valley bottom slope and this would be valley bottom centerline length divided by the difference in "nearest min elevation in a window (e.g. 10 x 10) at top and bottom.

Sinuosity Measures

I like what you have proposed for sinuosity measures. I would suggest NOT presupposing the sinuosity segmentation of channel netowork to 500 m to start. Start with existing NHD +HR segmentation and see how bad or different it is. I think for bigger rivers we may want to have longer segmentation. @KellyMWhitehead please lay out as you start doing this exactly what the choices are in making these calculations and or give us some different examples so we can make some informed decisions. The thing I don't think we're considering here yet is the discrepancy between NHD+HR channel segmentation (default) or our additional segmentaton, and our VBET cl segmentation. Any insights @KellyMWhitehead ?

Please also see @lauren-herbine what @kbartelt came up with in RIM for her Relative Flow Length metric, which was very good for multi-threaded systems (i.e. total channel(s) CL or thalweg length divided by riverscape length).

[joewheaton]

Thinking about @philipbaileynar :

@philipbaileynar let's start segmenting at 500m. I think level IV ecoregions would be too long in some cases and could "pinch" a reach in others:

and

The length over which we measure sinuosity affects the result. How do we pick the right length if our reaches are so fragmented that they are nearly all straight?

So another approach is with @KellyMWhitehead and Carol stole from the Fluvial Corridor Tool, in which some metrics can be calculated at length scales other than their segmentation with a moving window (also facilitates calculating at multiple scales). In other words, it doesn't matter what your segmentation length scale is, you choose a "metric length" calculation scale appropriate for the calculation (e.g. a fixed length window of VB for example). This is the approach laid out in Notebaert & Piegay (2013) of using DGOs and AGOs...

• Notebaert B, Piegay H. 2013. Multi-scale factors controlling the pattern of floodplain width at a network scale: The case of the Rhone basin, France. Geomorphology 200 : 155–171. DOI: 10.1016/j.geomorph.2013.03.014

[lauren-herbine]

Reigniting. We have a lot of good discussion on here, I will attempt to narrow everything down into what we actually need, starting with Channel Gradient

Stream Gradient

Attribute	Tool?	Line(s) of evidence needed	Measurements needed from LOE	Calculation	Alpha (ready today)	Beta v 0.1	V 1 (i.e. thesis)	V2 (i.e. later)
Stream Gradient	ChannelArea + VBET	DEM + ChannelArea polygons + Channel centerline OR Stream Network	ELEV_START - elevation (in native units of input raster dataset) of ‘From’ node, ELEV_END - elevation of ‘To’ node, LENGTH_M - length in meters per stream feature	((𝑡𝑜𝑝 𝑜𝑓 𝑟𝑒𝑎𝑐ℎ 𝑒𝑙𝑒𝑣𝑎𝑡𝑖𝑜𝑛 −𝑏𝑜𝑡𝑡𝑜𝑚 𝑜𝑓 𝑟𝑒𝑎𝑐ℎ 𝑒𝑙𝑒𝑣𝑎𝑡𝑖𝑜𝑛))/(channel centerline length (m) OR NHD stream segment polyline length)		TRUE

Elevation

If polygon is wider than DEM raster cell size, find the mid-point of the ChannelArea polygon segment at the top and bottom of the segment.

If not, use DEM cells that the top and bottom of the ChannelArea polygon fall into.

Length

Three approaches possible:

1. Use EITHER
   1. NHD polyline length if the ChannelArea polygon is just buffered off of the NHD polyline
   2. OR channel centerline off of ChannelArea polygons that were already established by the NHD (ie flowareas, waterbodies)
2. Just use NHD polyline lengths off of the stream network
3. Just use channel centerlines derived from the ChannelArea polygons

@philipbaileynar @joewheaton @KellyMWhitehead Please let me know if these are making sense or if you think they should be further discussed!

[joewheaton]

Thanks @lauren-herbine for this. What you're saying makes sense, but we need to really unpack this further and discuss. This discussion of valley bottom slope and sinuosity is actually pretty easy assuming a single center-line and such. However, channels are a mess. We're now dealing with Channel Areas. This becomes a more difficult problem than it is for NHD network.

On slope specifically, I think we should look at what we used to do in BRAT. Sometimes the precise start and end points of lines are positionally garbage. As such, we used to search in some sort of window (can't remember what it was maybe 7 x 7 or 9 x 9) for the minimum elevation point for top and do same for bottom to increase chances of getting the DEMs low point. Then for reaches where we ended up with negative slope, I think we just set slope to zero.

I know that NHD+HR has a slope estimate, but we can't just use that here as we'll have new segments defined by the valley bottom segmentation. The old version of GNAT this is easier because it is just all a function of your channel network segmentation. Our new GNAT might think about having channel network attributes vs riverscape network attributes. What is it you want/need to river style based on.

We need to think about what we actually want from a riverscapes perspective. This is easy for a riverscape where channel centerline exists (and lets not forget many do for NHD+HR). We stil have not defined what we want though for a riverscape with multiple channels and channel segments within the riverscape segment. Do we want:

• Attributes of primary anabranch
• Attributes of primary vs. secondary anabranch
• Average of all anabranches
• Average of all channels (anabranches + tributaries + yazoos)
• Max/Min of any of above?

River Styles itself is rather vague in this regard. However, one thing that we're assuming in our discussions here (implicitly) is that we need each attribute calculated for the entire network. River Styles DOES NOT operate that way. It needs certain attributes for certain valley settings... We can play that game too and maybe side-step the ones with difficult or ambiguous calculations as they have already been flagged as something totally different.

[lauren-herbine]

Hi @joewheaton thanks for starting this discussion- I'm going to move part of it over to a new discussion board (#565) as it is more broad and doesn't pertain specifically to these attributes.

[lauren-herbine]

Thanks @lauren-herbine for this. What you're saying makes sense, but we need to really unpack this further and discuss. This discussion of valley bottom slope and sinuosity is actually pretty easy assuming a single center-line and such. However, channels are a mess. We're now dealing with Channel Areas. This becomes a more difficult problem than it is for NHD network.

Just to clarify- Channel Areas that represent a single channel off of the NHD network should be simple, yes? We are mostly worried about/need more discussion on Channel Areas with multiple channels from the NHD network, yes?

On slope specifically, I think we should look at what we used to do in BRAT. Sometimes the precise start and end points of lines are positionally garbage. As such, we used to search in some sort of window (can't remember what it was maybe 7 x 7 or 9 x 9) for the minimum elevation point for top and do same for bottom to increase chances of getting the DEMs low point. Then for reaches where we ended up with negative slope, I think we just set slope to zero.

I am struggling to find where in the BRAT GitHub scripts I would figure this out- @MattReimer should I be looking in build_topography for this type of thing (tl;dr: how BRAT calculates slope using a search window for DEM high and low points)?

I know that NHD+HR has a slope estimate, but we can't just use that here as we'll have new segments defined by the valley bottom segmentation.

Yes, I agree- what I was suggesting was using the NHD polyline (segmented by our VB segments) to calculate a new slope. We will need some sort of line or something (window?) that tells us where in the VB segment we should sample the DEM at the top and bottom- so my numbered choices were more aimed at which methods we want to follow to get to those 2 points.

Does this make sense? Are we saying the same things?

[lauren-herbine]

Valley Gradient

Attribute	Tool?	Line(s) of evidence needed	Measurements needed from LOE	Calculation	Alpha (ready today)	Beta v 0.1	V 1 (i.e. thesis)	V2 (i.e. later)
Valley Gradient	VBET & DEM	Valley bottom segment centerline + DEM	elevation measurements from top and bottom of valley bottom segment (m), length of valley bottom centerline (m)	((𝑡𝑜𝑝 𝑜𝑓 𝑠𝑒𝑔𝑚𝑒𝑛𝑡 𝑒𝑙𝑒𝑣. −𝑏𝑜𝑡𝑡𝑜𝑚 𝑜𝑓 𝑠𝑒𝑔𝑚𝑒𝑛𝑡 𝑒𝑙𝑒𝑣.))/(𝑣𝑎𝑙𝑙𝑒𝑦 𝑏𝑜𝑡𝑡𝑜𝑚 𝑐𝑒𝑛𝑡𝑒𝑟𝑙𝑖𝑛𝑒)			TRUE

Some previous discussion points:

Is valley gradient stupid to base off of top of polygon/bottom of polygon?

Top of polygon/bottom of polygon would just find the center of the valley bottom polygons at both ends, and base the gradient calculations off of these elevations.

In the figure above, the pink dots are the valley bottom polygon elevations if these were to be chosen where the channel intersects the VB polygon. The teal dots are if the VB polygon elevations were taken from where the VB centerline intersects with the polygon.

I am tempted to lean towards the channel intersection method, as there are instances where the VB centerline could intersect at elevations that would artificially raise or lower the VB gradient (ie one elevation is taken from inactive floodplain, the other happens to be where the channel crosses).

After some discussion with @joewheaton : Sample across the top of the valley bottom polygon for the lowest elevation (do same for bottom of polygon) to calculate valley bottom gradient.

Do I even want valley gradient?

In most cases, valley gradient will likely be redundant with the combination of stream order and bankful width. I thought that, however, it could still be useful in some cases. I flew around to places I thought might be "odd" (a mainstem Mississippi HUC; Flambeau, WI; Kansas River). I was looking for places where the steeper part of basins were not co-located with the headwaters. I thought perhaps incision through terraces/legacy sediment might provide that.

Flambeau, WI. Very little relief across the entire basin. Headwaters are low slope, but bankful width is also still lowest here.

I thought I'd try where Mississippi tributaries come off the higher hills and enter the floodplain. Perhaps here we'd see a drastic change in slope mid-longitudinal profile? Here, from the headwaters of Knob Creek (the perennial stream in the above diagram) to the point where it spills onto the Mississippi floodplain, the valley bottom gradient is 0.7% and the channel gradient is 0.6%. The 500 meter segment (outlined in black) right where this transition occurs has a 0.5% valley bottom gradient. So- another strike out in trying to find places where VB gradient would be useful. Any thoughts?

After some discussion with @joewheaton: Minnesota- where multiple order tribs come into the mainstem with a range of bankfuls at steep slopes through the bluffs

[lauren-herbine]

Channel Sinuosity

Attribute	Tool?	Line(s) of evidence needed	Measurements needed from LOE	Calculation	Alpha (ready today)	Beta v 0.1	V 1 (i.e. thesis)	V2 (i.e. later)
Channel Sinuosity	GNAT	Stream Network + Centerlines	straight line distance between end points, segment length	(𝑛𝑒𝑡𝑤𝑜𝑟𝑘 𝑙𝑒𝑛𝑔𝑡ℎ)/(𝑠𝑡𝑟𝑎𝑖𝑔ℎ𝑡−𝑙𝑖𝑛𝑒 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒)		TRUE

This is currently relevant to single channel systems. See #565

Similar to stream gradient we have a couple of options when it comes to how we measure the channel length:

Length

Three approaches possible:

1. Use EITHER
   1. NHD polyline length if the ChannelArea polygon is just buffered off of the NHD polyline
   2. OR channel centerline off of ChannelArea polygons that were already established by the NHD (ie flowareas, waterbodies)
2. Just use NHD polylines off of the stream network (segmented at our VB sausage breaks)
3. Just use channel centerlines derived from the ChannelArea polygons

Straight line distance

Can we calculate distances from nodes? I've noticed whenever we are discussing topology, we reference nodes that will result at segment breaks. These seem like good places to measure our straight line distance from. Would these only occur where centerlines break/overlap (like in the figure below)? Or will nodes also be produced were the stream network and/or ChannelArea polygons break/overlap?

For example, from @joewheaton in #419:

[lauren-herbine]

An update on initial GNAT Stream Gradient:

Stream Gradient QA/QC

An initial look at GNAT stream gradient calculations.

HUC: 17060304

I asked @KellyMWhitehead to give me a range of stream orders to look at how everything is working at this point. Mainstems were still a struggle at this point, so those aren't included. I classified on a continuous scale with 10 bins for visualization.

Doing some data wrangling in R, it looks like there are 990 observations.

Some visual QA/QC: Looking at these different types of regions, things are lining up the way I would expect them to be.

• We have very steep slopes coming off a high mesa-like formation (2). This is where we are seeing the highest gradient- this lines up with our integrated width vs. stream gradient- our highest stream gradients occur in our smallest stream size (purple) and low widths.

• We have mid-range to low slopes at our source waters where our valley bottom width expands (1 and 3). These are the small stream size (purple) clusters we see at the mid to lower end of the stream gradient and higher end of the integrated width.

• Seeing medium sized streams (teal) fall in the mid range of integrated width and stream gradient makes sense, as those size streams would (typically) not be found in the headwaters; in narrow, steep, mesa-side channels; nor in mainstems.

• Note that I filtered for integrated widths less than 25000 (filter(integrated_width < 25000)) as there were many wide valley bottoms making visualization difficult.

  Manual testing

I did manual testing at a few places (using the closest DEM pixel value and measuring distances with the measure tool), and the values I am finding are falling within a reasonable range of the GNAT calculated values!

[KellyMWhitehead]

@philipbaileynar @lauren-herbine

gnat_17060304_20220816.zip

Metrics:

• Stream Gradient
• Valley Gradient
• Stream Order

Improvements

• centerline clipping and segment point distancing
• removed strict error checking for calculating metrics on multipart linestrings for flowlines

Known issues

• Some centerlines have extra segments that form polygons. This prevents the formation of a single linestring.
• Some metrics along flowlines at mainstem confluences are missing due to bug in clipping logic.

[lauren-herbine]

@KellyMWhitehead is stream order in "metrics"? When I open the attribute table, I am only seeing Stream Gradient and Valley Gradient. And then when I open "metric_values", I am only seeing values of metric_id 1 (STRMGRAD). Am I missing something?

[lauren-herbine]

Resolved offline

[KellyMWhitehead]

@philipbaileynar @lauren-herbine
GNAT_17060304_20220817.zip

This release includes a total of 9 metrics.

Additional Metrics:
-Headwater
-Stream Type
-Active Floodplain Area
-Active Channel Area
-Integrated Width
-Active Channel Ratio

Improvements

• Spatial view vx_metric_points is now working in qgis. This is the preferred way of displaying metric values.
• The zip includes the full gnat project, without the dem due to attachment size limits.

Known Issues (includes previously mentioned issues)

• Level paths with repeated values of seg_distance (due segmenting level paths with multipart features) are likely getting incorporated in the same window.
• Metric VALGRAD is incorrectly using the elevations found at the centerline endpoints instead of the stream endpoints.

[lauren-herbine]

This is great, thanks @KellyMWhitehead ! I'm doing some analysis window analysis that we talked about with @philipbaileynar yesterday- hopefully I can post that here by the end of the day