Add Zhou depression filling, flat resolution and stencil flow directions by evetion · Pull Request #90 · Deltares/Geomorphometry.jl

evetion · 2026-11-09T20:44:56Z

First of four, stacked on #83. Three serial pieces the package was missing, each added as a method alongside what is there rather than a replacement.

filldepressions takes a method. The existing PriorityFlood, or Zhou (Zhou et al. 2016), which reaches the same surface 1.1–1.4x faster by region-growing slope cells on a plain queue instead of a priority queue. There is a test that the two agree exactly, including on quantized DEMs where large flats make divergence most likely.

resolveflats gives flats a drainage gradient (Barnes et al. 2014b): two breadth-first gradients, away from higher terrain and towards lower, the latter at twice the weight so no iteration is needed. Filling leaves flats everywhere, and no flow direction is defined over one, so filldepressions now calls this by default; flats=false gives the filled surface alone. Cells are raised by single nextfloat steps, so the surface is unchanged at any meaningful precision — one existing assertion in runtests.jl moves from ==(10) to ≈(10) for that reason, with the exact form kept under flats=false.

flowdirection derives directions from the local neighbourhood alone, for D8, DInf and FD8, returning the share each receiver takes beside them. D8 has a single receiver that takes everything, so its weights are a Fill and cost nothing whatever the size of the grid. The trailing weights axis is only as long as the method can need: 1, 2 or 8.

The three methods share one stencil, and the partitioning rules are shared with flowaccumulation through _dinfreceivers and _fd8weights!, so the carved and derived routes cannot drift apart.

Two things worth a reviewer's eye:

_dinfreceivers guards a non-finite aspect. It is undefined over a flat, and bracketing a NaN picks the cell itself as both of its receivers.
FD8 caches the neighbourhood geometry per cell instead of recomputing every bearing per neighbour pair (5.6x faster), and clamps angular gaps to 90°. At a border a cell is missing neighbours; taken literally the contour length runs to tand(90) = Inf and that neighbour takes the cell's whole flow.

Weights live here rather than in the accumulation PR on purpose: they are what flowdirection returns, and splitting them would have meant writing a weightless version and then rewriting it.

🤖 Generated with Claude Code

https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf

evetion · 2026-11-09T21:02:16Z

Added one commit after the initial push, which also touches on #85.

flowdirection allocated its direction grid with similar(dem, FlowDirection{...}). On a Raster carrying a missingval that throws — Rasters converts the fill through typemax, and a direction eltype has no such method. flowaccumulation has the same bug on this base, so it is not a regression, but worth not extending:

flowaccumulation D8   THROWS: MethodError: no method matching typemax(::Type{FlowDirection{LDD, UInt8}})
flowdirection D8      THROWS: MethodError: no method matching typemax(::Type{FlowDirection{LDD, UInt8}})

Both now go through _directiongrid, taken from #85 unchanged in name and signature so the two converge on merge rather than collide, with the Rasters extension overriding it to allocate without a missingval.

The weights grid had a related problem: it used bare zeros, which drops a wrapper type's dimensions. It now reuses _alloc_directions, the hook horizon_angle already uses for exactly this (rows, cols, directions) shape and which the Rasters and GeoArrays extensions already implement. So FD8 weights come back as a Raster with a Band dimension rather than a bare Array.

After the change, on a Raster with a missingval:

flowaccumulation D8   ok  -> Raster (20, 20)
flowdirection D8      ok  -> Fill (20, 20, 1)      # D8 stores no weights
flowdirection FD8     ok  -> Raster (20, 20, 8)
watersheds            ok  -> Raster (20, 20)

🤖 Generated with Claude Code

https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf

evetion · 2026-11-10T09:35:52Z

+    acell, bcell, afraction, bfraction = flowneighbors
+
+    # A receiver the partition names but that is not lower takes nothing, so the shares of
+    # those that remain are renormalised before being recorded.


This is unreadable.

Rewritten as one loop over the two receivers, instead of two parallel sets of takesa/abit flags. Same logic, and I benchmarked it to be sure the loop cost nothing: DInf over 1500x1500 is 0.437 s against 0.431 s before, inside noise.

evetion · 2026-11-10T09:36:57Z

