Note

This page was generated from a jupyter notebook.

Comparison of FlowDirectors

Introduction

Landlab’s topographic flow-routing capability directs flow and accumulates it using two types of components:

FlowDirectors use the topography to determine how flow moves between adjacent nodes. For every node in the grid it determines the node(s) to receive flow, and the proportion of flow to send from one node to its receiver(s).

The FlowAccumulator uses the direction and proportion of flow moving between each node and (optionally) water runoff to calculate drainage area and discharge.

The FlowDirectors are method-specific. Presently landlab supports four different methods for determining flow direction:

1. FlowDirectorSteepest Flow is routed to only one node. The algorithm considers the link slopes leaving from each node and chooses the steepest downhill link to route flow along. In the case of a raster grid, only the links are considered (Landlab differentiates between links, which never cross and are located at North, South, East, and West on a raster grid, and diagonals which cross and are located at North East, North West, South East, and South West). For raster grids, this method is also known as D4 flow routing. In the case of irregular grids, all links originating from a node are consideded.

2. FlowDirectorD8 (raster only) Flow is only routed to one node but diagonals are also considered.

3. FlowDirectorMFD Flow is directed to all nodes that are located downhill of the source node. In the case of a raster grid, diagonals can be included using the keyword diagonals=True. Flow is partitioned between receiver nodes based on the relative slope along the links leading to the receiver nodes. The default method for partitioning is based on the sum of receiver slopes (partition_method='slope'). Partitioning can also be done on the basis of the square root of slope, which gives the result of a steady kinematic wave(partition_method='square_root_of_slope').

4. FlowDirectorDINF (raster only) Flow is directed to two cells based on the slope of the triangular facets that can be defined between a node and its neighbors. The steepest downhill facet is chosen and then flow is partitioned between the receiver nodes at the bottom of that facet based on the relative slopes along the facet-bounding links. (The method, known as “D-infinity”, is described by Tarboton (1997, Water Resources Research, 33(2), 309-319)).

In this tutorial we will go over more detailed examples that contrast the differences between each flow-direction algorithm. For information about how to initialize and run a FlowDirector or the FlowAccumulator, refer to the other tutorials in this section.

First, we import the necessary python modules and make a small plotting routine.

import matplotlib.pyplot as plt
from matplotlib import cm

from landlab import RasterModelGrid
from landlab.components import (
    FlowAccumulator,
    FlowDirectorD8,
    FlowDirectorDINF,
    FlowDirectorMFD,
    FlowDirectorSteepest,
)
from landlab.plot.drainage_plot import drainage_plot


# create a plotting routine to make a 3d plot of our surface.
def surf_plot(mg, surface="topographic__elevation", title="Surface plot of topography"):
    plt.figure()
    ax = plt.axes(projection="3d")

    # Plot the surface.
    Z = mg.at_node[surface].reshape(mg.shape)
    color = cm.gray((Z - Z.min()) / (Z.max() - Z.min()))
    ax.plot_surface(
        mg.x_of_node.reshape(mg.shape),
        mg.y_of_node.reshape(mg.shape),
        Z,
        rstride=1,
        cstride=1,
        facecolors=color,
        linewidth=0.0,
        antialiased=False,
    )
    ax.view_init(elev=35, azim=-120)
    ax.set_xlabel("X axis")
    ax.set_ylabel("Y axis")
    ax.set_zlabel("Elevation")
    plt.title(title)
    plt.show()

3 topographic grids

For this tutorial we will consider three different topographic surfaces that highlight the difference between each of the flow direction algorithms.

mg1 = RasterModelGrid((10, 10))
_ = mg1.add_field("topographic__elevation", mg1.y_of_node, at="node")
surf_plot(mg1, title="Grid 1: A basic ramp")

mg2 = RasterModelGrid((10, 10))
_ = mg2.add_field(
    "topographic__elevation", mg2.x_of_node + 2.0 * mg2.y_of_node, at="node"
)
surf_plot(mg2, title="Grid 2: A ramp inclined in X and in Y")

mg3 = RasterModelGrid((10, 10))
_ = mg3.add_field(
    "topographic__elevation",
    mg3.x_of_node**2 + mg3.y_of_node**2 + mg3.y_of_node,
    at="node",
)
surf_plot(mg3, title="Grid 3: A more complicated surface")

Comparing the different methods for each grid

We can illustrate the receiver node FlowDirectionSteepest has assigned to each donor node using a plotting function in Landlab called drainage_plot. We will see many of these plots in this tutorial so let’s take a moment to walk through the plot and what it contains.

• The background image (white to black) shows the values of topographic elevation of the underlying surface.

