Obtaining stream flow direction

5283
21
Jump to solution
02-20-2014 06:28 AM
MarkBrown
New Contributor
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
0 Kudos
1 Solution

Accepted Solutions
XanderBakker
Esri Esteemed Contributor
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

View solution in original post

0 Kudos
21 Replies
XanderBakker
Esri Esteemed Contributor
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
0 Kudos
MarkBrown
New Contributor
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
0 Kudos
XanderBakker
Esri Esteemed Contributor
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
0 Kudos
DuncanHornby
MVP Notable Contributor
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]
0 Kudos
MarkBrown
New Contributor
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
0 Kudos
MarkBrown
New Contributor
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
0 Kudos
MarkBrown
New Contributor
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
0 Kudos
DuncanHornby
MVP Notable Contributor
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
0 Kudos
MarkBrown
New Contributor
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
0 Kudos