+### Filling and flow directions
+
+[`filldepressions`](@ref) takes a `method`: the default [`PriorityFlood`](@ref) of
+[barnesPriorityFloodOptimalDepressionFilling2014](@citet), or [`Zhou`](@ref) after
+[zhouEfficientVariantPriorityFlood2016](@citet), which reaches the same surface 1.1-1.4x
+faster by region-growing slope cells on a plain queue instead of a priority queue.
+
+Filling leaves flats, over which no flow direction is defined, so by default it also gives
+them a drainage gradient with [`resolveflats`](@ref), after
+[barnesEfficientAssignmentDrainage2014](@citet). Cells are raised by single `nextfloat`
+steps, so the surface is unchanged at any meaningful precision. Pass `flats=false` for the
+filled surface alone.
+
+[`flowdirection`](@ref) then reads directions off that surface, looking only at each cell's
+neighbours, and returns them with the share each receiver takes:
+
+```@example plots
+z = filldepressions(cdtm)                            # fill depressions, resolve flats
+directions, weights = flowdirection(z; method=D8())  # directions off the filled surface
+nothing # hide
+```
+
+[`D8`](@ref) has a single receiver that takes everything, so its weights are a `Fill` and
+cost nothing whatever the size of the grid. [`DInf`](@ref) and [`FD8`](@ref) store a share
+per receiver, two and eight respectively, which `weighttype=Float16` halves.
+


This is way too much text and the example doesn't plot anything.

Section is 22 lines now from 48, and the example shows the arrow encoding (parent(directions[400:800, 300:700])) rather than assigning to nothing. My first attempt plotted sum(weights; dims=3), which throws on a Raster with a Band lookup, so it is the arrows instead.

evetion · 2026-11-10T09:40:19Z

+_convention(::D8) = LDD
+_convention(::DInf) = D8D
+_convention(::FD8) = D8D


Why do these have underscores, and the rest not? Stuff is only public if we export it explicitely, so feel free not to use _.

Dropped. direction.jl and flats.jl are new files that mixed both styles, so their helpers are plain now: convention, directionstate, setshare!, steepestdescent, celldirection!, flatmask, growflat!. _weights became allocweights because weights is a local name throughout that file. I left hydrology.jl alone deliberately: it is underscored throughout, so renaming only my additions there would have made it less consistent, not more. Say the word if you want that file swept too.

evetion · 2026-11-10T09:44:05Z

+@inline _setshare!(weights::AbstractArray{T}, cell, slot, fraction) where {T} =
+    @inbounds weights[cell, slot] = T(fraction)
+
+# Per-method state cheaper to derive once for the whole DEM than per cell.


What does that even mean? Is this a cache?

Fair. It is not a cache, it is per-method setup done once for the whole grid rather than per cell: a grid of aspects for DInf off a raster, and the reused scratch buffers for FD8. Methods needing neither get nothing. That is now the docstring.

Copilot

🟡 Changes recommended

Unresolved critical defects allow flat resolution to create new pits, produce inconsistent FD8 weights, and cause DInf to misroute flow on negative-spacing rasters.

5 open issues

What changed in this PR

Extends the generic hydrology APIs introduced in #83 with depression-filling options, flat resolution, and local flow directions with receiver weights.

Changes:

Adds Zhou filling and enables flat resolution by default, with flats=false preserving fill-only behavior.
Adds D8, DInf, and FD8 direction stencils and receiver weights, plus shared routing helpers.
Adds hydrology tests, raster allocation support, usage examples, and algorithm references.

File	Description
`test/runtests.jl`	Includes hydrology tests and adjusts filling assertions.
`test/hydrology.jl`	Tests filling, flat drainage, directions, and weights.
`src/hydrology.jl`	Adds filling methods and refactors routing helpers.
`src/Geomorphometry.jl`	Includes and exports the new APIs.
`src/flats.jl`	Implements drainage gradients across flats.
`src/direction.jl`	Implements local flow directions and receiver weights.
`Project.toml`	Adds Random to test dependencies.
`ext/GeomorphometryRastersExt.jl`	Allocates raster direction grids without a missing-value sentinel.
`docs/src/usage.md`	Documents filling options and direction derivation.
`docs/src/refs.bib`	Adds Barnes flat-resolution and Zhou filling references.

_{Review effort: Balanced}

+
+# Per-method state cheaper to derive once for the whole DEM than per cell.
+_directionstate(::FlowDirectionMethod, dem, cellsize) = nothing
+_directionstate(::DInf, dem, cellsize) = _dinf_aspects(dem, cellsize)


