10

Does anyone know what type of algorithm is used in the ArcGIS Watershed tool (in the Spatial Analyst package)?

Very little information given on Esri's website ... but I suspect it may be some kind of depth/breadth search.

I've looked at these ArcGIS Online Help pages:

So yes, it uses the flow direction raster, but what algorithm is it using to traverse the raster?

Please note, I'm not looking for answers along the lines of 'it uses D8..'...D8 is not really an algorithm, but a model to help define the algorithm you would use. I.E you could implement the D8 scheme within a depth-first search algorithm and/or a breadth-first search algorithm

  • James, I'm trying to do The same thing, i.e., create a app which take a determined coordinate and give us a watershed delineation. I'm using python. Let's talk about our progress. – user28433 Mar 24 '14 at 19:21
  • I'm also using Python. I'm starting with the simpler problem of computing a flow direction grid and moving on from there. – James Mar 25 '14 at 9:57
6

The method that I've implemented in a couple of languages and believe that ESRI uses (sorry, no references other than Jenson and Domingue cited elsewhere in this page) is to start at a user-supplied "pour-point" cell or a cell at the edge of the flow direction grid (fdr), examine its eight neighbors to find which of those direct flow into the current cell, and assign those cells to the current "watershed" in the output grid. Then the function recursively calls itself once for each of the inflowing neighbors. This process repeats until all inflowing cells are exhausted for a pour-point, and will then repeat for all pour-points.

The recursive algorithm design can be pretty expensive because it can end up trying to hold lots of data in memory, having to swap/page to the disk, and therefore generally suffering i/o slow downs.

(see whuber's comment below about different methods of recursion, if you're gonna RYO)

_____________ EDIT _____________