• The color of the dots inside of each pixel show the locations of the nodes and the type of node.

• The arrows show the direction of flow, and the color shows the proportion of flow that travels along that link.

• An X on top of a node indicates that node is a local sink and flows to itself.

Note that in Landlab boundary nodes, or nodes that are on the edge of a grid do not have area and do not contribute flow to nodes. These nodes can either be Fixed Gradient Nodes, Fixed Value Nodes, or Closed Nodes. With the exception of Closed Nodes the boundary nodes can receive flow.

An important step in all flow direction and accumulation is setting the proper boundary condition. Refer to the boundary condition tutorial for more information.

Grid 1: Basic Ramp

As with the Introduction to Flow Director tutorial, let’s start with the basic ramp.

mg1a = RasterModelGrid((10, 10))
_ = mg1a.add_field("topographic__elevation", mg1a.y_of_node, at="node")

fd1a = FlowDirectorSteepest(mg1a, "topographic__elevation")
fd1a.run_one_step()
plt.figure()
drainage_plot(mg1a, title="Basic Ramp using FlowDirectorSteepest")

Reassuringly we can see that the flow is being sent from high elevations at the top of the grid to low elevations at the bottom of the grid. We can also see that all of the arrows are yellow, and thus all of the flow is traveling on these links.

Now let’s see how the other FlowDirectors direct the flow on this simple grid. We don’t need to specify the surface so long as it is the field 'topographic__elevation'.

mg1b = RasterModelGrid((10, 10))
_ = mg1b.add_field("topographic__elevation", mg1b.y_of_node, at="node")

fd1b = FlowDirectorD8(mg1b)
fd1b.run_one_step()
plt.figure()
drainage_plot(mg1b, title="Basic Ramp using FlowDirectorD8")

For this ramp, the steepest slope is down a link, and not a diagonal, so FlowDirectorD8 gives the same result as FlowDirectorSteepest.

mg1c = RasterModelGrid((10, 10))
_ = mg1c.add_field("topographic__elevation", mg1c.y_of_node, at="node")

fd1c = FlowDirectorMFD(mg1c, diagonals=False)  # diagonals=False is the default option
fd1c.run_one_step()
plt.figure()
drainage_plot(mg1c, title="Basic Ramp using FlowDirectorMFD without diagonals")

Similarly, while there is more than one node below each core node, there is only one node that is connected by a link and not a diagonal. Thus FlowDirectorMFD with the keyword diagonals set to True provides the same results as FlowDirectorSteepest and FlowDirectorD8

mg1d = RasterModelGrid((10, 10))
_ = mg1d.add_field("topographic__elevation", mg1d.y_of_node, at="node")

fd1d = FlowDirectorMFD(mg1d, diagonals=True)
fd1d.run_one_step()
plt.figure()
drainage_plot(mg1d, title="Basic Ramp using FlowDirectorMFD with diagonals")

When we permit flow along diagonal connections between nodes and flow to all downhill nodes, we see a difference in the directing pattern on this simple ramp. The flow is partitioned between the three downhill nodes, and there is more flow being sent to along the link as compared with the diagonals (the links are a lighter color blue than the diagonals).

One issue we might have with the results from FlowDirectorMFD in this case is that the flow on the diagonals crosses. This is one of the problems with using diagonal connections between nodes.

mg1e = RasterModelGrid((10, 10))
_ = mg1e.add_field("topographic__elevation", mg1e.y_of_node, at="node")

fd1e = FlowDirectorDINF(mg1e)
fd1e.run_one_step()
plt.figure()
drainage_plot(mg1e, title="Basic Ramp using FlowDirectorDINF")

In FlowDirectorDINF flow is partitioned to two nodes based on steepness of the eight triangular facets surrounding each node. The partitioning is based on the relation between the link and diagonal slope that form the edge of the facet and the slope of the facet itself. When one of the facet edges has the same slope as the facet, as is the case in this ramp example, all of the flow is partitioned along that edge.

Grid 2: Inclined plane in two dimentions

Next let’s look at all the flow directors but with the inclined plane. Recall that this plane is tilted in both X and Y axes, and that is tilted more steeply in the Y direction.

mg2a = RasterModelGrid((10, 10))
_ = mg2a.add_field(
    "topographic__elevation", mg2a.x_of_node + 2.0 * mg2a.y_of_node, at="node"
)

fd2a = FlowDirectorSteepest(mg2a, "topographic__elevation")
fd2a.run_one_step()
plt.figure()
drainage_plot(mg2a, title="Grid 2 using FlowDirectorSteepest")

