Note

This page was generated from a jupyter notebook.

Using the HackCalculator Component

Background

Hack’s Law refers to a commonly observed power-law relationship between the length of the longest stream in a drainage basin, \(L\), and the area that it drains, \(A\):

\[L = C A^h\]

where \(h\) is commonly in the neighborhood of 0.5 or 0.6. It is named for the American geomorphologist John Hack (see Hack, 1957). A value of 0.5 represents perfect geometry similarity between the areal measure, \(A\), and the embedded length \(L\). Montgomery and Dietrich (1992) noted a useful “rule of thumb” empirical relationship:

\[A \approx \frac{1}{3} L^2\]

which says that, roughly speaking, the surface area of a typical drainage basin is on the order of 1/3 times the square of its longest stream length. But individual drainage basins often deviate somewhat from this idealized behavior, and often it is useful to extract the values of \(C\) and \(h\) for particular basins.

Component overview

The HackCalculator provides estimates of the best-fit \(C\) and \(h\) values for one or more drainage basins within a given area represented by a digital elevation model, which is contained in a grid field called topographic__elevation. The component requires a grid that includes the following fields:

  • topographic__elevation (node) elevation data for the region of interest

  • drainage_area (node) contributing drainage area for each node

  • flow__receiver_node (node) ID of the receiver (node that receives flow) for each node

  • flow__link_to_receiver_node (node) ID of the link connecting node and receiver for each node

  • flow__upstream_node_order (node) array of node IDs in downstream-to-upstream order

Apart from topographic__elevation, each of these fields are created and calculated by running the FlowAccumulator component.

The component uses the ChannelProfiler component to calculate cumulative downstream length for one or more channels in the terrain. Parameters for the ChannelProfiler may be given as arguments to HackCalculator.

When run (using the method calculate_hack_parameters), the component creates a Pandas DataFrame called hack_coefficient_dataframe. It is a pandas dataframe with one row for each basin for which Hack parameters are calculated. Thus, there are as many rows as the number of watersheds identified by the ChannelProfiler.

The dataframe has the following index and columns:

  • Index

    • basin_outlet_id: The node ID of the watershed outlet where each set of Hack parameters was estimated.

  • Columns

    • A_max: The drainage area of the watershed outlet.

    • C: The Hack coefficient as defined in the equations above.

    • h: The Hack exponent as defined in the equations above.

If you pass the argument save_full_df=True, HackCalculator will generate an additional DataFrame called full_hack_dataframe. It is pandas dataframe with a row for every model grid cell used to estimate the Hack parameters. It has the following index and columns:

  • Index

    • node_id: The node ID of the model grid cell.

  • Columns

    • basin_outlet_id: The node IDs of watershed outlet

    • A: The drainage are of the model grid cell.

    • L_obs: The observed distance to the divide.

    • L_est: The predicted distance to divide based on the Hack coefficient fit.

Imports and inline docs

First, import what we’ll need:

[ ]:

from landlab import RasterModelGrid
from landlab.components import FlowAccumulator, HackCalculator
from landlab.io import esri_ascii

The docstring describes the component and provides some simple examples:

[ ]:

print(HackCalculator.__doc__)

The __init__ docstring lists the parameters:

[ ]:

print(HackCalculator.__init__.__doc__)

Example 1

In this example, we read in a small digital elevation model (DEM) from NASADEM for an area on the Colorado high plains (USA) that includes a portion of an escarpment along the west side of a drainage known as West Bijou Creek (see Rengers & Tucker, 2014).

The DEM file is in ESRI Ascii format, but is in a geographic projection, with horizontal units of decimal degrees. To calculate slope gradients properly, we’ll first read the DEM into a Landlab grid object that has this geographic projection. Then we’ll create a second grid with 30 m cell spacing (approximately equal to the NASADEM’s resolution), and copy the elevation field from the geographic DEM. This isn’t a proper projection of course, but it will do for purposes of this example.

We use the lazy_load function rather than load function because we just want to just read the metadata from the data file needed to create a new grid. We will use this metadata to create a new grid for each of the examples below.

[ ]:

with open("west_bijou_escarpment_snippet.asc") as fp:
    grid_info, data = esri_ascii.lazy_load(fp, name="topographic__elevation", at="node")

[ ]:

grid = RasterModelGrid(grid_info.shape, xy_spacing=30.0)
grid.at_node["topographic__elevation"] = data

[ ]:

grid.imshow("topographic__elevation", colorbar_label="Elevation (m)")

Next, we run the HackCalculator on this DEM. First we need to instantiate and run FlowAccumulator to calculate flow directions and drainage area.

[ ]:

# instatiated and run flow accumulator
fa = FlowAccumulator(
    grid,
    flow_director="FlowDirectorD8",  # use D8 routing
    depression_finder="LakeMapperBarnes",  # pit filler
    method="D8",  # pit filler use D8 too
    redirect_flow_steepest_descent=True,  # re-calculate flow dirs
    reaccumulate_flow=True,  # re-calculate drainagea area
)
fa.run_one_step()  # run the flow accumulator

