Flow Routing

Flow graph

Defined in fastscapelib/flow/flow_graph.hpp

template<class G, class S = typename G::container_selector, class Tag = flow_graph_fixed_array_tag>
class flow_graph

Main class used to compute or follow flow routes on the topographic surface.

Template Parameters:
  • G – The grid type.

  • S – The xtensor container selector for data array members.

  • Tag – The flow graph implementation tag.

Public Types

using self_type = flow_graph<G, S, Tag>
using grid_type = G
using container_selector = S
using impl_type = detail::flow_graph_impl<G, S, Tag>
using operators_type = flow_operator_sequence<impl_type>
using size_type = typename grid_type::size_type
using data_type = typename grid_type::grid_data_type
using data_array_type = dynamic_shape_container_t<container_selector, data_type>
using shape_type = typename data_array_type::shape_type
using data_array_size_type = dynamic_shape_container_t<container_selector, size_type>
using graph_map = std::map<std::string, std::unique_ptr<self_type>>
using graph_impl_map = std::map<std::string, std::shared_ptr<impl_type>>
using elevation_map = std::map<std::string, std::unique_ptr<data_array_type>>

Public Functions

flow_graph(G &grid, operators_type operators)

Flow graph constructor.

Parameters:
  • grid – The grid object.

  • operators – The sequence of flow operators to use for updating the graph.

const std::vector<const flow_operator*> &operators() const

Returns a vector of pointers to the flow operators objects used in this graph.

Note: the flow_operator base class has very limited polymorphism (i.e., only its name method is virtual). It is still possible to downcast the returned pointers (e.g., using the visitor pattern).

bool single_flow() const

Returns true if the graph has single flow directions.

const std::vector<std::string> &graph_snapshot_keys() const
self_type &graph_snapshot(std::string name) const

Graph snapshot getter.

Parameters:

name – Name of the snapshot.

const std::vector<std::string> &elevation_snapshot_keys() const

Returns the names of the elevation snapshots.

const data_array_type &elevation_snapshot(std::string name) const

Elevation snapshot getter.

Parameters:

name – Name of the snapshot.

const data_array_type &update_routes(const data_array_type &elevation)

Update flow routes from the input topographic surface.

This applies in chain the flow operators and takes snapshots (if any).

Parameters:

elevation – The input topographic surface elevation.

Returns:

Either the input elevation unchanged or the elevation that has been updated in order to route flow accross closed depressions.

grid_type &grid() const

Returns a reference to the corresponding grid object.

size_type size() const

Returns the total number of nodes in the graph (equals to the size of the grid).

shape_type grid_shape() const

Returns the shape of the grid node arrays.

const impl_type &impl() const

Returns a reference to the flow graph implementation.

std::shared_ptr<impl_type> impl_ptr() const

Returns a shared pointer to the flow graph implementation.

While this might be useful in some cases (e.g., Python bindings), in general accessing the implementation through impl() (const reference) should be preferred.

std::vector<size_type> base_levels() const

Return the indices of the base level nodes.

template<class C>
void set_base_levels(C &&levels)

Clear all existing base level nodes and set new ones.

Template Parameters:

C – Any stl-compatible container type.

Parameters:

levels – The indices of the new base level nodes.

dynamic_shape_container_t<container_selector, bool> mask() const

Return a mask of where elements with a value true correspond to grid nodes that are not included in the flow graph.

template<class C>
void set_mask(C &&mask)

Set a new grid mask.

Template Parameters:

C – Any xtensor-compatible container or expression type.

Parameters:

mask – The new mask where elements with a value true correspond to grid nodes that are not included in the flow graph.

data_array_type accumulate(const data_array_type &src) const

Traverse the flow graph in the top->down direction and accumulate locally produced quantities or fluxes.

The local quantitites or fluxes (i.e., source) may for example correspond to precipitation, surface water runoff, sediment flux, etc.

The accumulated values represent at each node of the graph the integral of the source over the node upslope contributing area.

For example, if the source units are meters (height), the units of the output accumulated values are cubic meters (volume).

Parameters:

src – The source, must be given per area unit and have the same shape than grid_shape.

Returns:

The output accumulated values (array of same shape than grid_shape).

void accumulate(data_array_type &acc, const data_array_type &src) const

