Skip to content

Commit

Permalink
Flow graph flexible base levels and mask (#140)
Browse files Browse the repository at this point in the history 
* implement flexible base levels

Use an internal set for storing base level indices (in flow graph
implementation).

Refactor code that was relying (hard-coded) on
`node_status::fixed_value` for identifying base level nodes.

* implement grid mask for flow graph

Update flow operators so that masked grid node are ignored.

* add mask and base levels public API (C++)

* set_mask: test shape

* add C++ tests

* add python API

* disable python API sink algos

algo C++ API and bindings will be removed in a follow-up PR.

* add python tests

* doc: update user guide

* fix np.bool8 deprecation

* fix mypy (np.bool_ again)

* update basin ids computation

masked grid nodes have all the same assigned id (max integer limit).
  • Loading branch information
@benbovy
benbovy committed Oct 6, 2023
1 parent 26b6f6e commit f6d1a26
Show file tree
Hide file tree
Showing 20 changed files with 642 additions and 148 deletions.
27 changes: 25 additions & 2 deletions doc/source/guide_flow.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,9 +90,32 @@ These specific nodes may collect input flow but cannot release any output flow
(it could also mean that the flow is leaving the modeled domain or sub-domain,
e.g., surface land water entering the sea).

By default, all grid nodes with status
When creating a new flow graph, all grid nodes with status
{py:attr}`~fastscapelib.NodeStatus.FIXED_VALUE` are set as base
level nodes when creating a flow graph.
level nodes by default.

It is possible to set alternative base level nodes via the
{cpp:func}`~fastscapelib::flow_graph::set_base_levels` method (C++) or the
{py:attr}`~fastscapelib.FlowGraph.base_levels` property (Python). This can be
done each time prior to {ref}`updating the flow paths
<guide-compute-flow-paths>` in the cases where the base levels are evolving over
time (e.g., inner lake, sea or ocean dynamics, planet craters).

### Mask

A binary mask can be set via the {cpp:func}`~fastscapelib::flow_graph::set_mask`
method (C++) or the {py:attr}`~fastscapelib.FlowGraph.mask` property (Python) to
exclude grid nodes from the computation of the flow graph. By default, all grid
nodes are included in the flow graph (no mask).

Setting a mask is useful if the modeled domain encompass multiple sub-domains
such as land vs. ocean. Those sub-domains may each have their own flow graph
with different {ref}`flow routing strategies <guide-flow-routing-strategies>`,
{ref}`base level nodes <guide-base-level-nodes>` and masks.

Like for the base level nodes, a flow graph's mask can be set or updated each
time prior to {ref}`updating the flow paths <guide-compute-flow-paths>`, e.g.,
according to sea level change over time.

(guide-compute-flow-paths)=
### Initializing / Updating the Flow Paths
Expand Down
95 changes: 48 additions & 47 deletions include/fastscapelib/algo/pflood.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,10 +33,10 @@ namespace fastscapelib
* The main purpose of this container is for using with
* priority-flood algorithms.
*/
template <class G, class T>
template <class FG, class T>
struct pflood_node
{
using size_type = typename G::size_type;
using size_type = typename FG::size_type;

size_type m_idx;
T m_elevation;
Expand All @@ -50,21 +50,21 @@ namespace fastscapelib
{
}

bool operator>(const pflood_node<G, T>& other) const
bool operator>(const pflood_node<FG, T>& other) const
{
return m_elevation > other.m_elevation;
}
};


template <class G, class T>
using pflood_pr_queue = std::priority_queue<pflood_node<G, T>,
std::vector<pflood_node<G, T>>,
std::greater<pflood_node<G, T>>>;
template <class FG, class T>
using pflood_pr_queue = std::priority_queue<pflood_node<FG, T>,
std::vector<pflood_node<FG, T>>,
std::greater<pflood_node<FG, T>>>;


template <class G, class T>
using pflood_queue = std::queue<pflood_node<G, T>>;
template <class FG, class T>
using pflood_queue = std::queue<pflood_node<FG, T>>;


/**
Expand All @@ -73,60 +73,59 @@ namespace fastscapelib
* Add fixed value grid nodes to the priority queue and mark them as
* resolved.
*/
template <class G, class E, class elev_t = typename std::decay_t<E>::value_type>
void init_pflood(G& grid,
template <class FG, class E, class elev_t = typename std::decay_t<E>::value_type>
void init_pflood(FG& graph_impl,
E&& elevation,
xt::xtensor<bool, 1>& closed,
pflood_pr_queue<G, elev_t>& open)
pflood_pr_queue<FG, elev_t>& open)
{
using size_type = typename G::size_type;
using size_type = typename FG::size_type;

// TODO: assert elevation shape match grid shape

const auto elevation_flat = xt::flatten(elevation);

for (size_type idx = 0; idx < grid.size(); ++idx)
for (size_type idx : graph_impl.base_levels())
{
if (grid.nodes_status()[idx] == fastscapelib::node_status::fixed_value)
{
open.emplace(pflood_node<G, elev_t>(idx, elevation_flat(idx)));
closed(idx) = true;
}
open.emplace(pflood_node<FG, elev_t>(idx, elevation_flat(idx)));
closed(idx) = true;
}
}


/**
* fill_sinks_flat implementation.
*/
template <class G, class E>
void fill_sinks_flat_impl(G& grid, E&& elevation)
template <class FG, class E>
void fill_sinks_flat_impl(FG& graph_impl, E&& elevation)
{
using neighbors_indices_type = typename G::neighbors_indices_type;
using neighbors_indices_type = typename FG::grid_type::neighbors_indices_type;
using elev_t = typename std::decay_t<E>::value_type;

auto elevation_flat = xt::flatten(elevation);
neighbors_indices_type neighbors_indices;

pflood_pr_queue<G, elev_t> open;
xt::xtensor<bool, 1> closed = xt::zeros<bool>({ grid.size() });
pflood_pr_queue<FG, elev_t> open;
xt::xtensor<bool, 1> closed = xt::zeros<bool>({ graph_impl.size() });

init_pflood(grid, elevation, closed, open);
auto& grid = graph_impl.grid();

init_pflood(graph_impl, elevation, closed, open);

while (open.size() > 0)
{
pflood_node<G, elev_t> inode = open.top();
pflood_node<FG, elev_t> inode = open.top();
open.pop();

for (auto n_idx : grid.neighbors_indices(inode.m_idx, neighbors_indices))
{
if (closed(n_idx))
if (graph_impl.is_masked(n_idx) || closed(n_idx))
{
continue;
}

elevation_flat(n_idx) = std::max(elevation_flat(n_idx), inode.m_elevation);
open.emplace(pflood_node<G, elev_t>(n_idx, elevation_flat(n_idx)));
open.emplace(pflood_node<FG, elev_t>(n_idx, elevation_flat(n_idx)));
closed(n_idx) = true;
}
}
Expand All @@ -136,23 +135,25 @@ namespace fastscapelib
/**
* fill_sinks_sloped implementation.
*/
template <class G, class E>
void fill_sinks_sloped_impl(G& grid, E&& elevation)
template <class FG, class E>
void fill_sinks_sloped_impl(FG& graph_impl, E&& elevation)
{
using neighbors_indices_type = typename G::neighbors_indices_type;
using neighbors_indices_type = typename FG::grid_type::neighbors_indices_type;
using elev_t = typename std::decay_t<E>::value_type;

neighbors_indices_type neighbors_indices;

pflood_pr_queue<G, elev_t> open;
pflood_queue<G, elev_t> pit;
xt::xtensor<bool, 1> closed = xt::zeros<bool>({ grid.size() });
pflood_pr_queue<FG, elev_t> open;
pflood_queue<FG, elev_t> pit;
xt::xtensor<bool, 1> closed = xt::zeros<bool>({ graph_impl.size() });

auto& grid = graph_impl.grid();

init_pflood(grid, elevation, closed, open);
init_pflood(graph_impl, elevation, closed, open);

while (!open.empty() || !pit.empty())
{
pflood_node<G, elev_t> inode, knode;
pflood_node<FG, elev_t> inode, knode;

if (!pit.empty()
&& (open.empty() || open.top().m_elevation == pit.front().m_elevation))
Expand All @@ -171,20 +172,20 @@ namespace fastscapelib

for (auto n_idx : grid.neighbors_indices(inode.m_idx, neighbors_indices))
{
if (closed(n_idx))
if (graph_impl.is_masked(n_idx) || closed(n_idx))
{
continue;
}

if (elevation.flat(n_idx) <= elev_tiny_step)
{
elevation.flat(n_idx) = elev_tiny_step;
knode = pflood_node<G, elev_t>(n_idx, elevation.flat(n_idx));
knode = pflood_node<FG, elev_t>(n_idx, elevation.flat(n_idx));
pit.emplace(knode);
}
else
{
knode = pflood_node<G, elev_t>(n_idx, elevation.flat(n_idx));
knode = pflood_node<FG, elev_t>(n_idx, elevation.flat(n_idx));
open.emplace(knode);
}

Expand All @@ -206,13 +207,13 @@ namespace fastscapelib
* The algorithm is based on priority queues and is detailed
* in Barnes (2014). It fills sinks with flat areas.
*
* @param grid Grid object
* @param graph_impl Graph implementation object
* @param elevation Elevation values at grid nodes
*/
template <class G, class E>
void fill_sinks_flat(G& grid, xtensor_t<E>& elevation)
template <class FG, class E>
void fill_sinks_flat(FG& graph_impl, xtensor_t<E>& elevation)
{
detail::fill_sinks_flat_impl(grid, elevation.derived_cast());
detail::fill_sinks_flat_impl(graph_impl, elevation.derived_cast());
}


Expand All @@ -228,15 +229,15 @@ namespace fastscapelib
* (i.e. elevation is increased by small amount) so that there is no
* drainage singularities.
*
* @param grid Grid object
* @param graph_impl Graph implementation object
* @param elevation Elevation values at grid nodes
*
* @sa fill_sinks_flat
*/
template <class G, class E>
void fill_sinks_sloped(G& grid, xtensor_t<E>& elevation)
template <class FG, class E>
void fill_sinks_sloped(FG& graph_impl, xtensor_t<E>& elevation)
{
detail::fill_sinks_sloped_impl(grid, elevation.derived_cast());
detail::fill_sinks_sloped_impl(graph_impl, elevation.derived_cast());
}
} // namespace fastscapelib

Expand Down
18 changes: 13 additions & 5 deletions include/fastscapelib/flow/basin_graph.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,7 +314,6 @@ namespace fastscapelib
const auto& dfs_indices = m_flow_graph_impl.dfs_indices();

auto& grid = m_flow_graph_impl.grid();
const auto& nodes_status = grid.nodes_status();

// used to (re)initalize container index / position
size_type init_idx = static_cast<size_type>(-1);
Expand All @@ -339,13 +338,18 @@ namespace fastscapelib

for (const auto idfs : dfs_indices)
{
if (m_flow_graph_impl.is_masked(idfs))
{
continue;
}

const auto irec = receivers(idfs, 0);

// any new basin visited
// any new basin visited (inner or outer)
if (irec == idfs)
{
ibasin = basins(idfs);
is_inner_basin = nodes_status.flat(idfs) != node_status::fixed_value;
is_inner_basin = !m_flow_graph_impl.is_base_level(idfs);

if (!is_inner_basin)
{
Expand All @@ -368,14 +372,18 @@ namespace fastscapelib

for (auto n : grid.neighbors(idfs, neighbors))
{
if (m_flow_graph_impl.is_masked(n.idx))
{
continue;
}

const size_type nbasin = basins(n.idx);

// skip if neighbor node is in the same basin or in an
// already connected adjacent basin unless the latter is an
// outer basin
bool skip = ibasin >= nbasin;
node_status nstatus = nodes_status.flat(outlets()[nbasin]);
bool is_inner_nbasin = nstatus != node_status::fixed_value;
bool is_inner_nbasin = !m_flow_graph_impl.is_base_level(outlets()[nbasin]);
if (skip && is_inner_nbasin)
{
continue;
Expand Down
67 changes: 66 additions & 1 deletion include/fastscapelib/flow/flow_graph.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,9 @@ namespace fastscapelib
"must have at least one operator that defines the output flow direction type");
}

// initialize default base levels at fixed value nodes
m_impl.set_base_levels(m_grid.nodes_indices(node_status::fixed_value));

// pre-allocate graph and elevation snapshots
for (const auto& key : m_operators.graph_snapshot_keys())
{
Expand Down Expand Up @@ -228,6 +231,64 @@ namespace fastscapelib
return m_impl;
}

/**
* Return the indices of the base level nodes.
*/
std::vector<size_type> base_levels() const
{
const auto& impl_levels = m_impl.base_levels();
std::vector<size_type> indices(impl_levels.begin(), impl_levels.end());
return indices;
}

/**
* Clear all existing base level nodes and set new ones.
*
* @tparam C Any stl-compatible container type.
* @param levels The indices of the new base level nodes.
*/
template <class C>
void set_base_levels(C&& levels)
{
if (!m_writeable)
{
throw std::runtime_error("cannot set base levels (graph is read-only)");
}

m_impl.set_base_levels(levels);
}

/**
* Return a mask of where elements with a value ``true`` correspond
* to grid nodes that are not included in the flow graph.
*/
xt_array_t<xt_selector, bool> mask() const
{
return m_impl.mask();
}

/**
* Set a new grid mask.
*
* @tparam C Any xtensor-compatible container or expression type.
* @param mask The new mask 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)
{
if (!m_writeable)
{
throw std::runtime_error("cannot set mask (graph is read-only)");
}
if (!xt::same_shape(mask.shape(), m_grid.shape()))
{
throw std::runtime_error("cannot set mask (shape mismatch with grid shape)");
}

m_impl.set_mask(std::forward<C>(mask));
}

/**
* Traverse the flow graph in the top->down direction and accumulate
* locally produced quantities or fluxes.
Expand Down Expand Up @@ -303,7 +364,11 @@ namespace fastscapelib
*
* @return Catchment ids (array of same shape than ``grid_shape``).
*
* Note: results may be cached.
* @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.
*/
data_array_size_type basins()
{
Expand Down

0 comments on commit f6d1a26

Please sign in to comment.