enter image description here

I have a series of points representing a river network (See Image Above). Each point contains attributes on Latitude, Longitude, Z value and Length calculated using flow accumulation. The river network was created from Digital Elevation Models using hydrology tools in ArcGIS (Fill, Flow Direction, Flow Accumulation, Strahler Order > 4). I also have a polyline shapefile that also runs along this river network.

Starting at the first point at the beginning of each network I would like to determine another point "downstream" that is approximately 2000 metres away (Or closest to 2000 metres). I would then like to calculate channel slope based on the difference in elevation divided by the difference in length (2000m) using the two points.

I am aware that at the junctions there will be an issue, any suggestions on accomplishing this task in GIS to determine the point 2000m downstream and continuing downstream at junctions?

I am using ArcGIS 10.1 with an ArcInfo license and also have access to QGIS and GrassGIS. I am familiar with using python scripts as well.

Edit: Progress

The stream polyline has the following attributes:

x1: x coordinate of start point (Double)

y1: y coordinate of start point (Double)

x2: x coordinate of end point (Double)

y2: y coordinate of end point (Double)

length: length of polyline

shape: shape of polyline

startxy: [x1]&", "&[y1] (Text)

endxy: [x2]&", "&[y2] (Text)

I am using the tool to make arrays

start = arcpy.da.FeatureClassToNumPyArray("polyline",["x1","y1"])


G ={}
>>> for (st,sh,le,en) in zip (start,shape,length,end):
...     G[st]=(st,le,en)


def Downstream (p, x, G):
...     e = G[p]
...     if (IsNull(e)):
...         return ("Not Found")
...         if (x < 0):
...             return ("Invalid")
...             if (x < e.length):
...                 return (Along (e.shape, x))
...                 return (Downstream (e.end, x-e.length,G))

The Downstream function is receiving a runtime error where there is a key error. I have used G.keys() to see the keys, but they are not working when inputted in the function.

  • Does your polyline layer have features representing the entire flow path (start point to end point) or just the tributary (start point to junction)? I would first get features for each flow path. From there, there are many ways to get the location along river segment, such as splitting at the 2000m distance. You can then get the points of the itersections (where the splits meet), spatial join the values you need from the river segment to the points, and then do your slope calculations and what not.
    – evv_gis
    Feb 25, 2014 at 16:27
  • RivEX could do this. Use RivEX to extract your sources, create a distance field and set that to 2Km, attribute your network with RivEX then run the create downstream reach tool to create the downstream polylines. Not only do you have a polyline to do spatial selections but you could extract the end node and use that to select your nearest point.
    – Hornbydd
    Feb 25, 2014 at 23:02

1 Answer 1


This is an algorithm on a directed graph with suitably weighted edges. It is accomplished in two phases:

  1. Convert the GIS data into a suitable data structure.

  2. Run the algorithm.

Phase I: Create a data structure

The GIS capabilities you need are:

  1. "Explode" the polyline representation of the stream network into contiguous segments. (A segment merely has no branches. It does not necessarily need to start or end at a confluence or outlet.)

  2. Measure the length of each segment.

  3. Provide coordinates for the endpoints of each segment.

  4. Read the flow accumulations at the termini of each segment.

Applying these operations results in a collection of objects that eventually will be represented as directed edges in the graph. A generic object a tuple (start, end, length, shape), where

  • start is the location of the endpoint with the smaller flow accumulation,

  • end is the location of the other endpoint (with the larger flow accumulation),

  • length is the segment length, and

  • shape is its actual shape, oriented to run from start to end.

Store these in an associative array ("dictionary") keyed by start. Let's call this G.

Phase II: The algorithm

The GIS capability you need is the ability to obtain the coordinates of a point that is a distance y along a polyline q provided y is less than or equal to the length of q. Let's write this function. as Along(q, y).

Given any location p and a desired downstream distance x, the function Downstream(p, x, G) returns the coordinates of a point that is distance x downstream of p in the network G, provided such a point exists (and otherwise returns a suitable code indicating such a point cannot be found). Here is pseudocode for a recursive implementation. It performs one step in the process of moving downstream along the river system, keeping track of the distance traveled along the way.

Downstream(p, x, G) {
    # Find an edge beginning at point p.
    e = G[p]
    if (IsNull(e)) return ("Point not found in the network")
    # Travel along e as far as possible or necessary.
    if (x < 0) return ("Invalid distance")
    if (x < e.length) return (Along(e.shape, x))
    # Move to the end of e and iterate.
    return (Downstream(e.end, x - e.length, G))

You can implement it as a loop if you like by mechanically eliminating the tail recursion.

Furthermore, you can automatically find all the upstream start points (as indicated in the figure) by taking the set of all edges in G and removing those that can be fetched using the end points of some edge: this leaves those whose start points correspond to no end points; these are the sources of the uppermost reaches in the system.

  • This is a great answer, but it seems to go beyond my experience with python. I've organized the data from polyline representation to include attributes for the length of each segment, the endpoints of each segment and the flow accumulation of the termini of each segment. I have not used objects, tuples or dict in python before, so this is where I have hit a roadblock. What would be the next step? or Should I post another question on the forum with regards to this issue.
    – JC11
    Feb 26, 2014 at 20:30
  • 1
    You've already done the hard part of extracting the relevant data from the GIS. You don't need a formal "object"; that's just a convenient way of conceptualizing what the algorithm does. For instance, you could maintain the four attributes of each edge in four parallel arrays. You definitely want to learn how to use a dictionary: as you can see, they make short work of looking things up, which makes them powerful tools. For basic Python coding questions go to Stack Overflow--but first do some research online; there's lots of help out there.
    – whuber
    Feb 26, 2014 at 21:47
  • Just to make sure I'm following this correctly. I used arcpy.da.featureclasstoNumPyArray to make start, end, length and shape into arrays. I then create a dictionary where start is the key. Are end, length and shape 'values' associated with the key? From learning how to use a dictionary, there can only be one value per key. I'm not to sure how to use a dictionary in the way you describe.
    – JC11
    Mar 3, 2014 at 18:58
  • Let the value associated with a key be an index into each of the arrays.
    – whuber
    Mar 3, 2014 at 19:04
  • I have made significant progress and am just trying to get the algorithm to function. I have edited the above question describing my accomplishments so far. The only issue is a keyerror now.
    – JC11
    Mar 4, 2014 at 19:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.