<jupyter_start><jupyter_text>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.<jupyter_code># 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
from landlab.components import PriorityFloodFlowRouter

# 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>Topographic gridsFor this tutorial we will consider one topographic surface. Here it is plotted in three dimensions.<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Initalizing and running the FlowAccumulatorTo 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.htmlSome 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.<jupyter_code>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",
)<jupyter_output><empty_output><jupyter_text>The **FlowAccumulator** has one public methods: ```run_one_step()```.<jupyter_code>fa.run_one_step()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>plt.figure()
drainage_plot(mg)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>plt.figure()
drainage_plot(mg, "drainage_area")<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>print(mg.at_node["drainage_area"].reshape(mg.shape))<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>print(mg.number_of_core_nodes)<jupyter_output><empty_output><jupyter_text>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```.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Next, we re-run the **FlowAccumulator** and plot the discharge.<jupyter_code>fa.run_one_step()
plt.figure()
drainage_plot(mg, "surface_water__discharge", title="Discharge")<jupyter_output><empty_output><jupyter_text>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<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>* Rho8<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>* Rho4<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>* Quinn<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>* Freeman<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>* Holmgren<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>* Dinf<jupyter_code>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")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09334971047552203, 'fraction_numerical': 0.01052816283558519, 'mean_word_length': 5.6544074722708695, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 1, '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': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '281165', 'license': 'MIT', 'n_tokens_mistral': 3662, 'n_tokens_neox': 3430, 'n_words': 1181, 'path': 'notebooks/tutorials/flow_direction_and_accumulation/the_Flow_Director_Accumulator_PriorityFlood.ipynb', 'repo_name': 'Joseph-A-Mason/landlab', 'size': '11398'}