Dug out my old C code as an example (note: Although most pythoners may want to run from the room, shouldn't be too bad). Thought it might be of interest to illustrate. Although I'm only now superficially familiar w/breadth-first vs depth-first recursion, I'm thinking that my routine is indeed depth-first (and that my natural language description above was misleading) based on this stackoverflow posting (hopefully @whuber or another person smarter than me can confirm/deny).

Code: explanation: idir is the raster of flow direction values. offset refers to the center cell that is currently being analyzed, and off checks each of that cell's neighbors. This calls another function, does_it_flow_into_me, which returns a boolean as to whether the neighboring cell's flowdir points to the current cell. If true for a neighbor, then recurse to that location.

void shed(int init_x, int init_y, int basin_id){

int i, j, offset, off, flow_dir;

offset = ((init_y - 1) * nc) + (init_x - 1);
*(basin + offset) = basin_id;


/* kernel analysis */
for (i = -1; i <  2; i++) {
    for (j = -1; j <  2; j++) {
        if ((i) || (j)) {

            off = offset + (j * nc +  i);
            flow_dir = *(idir + off);


            if (does_it_flow_into_me(i,j,flow_dir)){
                shed(init_x+i, init_y+j,basin_id);
            }
        } /*not center */
    } /* do - j */
} /* do - i */
}
  • You describe breadth-first recursion. By means of a small stack you can implement efficient depth-first recursion, which requires little memory. The main performance issue would concern large watersheds where tiles of the grid might have to be swapped in and out of RAM repeatedly. As discussed in comments to other answers, though, the real issue concerns coping with cells where there is no uniquely determined D8 direction, especially cells lying within extensive flat horizontal patches (such as those created by preliminary sink-filling routines). – whuber Jan 27 '14 at 18:55
  • Definitely a garbage in-garbage out issue. What I and most GIss do will not clean up the input! Sounds like I need to go look up depth-first recursion to put some polish on my hack. – Roland Jan 27 '14 at 19:04
  • I don't think this is garbage in--remember, regardless of how the implementation is broken down, the original input is the DEM itself rather than somebody's D8 coding--but it's definitely a challenge. The real world has many places that are so flat that the flow direction is difficult to determine: any static waterbody is a good example. In effect you need to find outlets of lakes and other flat areas and you need to cope with flat areas that have multiple outlets. This requires non-local searches, which are hard to do. – whuber Jan 27 '14 at 19:10
  • I'm probably confused then. I'm thinking we're discussing help.arcgis.com/en/arcgisdesktop/10.0/help../index.html#//…, which takes flowdir as input. Don't want to pull us into the weeds if I've not read the rest closely enough! – Roland Jan 27 '14 at 19:27
  • No, I think you're right: as I re-read the question, I see it focuses specifically on processing the flow direction raster as input, rather than on the more general situation I was envisioning. So +1 to your answer for addressing it directly and with insight and helpful pointers. – whuber Jan 27 '14 at 20:16
4

The ArcGIS help says:

Watersheds can be delineated from a DEM by computing the flow direction and using it in the Watershed tool. To determine the contributing area, a raster representing the direction of flow must first be created with the Flow Direction tool.

The Flow Direction is calculated from the DEM using the D8 method, Where the flow is abstracted by calculating for each cell, which of it's 8 neighbors, the water from this cell will flow to.

There are many alternatives to D8, such as Rho8, Froh8 & Stream Tubes, but most GIS Software including ArcGIS tend to use D8, since it is simpler, and less computationally intensive than others.


A few years ago, I was working on a Watershed Delineation project, and we were facing several issues due to ArcGIS using the D8 method. The two main problems were

  • D8 allows only Uni Directional Flow. Water can flow out only in one direction from one cell.
  • The Stream flows generated had a huge bias along the diagonal axis. This gave rise to strange looking streams.

From Our data, we knew that these two issues were big problems, so I had developed some tools to generate flow directions using hybrid methods.

One of my earliest tasks was to reverse engineer the Catchment calculation tool. I found that it was logically quite simple. If you wish to find the catchment for given point (also called the pour point), you first find the cell in which it belongs. Often you will try to snap it to the point with the highest flow accumulation in a given tolerance.

For this cell you will find all the cells in the neighborhood that contribute to it. For each of these neighborhood cells, you find the cells that contribute to them and so on. You continue this iterative process till you find no new cells. That's when you have reached the ridge lines or the watershed boundary.

I found that my simple code which did this for ASCII rasters, gave almost similar output when compared to ArcGIS's Watershed tool. Sometimes there used to be a difference of a few cells on the boundary, so I'm convinced that ArcGIS follows an unmodified D8 algorithm.

  • Thank you for the elaboration. But what is the algorithm for using the D8 directions to find watersheds? Please see the comments following dmahr's answer. – whuber Jan 27 '14 at 15:31
  • Hi, thanks but this doesn't really answer the question. You hit upon it with the sentence "For this cell you will find all the cells in the neighborhood that contribute to it. For each of these neighborhood cells, you find the cells that contribute to them and so on". There are many different algorithms to implement that search. The question is which one – James Jan 28 '14 at 13:07
4

This has been asked before, though perhaps in a slightly different context. All of the geoprocessing tools in the Hydrological toolset of Spatial Analyst use the D8 flow direction model, as stated in the How Flow Direction Works page:

There are eight valid output directions relating to the eight adjacent cells into which flow could travel. This approach is commonly referred to as an eight-direction (D8) flow model and follows an approach presented in Jenson and Domingue (1988).

A copy of the Jenson and Domingue (1988) paper is available here.

All of the tools that use Flow Direction rasters as input utilize this flow direction model by association. This includes including Watershed, Flow Accumulation, Flow Length, Fill, etc.

  • So I suppose a follow on question would be, how is that algorithm amended to return the catchment? – James Jan 27 '14 at 14:52
  • The Watershed tool navigates up the flow direction raster from the pour points. It's the reverse of the Flow Accumulation tool, except that instead of the output raster describing the number of cells, it reports the ID of the pour point. – dmahr Jan 27 '14 at 15:01
  • 1
    Ok, I guess I need to be a bit more specific. I know the concept of what it does. I don't know what algorithm is implemented. I.e. I assume it is some kind of search algorithm, but it could still be; breadth-first, depth-first, iterative-deepening depth-first etc... – James Jan 27 '14 at 15:09
  • 1
    thanks dmahr. @whuber: As far as I know, different search algorithms could give slightly different results? And yes, finding a generic algorithm is not a problem, but learning how ESRI handles watershed-specific areas (such as flat parts of a DTM) is useful. – James Jan 27 '14 at 15:26
  • 1
    James Please edit your question to clarify that last point, so that this thread stops collecting otherwise useless "it's D8" answers. (What is helpful about the D8 comments is that if you accept that D8 leads to a unique flow direction graph, then there is a unique correct solution to the watershed delineation problem, because the watersheds are properties of the graph itself. Thus if there are any ambiguities they must lie in (a) the definition of "watershed," (b) how the D8 directions are computed, or (c) how horizontal cells (i.e., without a unique D8 direction) are handled.) – whuber Jan 27 '14 at 15:30
0

To give more thought to this question, I ran a watershed analysis in arc: I took a (filled) DEM, calculated the flow direction and placed a few points that corresponded to locations on a previously calculated stream network. I ran the 'watershed' tool and it gave me a few nice basins, pretty much covering most of the remaining area 'upstream' (as you would expect):

watershed image

I then coded up a quick search algorithm in Python (like the answer above), which inspects the flow direction grid and 'follows' the flow paths. For each node, I inspect the 8 neighbours and if a neighbour flows into the current node, I call the same function recursively with the neighbour node as the input.

Pseudo(ish) code:

class d8():
    def __init__(self, arr):
       self.catchment = set()
       self.arr = arr

    def search(self, node):
        """ Searches all neighbouring nodes to find flow paths """ 

        # add the current node to the catchment
        self.catchment.add(node)

        # search the neighbours, ignore ones we already visited
        for each_neighbour:
            if neighbour is in self.catchment:
               do nothing

            # if the neighbour flows into the current node, visit that neighbour
            elif neighbour_flows_into_me:
               self.search(neighbour)

I ran that function using the same flow direction input grid and one of the same points. The problem is, where arc returns a catchment of some 40000 cells for that point, my algorithm just returns 72 cells.

Anyone know what I am doing wrong?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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