Add adaptive Dinf-MFD hybrid flow direction

[brendancol]

Problem

Adaptive MFD (#980) transitions from dispersive routing on hillslopes to D8-like concentrated routing in channels by raising the partition exponent. But at high exponent values, flow collapses to the single steepest neighbor (D8), losing angular resolution. A channel running at 37 degrees from east dumps all flow into one cardinal neighbor instead of splitting it proportionally between E and NE.

Dinf solves the angular resolution problem but applies the same single-direction logic everywhere, which is wrong on hillslopes where flow should disperse.

Proposed algorithm

Blend MFD and Dinf based on upstream accumulation:

1. First pass: MFD accumulation with low exponent (e.g. p=1.0) to get an accumulation grid.
2. Second pass: at each cell, compute both MFD fractional weights and the Dinf angle. Blend them using a transition parameter derived from accumulation:

t = min(a / a_threshold, 1.0)
final_fractions = (1 - t) * mfd_fractions + t * dinf_fractions

Where dinf_fractions is an 8-band array with nonzero values only in the two neighbors bracketing the Dinf angle, split by angular proximity (the same decomposition flow_accumulation already uses for Dinf inputs).

• Low accumulation (hillslopes): t near 0, MFD dispersion dominates
• High accumulation (channels): t near 1, Dinf angular precision dominates

Why not just use adaptive MFD (#980)?

At high exponent, adaptive MFD approximates D8, which has 45-degree angular resolution. Dinf gives continuous angular resolution. For channels that don't align with cardinal/diagonal directions, this produces more accurate flow paths and better-shaped contributing areas.

Whether the angular precision matters in practice is an open question. D8 and Dinf give similar accumulation patterns in steep convergent terrain. This may only matter for applications sensitive to exact flow path geometry (e.g. channel network extraction, specific catchment area).

Function signature

def flow_direction_adaptive_dinf(
    agg: xr.DataArray,
    a_threshold: float = 100.0,
    p_hillslope: float = 1.0,
    name: str = 'flow_direction_adaptive_dinf',
    boundary: str = 'nan',
) -> xr.DataArray

Returns the same 3D (8, H, W) fractional flow array as flow_direction_mfd.

References

• Tarboton, D.G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research, 33(2), 309-319.
• Qin, C. et al. (2007). An adaptive approach to selecting a flow-partition exponent for a multiple-flow-direction algorithm. International Journal of Geographical Information Science, 21(4), 443-458.