Select points approx. 2000 metres from another point along a river?

Question

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"])

Dictionary

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

Algorithm

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.

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.