Note
This page was generated from a jupyter notebook.
Using the DrainageDensity Component¶
Overview¶
Drainage density is defined as the total map-view length of stream channels, \(\Lambda\), within a region with map-view surface area \(A\), divided by that area:
The measure has dimensions of inverse length. The traditional method for measuring drainage density was to measure \(\Lambda\) on a paper map by tracing out each stream. An alternative method, which lends itself to automated calculation from digital elevation models (DEMs), is to derive drainage density from a digital map that depicts the flow-path distance from each grid node to the nearest channel node, \(L\) (Tucker et al., 2001). If the average flow-path distance to channels is \(\overline{L}\), then the corresponding average drainage density is:
An advantage of this alternative approach is that \(L\) can be mapped and analyzed statistically to reveal spatial variations, correlations with other geospatial attributes, and so on.
The DrainageDensity component is designed to calculate \(L\), and then derive \(D_d\) from it using the second equation above. Given a grid with drainage directions and drainage area, along with either a grid of channel locations or a threshold from which to generate channel locations, DrainageDensity component calculates the flow-path distance to the nearest channel node for each node in the grid. The values of \(L\) are stored in a new at-node field called
surface_to_channel__minimum_distance.
The component assumes that drainage directions and drainage area have already been calculated and the results stored in the following node fields:
flow__receiver_node: ID of the neighboring node to which each node sends flow (its “receiver”)flow__link_to_receiver_node: ID of the link along which each node sends flow to its receiverflow__upstream_node_order: downstream-to-upstream ordered array of node IDstopographic__steepest_slope: gradient from each node to its receiver
The FlowAccumulator generates all four of these fields, and should normally be run before DrainageDensity.
Identifying channels¶
The DrainageDensity component is NOT very sophisticated about identifying channels. There are (currently) two options for handling channel identification:
specify the parameters of an area-slope channelization threshold, or
map the channels separately, and pass the result to
DrainageDensityas a “channel mask” array
Area-slope channel threshold¶
This option identifies a channel as occurring at any grid node where the actual drainage area, represented by the field drainage_area, exceeds a threshold, \(T_c\):
Here \(A\) is drainage_area, \(S\) is topographic__steepest_slope, and \(C_A\), \(C_s\), \(m_r\), and \(n_r\) are parameters. For example, to create a channel mask in which nodes with a drainage area greater than \(10^5\) m\(^2\) are identified as channels, the DrainageDensity component would be initialized as:
dd = DrainageDensity(
grid,
area_coefficient=1.0,
slope_coefficient=1.0,
area_exponent=1.0,
slope_exponent=0.0,
channelization_threshold=1.0e5,
)
Channel mask¶
This option involves creating a number-of-nodes-long array, of type np.uint8, containing a 1 for channel nodes and a 0 for others.
Imports and inline docs¶
First, import what we’ll need:
[ ]:
import numpy as np
from landlab import RasterModelGrid
from landlab.components import DrainageDensity, FlowAccumulator
from landlab.io import esri_ascii
The docstring describes the component and provides some simple examples:
[ ]:
print(DrainageDensity.__doc__)
The __init__ docstring lists the parameters:
[ ]:
print(DrainageDensity.__init__.__doc__)
Example 1: channelization threshold¶
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.
[ ]:
# read the DEM
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.add_field("topographic__elevation", data, at="node")
[ ]:
grid.imshow(grid.at_node["topographic__elevation"], colorbar_label="Elevation (m)")
To use DrainageDensity, we need to have drainage directions and areas pre-calculated. We’ll do that with the FlowAccumulator component. We’ll have the component do D8 flow routing (each DEM cell drains to whichever of its 8 neighbors lies in the steepest downslope direction), and fill pits (depressions in the DEM that would otherwise block the flow) using the LakeMapperBarnes component. The latter two arguments below tell the lake mapper to update the flow directions and drainage
areas after filling the pits.
[ ]:
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
grid.imshow(
np.log10(grid.at_node["drainage_area"] + 1.0), # log10 helps show drainage
colorbar_label="Log10(drainage area (m2))",
)
Now run DrainageDensity and display the map of \(L\) values:
[ ]:
dd = DrainageDensity(
grid,
area_coefficient=1.0,
slope_coefficient=1.0,
area_exponent=1.0,
slope_exponent=0.0,
channelization_threshold=2.0e4,
)
ddens = dd.calculate_drainage_density()
grid.imshow(
grid.at_node["surface_to_channel__minimum_distance"],
cmap="viridis",
colorbar_label="Distance to channel (m)",
)
print(f"Drainage density = {ddens} m/m2")
Display the channel mask:
[ ]:
grid.imshow(grid.at_node["channel__mask"], colorbar_label="Channel present (1 = yes)")
Example 2: calculating from an independently derived channel mask¶
This example demonstrates how to run the component with an independently derived channel mask. For the sake of illustration, we will just use the channel mask from the previous example, in which case the \(L\) field should look identical.
[ ]:
# make a copy of the mask from the previous example
chanmask = grid.at_node["channel__mask"].copy()
# re-make the grid (this will remove all the previously created fields)
grid = RasterModelGrid(grid_info.shape, xy_spacing=30.0)
grid.add_field("topographic__elevation", data, at="node")
[ ]:
# 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 DrainageDensity component
dd = DrainageDensity(grid, channel__mask=chanmask)
ddens = dd.calculate_drainage_density()
# display distance-to-channel
grid.imshow(
grid.at_node["surface_to_channel__minimum_distance"],
cmap="viridis",
colorbar_label="Distance to channel (m)",
)
print(f"Drainage density = {ddens} m/m2")
References¶
Rengers, F. K., & Tucker, G. E. (2014). Analysis and modeling of gully headcut dynamics, North American high plains. Journal of Geophysical Research: Earth Surface, 119(5), 983-1003. https://doi.org/10.1002/2013JF002962
Tucker, G. E., Catani, F., Rinaldo, A., & Bras, R. L. (2001). Statistical analysis of drainage density from digital terrain data. Geomorphology, 36(3-4), 187-202, https://doi.org/10.1016/S0169-555X(00)00056-8.