+function _flowdirection!(weights, fd8::FD8, scratch, dem, cell, cellsize)
+    z = dem[cell]
+    isfinite(z) || return zero(UInt8)
+    total = _fd8weights!(scratch, fd8, dem, cell, _Unvisited(), cellsize)


+    increments = _flatmask(dem)
+    @inbounds for cell in eachcell(dem)
+        increment = increments[cell]
+        iszero(increment) || (dem[cell] = nextfloat(dem[cell], increment))


+    flats::Bool = true,
+)
+    _filldepressions!(method, dem, queued)
+    flats && resolveflats!(dem)


+    flats::Bool = true,
+)
+    _filldepressions!(method, dem, queued)
+    flats && resolveflats!(dem)


Three serial pieces the package was missing, each as a method alongside what is already there rather than a replacement. `filldepressions` takes a `method`: the existing `PriorityFlood` or `Zhou`, after Zhou et al. (2016), which reaches the same surface 1.1-1.4x faster by region-growing slope cells on a plain queue instead of a priority queue. `resolveflats` gives flats a drainage gradient, after Barnes et al. (2014b): two breadth-first gradients, away from higher terrain and towards lower, the latter at twice the weight so no iteration is needed. `filldepressions` calls it by default, since a filled surface is otherwise full of flats that have no flow direction; `flats=false` gives the filled surface alone. `flowdirection` derives directions from the local neighbourhood alone, for D8, DInf and FD8, returning the share each receiver takes beside them. D8 has one receiver that takes everything, so its weights are a `Fill` and cost nothing. The three share one stencil and the partitioning rules are shared with `flowaccumulation` through `_dinfreceivers` and `_fd8weights!`, so the carved and derived routes cannot drift apart. `_dinfreceivers` also guards a non-finite aspect, which is undefined over a flat and would otherwise bracket to the cell itself as both receivers. FD8 caches the neighbourhood geometry per cell rather than recomputing every bearing per neighbour pair, which is 5.6x faster, and clamps angular gaps to 90 degrees: at a border a cell is missing neighbours, and the contour length otherwise runs to tand(90) = Inf and takes the cell's whole flow. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf

…heir type `flowdirection` allocated its direction grid with `similar(dem, FlowDirection{...})`, which throws on a `Raster` carrying a `missingval`: Rasters converts the fill through `typemax`, and a direction eltype has no such method. `flowaccumulation` has the same bug on this base. Both now go through `_directiongrid`, which the Rasters extension overrides to allocate without a missingval. The weights grid used bare `zeros`, which drops the dimensions of a wrapper type. It now reuses `_alloc_directions`, the hook `horizon_angle` already uses for exactly this (rows, cols, directions) shape and which the Rasters and GeoArrays extensions already implement. `_directiongrid` is taken from #85 unchanged, name and signature, so the two converge rather than collide. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf

Four defects, all from making flat resolution the default: - An integer DEM reached `nextfloat(::Int, ::Int32)` and threw. `flats` now defaults to true only for a floating point DEM, since integer elevations cannot express a gradient of single float steps, so integer input behaves as it did before. - `resolveflats!` ignored the caller's exclusion mask, so cells documented as excluded were modified. The flood overwrites `queued`, so the mask is kept and handed to `resolveflats!`, which treats those cells as already draining. - The increments assumed the surrounding relief dwarfs a few float steps. Where it does not, a raised cell could overtake a neighbour that was strictly higher and take away that neighbour's drainage. Each increment is now capped below its lowest higher neighbour. It costs 0.08 s on a quantized 1000x1000 DEM, since ample headroom is one comparison. - `_dinf_aspects` dropped the signs of `cellsize` while `_cellbearing` keeps them, leaving the aspect and the receiver bearings in different frames. On a negative-spacing raster DInf picked a neighbour of gradient 0.707 where 2.121 was available. FD8 contour lengths now come from one place, `_rastercontour`, for both `flowdirection` and `flowaccumulation`. They had separate formulas, so the two disagreed on the same raster: a symmetric peak split eight equal ways through the stencil and 2:1 cardinal-to-diagonal through accumulation. The Moore neighbourhood is not a Voronoi one, so the angular construction degenerates on it and a square grid needs the fixed table. Also, on review comments: `direction.jl` and `flats.jl` are new files that mixed underscore and plain names, so the helpers they add lose the underscore (`hydrology.jl` keeps its own convention, which is underscored throughout). `directionstate` says what it holds rather than that it is cheaper. The DInf branch reads as one loop over the two receivers instead of two parallel sets of flags, at no cost in speed. The manual section is a third shorter and its example now shows something. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf

