Implement D-infinity flow routing

[wkearn]

This PR implements the D-infinity flow routing algorithm of Tarboton (1997). D-infinity routes flow to up to two downstream neighbors of each pixel, with flow proportioned according to the angle of the normal vector of 8 triangular facets surrounding the center pixel.

This PR also introduces a major redesign of the flow routing algorithms to accommodate better multiple flow direction algorithms and to separate flow routing from flat resolution. This design is used in the D-infinity implementation and several helper functions are added. The existing D8 implementation has not yet been modified.

The new flow routing procedure is:

1. Use a flow routing algorithm to determine flow directions over the DEM. These are stored in a bit array of type uint8_t where bit k is set in pixel (i, j) if there is an outgoing edge from that pixel in the kth direction. The directions are numbered starting from 1 to the right and moving counterclockwise.

Any unresolved pixels have a direction of 0. These could be missing data in the DEM, boundary pixels or sinks.

2. Count the number of edges in the flow direction bitmap. We provide a function edgeset_count to do this, but it could also be done externally. It is the number of set bits in the direction bitmap.

3. Allocate a weight array of type float containing the counted number of edges. Fill it up according to the weighting scheme of the flow routing algorithm. These weights are packed into the array to minimize space usage. This is why the number of edges need to be counted prior to computing the weights. The flow_routing_dinf_weights function computes and stores them in the correct order.

4. If you don't want to resolve flats, skip to step 7.

5. Apply a flat resolution algorithm to create a new direction bitmap and weight array for the unresolved pixels. This process looks a lot like the regular flow routing: fill up the bitmap, count the edges, allocate and fill up the weights.

6. Count the total number of edges in the originally resolved direction bitmap and the bitmap of directions on
   flats. edgeset_count_merged is provided for this purpose. Allocate a weight array and merge the originally resolved directions and the resolved directions on flats using the merge_edges function.

You have to do this after flat resolution because of the way the edges are sorted in the weight array. There might be a better way to do it that would allow you to skip this step and merge on the fly.

7. Now you have a bitmap of directions and float array of weights for all the edges you want in your Flow object. Use the edgeset_scan function to prepare an array of type ptrdiff_t with the offsets of each pixels edges in the weight array. This is done for you by edgeset_merge if you used that. edgeset_count, edgeset_scan and edgeset_merge all return the total number of edges in the edge set.

8. Allocate your Flow object arrays:

• ptrdiff_t *stream: the topologically sorted list of nodes. Size is the number of pixels in your DEM
• ptrdiff_t *source: the topologically sorted list of edge sources. Size is the number of edges
• ptrdiff_t *target: the topologically sorted list of edge targets. Size is the number of edges.
• float *weight: the topologically sorted list of edge weights. Size is the number of edges.

9. Use direction_tsort to traverse the flow network and topologically sort the edges. The node, source, target and weight arrays are filled up in the right order.

The bitmap edge set provides a uniform interface for flow routing and flat resolution algorithms to exchange data. This could also be done with topologically sorted edge lists as stored in the Flow objects. However, it is harder to merge two topological sorts than it is to merge two edge sets. The topological sort has to be redone with another graph traversal, and edge lists don't let you compute traversals from an arbitrary pixel so easily.

This design lets you interchange flow routing and flat resolution algorithms as desired. Each algorithm returns an edge set that can be merged and then topologically sorted. We can, for example, decouple the D8 flow routing from the least cost auxiliary topography carving that we currently do. And you could conceivably use the LCAT carving with Dinf if you wanted.

The D-inf implementation consists of two functions in src/dinf.c.

• flow_routing_dinf_directions computes the direction bitmap and an array of outgoing flow proportions for each pixel

• flow_routing_dinf_weights computes the edge weights from the array of flow proportions.

A new function flow_routing_tsort is added to src/flow_routing.c to topologically sort bitmap edge sets. This is intended to replace the topological sorting functionality of flow_routing_d8_carve.

A shortest-path flat resolution algorithm is added to src/flow_routing.c with two functions

• resolve_flats_shortest_path computes the direction bitmap using Dijkstra's algorithm to route flow along flats following the shortest distance to a sill
• resolve_flats_shortest_path_weights computes the flow weights for the shortest-path flat resolution. These are all 1s, but this function is provided to maintain the direction/weights interface.

The shortest path resolution has been added because it is what is described by Tarboton (1997). He uses an iterative algorithm, but it amounts to the same thing. This is not what is used by default in MATLAB TopoToolbox or, as far as I can tell, in Tarboton's TauDEM.

Moreover this implementation has been done entirely from scratch based on the description in Tarboton (1997) without reliance on the D-inf implementation in MATLAB TopoToolbox or TauDEM. It probably does not behave exactly as in MATLAB TopoToolbox. The flat resolution approach is certainly different. MATLAB TopoToolbox's D-inf implementation also fails on some DEMs with NaNs such as the Taiwan example. I think the problem there is when there are embedded NaNs that are not contiguous with the border.

libtopotoolbox handles embedded NaNs as missing data. Because of the way D-inf computes flow angles based on triangular facets, this can lead to some surprising results, such as the following situation:

 8   9   7
 5   5   6
NaN  1  NaN

The center pixel should probably flow down, but because both facets containing the downward neighbor have NaNs, the slope of each facet is NaN, and the flow is unresolved. The flat resolution algorithm will probably flow to the neighbor of equal elevation to the left and from there to the bottom right. This actually occurs in the Kunashiri example DEM. It could be possible to fix this by special casing the situation where one of the two neighboring pixels for a facet is NaN.

Functions for dealing with the edge sets have been added to src/helpers/edgeset.c. These could conceivably also be done in higher level languages, but making sure you get all the ordering right is a bit tricky.

Test cases have been added to test/dinf.c. These tests use the Grid and Flow object interface added in #221 and test on a variety of small test cases, including those from the Tarboton paper, the snapshot DEMs and randomly generated ones. Memory order tests are also conducted. I have also taken this opportunity to prototype a new testing harness that provides TAP output.