Flow is directed down parallel to to the the Y-axis of the plane. This makes sense in the context of the FlowDirectorSteepest algorithm; it only sends flow to one node, so it an idealized geometry such as the plane in this example, it provides flow direction that is non-realistic.

As we will discuss throughout this tutorial, there are benefits and drawbacks to each FlowDirector algorithm.

mg2b = RasterModelGrid((10, 10))
_ = mg2b.add_field(
    "topographic__elevation", mg2b.x_of_node + 2.0 * mg2b.y_of_node, at="node"
)

fd2b = FlowDirectorD8(mg2b)
fd2b.run_one_step()
plt.figure()
drainage_plot(mg2b, title="Grid 2 using FlowDirectorD8")

FlowDirectorD8 consideres the diagonal connections between nodes. As the plane is inclined to the southwest the flow direction looks better here, though as we will see later, sometimes FlowDirectorD8 does non-realistic directing too.

mg2c = RasterModelGrid((10, 10))
_ = mg2c.add_field(
    "topographic__elevation", mg2c.x_of_node + 2.0 * mg2c.y_of_node, at="node"
)

fd2c = FlowDirectorMFD(mg2c, diagonals=False)  # diagonals=False is the default option
fd2c.run_one_step()
plt.figure()
drainage_plot(mg2c, title="Grid 2 using FlowDirectorMFD without diagonals")

As FlowDirectorMFD can send flow to all the nodes downhill it doesn’t have the same problem that FlowDirectorSteepest had. Because the plane is tilted down more steeply to the south than to the east, it sends more flow on the steeper link.

mg2d = RasterModelGrid((10, 10))
_ = mg2d.add_field(
    "topographic__elevation", mg2d.x_of_node + 2.0 * mg2d.y_of_node, at="node"
)

fd2d = FlowDirectorMFD(mg2d, diagonals=True)
fd2d.run_one_step()
plt.figure()
drainage_plot(mg2d, title="Grid 2 using FlowDirectorMFD with diagonals")

When FlowDirectorMFD considers diagonals in addition to links, we see that it sends the flow to four nodes instead of three. While all of the receiver nodes are downhill from their donor nodes, we see again that using diagonals permits flow to cross itself. We also see that the most flow is routed to the south and the south east, which makes sense based on how the plane is tilted.

mg2e = RasterModelGrid((10, 10))
_ = mg2e.add_field(
    "topographic__elevation", mg2e.x_of_node + 2.0 * mg2e.y_of_node, at="node"
)

fd2e = FlowDirectorDINF(mg2e)
fd2e.run_one_step()
plt.figure()
drainage_plot(mg2e, title="Basic Ramp using FlowDirectorDINF")

Here FlowDirectorDINF routes flow in two directions, to the south and southeast. The plane is steeper to from north to south than from east to west and so more flow is directed on the diagonal to the southeast.

Grid 3: Curved surface

Finally, let’s consider our curved surface.

mg3a = RasterModelGrid((10, 10))
_ = mg3a.add_field(
    "topographic__elevation",
    mg3a.x_of_node**2 + mg3a.y_of_node**2 + mg3a.y_of_node,
    at="node",
)

fd3a = FlowDirectorSteepest(mg3a, "topographic__elevation")
fd3a.run_one_step()
plt.figure()
drainage_plot(mg3a, title="Grid 3 using FlowDirectorSteepest")

Flow on this surface using FlowDirectorSteepest looks realistic, as flow is routed down into the bottom of the curved surface.

mg3b = RasterModelGrid((10, 10))
_ = mg3b.add_field(
    "topographic__elevation",
    mg3b.x_of_node**2 + mg3b.y_of_node**2 + mg3b.y_of_node,
    at="node",
)

fd3b = FlowDirectorD8(mg3b)
fd3b.run_one_step()
plt.figure()
drainage_plot(mg3b, title="Grid 3 using FlowDirectorD8")

Near the bottom left of the grid, the steepest descent is on a diagonal, so using FlowDirectorD8 gives a different drainage pattern.

mg3c = RasterModelGrid((10, 10))
_ = mg3c.add_field(
    "topographic__elevation",
    mg3c.x_of_node**2 + mg3c.y_of_node**2 + mg3c.y_of_node,
    at="node",
)

fd3c = FlowDirectorMFD(mg3c, diagonals=False)  # diagonals=False is the default option
fd3c.run_one_step()
plt.figure()
drainage_plot(mg3c, title="Grid 3 using FlowDirectorMFD without diagonals")

Permitting multiple receivers with and without diagonals give an additional two different drainage patterns.

