SinkFiller: Permanently fill pits in a topography#

fill_sinks_barnes.py.

Fill sinks in a landscape to the brim, following the Barnes et al. (2014) algorithms.

class SinkFillerBarnes(*args, **kwds)[source]#

Bases: LakeMapperBarnes

Uses the Barnes et al (2014) algorithms to replace pits in a topography with flats, or optionally with very shallow gradient surfaces to allow continued draining.

This component is NOT intended for use iteratively as a model runs; rather, it is to fill in an initial topography. If you want to repeatedly fill pits as a landscape develops, you are after the LakeMapperBarnes component. If you want flow paths on your filled landscape, manually run a FlowDirector and FlowAccumulator for yourself.

The locations and depths etc. of the fills will be tracked, and properties are provided to access this information.

References

Required Software Citation(s) Specific to this Component

Barnes, R., Lehman, C., Mulla, D. (2014). Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers and Geosciences 62(C), 117 - 127. https://dx.doi.org/10.1016/j.cageo.2013.04.024

Additional References

None Listed

Initialise the component.

Parameters:

  • grid (ModelGrid) – A grid.

  • surface (field name at node or array of length node) – The surface to fill.

  • method ({'Steepest', 'D8'}) – Whether or not to recognise diagonals as valid flow paths, if a raster. Otherwise, no effect.

  • fill_flat (bool) – If True, pits will be filled to perfectly horizontal. If False, the new surface will be slightly inclined (at machine precision) to give steepest descent flow paths to the outlet, once they are calculated.

  • ignore_overfill (bool) – If True, suppresses the Error that would normally be raised during creation of a gentle incline on a fill surface (i.e., if not fill_flat). Typically this would happen on a synthetic DEM where more than one outlet is possible at the same elevation. If True, the was_there_overfill property can still be used to see if this has occurred.

__init__(grid, surface='topographic__elevation', method='D8', fill_flat=False, ignore_overfill=False)[source]#

Initialise the component.

Parameters:

  • grid (ModelGrid) – A grid.

  • surface (field name at node or array of length node) – The surface to fill.

  • method ({'Steepest', 'D8'}) – Whether or not to recognise diagonals as valid flow paths, if a raster. Otherwise, no effect.

  • fill_flat (bool) – If True, pits will be filled to perfectly horizontal. If False, the new surface will be slightly inclined (at machine precision) to give steepest descent flow paths to the outlet, once they are calculated.

  • ignore_overfill (bool) – If True, suppresses the Error that would normally be raised during creation of a gentle incline on a fill surface (i.e., if not fill_flat). Typically this would happen on a synthetic DEM where more than one outlet is possible at the same elevation. If True, the was_there_overfill property can still be used to see if this has occurred.

property fill_areas#

A nlakes-long array of the area of each fill.

The order is the same as that of the keys in fill_dict, and of fill_outlets.

property fill_at_node#

Return a boolean array, True if the node is filled, False otherwise.

property fill_depths#

Return the change in surface elevation at each node this step.

property fill_dict#

Return a dictionary where the keys are the outlet nodes of each filled area, and the values are deques of nodes within each.

Items are not returned in ID order.

property fill_map#

Return an array of ints, where each filled node is labelled with its outlet node ID.

Nodes not in a filled area are labelled with BAD_INDEX_VALUE (default -1).

property fill_outlets#

Returns the outlet for each filled area, not necessarily in ID order.

property fill_volumes#

A nlakes-long array of the volume of each fill.

The order is the same as that of the keys in fill_dict, and of fill_outlets.

property number_of_fills#

Return the number of individual filled areas.

run_one_step()[source]#

