<jupyter_start><jupyter_text>Comparison of FlowDirectors IntroductionLandlab'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.<jupyter_code>%matplotlib inline
# import plotting tools
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib as mpl

# import numpy
import numpy as np

# import necessary landlab components
from landlab import RasterModelGrid, HexModelGrid
from landlab.components import FlowAccumulator
from landlab.components import (
    FlowDirectorD8,
    FlowDirectorDINF,
    FlowDirectorMFD,
    FlowDirectorSteepest,
)

# import landlab plotting functionality
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"):

    fig = plt.figure()
    ax = fig.gca(projection="3d")

    # Plot the surface.
    Z = mg.at_node[surface].reshape(mg.shape)
    color = cm.gray((Z - Z.min()) / (Z.max() - Z.min()))
    surf = 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()<jupyter_output><empty_output><jupyter_text>3 topographic gridsFor this tutorial we will consider three different topographic surfaces that highlight the difference between each of the flow direction algorithms.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Comparing the different methods for each gridWe 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 RampAs with the Introduction to Flow Director tutorial, let's start with the basic ramp.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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'```.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>For this ramp, the steepest slope is down a link, and not a diagonal, so **FlowDirectorD8** gives the same result as **FlowDirectorSteepest**.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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**<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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 dimentionsNext 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.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>**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.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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 surfaceFinally, let's consider our curved surface.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Flow on this surface using **FlowDirectorSteepest** looks realistic, as flow is routed down into the bottom of the curved surface.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Near the bottom left of the grid, the steepest descent is on a diagonal, so using **FlowDirectorD8** gives a different drainage pattern.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Permitting multiple receivers with and without diagonals give an additional two different drainage patterns.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Again we see flow paths crossing when we permit consideration of flow along the diagonals.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>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 AreaBefore 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**.<jupyter_code>from landlab.components import FlowAccumulator

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"
)<jupyter_output><empty_output><jupyter_text>Here we see that flow has accumulated into one channel in the bottom of the curved surface.<jupyter_code>fa = FlowAccumulator(mg3, "topographic__elevation", flow_director="D8")
fa.run_one_step()
plt.figure()
drainage_plot(mg3, "drainage_area", title="Flow Accumulation using FlowDirectorD8")<jupyter_output><empty_output><jupyter_text>When diagonals are considered, as in **FlowDirectorD8**, the drainage patter looks very diferent. Instead of one channel we have two smaller channels.<jupyter_code>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",
)<jupyter_output><empty_output><jupyter_text>Flow is distributed much more when we use **FlowDirectorMFD**.<jupyter_code>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"
)<jupyter_output><empty_output><jupyter_text>Adding diagonals to **FlowDirectorMFD** gives a channel somewhat similar to the one created by **FlowDirectorSteepest** but much more distributed.<jupyter_code>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")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08713626114034145, 'fraction_numerical': 0.01868356061419521, 'mean_word_length': 5.631185475258099, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 89, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '74058', 'license': 'MIT', 'n_tokens_mistral': 6208, 'n_tokens_neox': 5830, 'n_words': 2056, 'path': 'notebooks/tutorials/flow_direction_and_accumulation/compare_FlowDirectors.ipynb', 'repo_name': 'Joseph-A-Mason/landlab', 'size': '18626'}