nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Obtaining stream flow direction

MarkBrown

Greetings,

I have a polyline stream layer (shapefile) and I want to know if there is a way to find the cardinal flow direction (N, NE, E, SE, S...) of each individual stream segment? I have a DEM, but I am not interested in the terrain flow direction that water would normally flow. I am looking to find the cardinal flow direction of the stream(s) segments themselves. My ultimate goal is to determine the average stream flow direction within a polygon layer using a spatial join.

I am not sure if this matters, but I have From_Nodes and To_Nodes, reach codes and distances in my attribute table.

I hope this makes sense.

Any and all information, tools, tips, tricks is greatly appreciated.

Thanks....

Mark

Find more posts tagged with

Accepted answers

XanderBakker

Hi Mark,

Yes you can. In the Python script locate the following section of code:

def getCardinal(angle):     lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"],                    [157.5, "NW"], [202.5, "W"], [247.5, "SW"],                    [292.5, "S"], [337.5, "SE"], [360, "E"]]

This is a nested list, where each element (e.g. "[22.5,"E"]") in the list holds information on the upper limit of the angles class (e.g. 22.5°) and the cardinal flow direction (e.g. "E"). You can expand the list with as many items as you want, as long as the angles are in the right order. So:

def getCardinal(angle):     lstCardinal = [[11.25,"E"], [33.75,"ENE"],                     [56.25, "NE"], [78.75, "NNE"],                     [101.25, "N"], [123.75, "NNW"],                      etc...                    [348.75, "ESE"], [360, "E"]]

Just remember that the angle is arithmetic (not geographic). 0° is East and the angle increments counter clockwise...

BTW; if you consider your question answered, please check the mark (turns green when you click on it) next to the post that answered you question most.

Kind regards,

Xander

All comments

XanderBakker

Hi Mark,

Yes you can. In the Python script locate the following section of code:

def getCardinal(angle):     lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"],                    [157.5, "NW"], [202.5, "W"], [247.5, "SW"],                    [292.5, "S"], [337.5, "SE"], [360, "E"]]

def getCardinal(angle):     lstCardinal = [[11.25,"E"], [33.75,"ENE"],                     [56.25, "NE"], [78.75, "NNE"],                     [101.25, "N"], [123.75, "NNW"],                      etc...                    [348.75, "ESE"], [360, "E"]]

MarkBrown

Hi Mark,

I had some problems with the "custom" coordinate system, but I have attached a toolbox with the functionality and also an updated file geodatabase. Unzip the toolbox ("toolshare2.zip") and navigate with the ArcCatalog window to the unzip location. Try the tools inside. Hope they work for you.

The SL distance (length) that is returned follows the linear unit of the input featureclass. In your case it is US_Foot.

Kind regards,

Xander

Xander...

Worked like a charm. I tried it on a couple of other datasets and there does not seem to be any issues. The data is consistent to what I was looking for. Much more manageable as well.

One other question.....if in the future I wanted to further divide the cardinal flow directions from the earlier script into even smaller units of degrees (i.e. 0-11.25 degrees N, 11.25 - 22.5 NNE, etc...) I could modify that script to do calculate those flow directions/angles, correct?

You and Duncan have been of great help. Thanks for all your time and patience with me.

Greatly appreciated...

v/r

Mark

XanderBakker

Hi Mark,

I had some problems with the "custom" coordinate system, but I have attached a toolbox with the functionality and also an updated file geodatabase. Unzip the toolbox ("toolshare2.zip") and navigate with the ArcCatalog window to the unzip location. Try the tools inside. Hope they work for you.

The SL distance (length) that is returned follows the linear unit of the input featureclass. In your case it is US_Foot.

Kind regards,

Xander

toolshare2.zip

Export_Streams_Test_results.gdb.zip

MarkBrown

Hi Mark,

When I was asking for a few streams, I was refering to a featureclass. Maybe you can attach a few streams as featureclass in a zipped fgdb or as a zipped shapefile (although some fieldnames will be renamed when exporting to shapefile). This way I can look at the geometry itself and the spatial reference and test my script.

Kind regards,

Xander

Xander...

Please find the attached zipped FGDB. I hope this works for you.

thanks...

Mark

Export.zip

XanderBakker

Hi Mark,

When I was asking for a few streams, I was refering to a featureclass. Maybe you can attach a few streams as featureclass in a zipped fgdb or as a zipped shapefile (although some fieldnames will be renamed when exporting to shapefile). This way I can look at the geometry itself and the spatial reference and test my script.

Kind regards,

Xander

MarkBrown

Hi Mark,

If your data is located somewhere at the border of Pennsylvania and West Virginia, the coordinates are in decimal degrees (not miles). The coordinates seem to correspond to a stream, just next to the Golden Oaks rd (PA 18).

Using decimal degrees to calculate the angle and the straight distance in the scrips will not give correct answers. The coordinates should be projected first.

I don't know what the reason is for the calculation to trow a 99999 error. If you are willing to share a small part of your data (streams) I can have a look.

Kind regards,

Xander

Hi Xander...

Indeed, I should have originally provided you with a snippit of my data.

In the first image (original data)....this was my data as I started with it. As you can see, I have multiple reach codes for individual segments. In my From and To nodes, there are decimal units....which are likely in miles because the stream line distances match up between the nodes. I did not originally have Lat/Long coordinates and used the calculate geometry to calculate those fields. The FlowAngle and CardFD fields were calculated based on the previous python scripts. BEFORE I HAD EVEN OBTAINED THE LAT/LONG COORDINATES. And the script appeared to do what you had said it would.

But for my flow direction data to be more representative of my polygon layer, my thinking is that it would be best to combine the streams based on the same reach codes. Thus the combined stream would have a single From and To node. This meant DISSOLVING my data or using the unsplit line tool to accomplish that task (Data Dissolve). Both seemed to acquire the same reasults. After dissolviing my data, I again calculated the Lat/Long coordinates of the From and To nodes. From there I believed I should be able to get the actual stream length, the straight line stream length and the flow angle or flow bearing of the straight line distance.....my original intended goal.

When I again tried running the scripts, I acquired the 99999 error code. That is where I stand at this point.

I hope this makes sense...

thanks for your time.

v/r

Mark

Original_data.jpg

Data_Dissolve.JPG

XanderBakker

Hi Mark,

If your data is located somewhere at the border of Pennsylvania and West Virginia, the coordinates are in decimal degrees (not miles). The coordinates seem to correspond to a stream, just next to the Golden Oaks rd (PA 18).

Using decimal degrees to calculate the angle and the straight distance in the scrips will not give correct answers. The coordinates should be projected first.

I don't know what the reason is for the calculation to trow a 99999 error. If you are willing to share a small part of your data (streams) I can have a look.

Kind regards,

Xander

MarkBrown

Hi Mark,

This sounds a lot like it has to do with the coordinate systems. Can you specify the following:
- coordinate system of your streams dataset
- coordinate system of your data frame

If both are geographic, can you provide me a projected coordinate that suits your data (for the length calculation)? Probably the straight line geometry has to be created with the spatial reference of the source, and if that is geographic, it needs to be projected to obtain the proper length. This can all be done in the field calculator script.

Just to be sure can you also provide an example:
- start coordinate
- end coordinate
- distance measured manually (the correct distance)
- distance calculated by the script

Kind regards,

Xander

Hi Xander,

Let me back up just a bit and give you an idea of my data set. I have NHD data with thousands of individual stream segments. If you recall my initial pursuit was to find the stream length and flow direction of the individual segments. The information from you and Duncan proved fruitful in that effort. But as you recall, my combined flow direction was not matching up to my needs. I then asked how I could get a straight line distance between my starting and ending nodes.

My nodes were in some kind of decimal units....(I am guessing miles) because when I calculated the stream length in miles (using field calculator) everything pretty much agreed, but my straight line distances from the above python script did not. As a novice to all this, I discovered that I probably needed coordinates for my nodes in order to get a straight line distance between the first and last nodes as well as a bearing.

So....in order to get to a straight line distance, I had to Dissolve my dataset based on my reachcodes....meaning I had multilple segments that made up one complete stream segment. That reduced the number of individual segments I had to deal with, which did not carry over the original decimal node values.