# instantiate and run HackCalculator component
hc = HackCalculator(grid)
hc.calculate_hack_parameters()

Here is the resulting DataFrame, containing the area of the largest drainage basin and its corresponding \(C\) and \(h\) values:

[ ]:

hc.hack_coefficient_dataframe

You can access the embedded ChannelProfiler via HackCalculator.profiler:

The ChannelProfiler data is an ordered dict, in this case containing data for one watershed: the one that drains to node 6576 (for details see the reference documentation and tutorial resources for ChannelProfiler).

For this example, we might wish to visualize the main channel for which the Hack coefficient and exponent were calculated. We can do that with the profiler’s plot_profiles_in_map_view method:

[ ]:

hc.profiler.plot_profiles_in_map_view(colorbar_label="Elevation (m)")

The component also provides a node field called distance_to_divide that, as the name implies, contains the streamwise distance between a node and its source at a drainage divide:

[ ]:

grid.imshow("distance_to_divide", colorbar_label="Distance from drainage divide (m)")

Example 2: full data frame

The next example is the same as the first, but here we request and examine the “full dataframe”:

[ ]:

# create grid and copy DEM into it
grid = RasterModelGrid(grid_info.shape, xy_spacing=30.0)
grid.at_node["topographic__elevation"] = data

# instatiated and run flow accumulator
fa = FlowAccumulator(
    grid,
    flow_director="FlowDirectorD8",  # use D8 routing
    depression_finder="LakeMapperBarnes",  # pit filler
    method="D8",  # pit filler use D8 too
    redirect_flow_steepest_descent=True,  # re-calculate flow dirs
    reaccumulate_flow=True,  # re-calculate drainagea area
)
fa.run_one_step()  # run the flow accumulator

# instantiate and run HackCalculator component
hc = HackCalculator(grid, save_full_df=True)
hc.calculate_hack_parameters()

[ ]:

hc.hack_coefficient_dataframe

[ ]:

hc.full_hack_dataframe

Example 3: multiple watersheds

By default, the ChannelProfiler extracts data from just one watershed, which is why the above example reports Hack parameters for just one basin. Here we re-run the analysis with five basins.

[ ]:

# create grid and copy DEM into it
grid = RasterModelGrid(grid_info.shape, xy_spacing=30.0)
grid.at_node["topographic__elevation"] = data

# instatiated and run flow accumulator
fa = FlowAccumulator(
    grid,
    flow_director="FlowDirectorD8",  # use D8 routing
    depression_finder="LakeMapperBarnes",  # pit filler
    method="D8",  # pit filler use D8 too
    redirect_flow_steepest_descent=True,  # re-calculate flow dirs
    reaccumulate_flow=True,  # re-calculate drainagea area
)
fa.run_one_step()  # run the flow accumulator

# instantiate and run HackCalculator component
hc = HackCalculator(grid, number_of_watersheds=5)
hc.calculate_hack_parameters()

[ ]:

hc.hack_coefficient_dataframe

[ ]:

hc.profiler.plot_profiles_in_map_view(colorbar_label="Elevation (m)")

Example 4: multiple channels per basin

So far, we have only performed the calculation on the main channel in each drainage basin. We can operate on all the channels in each basin by setting the ChannelProfiler parameter main_channel_only to False. While we’re at it, we will also specify a drainage area threshold for channels of 20,000 m\(^2\).

[ ]:

# create grid and copy DEM into it
grid = RasterModelGrid(grid_info.shape, xy_spacing=30.0)
grid.at_node["topographic__elevation"] = data

# instatiated and run flow accumulator
fa = FlowAccumulator(
    grid,
    flow_director="FlowDirectorD8",  # use D8 routing
    depression_finder="LakeMapperBarnes",  # pit filler
    method="D8",  # pit filler use D8 too
    redirect_flow_steepest_descent=True,  # re-calculate flow dirs
    reaccumulate_flow=True,  # re-calculate drainagea area
)
fa.run_one_step()  # run the flow accumulator

# instantiate and run HackCalculator component
hc = HackCalculator(
    grid,
    number_of_watersheds=5,
    main_channel_only=False,
    minimum_channel_threshold=2.0e4,
)
hc.calculate_hack_parameters()

[ ]:

hc.hack_coefficient_dataframe

[ ]:

hc.profiler.plot_profiles_in_map_view(colorbar_label="Elevation (m)")

References

Hack, J. T. (1957). Studies of longitudinal stream profiles in Virginia and Maryland. Geological Survey Professional Paper 294-B. US Government Printing Office.

Montgomery, D. R., & Dietrich, W. E. (1992). Channel initiation and the problem of landscape scale. Science, 255(5046), 826-830.

Generated by nbsphinx from a Jupyter notebook.