-
Notifications
You must be signed in to change notification settings - Fork 6
Design notes: flow routing
There are many ways to compute flow routing from digital elevation data. Each algorithm differs by:
- The representation of the flow paths, i.e., convergent (D8, D-Inf) vs divergent (MDF) flow.
- The grid(s) supported or for which the algorithm is optimal
- Parallel vs. sequential execution
A sub-problem to flow routing is resolving the flow paths that are trapped in local depressions (sinks) of the topographic surface, which don't always represent real features on the surface but could rather be artifacts, e.g., resulting from the generation of digital elevation data. Again here, there exists a variety of strategies/algorithms/variants/combinations.
Fastscapelib aims to expose through a common API various, pugglable options for the two problems described above. Ideally, it should be extensible and allow plugin-in custom options implemented in third-party libraries or applications.
The proposed design is based on one main class flow_graph and two entity semantic base classes flow_router and sink_resolver. These three classes are each detailled below.
Flow paths computed from a DEM can be represented as a directed acyclic graph that connects each grid node to one or several of its downslope, flow receiver neighbor(s). This class holds the graph data and provides the API for graph initialization, update and traversal. Its (incomplete) declaration would look like:
/**
* Main class used to compute or follow flow routes on
* the topographic surface.
*
* @tparam G The grid type.
* @tparam elev_t Type used to store elevation values.
*/
template <class G, class elev_t>
class flow_graph
{
public:
using self_type = flow_graph<G, elev_t>;
using elevation_type = xt_container_t<typename G::xt_selector, elev_t, G::xt_ndims>;
using flow_router_ptr = std::unique_ptr<flow_router<self_type>>;
using sink_resolver_ptr = std::unique_ptr<sink_resolver<self_type>>;
flow_graph(const G& grid,
flow_router_ptr f_router,
sink_resolver_ptr s_resolver)
: p_flow_router(f_router), p_sink_resolver(s_resolver)
{
}
void update(const elevation_type& elevation)
{
p_flow_router->route(elevation, *this, *p_sink_resolver);
}
private:
flow_router_ptr p_flow_router;
sink_resolver_ptr p_sink_resolver;
};A number of types, like the xtensor container types used to hold graph data, can be retrieved directly from the type of grid used (G). See PR#42.
flow_router subclasses provide an implementation for the computation of flow paths. The route method defined in the base class is called from flow_graph. This method also calls sink resolving methods (see below). The declaration of the base class would look like:
/**
* Base class for the implementation of flow routing
* methods.
*
* All derived classes must implement ``route_impl``.
*
* @tparam FG The flow_graph class.
*/
template <class FG>
class flow_router
{
public:
using elevation_type = typename FG::elevation_type;
// Entity semantic
virtual ~flow_router() = default;
flow_router(const flow_router&) = delete;
flow_router(flow_router&&) = delete;
flow_router& operator=(const flow_router&) = delete;
flow_router& operator=(flow_router&&) = delete;
void route(const elevation_type& elevation,
FG& fgraph,
sink_resolver<FG>& sresolver)
{
sresolver.resolve_before_route(elevation, fgraph);
route_impl(elevation, fgraph);
sresolver.resolve_after_route(elevation, fgraph);
}
private:
virtual void route_impl(const elevation_type& elevation, FG& fgraph) = 0;
protected:
flow_router() = default;
};An example of declaration of a derived class for straightforward flow routing on a 1-d channel grid would be like:
template <class FG>
class one_channel_router : public flow_router<FG>
{
public:
using base_type = flow_router<FG>;
using elevation_type = typename base_type::elevation_type;
one_channel_router() = default;
virtual ~one_channel_router() = default;
private:
void route_impl(const elevation_type& elevation, FG& fgraph) override;
};To support both D4 and D8 cases on raster grids, we could add a template parameter to the derived class that implements single flow routing (this parameter would be ignored for non-raster grids). For example:
template <class FG, raster_connect RC = raster_connect::queen>
class single_flow_router : public flow_router<FG>
{
public:
using base_type = flow_router<FG>;
using elevation_type = typename base_type::elevation_type;
single_flow_router() = default;
virtual ~single_flow_router() = default;
private:
void route_impl(const elevation_type& elevation, FG& flow_graph) override;
};Some routing methods like MFD have parameters that we need to add in the constructor of the derived class, e.g.,
template <class FG, raster_connect RC = raster_connect::queen>
class multiple_flow_router : public flow_router<FG>
{
public:
using base_type = flow_router<FG>;
using elevation_type = typename base_type::elevation_type;
multiple_flow_router(double slope_exp, double dist_exp);
virtual ~multiple_flow_router() = default;
private:
void route_impl(const elevation_type& elevation, FG& flow_graph) override;
};sink_resolver subclasses provide an implementation of a given sink resolving method before and/or after computing the flow paths. Some methods like priority-flood variants consists of updating (a copy) of the DEM values -- e.g., by filling/carving the sinks -- before computing the flow paths. Other methods like this one update the flow paths within the sinks after their computation. Advanced strategies might do some computation before and after computing the flow paths. Therefore, sink_resolver subclasses must provide an implementation (that could just be void) for the two corresponding methods. The declaration of the base class would look like:
template <class FG>
class sink_resolver
{
public:
using elevation_type = typename FG::elevation_type;
// Entity semantic
virtual ~sink_resolver() = default;
sink_resolver(const sink_resolver&) = delete;
sink_resolver(sink_resolver&&) = delete;
sink_resolver& operator=(const sink_resolver&) = delete;
sink_resolver& operator=(sink_resolver&&) = delete;
virtual void resolve_before_route(const elevation_type& elevation, FG& fgraph) = 0;
virtual void resolve_after_route(const elevation_type& elevation, FG& fgraph) = 0;
protected:
sink_resolver() = default;
};The most basic derived class would correspond to the case where we don't want to do anything special with sinks, e.g.,
template <class FG>
class dummy_resolver : public sink_resolver<FG>
{
public:
using base_type = sink_resolver<FG>;
using elevation_type = typename base_type::elevation_type;
dummy_resolver() = default;
virtual ~dummy_resolver() = default;
void resolve_before_route(const elevation_type& elevation, FG& fgraph) override
{
}
void resolve_after_route(const elevation_type& elevation, FG& fgraph) override
{
}
};Another derived class implementing a priority flood algorithm would let resolve_after_route empty. We would also need the additional template parameter for raster connectivity here:
template <class FG, raster_connect RC>
class pflood_resolver : public sink_resolver<FG>
{
public:
using base_type = sink_resolver<FG>;
using elevation_type = typename base_type::elevation_type;
pflood_resolver() = default;
virtual ~pflood_resolver() = default;
void resolve_before_route(const elevation_type& elevation, FG& fgraph) override
{
// run priority flood
}
void resolve_after_route(const elevation_type& elevation, FG& fgraph) override
{
}
};Note that for the later case we still need to figure out how the updated DEM will be passed and used in flow_router.route_impl().
Most of the time, users would build flow_graph objects with a combination of flow_router / sink_resolver derived classes that are chosen among the few ones already available in the library. It would be a bit cumbersome to manually create pointers to those classes and pass them to the flow_graph constructor above. For convenience, we could provide some API sugar using enums, makers and overloaded constructors.
Example of a selector for sink resolving algorithms:
enum class sink_resolve_algo
{
none = 0,
fill_pflood,
fill_mst_kruskal,
fill_mst_boruvka,
fill_auto, // either pflood or mst depending on number of sinks
carve_mst_kruskal,
carve_mst_boruvka,
};The selector for flow routing algorithms is a bit more complex as we might need to pass some parameters depending on the algorithm. A solution might be:
enum class flow_route_algo
{
one_channel = 0,
single,
multiple, // there are here parameters (floating point exponents)
single_parallel
};
struct flow_route_method
{
flow_route_algo algo;
double param1;
double param2;
flow_route_method(flow_route_algo algo, double p1 = 0., double p2 = 0.)
: algo(algo), param1(p1), param2(p2)
{
}
};The default values set for the parameters p1 and p2 allow using flow_route_algo directly (implicit flow_route_method constructor) when those parameters are not required.
We then need something that returns from those selectors new pointers to the corresponding flow_router / sink_resolver derived classes. The tricky part is to handle the raster connectivity template parameter in both cases. A solution could look like, e.g., for sink_resolver:
template<class FG, raster_connect RC>
sink_resolver<FG>* make_sink_resolver_impl(sink_resolve_algo sink_algo)
{
switch(sink_algo)
{
case sink_resolve_algo::none:
default:
return new dummy_resolver<FG>();
case sink_resolve_algo::fill_pflood:
return new pflood_resolver<FG, RC>();
}
}
template<class FG>
sink_resolver<FG>* make_sink_resolver(sink_resolve_algo sink_algo,
raster_connect rconnect)
{
switch(rconnect)
{
case raster_connect::queen:
default:
return make_sink_resolver_impl<FG, raster_connect::queen>(sink_algo);
case raster_connect::rook:
return make_sink_resolver_impl<FG, raster_connect::rook>(sink_algo);
}
}The same applies for flow_router (with additional handling of algorithm parameters):
template<class FG, raster_connect RC>
flow_router<FG>* make_flow_router_impl(flow_route_method route_method)
{
switch(route_method.algo)
{
case flow_route_algo::one_channel:
return new one_channel_router<FG>();
case flow_route_algo::single:
default:
return new single_flow_router<FG, RC>();
case flow_route_algo::multiple:
return new multiple_flow_router<FG, RC>(fmethod.param1, fmethod.param2);
}
}
template<class FG>
flow_router<FG>* make_flow_router(flow_route_method route_method,
raster_connect rconnect)
{
switch(rconnect)
{
case raster_connect::queen:
default:
return make_flow_router_impl<FG, raster_connectivity::queen>(route_method);
case raster_connect::rook:
return make_flow_router_impl<FG, raster_connectivity::rook>(route_method);
}
}Finally, we add a flow_graph constructor for easy selection of built-in methods:
template <class G, class elev_t>
flow_graph::flow_graph(const G& grid,
flow_route_method route_method,
sink_resolve_algo sink_algo = sink_resolve_algo::none,
raster_connect rconnect = raster_connect::queen)
: p_flow_router(make_flow_router<self_type>(route_method, rconnect)),
p_sink_resolver(make_sink_resolver<self_type>(sink_algo, rconnect))
{
}Additional constructors with SFINAE might be useful in cases where some combination of flow_router / sink_resolver don't make sense or where we don't need to choose a flow routing method, like the straightforward case of 1-d channel grids:
template <class G, class elev_t>
template<class PG = G, std::enable_if_t<std::is_same<PG, profile_grid_xt<PG::xt_selector>>::value, int> = 0>
flow_graph::flow_graph(const PG& channel_grid,
sink_resolve_algo sink_algo = sink_resolve_algo::none)
: flow_graph(channel_grid, flow_route_algo::one_channel, sink_algo)
{
// ...
}- Maybe it would be more consistent to move the
elev_ttemplate parameter fromflow_graphto grids and re-use it inflow_graph<G>? This would also make sense if we later need it in other classes that are independent offlow_graph, in a totally different context.