Traverse the flow graph in the top->down direction and accumulate locally produced quantities or fluxes.

This version allows reusing a pre-allocated array to store the output.

Both acc and src must of the same shape then grid_shape.

Parameters:
  • acc – The output accumulated values.

  • src – The source, must be given per area unit.

data_array_type accumulate(data_type src) const

Traverse the flow graph in the top->down direction and accumulate locally produced quantities or fluxes.

Parameters:

src – The spatially uniform source (scalar).

Returns:

The output accumulated values (array of same shape than grid_shape).

void accumulate(data_array_type &acc, data_type src) const

Traverse the flow graph in the top->down direction and accumulate locally produced quantities or fluxes.

Both acc and src must of the same shape then grid_shape.

Parameters:
  • acc – The output accumulated values.

  • src – The spatially uniform source (scalar).

data_array_size_type basins()

Delineate catchments (or basins).

A catchment is defined by all adjacent nodes that flow towards the same outlet (or pit) graph node.

Note

Results may be cached. All masked grid nodes have the same assigned catchment id set by the maximum limit of the integer value range. It is therefore preferable to mask the results prior to, e.g., plotting it.

Returns:

Catchment ids (array of same shape than grid_shape).

template<class FK, class FKD>
int apply_kernel(FK &kernel, FKD &data)

Apply a given kernel along the flow graph.

Visit the graph nodes in the direction and order given in the kernel object, call the kernel function and fill the output variables referenced in the kernel data.

Template Parameters:
  • FK – The flow kernel type.

  • FKD – The flow kernel data type.

Parameters:
  • kernel – The flow kernel object to apply along the graph.

  • data – The object holding or referencing input and output data used by the flow kernel.

Implementations

Defined in fastscapelib/flow/flow_graph_impl.hpp

struct flow_graph_fixed_array_tag

Fixed array flow graph implementation.

This implementation uses fixed shape arrays to store the graph topology (node receivers, donors, etc.).


Flow operators

Defined in fastscapelib/flow/flow_operator.hpp

enum class fastscapelib::flow_direction

Flow direction enum class.

Used as a tag to specify the kind of flow graph (i.e., single, multiple or undefined) expected as input or computed as output of a flow operator.

A single direction flow (directed acyclic) graph is a tree (or a forest of trees) where the total amount of matter (e.g., water, sediment) at each grid node is propagated to a unique (downslope) node neighbor. In a multiple direction flow graph, that amount is partionned among one or more node neighbors.

Use undefined for a flow operator that works on both single and multiple flow direction graphs or that do not use the graph (e.g., sink resolver only correcting the topographic elevation).

Values:

enumerator undefined
enumerator single
enumerator multi
class flow_operator

Flow operator.

It represents a logical unit that can read and/or modify in-place the flow graph and topographic elevation. It is a common type for flow routers, sink resolvers and snapshots (save intermediate graph and/or elevation states).

A flow operator may have one or more implementations, each relative to a specific flow graph implementation. Note: this class and its derived classes only contain the operators’ (static) properties and parameters (implementations are in separate classes).

Do not use this class directly (it has no implementation). Use instead its derived classes.

Subclassed by fastscapelib::flow_snapshot, fastscapelib::mst_sink_resolver, fastscapelib::multi_flow_router, fastscapelib::pflood_sink_resolver, fastscapelib::single_flow_router

Public Functions

inline virtual std::string name() const noexcept = 0

Returns the name of the operator.

Public Static Attributes

static constexpr bool elevation_updated = false

Whether the operator updates topographic elevation.

static constexpr bool graph_updated = false

Whether the operator updates the flow graph.

static constexpr flow_direction in_flowdir = flow_direction::undefined

The type of flow direction required as input of this operator.

static constexpr flow_direction out_flowdir = flow_direction::undefined

The output flow direction type of this operator.

template<class FG>
class flow_operator_sequence

Immutable container of flow operators (e.g., flow routers, sink resolvers, flow snapshots) that are applied in chain when updating a flow graph.

This class is not intended to be used as a stand-alone container. It is used internally as an entity of flow_graph and can (should) be created implicitly in the flow_graph constructor.

Template Parameters:

FG – The flow graph implementation type.

Public Functions

inline iterator_type begin()

STL-compatible iterator pointing to the 1st operator.

inline iterator_type end()