I then created separate fields to calculate the lat/long coordinates of my TO and FROM nodes. It appears that I can calculate the straight line distance. It will calculate it in miles, but I can do a conversion with field calculator. If I calculate geometry, it will calculate the "actual stream length". So, what I need is to calculate the flow angle or bearing from that line, but I get a 99999 error It will not calculate the field based on the python script provided earlier nor will it calculate my flow angles or cardinal bearing.

Here is an example in my combined streams layer.....

Streams layer and data frame: Projected Coordinate System:NAD_1983_Lambert_Conformal_Conic
Geographic Coordinate System: GCS_North_American_1983

In this particular example it is the result of two stream segments that I combined through DISSOLVE based on reach code.

From_Node_X (-80.42511)
From_Node_Y (39.72707)

To_Node_X (-80.43147)
To_Node_Y (39.72133)

Actual stream length (calculate geometry) 0.57185 miles
Straight line distance: 2752.176 feet

Long story short....from my new dissolved data, I need to once again find the resultant flow angle or bearing from my start and end points of my nodes.

Thanks for your time...

v/r

Mark

XanderBakker

Hi Mark,

This sounds a lot like it has to do with the coordinate systems. Can you specify the following:
- coordinate system of your streams dataset
- coordinate system of your data frame

If both are geographic, can you provide me a projected coordinate that suits your data (for the length calculation)? Probably the straight line geometry has to be created with the spatial reference of the source, and if that is geographic, it needs to be projected to obtain the proper length. This can all be done in the field calculator script.

Just to be sure can you also provide an example:
- start coordinate
- end coordinate
- distance measured manually (the correct distance)
- distance calculated by the script

Kind regards,

Xander

MarkBrown

Hi Mark,

I see Duncan has already provided you with a good way to calculate the intra-nodal distance. I was wondering however, since your goal is to "determine the average stream flow direction within a polygon layer" and you don't want to use your DEM (aspect) for this purpose, you could do something else.

What if you would create fields for each cardinal flow direction and each field would be filled with the sum of the lengths of coordinate pairs within each stream in that cardinal flow direction. My guess is that this would provide a better understanding of the flow direction once you summarize it within a polygon.

Kind regards,

Xander

Hi Xander and Duncan...

Before I saw Duncan's latest post I was able to figure out how to run the script. I am getting a calculation, but for some reason the straight line distances that are calculated aren't what they should be. They are much smaller (way smaller) than the actual stream lengths. My projections and coordinate systems agree. I'm not sure if the issue resides in the from and to nodes in the attribute table or what. You folks have been a tremendous help thus far, but I'm still perplexed in trying to obtain the straight line distances.

Thanks....

Mark

XanderBakker

Hi Mark,

I see Duncan has already provided you with a good way to calculate the intra-nodal distance. I was wondering however, since your goal is to "determine the average stream flow direction within a polygon layer" and you don't want to use your DEM (aspect) for this purpose, you could do something else.

What if you would create fields for each cardinal flow direction and each field would be filled with the sum of the lengths of coordinate pairs within each stream in that cardinal flow direction. My guess is that this would provide a better understanding of the flow direction once you summarize it within a polygon.

Kind regards,

Xander

Hornbydd

Mark,

Having created a field of type double, you right click on field header run field calculator and complete as shown below.

[ATTACH=CONFIG]31854[/ATTACH]

Duncan

Field Calc.png

MarkBrown

Mark,

You can compute the straight line distance of any polyline with the following python field calculate. First create a field of type double then run a calculate with the python expression getSLDist(!shape!) and the pre-logic script code is:
import arcpy
def getSLDist(geom):
  # Get the end points of the polyline
  fp = geom.firstPoint
  tp = geom.lastPoint

  # Create an Array object and add points to it
  array = arcpy.Array()
  array.add(fp)
  array.add(tp)

  # Create a polyline from array, as it only has 2 points in it
  # it must be a straight line, then return its length
  line = arcpy.Polyline(array)
  dist = line.length
  return dist
