Add river-constrained D8 flow direction

[kongdd]

Summary

Build hydrologically consistent D8/LDD flow direction from a coarse DEM while enforcing trusted high-resolution river centerlines.

The design deliberately separates the information sources:

• the DEM controls hillslope drainage;
• trusted vector hydrography controls channel location and downstream routing;
• a routing-only conditioned DEM connects noisy/flat terrain to the trusted river network.

Full pipeline

hydro_enforced_flowdir performs:

1. vector river -> strictly 8-connected raster flow paths;
2. optional upstream/downstream orientation from the full DEM profile;
3. conservative raster confluence repair (junction_radius);
4. river-corridor DEM conditioning;
5. topological monotonic channel-profile enforcement through confluences;
6. river-seeded Priority-Flood to resolve remaining depressions/flats;
7. D8 calculation for non-river terrain;
8. forced river-cell LDD from trusted flow paths;
9. O(N), low-memory cycle validation of the final D8 functional graph.

The low-memory validator avoids constructing a full Graphs.jl adjacency graph just to validate large 30 m routing rasters. RiverGraph remains available for downstream network analysis and detailed QC.

GIS vector input

ArchGDAL is now an optional Julia package extension rather than a hard dependency. It reads LineString and MultiLineString layers directly from any GDAL-supported vector source; MultiLineString parts are flattened to independent river lines. No implicit reprojection is performed.

using RiverGraphs, ArchGDAL

state = hydro_enforced_flowdir(
    dem,
    "rivers.gpkg",
    lon,
    lat;
    cellsize=(30.0, 30.0),
    layer=0,
    direction=:geometry,
    junction_radius=1,
    burn_depth=0.0,
    bank_drop=0.01,
    min_slope=1e-4,
    boundary_outlets=true,
    return_dem=true,
)

flowdir = state.flowdir
routing_dem = state.dem
river_paths = state.paths
qc = flowdir_qc(flowdir, river_paths)

If the source river lines are not reliably digitized upstream -> downstream, use direction=:dem. If vector topology is already clean, keep junction_radius=0; one-cell snapping is intentionally conservative and should only repair rasterization gaps.

Priority-Flood behavior

Priority-Flood can use all trusted river cells as drainage seeds. It uses a small constant vertical epsilon per 8-neighbour step to remove flats while preserving monotone heap ordering. The implementation uses separate Float64 priority and linear-Int index heaps, avoiding a DataStructures.jl dependency.

• boundary_outlets=true: for a regional rectangular DEM where drainage may legitimately leave the raster boundary.
• boundary_outlets=false: for a clipped basin/domain that should drain only to the supplied river network; disconnected valid DEM components become explicit errors instead of silent pits.

Quality control

flowdir_qc reports channel-edge agreement and classifies final sinks into trusted river sinks, active-domain boundary sinks, and interior non-river sinks. A typical accepted result should have:

• river_match == 1.0;
• interior_nonriver_sinks == 0;
• expected river outlet count/topology.

Main API

• hydro_enforced_flowdir
• read_river_lines (ArchGDAL extension)
• priority_flood_dem
• condition_river_dem
• river_constrained_flowdir
• d8_flowdir
• rasterize_flowpath / rasterize_flowpaths
• orient_flowpath
• snap_flowpath_junctions
• force_flowpaths!
• validate_flowdir
• flowdir_qc
• ldd_code

Tests

Tests cover D8 coding, vector -> 8-connected raster paths, DEM-based orientation, confluence snapping, divergence rejection, monotonic channel conditioning, tributary confluences, river-seeded Priority-Flood, full hydro enforcement, ArchGDAL LineString/MultiLineString ingestion, direct vector-file routing, low-memory cycle validation, and QC metrics.

Remaining before merge

• real-data validation on a representative 30 m DEM + high-resolution river layer;
• large-raster runtime/peak-memory benchmark;
• preserve Float32 through river-conditioning to reduce peak memory on regional products;
• optional explicit reach-to-reach topology input, preferable to heuristic junction snapping when available.

This PR remains draft until real-data validation and large-raster benchmarking are done.

[codecov]

Welcome to Codecov 🎉

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

Thanks for integrating Codecov - We've got you covered ☂️