What Algorithm is used by ArcGIS Watershed tool?

Question

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:

• Watershed (Spatial Analyst)
• How Watershed works

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

Answer 1

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 */
}

Answer 2

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.

Answer 3

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.

Answer 4

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):

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?