Duncan

Hi Duncan...

Thanks for the reply. What am I supposed to add to the CodeBlock box (i.e. Straight_Length

? I placed the script code into the pre-logic box, but my OK box is greyed out.

If I run the field calculator with the script in the codeblock, I receive a Python error: 000539 : Error message from Python.

Thanks...

Mark

Hornbydd

Mark,

You can compute the straight line distance of any polyline with the following python field calculate. First create a field of type double then run a calculate with the python expression getSLDist(!shape!) and the pre-logic script code is:

import arcpy
def getSLDist(geom):
  # Get the end points of the polyline
  fp = geom.firstPoint
  tp = geom.lastPoint

  # Create an Array object and add points to it
  array = arcpy.Array()
  array.add(fp)
  array.add(tp)

  # Create a polyline from array, as it only has 2 points in it
  # it must be a straight line, then return its length
  line = arcpy.Polyline(array)
  dist = line.length
  return dist

Duncan

MarkBrown

Hi Mark,

Although this type of analysis is normally done using the raster format, it can be done using a featureclass (shapefile). Assuming that you are interested in an angle between the start point and the end point of each feature you could do this:

open the attribute table of your shapefile
add a new field that will hold the cardinal flow direction, make sure it is a text field (let's call this field "CardFD")
right click on your new field "CardFD" and select "Field Calculator..." to open the Field Calculator
Choose Python as parser
Switch the "Show Codeblock" on
Paste the code below in the "Pre-Logic Script Code":
import arcpy
from math import atan2, pi

def getCardinal(angle):
    lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"],
                   [157.5, "NW"], [202.5, "W"], [247.5, "SW"],
                   [292.5, "S"], [337.5, "SE"], [360, "E"]]

    for item in lstCardinal:
        value = item[0]
        if angle < value:
            cardinal = item[1]
            break
    return cardinal

def calcGeomCardinality(polyline):
    pnt1 = polyline.firstPoint
    pnt2 = polyline.lastPoint
    angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi
    if angle_deg < 0:
        angle_deg = 360 + angle_deg
    return getCardinal(angle_deg)
Paste the line below in as formula (just below "CardFD =")
calcGeomCardinality( !SHAPE!)
Click "OK"
This will take each polyline feature and extract the from and to node from it. Next it will calculate the angle and parse it to cardinal flow direction.

You can also write the angle itself to a new field:

add a new field (double or float) that can hold the angle between the start and end point . Let's call this field "FlowAngle"
right click on your new field "FlowAngle" and select "Field Calculator..." to open the Field Calculator
it will probably show the code from cardinal flow direction. Replace the code in the "Pre-Logic Script Code" by
import arcpy
from math import atan2, pi

def getAngle(polyline):
    pnt1 = polyline.firstPoint
    pnt2 = polyline.lastPoint
    angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi
    if angle_deg < 0:
        angle_deg = 360 + angle_deg
    return angle_deg
Paste the line below in as formula (just below "FlowAngle=")
getAngle( !SHAPE!)
Click "OK"

If you are going to perform some analysis to determine the predominant angle, please be aware that the resulting angle does not have to be representative for the area.

Kind regards,

Xander

Hi Xander...

The steps you outlined were great and presented me with some other challenges and questions.

Forgive me as I am far from understanding Python script and programing. Although I have the actual stream meandering lengths between the nodes, is there a way to modify the script to calculate the intra-nodal distance as well as the direction? As you mentioned, "the resulting angle I received was not always representative for the area when taking into account all of my stream features."

The lengths of the features in the feature class are (almost) invariably going to be greater than the distance between nodes. I believe the stream-network statistics I am trying to obtain will be better served by the straight distance between nodes (as the crow flies) than they will by the actual length of the feature.

So, is there a way to modify the script or a tool I can use to calculate the inter nodal distance and obtain a flow angle and cardinal direction?

Again, any and all information is greatly appreciated.

v/r

Mark

Stream_example.pdf

MarkBrown

Hi Duncan,

I have not tried the Linear Direction Mean Tool. I will have a look at it and see how it works. I greatly appreciate the response.

v/r

Mark

MarkBrown

Hi Xander...

Thanks for the follow up. After I sent my post I realized I was thinking geographic. It was too early in the morning. Again thanks for the post. Everything seems to be in line and makes sense.

Mark

Hornbydd

Mark,

Have you considered the Linear Direction Mean tool? You can wrap this up in a very simple model. May not be the most speedy way of doing this but will generate a direction line for every segment. See the model below, knocked this up in 2 minutes...

You could run this on a subset of your data if you select some lines.

Duncan

[ATTACH=CONFIG]31729[/ATTACH]

Capture.PNG

XanderBakker

Hi Mark,

The error refers to something that went wrong with the field calculation. Since the code does not account for any type of error (polyline of 0 length, or polyline that has the same start and end point) there can be many reason for the error to occur.

The angle is arithmetic (not geographic). 0° is East and the angle increments counter clockwise:

< 22.5: "E"
< 67.5: "NE"
<112.5: "N"
<157.5: "NW"
<202.5: "W"
<247.5: "SW"
<292.5: "S"
<337.5: "SE"
<360.0: "E"

Kind regards,

Xander

MarkBrown

Hi Xander....

Thanks for your informative reply. I did just as you indicated and although I received an error (000539 : Error message from Python), the CardFD field did populate with a cardinal direction. I have spot checked a number of the individual segments and they appear to be in line with the cardinal direction.

I also received the same error for the FlowAngle field, but the field did populate with an angle. However, the FlowAngle does not always appear to be representative of the cardinal direction (i.e. FlowAngle 2.829 was populated for an East cardinal direction). In my thinking that should be North.

I greatly appreciate your help and time in responding to my inquery.

v/r

Mark

XanderBakker

Hi Mark,

Although this type of analysis is normally done using the raster format, it can be done using a featureclass (shapefile). Assuming that you are interested in an angle between the start point and the end point of each feature you could do this:

open the attribute table of your shapefile
add a new field that will hold the cardinal flow direction, make sure it is a text field (let's call this field "CardFD")
right click on your new field "CardFD" and select "Field Calculator..." to open the Field Calculator
Choose Python as parser
Switch the "Show Codeblock" on
Paste the code below in the "Pre-Logic Script Code":

import arcpy
from math import atan2, pi

def getCardinal(angle):
    lstCardinal = [[22.5,"E"], [67.5, "NE"], [112.5, "N"],
                   [157.5, "NW"], [202.5, "W"], [247.5, "SW"],
                   [292.5, "S"], [337.5, "SE"], [360, "E"]]

    for item in lstCardinal:
        value = item[0]
        if angle < value:
            cardinal = item[1]
            break
    return cardinal

def calcGeomCardinality(polyline):
    pnt1 = polyline.firstPoint
    pnt2 = polyline.lastPoint
    angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi
    if angle_deg < 0:
        angle_deg = 360 + angle_deg
    return getCardinal(angle_deg)

Paste the line below in as formula (just below "CardFD =")

calcGeomCardinality( !SHAPE!)

Click "OK"

This will take each polyline feature and extract the from and to node from it. Next it will calculate the angle and parse it to cardinal flow direction.

You can also write the angle itself to a new field:

add a new field (double or float) that can hold the angle between the start and end point . Let's call this field "FlowAngle"
right click on your new field "FlowAngle" and select "Field Calculator..." to open the Field Calculator
it will probably show the code from cardinal flow direction. Replace the code in the "Pre-Logic Script Code" by

import arcpy
from math import atan2, pi

def getAngle(polyline):
    pnt1 = polyline.firstPoint
    pnt2 = polyline.lastPoint
    angle_deg = (atan2(pnt2.Y - pnt1.Y, pnt2.X - pnt1.X)) * 180.0 / pi
    if angle_deg < 0:
        angle_deg = 360 + angle_deg
    return angle_deg

Paste the line below in as formula (just below "FlowAngle=")

getAngle( !SHAPE!)

Click "OK"

If you are going to perform some analysis to determine the predominant angle, please be aware that the resulting angle does not have to be representative for the area.

Kind regards,

Xander

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.