STL-compatible iterator pointing to the last operator.

inline bool elevation_updated() const

Returns true if at least one operator in the sequence updates topographic elevation.

inline bool graph_updated() const

Returns true if at least one operator in the sequence updates the flow graph.

inline flow_direction out_flowdir() const

Returns the flow direction type (single, multiple or undefined) of the final state of the flow graph, after having applied all operators.

inline bool all_single_flow() const

Returns true if all intermediate states of the flow graph have single flow directions.

inline bool snapshot_single_flow(std::string name)

Returns true if the graph snapshot given by name is single flow.

inline const std::vector<std::string> &graph_snapshot_keys() const

Returns the names of all flow graph snapshots to create.

inline const std::vector<std::string> &elevation_snapshot_keys() const

Returns the names of all topographic elevation snapshots to create.


Flow routers

Defined in fastscapelib/flow/flow_router.hpp


class single_flow_router : public fastscapelib::flow_operator

Single direction flow router operator.

This flow operator routes all the flow passing through a grid node towards its neighbors node of steepest slope.

On a raster grid with 8-node connectivity, this is equivalent to the so-called “D8” algorithm [OCallaghanM84].

Public Functions

single_flow_router() = default
inline single_flow_router(int n_threads)
inline virtual std::string name() const noexcept override

Returns the name of the operator.

inline int threads_count() const noexcept

Public Static Attributes

static constexpr bool graph_updated = true
static constexpr flow_direction out_flowdir = flow_direction::single
class multi_flow_router : public fastscapelib::flow_operator

Multiple direction flow router operator.

This flow operator partitions the flow passing through a grid node among its downslope neighbor nodes. Flow partitioning is proportional to the local slope between a node and its neighbors (power relationship with a fixed exponent parameter).

The fraction \(f_{i,j}\) of flow routed from node \(i\) to its neighbor \(j\) is given by

\[ f_{i,j} = \frac{\max (0, S_{i, j}^p)}{\sum_{k \in N} \max (0, S_{i, k}^p)} \]

where \(p\) is the slope exponent parameter, \(S_{i, j}\) is the slope between \(i\), \(j\) and \(N\) is the set of all neighbors of \(i\).

Depending on the value of the slope exponent parameter, this is equivalent to the methods described in Quinn et al. [QBCP91] or Holmgren [Hol94].

Public Functions

inline virtual std::string name() const noexcept override

Returns the name of the operator.

inline multi_flow_router(double slope_exp)

Create a new multi flow router operator.

Parameters:

slope_exp – The flow partition slope exponent.

Public Members

double m_slope_exp = 1.0

Public Static Attributes

static constexpr bool graph_updated = true
static constexpr flow_direction out_flowdir = flow_direction::multi

Sink resolvers

Defined in fastscapelib/flow/sink_resolver.hpp


struct pflood_sink_resolver : public fastscapelib::flow_operator

Priority-flood sink resolver operator.

This flow operator fills the closed depressions in the topographic surface using the priority flood algorithm +epsilon variant (Barnes et al., 2014). This variant prevents flat surfaces and hence ensure that the flow can be routed towards the outlets without disruption.

See Barnes et al. [BLM14] for a more detailed description of the algorithm.

Public Functions

inline virtual std::string name() const noexcept override

Returns the name of the operator.

Public Static Attributes

static constexpr bool elevation_updated = true
enum class fastscapelib::mst_route_method

Method used by mst_sink_resolver to route flow within each closed depressions.

The basic method is the most efficient one but does not result in a “realistic”, planar flow graph.

The carve method mimics carving a narrow canyon within the depression.

Values:

enumerator basic

Connect the pit node directly to the depression spill node.

enumerator carve

Revert the (unique) flow path between the spill and pit nodes.

class mst_sink_resolver : public fastscapelib::flow_operator

Minimum Spanning Tree (MST) sink resolver operator.

This flow operator re-routes the flow trapped in closed depressions towards their spill, using an efficient algorithm that explicitly computes a graph of (inner and outer) basins and reduces it as a tree (Cordonnier et al., 2019).

It requires a single flow graph as input.

This operator also use the updated routes in closed depressions to fill these with nearly flat surfaces (a tiny slope ensure natural flow routing for the operators applied after this one).

See Cordonnier et al. [CBB19] for a more detailed description of the algorithm.

