How to do “D4” pit-filling on a digital elevation model (DEM)¶

(Greg Tucker, July 2021)

Digital elevation models (DEMs) often contain “pits”: cells or groups of cells that represent depressions, and are surrounded by higher-elevation cells. Sometimes these pits represent real landscape features, such as sinkholes in a karst landscape, or “pools” along a dry river channel. In other cases, they are artifacts of errors in the original data or of data processing. Sometimes they reflect a truncated representation of elevation values. For example, some DEMs use integer values to represent elevation in meters. Adjacent cells that differ by less than a meter in elevation can register as having exactly the same elevation in an integer DEM, making the drainage pathways in ambiguous.

For some overland flow models, depressions in a DEM are no problem: water will simply pool in the depressions and overspill when the depression is full. But overland flow models that use the kinematic wave approximation, in which the ground-surface slope is used as a proxy for the hydraulic energy slope, can not do this. Therefore, 2D overland flow models, like the Landlab components KinwaveOverlandFlowModel and KinwaveImplicitOverlandFlow, require a DEM that has had its pits and flat areas digitally removed before they can be run reliably.

This tutorial describes how to do this using Landlab’s flow-routing and pit-filling components.

First, some imports:

[ ]:

from landlab.components import FlowAccumulator
from landlab.io import esri_ascii

Read the un-filled DEM, which happens to be in Arc/Info ASCII Grid format (a.k.a., ESRI ASCII). We will use the set_watershed_boundary_condition function to set all nodes with an elevation value equal to a “no data” code (default -9999) to closed boundaries, and any nodes with valid elevation values that lie on the grid’s perimeter to open (fixed value) boundaries. You can learn more about the raster version of this handy function here.

[ ]:

# Read a small DEM
with open("hugo_site.asc") as fp:
    grid = esri_ascii.load(fp, name="topographic__elevation", at="node")

# This sets nodes with a no-data code (default -9999) to closed boundary status
# (For a perimeter node to be considered )
grid.set_watershed_boundary_condition(grid.at_node["topographic__elevation"])

Having read the DEM and set its boundary conditions, we now instantiate and run FlowAccumulator. We will tell it to use “D4” routing, which is indicated by the parameter choice flow_director = FlowDirectorSteepest. (Note that the term “steepest” means “choose the steepest drop among all directly adjacent nodes”. In a RasterModelGrid, there are normally four neighboring nodes, and so when FlowDirectorSteepest is applied to a raster grid, it does ‘D4’ routing. When it is applied to a hex grid, for example, it chooses among six possible directions. This variation in the number of potential directions depending on grid type is why it this particular flow-direction method is called “steepest” rather than “D4”. For more on this, see the tutorials the_FlowDirectors and compare_FlowDirectors.)

First, we will run the flow accumulator without any depression handling, in order to visualize the drainage in this case:

[ ]:

fa_no_fill = FlowAccumulator(grid, flow_director="FlowDirectorSteepest")
fa_no_fill.run_one_step()

[ ]:

grid.imshow("drainage_area", colorbar_label="Drainage area (m2)")

Notice that there are interruptions in the drainage patterns: places where drainage area abruptly drops in the downstream direction. These spots represent drainages that terminate in pits.

Now we’ll repeat the process, but this time we will also tell the FlowAccumulator to use the LakeMapperBarnes component for depression handling. (Note that we could run LakeMapperBarnes separately, but because depression handling and flow accumulation work together, the FlowAccumulator provides a way to directly “embed” a depression handler and its arguments, which is the approach used here; you can learn more about this in the tutorial the_FlowAccumulator.) We can examine the parameters of LakeMapperBarnes:

[ ]:

from landlab.components import LakeMapperBarnes

print(LakeMapperBarnes.__init__.__doc__)

If we don’t wish to accept the built-in defaults, we can send any of these parameters as additional arguments to the FlowAccumulator constructor. For this application, we want:

surface = topographic__elevation: use the topography (this is the default)
method = 'Steepest': use D4 routing (the default)
fill_flat = False: we want a slight slope assigned to otherwise flat areas
fill_surface = 'topographic__elevation': fill the topography (the default)
redirect_flow_steepest_descent = True: so we can plot the revised drainage
reaccumulate_flow = True: so we can plot the revised drainage

[ ]:

fa = FlowAccumulator(
    grid,
    flow_director="FlowDirectorSteepest",
    depression_finder="LakeMapperBarnes",
    surface="topographic__elevation",
    method="Steepest",
    fill_flat=False,
    fill_surface="topographic__elevation",
    redirect_flow_steepest_descent=True,
    reaccumulate_flow=True,
)

[ ]:

original_elev = grid.at_node[
    "topographic__elevation"
].copy()  # keep a copy of the original
fa.run_one_step()
grid.imshow("drainage_area", colorbar_label="Drainage area (m2)")

Now we have a continuous drainage system that passes smoothly over the now-filled pits. Drainage from every cell in the watershed can reach the outlet. By plotting the difference between the original and modified topography, we can inspect the depth and spatial patterns of pit filling:

[ ]:

grid.imshow(
    grid.at_node["topographic__elevation"] - original_elev,
    colorbar_label="fill depth (m)",
)

Generated by nbsphinx from a Jupyter notebook.