mg3d = RasterModelGrid((10, 10))
_ = mg3d.add_field(
    "topographic__elevation",
    mg3d.x_of_node**2 + mg3d.y_of_node**2 + mg3d.y_of_node,
    at="node",
)

fd3d = FlowDirectorMFD(mg3d, diagonals=True)
fd3d.run_one_step()
plt.figure()
drainage_plot(mg3d, title="Grid 3 using FlowDirectorMFD with diagonals")

Again we see flow paths crossing when we permit consideration of flow along the diagonals.

mg3e = RasterModelGrid((10, 10))
_ = mg3e.add_field(
    "topographic__elevation",
    mg3e.x_of_node**2 + mg3e.y_of_node**2 + mg3e.y_of_node,
    at="node",
)

fd3e = FlowDirectorDINF(mg3e)
fd3e.run_one_step()
plt.figure()
drainage_plot(mg3e, title="Grid 3 using FlowDirectorDINF")

Finally we see yet a different drainage pattern when we use FlowDirectorDINF and flow is routed along an adjacent diagonal-link pair.

Comparison of Accumulated Area

Before concluding, let’s examine the accumulated drainage area using each of the FlowDirector methods and the third grid. For an introduction to creating and running a FlowAccumulator see the tutorial “Introduction to Flow Accumulators”.

Often we do flow routing and accumulation because we want to use the accumulated area as a proxy for the water discharge. So the details of how the flow is routed are important because they influence how the drainage area pattern evolves.

Lets begain with FlowDirectorSteepest.

mg3 = RasterModelGrid((10, 10))
_ = mg3.add_field(
    "topographic__elevation",
    mg3.x_of_node**2 + mg3.y_of_node**2 + mg3.y_of_node,
    at="node",
)

fa = FlowAccumulator(mg3, "topographic__elevation", flow_director="Steepest")
fa.run_one_step()
plt.figure()
drainage_plot(
    mg3, "drainage_area", title="Flow Accumulation using FlowDirectorSteepest"
)

Here we see that flow has accumulated into one channel in the bottom of the curved surface.

fa = FlowAccumulator(mg3, "topographic__elevation", flow_director="D8")
fa.run_one_step()
plt.figure()
drainage_plot(mg3, "drainage_area", title="Flow Accumulation using FlowDirectorD8")

When diagonals are considered, as in FlowDirectorD8, the drainage patter looks very diferent. Instead of one channel we have two smaller channels.

mg3 = RasterModelGrid((10, 10))
_ = mg3.add_field(
    "topographic__elevation",
    mg3.x_of_node**2 + mg3.y_of_node**2 + mg3.y_of_node,
    at="node",
)

fa = FlowAccumulator(mg3, "topographic__elevation", flow_director="MFD")
fa.run_one_step()
plt.figure()
drainage_plot(
    mg3,
    "drainage_area",
    title="Flow Accumulation using FlowDirectorMFD without diagonals",
)

Flow is distributed much more when we use FlowDirectorMFD.

mg3 = RasterModelGrid((10, 10))
_ = mg3.add_field(
    "topographic__elevation",
    mg3.x_of_node**2 + mg3.y_of_node**2 + mg3.y_of_node,
    at="node",
)

fa = FlowAccumulator(mg3, "topographic__elevation", flow_director="MFD", diagonals=True)
fa.run_one_step()
plt.figure()
drainage_plot(
    mg3, "drainage_area", title="Flow Accumulation using FlowDirectorMFD with diagonals"
)

Adding diagonals to FlowDirectorMFD gives a channel somewhat similar to the one created by FlowDirectorSteepest but much more distributed.

mg3 = RasterModelGrid((10, 10))
_ = mg3.add_field(
    "topographic__elevation",
    mg3.x_of_node**2 + mg3.y_of_node**2 + mg3.y_of_node,
    at="node",
)

fa = FlowAccumulator(mg3, "topographic__elevation", flow_director="DINF")
fa.run_one_step()
plt.figure()
drainage_plot(mg3, "drainage_area", title="Flow Accumulation using FlowDirectorDINF")

Finally, FlowDirectorDINF gives yet another pattern for the accumulation of drainage area.

Conclusion

This tutorial compared the different methods in more detail and over surfaces that are more complicated than a simple sloping ramp. It also described how these different FlowDirector methods change the patterns of accumulated drainage area.

Next consider one of two additional tutorials about directing and accumulating flow in Landlab.

1. Introduction to FlowDirector: A tutorial that goes over the different FlowDirectors present in Landlab and how to create and run a FlowDirector.

2. Introduction to FlowAccumulator: A tutorial that describes how to use the FlowAccumulator.

---

Generated by nbsphinx from a Jupyter notebook.