Public Functions

inline virtual std::string name() const noexcept override

Returns the name of the operator.

inline mst_sink_resolver(mst_method basin_method = mst_method::kruskal, mst_route_method route_method = mst_route_method::carve)

Public Members

mst_method m_basin_method = mst_method::kruskal
mst_route_method m_route_method = mst_route_method::carve

Public Static Attributes

static constexpr bool graph_updated = true
static constexpr bool elevation_updated = true
static constexpr flow_direction in_flowdir = flow_direction::single
static constexpr flow_direction out_flowdir = flow_direction::single

Flow snapshots

Defined in fastscapelib/flow/flow_snapshot.hpp


class flow_snapshot : public fastscapelib::flow_operator

Flow snapshot operator.

A special flow operator used to save intermediate states of the flow graph and/or topographic elevation values while applying the other operators in chain.

Those saved states are accessible from the flow_graph object.

Public Functions

inline virtual std::string name() const noexcept override

Returns the name of the operator.

inline flow_snapshot(std::string snapshot_name, bool save_graph = true, bool save_elevation = false)
inline const std::string &snapshot_name() const noexcept
inline bool save_graph() const noexcept
inline bool save_elevation() const noexcept

Basin graph

Defined in fastscapelib/flow/basin_graph.hpp

enum class fastscapelib::mst_method

Algorithm used to compute the reduced graph (tree) of basins (minimum spanning tree).

Kruskal’s algorithm is rather simple and has non-linear complexity \(O(E \log V)\) where \(E\) is the number of basin connections and \(V\) is the number of basins.

Boruvka’s algorithm is more advanced and has linear complexity for planar graphs (most flow graphs are planar except, e.g., those built on raster grids with connectivity including the diagonals).

The gain in performance of Boruvka’s algorithm over Kruskal’s one is significant only for (very) large graphs.

Values:

enumerator kruskal

Kruskal’s algorithm.

enumerator boruvka

Boruvka’s algorithm.

template<class FG>
class basin_graph

Represents a graph of adjacent basins (catchments) connected via a two “pass” nodes where the water flow spills from one basin into another.

Such graph is built and used by mst_sink_resolver to resolve flow routing accross closed depressions in the surface topography.

Template Parameters:

FG – The flow graph implementation type (only flow_graph_fixed_array_tag is supported).

Public Types

using flow_graph_impl_type = FG
using size_type = typename flow_graph_impl_type::size_type
using data_type = typename flow_graph_impl_type::data_type
using data_array_type = typename flow_graph_impl_type::data_array_type

Public Functions

inline basin_graph(const flow_graph_impl_type &flow_graph_impl, mst_method basin_method)

Create a new graph of basins from a flow graph.

Warning

A basin graph requires a flow graph with single direction flow paths as its current state. However, no validation is done here since the state of the same flow graph may later change from single to multiple direction after the basin graph has been (re)computed.

Parameters:
  • flow_graph_impl – A flow graph implementation instance (fixed array).

  • basin_method – The algorithm used to compute the reduced tree of basins.

inline size_type basins_count() const

Returns the total number of basins.

inline const std::vector<size_type> &outlets() const

Returns a vector of node indices that represent the outlets of each of the basins in this graph.

inline const std::vector<edge> &edges() const

Returns the edges of the graph of basins.

inline const std::vector<size_type> &tree() const

Returns the reduced graph (tree) of basins after applying the selected minimum spanning tree algorithm.

void update_routes(const data_array_type &elevation)

Update (re-compute from-scratch) the graph and its reduced tree from the given topographic surface.

Parameters:

elevation – The topographic surface elevation at each grid node.

inline size_t perf_boruvka() const

Boruvka’s algorithm performance diagnostic.

Friends

friend class testing::basin_graph_orient_edges_Test
struct edge

Represents an edge of the graph of basins.

Public Functions

inline bool operator==(const edge &other)

Public Members

size_type link[2]

Indices of the two connected basins.

size_type pass[2]

Indices of the grid nodes forming the pass accross the basins.

data_type pass_elevation

Maximum elevation among the two pass nodes.

data_type pass_length

Distance between the two pass nodes.

Public Static Functions

static inline edge make_edge(const size_type &from, const size_type &to)

Create a new edge with no assigned pass yet.