Fills the surface to remove all pits.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import SinkFillerBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ("left", "top", "bottom"):
...     mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
...
>>> mg.at_node["topographic__elevation"] = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 2.1, 1.1, 0.6, 1.6, 0.0],
...     [0.0, 2.0, 1.0, 0.5, 1.5, 0.0],
...     [0.0, 2.2, 1.2, 0.7, 1.7, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> z = mg.at_node["topographic__elevation"]
>>> z_init = z.copy()
>>> sfb = SinkFillerBarnes(mg, method="Steepest")  # , surface=z

TODO: once return_array_at_node is fixed, this example should also take surface… GIVE IT surface=z !!

>>> sfb.run_one_step()
>>> z_out = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 2.1, 1.5, 1.5, 1.6, 0.0],
...     [0.0, 2.0, 1.5, 1.5, 1.5, 0.0],
...     [0.0, 2.2, 1.5, 1.5, 1.7, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> np.allclose(z.reshape(mg.shape), z_out)
True
>>> fill_map = [
...     [-1, -1, -1, -1, -1, -1],
...     [-1, -1, 16, 16, -1, -1],
...     [-1, -1, 16, 16, -1, -1],
...     [-1, -1, 16, 16, -1, -1],
...     [-1, -1, -1, -1, -1, -1],
... ]
>>> np.all(sfb.fill_map.reshape(mg.shape) == fill_map)
True
>>> np.all(sfb.fill_at_node == (sfb.fill_map > -1))
True
>>> sfb.was_there_overfill  # everything fine with slope adding
False

>>> fr = FlowAccumulator(mg, flow_director="D4")  # routing now works
>>> fr.run_one_step()
>>> np.all(mg.at_node["flow__sink_flag"][mg.core_nodes] == 0)
True
>>> drainage_area = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 1.0, 2.0, 3.0, 1.0, 1.0],
...     [0.0, 1.0, 4.0, 9.0, 10.0, 10.0],
...     [0.0, 1.0, 2.0, 1.0, 1.0, 1.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> np.allclose(mg.at_node["drainage_area"].reshape(mg.shape), drainage_area)
True

Test two pits and a flat fill:

>>> from collections import deque
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1  # raise "guard" exit nodes
>>> z[7] = 2.0  # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6  # [9, 14, 15] is a lake
>>> z[22] = 0.9  # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> sfb = SinkFillerBarnes(mg, method="Steepest", fill_flat=True)
>>> sfb.run_one_step()
>>> sfb.fill_dict == {8: deque([7]), 16: deque([15, 9, 14, 22])}
True
>>> sfb.number_of_fills
2
>>> sfb.fill_outlets == [16, 8]
True
>>> np.allclose(sfb.fill_areas, np.array([4.0, 1.0]))  # same order
True

Unlike the LakeMapperBarnes equivalents, fill_depths and fill_volumes are always available through this component:

>>> fill_depths = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 1.0, 0.0, 0.5, 0.0, 0.0],
...     [0.0, 0.0, 0.4, 0.7, 0.0, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.1, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> np.allclose(sfb.fill_depths.reshape(mg.shape), fill_depths)
True
>>> np.allclose(sfb.fill_volumes, np.array([1.7, 1.0]))
True

Note that with a flat fill, we can’t drain the surface. The surface is completely flat, so each and every core node within the fill becomes a sink.

>>> fr.run_one_step()
>>> where_is_filled = np.where(sfb.fill_map > -1, 1, 0)
>>> np.all(
...     mg.at_node["flow__sink_flag"][mg.core_nodes]
...     == where_is_filled[mg.core_nodes]
... )
True

(Note that the fill_map does not think that the perimeter nodes are sinks, since they haven’t changed elevation. In contrast, the FlowAccumulator does think they are, because these nodes are where flow terminates.)

Created on Mon Oct 19.

@author: dejh

class SinkFiller(*args, **kwds)[source]#

Bases: Component

This component identifies depressions in a topographic surface, then fills them in in the topography. No attempt is made to conserve sediment mass. User may specify whether the holes should be filled to flat, or with a gradient downwards towards the depression outlet. The gradient can be spatially variable, and is chosen to not reverse any drainage directions at the perimeter of each lake.

The primary method of this class is ‘run_one_step’. ‘fill_pits’ is a synonym.

Constructor assigns a copy of the grid, and calls the initialize method.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, SinkFiller
>>> import numpy as np
>>> lake1 = np.array([34, 35, 36, 44, 45, 46, 54, 55, 56, 65, 74])
>>> lake2 = np.array([78, 87, 88])
>>> guard_nodes = np.array([23, 33, 53, 63, 73, 83])
>>> lake = np.concatenate((lake1, lake2))
>>> mg = RasterModelGrid((10, 10))
>>> z = np.ones(100, dtype=float)
>>> z += mg.node_x  # add a slope
>>> z[guard_nodes] += 0.001  # forces the flow out of a particular node
>>> z[lake] = 0.0
>>> field = mg.add_field(
...     "topographic__elevation",
...     z,
...     at="node",
...     units="-",
...     copy=True,
... )
>>> fr = FlowAccumulator(mg, flow_director="D8")
>>> fr.run_one_step()
>>> mg.at_node["flow__sink_flag"][mg.core_nodes].sum()
14
>>> hf = SinkFiller(mg, apply_slope=False)
>>> hf.run_one_step()
>>> np.allclose(mg.at_node["topographic__elevation"][lake1], 4.0)
True
>>> np.allclose(mg.at_node["topographic__elevation"][lake2], 7.0)
True

Now reset and demonstrate the adding of an inclined surface:

>>> field[:] = z
>>> hf = SinkFiller(mg, apply_slope=True)
>>> hf.run_one_step()
>>> hole1 = np.array(
...     [
...         4.00007692,
...         4.00015385,
...         4.00023077,
...         4.00030769,
...         4.00038462,
...         4.00046154,
...         4.00053846,
...         4.00061538,
...         4.00069231,
...         4.00076923,
...         4.00084615,
...     ]
... )
>>> hole2 = np.array([7.4, 7.2, 7.6])
>>> np.allclose(mg.at_node["topographic__elevation"][lake1], hole1)
True
>>> np.allclose(mg.at_node["topographic__elevation"][lake2], hole2)
True
>>> fr.run_one_step()
>>> mg.at_node["flow__sink_flag"][mg.core_nodes].sum()
0

References

Required Software Citation(s) Specific to this Component

None Listed

Additional References

Tucker, G., Lancaster, S., Gasparini, N., Bras, R., Rybarczyk, S. (2001). An object-oriented framework for distributed hydrologic and geomorphic modeling using triangulated irregular networks. Computers & Geosciences 27(8), 959-973. https://dx.doi.org/10.1016/s0098-3004(00)00134-5

Parameters:

  • grid (ModelGrid) – A landlab grid.

  • routing ({'D8', 'D4'} (optional)) – If grid is a raster type, controls whether fill connectivity can occur on diagonals (‘D8’, default), or only orthogonally (‘D4’). Has no effect if grid is not a raster.

  • apply_slope (bool) – If False (default), leave the top of the filled sink flat. If True, apply the slope fill_slope to the top surface to allow subsequent flow routing. A test is performed to ensure applying this slope will not alter the drainage structure at the edge of the filled region (i.e., that we are not accidentally reversing the flow direction far from the outlet.)

  • fill_slope (float (m/m)) – The slope added to the top surface of filled pits to allow flow routing across them, if apply_slope.

__init__(grid, routing='D8', apply_slope=False, fill_slope=1e-05)[source]#
Parameters:

  • grid (ModelGrid) – A landlab grid.

  • routing ({'D8', 'D4'} (optional)) – If grid is a raster type, controls whether fill connectivity can occur on diagonals (‘D8’, default), or only orthogonally (‘D4’). Has no effect if grid is not a raster.

  • apply_slope (bool) – If False (default), leave the top of the filled sink flat. If True, apply the slope fill_slope to the top surface to allow subsequent flow routing. A test is performed to ensure applying this slope will not alter the drainage structure at the edge of the filled region (i.e., that we are not accidentally reversing the flow direction far from the outlet.)

  • fill_slope (float (m/m)) – The slope added to the top surface of filled pits to allow flow routing across them, if apply_slope.

drainage_directions_change(lake_nodes, old_elevs, new_elevs)[source]#

True if the drainage structure at lake margin changes, False otherwise.

fill_pits()[source]#

This is a synonym for the main method run_one_step.

run_one_step()[source]#

This is the main method.

Call it to fill depressions in a starting topography.