Note

This page was generated from a jupyter notebook.

Introduction to PriorityFloodFlowRouter (priorityFlood filler, director and accumulator)

Landlab directs flow and accumulates it using FlowDirectors to determine how flow moves between adjacent nodes and a FlowAccumulator that uses the direction and proportion of flow moving between each node and (optionally) water runoff to calculate drainage area and discharge.

Here we showcase an alternative flow router and accumualtor that wraps an external python package (richdem). The PriorityFloodFlowRouter component combines filling, flow direction and accumualtion of water over structured grids in one component and significantly improves performaces of flow accumualtion and direction operations over large grids, especially when filling operations are required to route water over landscapes.

This notebook follows the structure of the FlowAccumulator and flow director notebook.
First, we import the necessary python modules and make a small plotting routine.
[ ]:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm

from landlab import RasterModelGrid
from landlab.components import PriorityFloodFlowRouter
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()

Topographic grids

For this tutorial we will consider one topographic surface. Here it is plotted in three dimensions.

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
surf_plot(mg, title="Grid 1")

Initalizing and running the FlowAccumulator

To instantiate the PriorityFloodFlowRouter, you must pass it the minimum of a model grid that has a field called 'topographic__elevation'. The default behavior to deal with depressions internally. Taking care of depressions is done using the ‘RichDEM’ filling or breaching algorithm. The depressionHandler variable can be set to ‘fill or ‘breach’ to use respectively a filling or breaching algorithm to handle depressions. If DEMs do not need to be filled, set updateFlowDepressions to False.

Similar to other landlab flow accumulators, PriorityFloodFlowRouter can take a constant or spatially variable input called runoff_rate, which it uses to calculate discharge. Alternatively, if there is an at_node field called water__unit_flux_in and no value is specified as the runoff_rate, PriorityFloodFlowRouter will use the values stored in water__unit_flux_in.

In addition to directing flow and accumulating it in one step, FlowAccumulator can also deal with depression finding internally. This can be done by passing a DepressionFinder to the keyword argument depression_finder. The default behavior is to not deal with depressions internally.

The flow_metric variable controlls the way in which water is routed over teh landscape: Various options possible: * D8 (O’Callaghan and Mark, 1984) {default} * Rho8 (Fairfield and Leymarie, 1991) * Quinn (1991) * Freeman (1991) * Holmgren (1994) * D∞ (Tarboton, 1997) For details and comparison, see https://richdem.readthedocs.io/en/latest/flow_metrics.html Some of these methods require an exponent, which can be provided throught the exponent variable. When routing flow, one challenge is to define water directions in flat regions. This is way an articifical small graident is introduced in flat areas. If you want to have such routing option in flat areas, set the epsilon variable to True {default is True}. This ensures that there is always a local gradient.

For some applications both single and multiple flow direction and accumulation is required. By calculating them in the same component, we can optimize procedures invovled with filling and breaching of DEMs. Set the seperate_Hill_Flow variable to True if you want a second flow accumualtor calcuatled on the same filled DEM data.

[ ]:

fa = PriorityFloodFlowRouter(mg)

# # this is the same as writing:
fa = PriorityFloodFlowRouter(
    mg,
    surface="topographic__elevation",
    flow_metric="D8",
    update_flow_depressions=True,
    runoff_rate=None,
    depression_handler="fill",
)

The FlowAccumulator has one public methods: run_one_step().

[ ]:

fa.run_one_step()

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 or any other at_node field we choose to plot.

  • The colors 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 tutorials for more information.

[ ]:

plt.figure()
drainage_plot(mg)

In this drainage plot, we can see that all of the flow is routed down the steepest link. A plot of the drainage area would illustrate how the flow would move. Next let’s make a similar plot except that instead of plotting the topographic elevation as the background, we will plot the drainage area.

[ ]:

plt.figure()
drainage_plot(mg, "drainage_area")

If we print out the drainage area, we can see that its maximum reaches 64, which is the total area of the interior of the grid.

[ ]:

print(mg.at_node["drainage_area"].reshape(mg.shape))

This is the same number as the number of core nodes. This makes sense becaue these are the only nodes in Landlab that have area, and in our model grid they each have an area of one.

[ ]:

print(mg.number_of_core_nodes)

We can rain on the surface, store that rain in the field water__unit_flux_in, and then re-run the FlowAccumulator. As an example, we will ‘rain’ a uniformly distributed random number between 0 and 1 on every node.

Since we already ran the FlowAccumulator, under the hood our grid already has a field called water__unit_flux_in and we need to set the clobber keyword to True.

[ ]:

rain = 1.0 + 5.0 * np.random.rand(mg.number_of_nodes)

plt.imshow(rain.reshape(mg.shape), origin="lower", cmap="PuBu", vmin=0)
plt.colorbar()
plt.show()

_ = mg.add_field("water__unit_flux_in", rain, at="node", clobber=True)

Next, we re-run the FlowAccumulator and plot the discharge.

[ ]:

fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

The basic pattern of drainage is the same but the values for the surface water discharge are different than for drainage area.

Various flow routing algorithms

  • D4

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(mg, surface="topographic__elevation", flow_metric="D4")
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

  • Rho8

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(mg, surface="topographic__elevation", flow_metric="Rho8")
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

  • Rho4

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(mg, surface="topographic__elevation", flow_metric="Rho4")
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

  • Quinn

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(mg, surface="topographic__elevation", flow_metric="Quinn")
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

  • Freeman

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(
    mg, surface="topographic__elevation", flow_metric="Freeman"
)
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

  • Holmgren

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(
    mg, surface="topographic__elevation", flow_metric="Holmgren"
)
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

  • Dinf

[ ]:

mg = RasterModelGrid((10, 10))
_ = mg.add_field(
    "topographic__elevation", 3.0 * mg.x_of_node**2 + mg.y_of_node**2, at="node"
)
fa = PriorityFloodFlowRouter(mg, surface="topographic__elevation", flow_metric="Dinf")
fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")

Generated by nbsphinx from a Jupyter notebook.