evetion · 2026-11-10T13:50:45Z

Rebased onto main now that #83 is merged, and addressed the review. I verified every Copilot finding against a repro before acting on it; all of them were real, four of them regressions I introduced with flats-by-default.

	finding	verdict
mask	flat resolution ignored the exclusion mask	real, mine
integers	`filldepressions(fill(5, 5, 5))` threw `nextfloat(::Int, ::Int32)`	real, mine
headroom	increments could overtake strictly higher terrain and remove its drainage	real, mine
FD8	stencil and accumulation used different contour lengths	real
signs	`_dinf_aspects` dropped the signs of `cellsize`	real, pre-existing

Two were wrong in their details while right in substance, which is worth knowing when reading the rest:

The FD8 one predicted 0.2 / 0.05 for accumulation. The actual split was 0.1667 / 0.0833; the disagreement was real, the numbers were not.
The cellsize sign one is inverted. It called (1, 0)/(1, -1) the bug and (1, 0)/(1, 1) correct, but those are the stored codes, which consumers re-orient. Resolved to index offsets, the buggy path picked a neighbour of gradient 0.707 where 2.121 was available, and the output it called correct is the correct one. The root cause and the file:line it named were right, so the fix stands; the symptom as described was backwards.

On the FD8 split — it directly contradicted this PR's own claim that the rules are shared, so there is now one _rastercontour that both flowdirection and flowaccumulation read. A symmetric peak gives the same shares either way, and a test asserts it by comparing the stencil's weights against what accumulation actually delivered. The Moore neighbourhood is not a Voronoi one — a diagonal shares only a corner, so there is no face to measure — which is why a square grid needs the fixed table and the angular construction is kept for grids where it is exact.

flats now defaults to true only for a floating point DEM. Integer elevations cannot express a gradient of single float steps, so integer input behaves exactly as it did before rather than throwing.

On your comments: direction.jl and flats.jl are new files that mixed both naming styles, so their helpers lose the underscore; hydrology.jl keeps its own, which is underscored throughout, since renaming only my additions there would make it less consistent rather than more. _weights became allocweights because weights is a local name all over that file. directionstate now says what it holds. The DInf branch is one loop over the two receivers instead of two parallel sets of flags, and I benchmarked it to be sure that cost nothing (0.437 s vs 0.431 s for DInf over 1500x1500, inside noise). The manual section is a third shorter and its example shows the arrow encoding.

🤖 Generated with Claude Code

https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf

evetion mentioned this pull request Nov 9, 2026

Fix three flow-direction defects: D8 steepest descent, DInf facets, FD8 contours #91

Open

evetion added this pull request to stack #94 November 9, 2026 21:27

evetion force-pushed the feat/zhou-flats-flowdirection branch 3 times, most recently from 654ea7a to 943c56a Compare November 10, 2026 07:06

Base automatically changed from feat/generic to main November 10, 2026 09:14

evetion force-pushed the feat/zhou-flats-flowdirection branch from 943c56a to ca78295 Compare November 10, 2026 09:14

evetion commented Nov 10, 2026

View reviewed changes

evetion requested a balanced review from Copilot November 10, 2026 12:10

Copilot started reviewing on behalf of evetion November 10, 2026 12:11

View session

Copilot AI reviewed Nov 10, 2026

View reviewed changes

evetion and others added 3 commits November 10, 2026 14:45

evetion force-pushed the feat/zhou-flats-flowdirection branch from ca78295 to 0dd6635 Compare November 10, 2026 13:49

Navigation Menu

Repository navigation

Add Zhou depression filling, flat resolution and stencil flow directions - #90

Add Zhou depression filling, flat resolution and stencil flow directions#90
evetion wants to merge 3 commits into
mainfrom
feat/zhou-flats-flowdirection

evetion commented Nov 9, 2026

Uh oh!

evetion commented Nov 9, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

evetion Nov 10, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

evetion commented Nov 10, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

evetion commented Nov 9, 2026

Uh oh!

evetion commented Nov 9, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🟡 Changes recommended

Uh oh!

evetion commented Nov 10, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants