From 0b38118b3d8f2e5cd80f48e276f628544d6a3856 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Mon, 9 Nov 2026 21:33:39 +0100 Subject: [PATCH 1/3] Add Zhou depression filling, flat resolution and stencil flow directions 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) Claude-Session: https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf --- Project.toml | 4 +- docs/src/refs.bib | 34 ++++- docs/src/usage.md | 26 ++++ src/Geomorphometry.jl | 4 + src/direction.jl | 179 +++++++++++++++++++++++ src/flats.jl | 175 +++++++++++++++++++++++ src/hydrology.jl | 326 +++++++++++++++++++++++++++++++++++------- test/hydrology.jl | 140 ++++++++++++++++++ test/runtests.jl | 6 +- 9 files changed, 836 insertions(+), 58 deletions(-) create mode 100644 src/direction.jl create mode 100644 src/flats.jl create mode 100644 test/hydrology.jl diff --git a/Project.toml b/Project.toml index bc3e231..4c180d0 100644 --- a/Project.toml +++ b/Project.toml @@ -47,6 +47,7 @@ LocalFilters = "1.2" OffsetArrays = "1.10" PaddedViews = "0.5" QuickHeaps = "0.2" +Random = "1" Rasters = "0.13, 0.14, 0.15" StaticArrays = "1" Statistics = "1" @@ -61,8 +62,9 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" Eikonal = "a6aab1ba-8f88-4217-b671-4d0788596809" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" GeoArrays = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "GeoArrays", "Rasters", "ArchGDAL", "Eikonal", "Aqua", "ExplicitImports"] +test = ["Test", "Random", "GeoArrays", "Rasters", "ArchGDAL", "Eikonal", "Aqua", "ExplicitImports"] diff --git a/docs/src/refs.bib b/docs/src/refs.bib index bb65f7e..4066b7a 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -17,7 +17,6 @@ @article{barnesPriorityFloodOptimalDepressionFilling2014 keywords = {Computer Science - Data Structures and Algorithms} } - @inproceedings{bartelsDTMGenerationLIDAR2006, title = {{{DTM Generation}} from {{LIDAR Data}} Using {{Skewness Balancing}}}, booktitle = {18th {{International Conference}} on {{Pattern Recognition}} ({{ICPR}}'06)}, @@ -34,7 +33,6 @@ @inproceedings{bartelsDTMGenerationLIDAR2006 ZSCC: 0000080} } - @article{bartelsThresholdfreeObjectGround2010, title = {Threshold-Free Object and Ground Point Separation in {{LIDAR}} Data}, author = {Bartels, Marc and Wei, Hong}, @@ -53,7 +51,6 @@ @article{bartelsThresholdfreeObjectGround2010 keywords = {_tablet,DSM,DTM,GIS,LIDAR,Remote sensing,Skewness Balancing} } - @book{burroughPrinciplesGeographicalInformation2015, title = {Principles of {{Geographical Information Systems}}}, author = {Burrough, P. A. and McDonnell, Rachael A. and Lloyd, Christopher D.}, @@ -66,7 +63,6 @@ @book{burroughPrinciplesGeographicalInformation2015 keywords = {Business & Economics / Economic History,Science / Earth Sciences / Geography,Technology & Engineering / Cartography,Technology & Engineering / Remote Sensing & Geographic Information Systems,Technology & Engineering / Surveying,Travel / General} } - @article{cavalliEffectivenessAirborneLidar2008, title = {The Effectiveness of Airborne {{LiDAR}} Data in the Recognition of Channel-Bed Morphology}, author = {Cavalli, Marco and Tarolli, Paolo and Marchi, Lorenzo and {Dalla Fontana}, Giancarlo}, @@ -106,7 +102,6 @@ @inproceedings{eastman1989pushbroom year = {1989} } - @article{hornHillShadingReflectance1981, title = {Hill Shading and the Reflectance Map}, author = {Horn, B.K.P.}, @@ -123,7 +118,6 @@ @article{hornHillShadingReflectance1981 keywords = {Automation,Graphics,History,Photography,Production,Reflectivity,Shadow mapping,Streaming media,Surface topography,Surface treatment} } - @article{ilichMultiscaleDTMOpensourcePackage2023, title = {{{MultiscaleDTM}}: {{An}} Open-Source {{R}} Package for Multiscale Geomorphometric Analysis}, shorttitle = {{{MultiscaleDTM}}}, @@ -475,3 +469,31 @@ @article{zhaoFastSweepingMethod2005 keywords = {characteristics,Courant-Friedrichs-Lewy (CFL) condition,Eikonal equation,Gauss-Seidel iteration,Godunov scheme,upwind difference} } +@article{barnesEfficientAssignmentDrainage2014, + title = {An Efficient Assignment of Drainage Direction over Flat Surfaces in Raster Digital Elevation Models}, + author = {Barnes, Richard and Lehman, Clarence and Mulla, David}, + year = {2014}, + month = jan, + journal = {Computers \& Geosciences}, + volume = {62}, + pages = {128--135}, + issn = {00983004}, + doi = {10.1016/j.cageo.2013.01.009}, + abstract = {Depression-filling leaves flat regions with no local relief, over which drainage directions are undefined. This paper presents an algorithm which assigns drainage directions over flats by combining a gradient away from higher terrain with a gradient towards lower terrain, computed with two breadth-first traversals. The result is a realistic, convergent drainage pattern in O(n) time.}, + keywords = {flats,flow direction,hydrology} +} + +@article{zhouEfficientVariantPriorityFlood2016, + title = {An Efficient Variant of the {{Priority-Flood}} Algorithm for Filling Depressions in Raster Digital Elevation Models}, + author = {Zhou, Guiyun and Sun, Zhongxuan and Fu, Suhua}, + year = {2016}, + month = may, + journal = {Computers \& Geosciences}, + volume = {90}, + pages = {87--96}, + issn = {00983004}, + doi = {10.1016/j.cageo.2016.02.021}, + abstract = {A variant of Priority-Flood that keeps the majority of cells out of the priority queue. Slope cells are processed by region growing on two plain queues, and only potential spill cells, where terrain turns back downhill, are pushed onto the priority queue. This reduces the complexity to O(m log m) with m much smaller than the number of cells, and runs 44.6\% faster than Barnes et al. (2014) in testing.}, + keywords = {depression filling,hydrology,priority flood} +} + diff --git a/docs/src/usage.md b/docs/src/usage.md index 68397ae..7197c2a 100644 --- a/docs/src/usage.md +++ b/docs/src/usage.md @@ -378,6 +378,32 @@ depression_volume(dtm; filled=fdtm) ::: +### 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. + Filling the depressions in a DEM is not necessary to calculate the flow accumulation. Here we use [`flowaccumulation`](@ref) to do so, and it automatically carves out depressions. Note that the local drainage direction is also returned. By default the FD8 algorithm is used, but you can also use the D∞ or D8 algorithm by setting the `method` keyword argument to `DInf()` or `D8()`. ::: warning diff --git a/src/Geomorphometry.jl b/src/Geomorphometry.jl index b8767e4..a616ad7 100644 --- a/src/Geomorphometry.jl +++ b/src/Geomorphometry.jl @@ -40,6 +40,8 @@ include("spread.jl") include("skew.jl") include("flowdir.jl") include("hydrology.jl") +include("flats.jl") +include("direction.jl") include("horizon.jl") include("geomorphon.jl") @@ -60,6 +62,8 @@ export roughness, export slope, aspect, laplacian, plan_curvature, profile_curvature, tangential_curvature export skewness_balancing +export PriorityFlood, Zhou +export resolveflats, resolveflats!, flowdirection export filldepressions, flowaccumulation, topographic_wetness_index, diff --git a/src/direction.jl b/src/direction.jl new file mode 100644 index 0000000..e538823 --- /dev/null +++ b/src/direction.jl @@ -0,0 +1,179 @@ +# Single-receiver methods use the same LDD encoding as `flowaccumulation`; methods that can +# send flow several ways need the bitwise D8D encoding. +_convention(::D8) = LDD +_convention(::DInf) = D8D +_convention(::FD8) = D8D + +# A method that can send flow several ways needs the share each receiver takes; a single +# receiver takes everything, so there is nothing to carry. +hasweights(::FlowDirectionMethod) = false +hasweights(::DInf) = true +hasweights(::FD8) = true + +"""Most receivers a method can ever give one cell.""" +maxreceivers(::D8) = 1 +maxreceivers(::DInf) = 2 +maxreceivers(::FD8) = 8 + +""" +Share of each cell's flow per receiver, one dimension larger than the DEM, the trailing +axis in [`D8D`](@ref) bit order. A single-receiver method stores nothing: its weights are a +[`Fill`](@extref FillArrays.Fill) and cost no memory. +""" +_weights(method::FlowDirectionMethod, dem, T) = + hasweights(method) ? zeros(T, size(dem)..., maxreceivers(method)) : + Fill(one(T), size(dem)..., maxreceivers(method)) + + +"""Record the share of a cell's `slot`-th receiver. Nothing to record without weights.""" +@inline _setshare!(::Fill, cell, slot, fraction) = fraction +@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. +_directionstate(::FlowDirectionMethod, dem, cellsize) = nothing +_directionstate(::DInf, dem, cellsize) = _dinf_aspects(dem, cellsize) +_directionstate(::FD8, dem, cellsize) = _fd8scratch(dem) + +""" + flowdirection(dem; method=D8(), cellsize=cellsize(dem), weighttype=Float32) + +Derive flow directions from `dem` by comparing each cell with its neighbours, and return +them with the share each receiver takes. The encoding matches [`flowaccumulation`](@ref): +[`LDD`](@ref) for the single receiver of [`D8`](@ref), [`D8D`](@ref) for the methods that +can send flow several ways at once. + +These depend only on the local neighbourhood, so a cell with no lower neighbour is a pit. +Pass a surface from [`filldepressions`](@ref) if every cell must drain. + +# Arguments +- `dem`: An indexed digital elevation model. +- `method::FlowDirectionMethod=D8()`: [`D8`](@ref), [`DInf`](@ref) or [`FD8`](@ref). +- `cellsize`: Cell size, used to scale horizontal distances. Derived from `dem` by default. +- `weighttype=Float32`: Element type of the weights. `Float16` halves the memory a + divergent method needs, at the cost of a little of the flow, as a cell's shares then no + longer sum to exactly one. Integer types are not usable. +""" +function flowdirection( + dem; + method::FlowDirectionMethod = D8(), + cellsize = cellsize(dem), + weighttype::Type{<:AbstractFloat} = Float32, +) + C = _convention(method) + directions = similar(dem, FlowDirection{C, UInt8}) + weights = _weights(method, dem, weighttype) + state = _directionstate(method, dem, cellsize) + @inbounds for cell in eachcell(dem) + code = _flowdirection!(weights, method, state, dem, cell, cellsize) + directions[cell] = FlowDirection{C}(code) + end + return directions, weights +end + + +""" +Offset towards `cell`'s steepest downslope neighbour, or zero when it has none: the +[`D8`](@ref) rule of +[Jenson and Domingue (1988)](@cite jensonExtractingTopographicStructure1988), and the +fallback wherever a method leaves a cell undecided. +""" +function _steepestdescent(dem, cell, cellsize) + offset = cell - cell + z = dem[cell] + isfinite(z) || return offset + steepest = zero(float(typeof(z))) + @inbounds for neighbor in neighbors(dem, cell) + nz = dem[neighbor] + (isfinite(nz) && nz < z) || continue + distance = celldistance(dem, cell, neighbor; cellsize) + iszero(distance) && continue + gradient = (z - nz) / distance + if gradient > steepest + steepest = gradient + offset = neighbor - cell + end + end + return offset +end + +_steepestcode(dem, cell, cellsize) = + _directioncode(dem, _steepestdescent(dem, cell, cellsize), cellsize) + +_flowdirection!(weights, ::D8, _, dem, cell, cellsize) = + _orient(_steepestdescent(dem, cell, cellsize), cellsize) + +# DInf: the two neighbours bounding the steepest facet, using the same `_dinfreceivers` +# partition as `flowaccumulation`. +function _flowdirection!(weights, ::DInf, state, dem, cell, cellsize) + z = dem[cell] + isfinite(z) || return zero(UInt8) + + flowneighbors = _dinfreceivers(dem, cell, state, cellsize) + isnothing(flowneighbors) && return _steepestshare!(weights, dem, cell, cellsize) + 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. + takesa = afraction > 0 && isfinite(dem[acell]) && dem[acell] < z + takesb = bfraction > 0 && isfinite(dem[bcell]) && dem[bcell] < z + total = (takesa ? afraction : 0.0) + (takesb ? bfraction : 0.0) + # The partition can name terrain that is not lower, most often beside a flat. Fall back + # rather than leave a cell that can drain without a direction. + total > 0 || return _steepestshare!(weights, dem, cell, cellsize) + + abit = takesa ? _directioncode(dem, acell - cell, cellsize) : zero(UInt8) + bbit = takesb ? _directioncode(dem, bcell - cell, cellsize) : zero(UInt8) + # Shares are stored in the order the directions come back out, which is by bit + if iszero(bbit) + _setshare!(weights, cell, 1, afraction / total) + elseif iszero(abit) + _setshare!(weights, cell, 1, bfraction / total) + elseif abit < bbit + _setshare!(weights, cell, 1, afraction / total) + _setshare!(weights, cell, 2, bfraction / total) + else + _setshare!(weights, cell, 1, bfraction / total) + _setshare!(weights, cell, 2, afraction / total) + end + return abit | bbit +end + +# FD8: every downslope neighbour receives flow, weighted by slope and contour length, using +# the same `_fd8weights!` as `flowaccumulation`. +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) + iszero(total) && return _steepestshare!(weights, dem, cell, cellsize) + + # Gather the shares against their direction bits, then emit them in bit order, which is + # the order the directions come back out in. + fill!(scratch.bitweights, 0.0) + code = zero(UInt8) + @inbounds for (index, neighbor) in enumerate(scratch.neighbors) + weight = scratch.weights[index] + iszero(weight) && continue + bit = _directioncode(dem, neighbor - cell, cellsize) + scratch.bitweights[trailing_zeros(Int(bit)) + 1] = weight / total + code |= bit + end + slot = 0 + @inbounds for bit in decompose(D8D, code) + iszero(bit) && continue + slot += 1 + _setshare!(weights, cell, slot, scratch.bitweights[trailing_zeros(Int(bit)) + 1]) + end + return code +end + +"""The steepest descent fallback: one receiver, taking everything.""" +function _steepestshare!(weights, dem, cell, cellsize) + code = _steepestcode(dem, cell, cellsize) + iszero(code) || _setshare!(weights, cell, 1, one(Float32)) + return code +end + +"""A stencil has no traversal order, so nothing is excluded for being already visited.""" +struct _Unvisited end +Base.getindex(::_Unvisited, _...) = false diff --git a/src/flats.jl b/src/flats.jl new file mode 100644 index 0000000..e6e8a55 --- /dev/null +++ b/src/flats.jl @@ -0,0 +1,175 @@ +""" + resolveflats(dem) + resolveflats!(dem) + +Give the flats of `dem` a drainage gradient, following +[Barnes et al. (2014)](@cite barnesEfficientAssignmentDrainage2014), and return the +elevation model. `resolveflats!` writes into `dem` instead of a copy. + +A flat is a connected region of equal elevation, over which no flow direction is defined. +The gradient is the superposition of one away from higher terrain, which gives convergent +flow, and one towards lower terrain, which guarantees the flat drains. The latter carries +twice the weight, so it settles cells where the two disagree without iterating. + +Cells are raised by whole `nextfloat` steps, so the surface is unchanged at any meaningful +precision but every flat cell gains a strictly lower neighbour. A flat with no outlet is a +depression, not a flat, and is left alone. [`filldepressions`](@ref) calls this by default. +""" +resolveflats(dem) = resolveflats!(copy(dem)) + +function resolveflats!(dem) + increments = _flatmask(dem) + @inbounds for cell in eachcell(dem) + increment = increments[cell] + iszero(increment) || (dem[cell] = nextfloat(dem[cell], increment)) + end + return dem +end + +""" +Per-cell increment counts for [`resolveflats`](@ref): zero wherever a cell already drains, +positive on cells that have no strictly lower neighbour and belong to a flat that drains. +""" +function _flatmask(dem) + R = eachcell(dem) + Cell = typeof(first(R)) + + drains = fill!(similar(dem, Bool), false) + labels = fill!(similar(dem, Int32), 0) + flatmask = fill!(similar(dem, Int32), 0) + + highedges = Cell[] + lowedges = Cell[] + + # Cells on the boundary drain out of the grid, so they never need resolving. + @inbounds for cell in outlets(dem) + drains[cell] = true + end + + # Step 1: a cell drains if it has a strictly lower neighbour; a high edge cell cannot + # drain but touches higher terrain. + @inbounds for cell in R + z = dem[cell] + if !isfinite(z) + drains[cell] = true + continue + end + drains[cell] && continue + lower = false + higher = false + for ncell in neighbors(dem, cell) + nz = dem[ncell] + isfinite(nz) || continue + nz < z && (lower = true) + nz > z && (higher = true) + end + drains[cell] = lower + (!lower && higher) && push!(highedges, cell) + end + + # A low edge cell drains, but borders a cell of its own elevation that does not: it is + # an outlet of that flat. + @inbounds for cell in R + drains[cell] || continue + isfinite(dem[cell]) || continue + for ncell in neighbors(dem, cell) + if dem[ncell] == dem[cell] && !drains[ncell] + push!(lowedges, cell) + break + end + end + end + isempty(lowedges) && return flatmask + + # Label each flat by flood filling its outlets over cells of equal elevation. Flats + # without an outlet keep label 0 and are left alone: they are depressions. + nlabels = zero(Int32) + queue = Cell[] + @inbounds for seed in lowedges + iszero(labels[seed]) || continue + nlabels += one(nlabels) + z = dem[seed] + labels[seed] = nlabels + empty!(queue) + push!(queue, seed) + head = 1 + while head <= length(queue) + cell = queue[head] + head += 1 + for ncell in neighbors(dem, cell) + if iszero(labels[ncell]) && dem[ncell] == z + labels[ncell] = nlabels + push!(queue, ncell) + end + end + end + end + filter!(cell -> !iszero(labels[cell]), highedges) + + # Step 2: gradient away from higher terrain, breadth-first inwards from the high + # edges. `flatheights` is how far each flat reaches, which Step 3 needs to invert it. + flatheights = zeros(Int32, nlabels) + seen = fill!(similar(dem, Bool), false) + _seed!(seen, highedges) + increment = one(Int32) + frontier = copy(highedges) + next = Cell[] + @inbounds while !isempty(frontier) + empty!(next) + for cell in frontier + flatmask[cell] = increment + flatheights[labels[cell]] = increment + _growflat!(next, seen, dem, labels, drains, cell) + end + frontier, next = next, frontier + increment += one(increment) + end + + # Step 3: gradient towards lower terrain, growing back from the outlets. Negating + # Step 2 and adding the flat's height inverts it into the gradient we want; every cell + # then gains twice its distance from the outlet, so draining dominates on disagreement. + @inbounds for cell in R + flatmask[cell] = -flatmask[cell] + end + fill!(seen, false) + _seed!(seen, lowedges) + increment = one(Int32) + frontier = copy(lowedges) + empty!(next) + @inbounds while !isempty(frontier) + empty!(next) + for cell in frontier + _growflat!(next, seen, dem, labels, drains, cell) + flatmask[cell] < 0 && (flatmask[cell] += flatheights[labels[cell]]) + flatmask[cell] += 2 * increment + end + frontier, next = next, frontier + increment += one(increment) + end + + # Cells that already drain keep their elevation: the gradient is only for cells that + # would otherwise have no flow direction. + @inbounds for cell in R + drains[cell] && (flatmask[cell] = 0) + end + return flatmask +end + +function _seed!(seen, cells) + for cell in cells + seen[cell] = true + end + return seen +end + +"""Queue the not-yet-visited cells of `cell`'s own flat that still need a gradient.""" +@inline function _growflat!(next, seen, dem, labels, drains, cell) + label = labels[cell] + @inbounds for ncell in neighbors(dem, cell) + if !seen[ncell] && labels[ncell] == label && !drains[ncell] + seen[ncell] = true + push!(next, ncell) + end + end + return next +end diff --git a/src/hydrology.jl b/src/hydrology.jl index 05ccf24..2fb1da0 100644 --- a/src/hydrology.jl +++ b/src/hydrology.jl @@ -12,17 +12,57 @@ end cellarea(dem::AbstractMatrix, cell; cellsize = cellsize(dem)) = abs(cellsize[1] * cellsize[2]) +abstract type DepressionFillingMethod end + +""" + PriorityFlood() + +Improved Priority-Flood, Algorithm 2 of +[Barnes et al. (2014)](@cite barnesPriorityFloodOptimalDepressionFilling2014). Cells are +flooded inwards from the outlets through a priority queue, and cells found below the +current flood level are filled through a plain queue. Runs in `O(m log m)` time, where +`m` is the number of cells that reach the priority queue. """ - filldepressions(dem, mask=fill!(similar(dem, Bool), false)) +struct PriorityFlood <: DepressionFillingMethod end -Performs the Priority-Flood algorithm [barnesPriorityFloodOptimalDepressionFilling2014](@cite) on the given digital elevation model (DEM) `dem` with an optional `mask`. +""" + Zhou() + +Priority-Flood variant of [Zhou et al. (2016)](@cite zhouEfficientVariantPriorityFlood2016). +Terrain that rises away from an already flooded cell needs no priority ordering, so slope +cells are region-grown on a plain queue and only *potential spill cells*, where the terrain +turns back downhill, reach the priority queue. Runs in `O(N + m log m)` time with `m` much +smaller than the number of cells `N`, and reported 44.6% faster than +[`PriorityFlood`](@ref). Measured here it is worth 1.1-1.4x, most on continuous elevations +with few ties and least on tiles with much nodata. +""" +struct Zhou <: DepressionFillingMethod end + +""" + filldepressions(dem, mask=fill!(similar(dem, Bool), false); method=PriorityFlood(), flats=true) + +Fill the depressions of a digital elevation model, so that every cell drains to an outlet. # Arguments - `dem`: An indexed digital elevation model. - `mask`: An optional boolean array representing the mask. Cells with `true` values are treated as already filled (queued) and excluded from the computation. Defaults to all `false`. +- `method::DepressionFillingMethod=PriorityFlood()`: Algorithm used, either + [`PriorityFlood`](@ref) or [`Zhou`](@ref). The two produce the same filled surface for a + DEM of finite elevations. They do **not** agree once the DEM carries nodata: neither + algorithm defines an ordering for `NaN`, so which cells a flood reaches, and in what + order, is left to the comparison rules rather than decided. Mask nodata out beforehand if + the choice of method must not change the result. +- `flats=true`: Also give the resulting flats a drainage gradient with + [`resolveflats`](@ref), so that every cell of the output has a strictly lower neighbour + and a defined flow direction. Pass `false` for the filled surface on its own. """ -function filldepressions(dem, mask = fill!(similar(dem, Bool), false)) - filldepressions!(copy(dem), mask) +function filldepressions( + dem, + mask = fill!(similar(dem, Bool), false); + method::DepressionFillingMethod = PriorityFlood(), + flats::Bool = true, +) + filldepressions!(copy(dem), mask; method, flats) end abstract type FlowDirectionMethod end @@ -44,8 +84,20 @@ Base.@kwdef struct FD8 <: FlowDirectionMethod p::Float32 = 1.1 end -function filldepressions!(dem, queued = fill!(similar(dem, Bool), false)) - first_cell = first(eachcell(dem)) +function filldepressions!( + dem, + queued = fill!(similar(dem, Bool), false); + method::DepressionFillingMethod = PriorityFlood(), + flats::Bool = true, +) + _filldepressions!(method, dem, queued) + flats && resolveflats!(dem) + return dem +end + +function _filldepressions!(::PriorityFlood, dem, queued) + R = eachcell(dem) + first_cell = first(R) open = PriorityQueue{typeof(first_cell), eltype(dem)}() pit = DataStructures.Queue{typeof(first_cell)}() @@ -69,6 +121,98 @@ function filldepressions!(dem, queued = fill!(similar(dem, Bool), false)) return dem end +# Zhou et al. (2016). Three queues: `open` holds only potential spill cells, while `trace` +# and `depression` region-grow slope and depression cells on plain queues. A cell reached +# by walking uphill from a processed cell drains back the way it came, so its elevation +# does not depend on visit order and it never needs the priority queue. +function _filldepressions!(::Zhou, dem, processed) + first_cell = first(eachcell(dem)) + Cell = typeof(first_cell) + open = PriorityQueue{Cell, eltype(dem)}() + trace = DataStructures.Queue{Cell}() + depression = DataStructures.Queue{Cell}() + + @inbounds for cell in outlets(dem) + enqueue!(open, cell, dem[cell]) + processed[cell] = true + end + + @inbounds while !isempty(open) + cell = dequeue!(open) + spill = dem[cell] + for ncell in neighbors(dem, cell) + processed[ncell] && continue + processed[ncell] = true + if dem[ncell] <= spill + dem[ncell] = spill + DataStructures.enqueue!(depression, ncell) + _processdepression!(dem, processed, depression, trace) + else + DataStructures.enqueue!(trace, ncell) + end + _processtrace!(dem, processed, trace, open) + end + end + return dem +end + +"""Region-grow a depression at the spill elevation, handing higher cells to `trace`.""" +function _processdepression!(dem, processed, depression, trace) + @inbounds while !isempty(depression) + cell = DataStructures.dequeue!(depression) + spill = dem[cell] + for ncell in neighbors(dem, cell) + processed[ncell] && continue + processed[ncell] = true + if dem[ncell] > spill + DataStructures.enqueue!(trace, ncell) + else + dem[ncell] = spill + DataStructures.enqueue!(depression, ncell) + end + end + end + return dem +end + +""" +Region-grow slope cells uphill on a plain queue. A cell is pushed onto the priority +queue only when it borders unprocessed terrain at or below its own elevation and no +processed neighbour of that terrain is lower, which makes it a potential spill cell. +""" +function _processtrace!(dem, processed, trace, open) + @inbounds while !isempty(trace) + cell = DataStructures.dequeue!(trace) + spill = dem[cell] + queued = false + for ncell in neighbors(dem, cell) + processed[ncell] && continue + nspill = dem[ncell] + if nspill <= spill + if !queued && _isspillcell(dem, processed, ncell, nspill) + enqueue!(open, cell, spill) + queued = true + end + continue + end + processed[ncell] = true + DataStructures.enqueue!(trace, ncell) + end + end + return dem +end + +@inline function _isspillcell(dem, processed, ncell, nspill) + @inbounds for jcell in neighbors(dem, ncell) + processed[jcell] && dem[jcell] < nspill && return false + end + return true +end + +# Non-diagonal neighbors come first: Algorithm 4 of +# [barnesPriorityFloodOptimalDepressionFilling2014](@cite) processes them first on its +# line 13, as they have the greatest centre-to-centre slopes and so take priority when +# several neighbors sit at the same elevation. const nbs = CartesianIndex.([(-1, 0), (0, -1), (0, 1), (1, 0), (-1, -1), (-1, 1), (1, -1), (1, 1)]) @@ -160,6 +304,21 @@ function _directioncode(dem, direction, cellsize) UInt8(Int(FlowDirection{D8D}(direction))) end +""" + _dinfreceivers(dem, cell, state, cellsize) + +The two neighbours [`DInf`](@ref) splits `cell`'s flow between and the fraction each takes, +or `nothing` when it has nowhere to send it. Shared by [`flowaccumulation`](@ref) and +[`flowdirection`](@ref) so the two cannot drift apart. +""" +function _dinfreceivers(dem, cell, state, cellsize) + aspect = _dinf_aspect(state, dem, cell, cellsize) + # Undefined over a flat or beside nodata, and bracketing a NaN would pick the cell + # itself as both of its receivers. + isfinite(aspect) || return nothing + return _flowneighbors(dem, cell, aspect, cellsize) +end + function _flowneighbors(dem, cell, aspect, cellsize) aspect = mod(aspect, 360) lower_cell = upper_cell = cell @@ -283,16 +442,7 @@ function _accumulate!(::DInf, acc, order, dir, R, dem, cellsize) for i in reverse(order) iszero(dir[i]) && continue cell = R[i] - aspect = _dinf_aspect(aspects, dem, cell, cellsize) - - if !isfinite(aspect) - acc[cell + dir[i]] += acc[i] - output[i] = _directioncode(dem, dir[i], cellsize) - visited[i] = true - continue - end - - flowneighbors = _flowneighbors(dem, cell, aspect, cellsize) + flowneighbors = _dinfreceivers(dem, cell, aspects, cellsize) if isnothing(flowneighbors) acc[cell + dir[i]] += acc[i] output[i] = _directioncode(dem, dir[i], cellsize) @@ -407,17 +557,16 @@ function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize) visited = fill!(similar(dem, Bool), false) output = similar(dem, FlowDirection{D8D, UInt8}) fill!(output, 0) + scratch = _fd8scratch(dem) for i in reverse(order) iszero(dir[i]) && continue cell = R[i] - total_weight = 0.0 downstream = cell + dir[i] + total_weight = _fd8weights!(scratch, fd8, dem, cell, visited, cellsize) downstream_weight = 0.0 - for neighbor in neighbors(dem, cell) - weight = _fd8weight(fd8, dem, cell, neighbor, visited, cellsize) - total_weight += weight - neighbor == downstream && (downstream_weight = weight) + for (index, neighbor) in enumerate(scratch.neighbors) + neighbor == downstream && (downstream_weight = scratch.weights[index]) end if iszero(total_weight) || iszero(downstream_weight) @@ -428,8 +577,8 @@ function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize) end directions = zero(UInt8) - for neighbor in neighbors(dem, cell) - weight = _fd8weight(fd8, dem, cell, neighbor, visited, cellsize) + for (index, neighbor) in enumerate(scratch.neighbors) + weight = scratch.weights[index] iszero(weight) && continue acc[neighbor] += acc[i] * (weight / total_weight) directions |= _directioncode(dem, neighbor - cell, cellsize) @@ -440,30 +589,107 @@ function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize) return output end -function _fd8weight(fd8, dem, cell, neighbor, visited, cellsize) - visited[neighbor] && return 0.0 - difference = dem[cell] - dem[neighbor] - (difference < 0 || isnan(difference)) && return 0.0 +""" +Scratch space for [`_fd8weights!`](@ref), reused across cells. + +Bearings are computed once per cell rather than once per neighbour pair, and the contour +lengths they feed depend only on the *shape* of the neighbourhood, not on elevation, so +they are recomputed only when that shape changes — on a raster, only at the edges. +""" +mutable struct FD8Scratch{C, O, T} + neighbors::Vector{C} + offsets::Vector{O} + bearings::Vector{T} + weights::Vector{T} + contours::Vector{T} + distances::Vector{T} + shape::Vector{O} + bitweights::Vector{T} # a cell's shares, gathered against their direction bit +end + +function _fd8scratch(dem) + cell = first(eachcell(dem)) + C, O = typeof(cell), typeof(cell - cell) + FD8Scratch(C[], O[], Float64[], Float64[], Float64[], Float64[], O[], zeros(Float64, 8)) +end + +""" +Recompute the contour length and distance of each neighbour for the neighbourhood shape now +in `scratch.offsets`, after [Quinn et al. (1991)](@cite quinnPredictionHillslopeFlow1991). +""" +function _contourlengths!(scratch::FD8Scratch, dem, cell, cellsize) + n = length(scratch.offsets) + resize!(scratch.contours, n) + resize!(scratch.distances, n) + bearings = scratch.bearings + resize!(bearings, n) + @inbounds for (index, neighbor) in enumerate(scratch.neighbors) + bearings[index] = _cellbearing(dem, cell, neighbor, cellsize) + scratch.distances[index] = celldistance(dem, cell, neighbor; cellsize) + end - distance = celldistance(dem, cell, neighbor; cellsize) - iszero(distance) && return 0.0 - bearing = _cellbearing(dem, cell, neighbor, cellsize) - previous_gap = next_gap = Inf - for other in neighbors(dem, cell) - other == neighbor && continue - other_bearing = _cellbearing(dem, cell, other, cellsize) - candidate_previous_gap = mod(bearing - other_bearing, 360) - candidate_next_gap = mod(other_bearing - bearing, 360) - iszero(candidate_previous_gap) || - (previous_gap = min(previous_gap, candidate_previous_gap)) - iszero(candidate_next_gap) || (next_gap = min(next_gap, candidate_next_gap)) + @inbounds for index in 1:n + bearing = bearings[index] + previous_gap = next_gap = Inf + for other in 1:n + other == index && continue + candidate_previous_gap = mod(bearing - bearings[other], 360) + candidate_next_gap = mod(bearings[other] - bearing, 360) + iszero(candidate_previous_gap) || + (previous_gap = min(previous_gap, candidate_previous_gap)) + iszero(candidate_next_gap) || (next_gap = min(next_gap, candidate_next_gap)) + end + # A gap wider than a right angle means neighbours are missing, at the edge of the + # DEM or beside nodata. Clamp it: the formula otherwise runs the contour out to + # tand(90) = Inf, and an infinite weight takes the whole cell's flow. + scratch.contours[index] = + (isfinite(previous_gap) && isfinite(next_gap)) ? + scratch.distances[index] / 2 * + (tand(min(previous_gap, 90) / 2) + tand(min(next_gap, 90) / 2)) : 0.0 end - (!isfinite(previous_gap) || !isfinite(next_gap)) && return 0.0 - contour_length = - distance / 2 * (tand(previous_gap / 2) + tand(next_gap / 2)) - contour_length > 0 || return 0.0 - return (difference / distance * contour_length)^fd8.p + resize!(scratch.shape, n) + copyto!(scratch.shape, scratch.offsets) + return scratch +end + +""" + _fd8weights!(scratch, fd8, dem, cell, visited, cellsize) + +Fill `scratch` with the neighbours of `cell` and their unnormalised FD8 weights, and return +the total. A neighbour weighs the downslope gradient towards it times the length of the +cell's contour facing it, raised to `fd8.p`, after +[Quinn et al. (1991)](@cite quinnPredictionHillslopeFlow1991). Neighbours that are not +lower, or that `visited` marks as already processed, weigh nothing. +""" +function _fd8weights!(scratch::FD8Scratch, fd8, dem, cell, visited, cellsize) + empty!(scratch.neighbors) + empty!(scratch.offsets) + @inbounds for neighbor in neighbors(dem, cell) + push!(scratch.neighbors, neighbor) + push!(scratch.offsets, neighbor - cell) + end + scratch.offsets == scratch.shape || _contourlengths!(scratch, dem, cell, cellsize) + resize!(scratch.weights, length(scratch.neighbors)) + + center = dem[cell] + total = 0.0 + @inbounds for (index, neighbor) in enumerate(scratch.neighbors) + weight = 0.0 + difference = center - dem[neighbor] + # No drop, no flow. Skipping here also keeps a zero difference away from a contour + # length of zero, which a truncated neighbourhood can produce. + if !visited[neighbor] && difference > 0 && !isnan(difference) + distance = scratch.distances[index] + contour_length = scratch.contours[index] + if distance > 0 && contour_length > 0 + weight = (difference / distance * contour_length)^fd8.p + total += weight + end + end + scratch.weights[index] = weight + end + return total end """ @@ -684,7 +910,7 @@ function _hand!(output, order, flowdirs, R, filled, stream_mask, cellsize) end """ - depression_depth(dem; filled=filldepressions(dem)) + depression_depth(dem; filled=filldepressions(dem; flats=false)) Computes the depth of each cell below the filled surface. @@ -695,13 +921,13 @@ depressions will have a depth of zero. This is useful for identifying potential cold air pooling zones and water retention areas. """ -function depression_depth(dem; filled = filldepressions(dem)) +function depression_depth(dem; filled = filldepressions(dem; flats=false)) filled .- dem end """ - depression_volume(dem::AbstractMatrix; filled=filldepressions(dem), cellsize=cellsize(dem)) - depression_volume(dem; filled=filldepressions(dem), cellsize=cellsize(dem)) + depression_volume(dem::AbstractMatrix; filled=filldepressions(dem; flats=false), cellsize=cellsize(dem)) + depression_volume(dem; filled=filldepressions(dem; flats=false), cellsize=cellsize(dem)) Computes the total volume of all depressions/basins in the DEM. @@ -709,14 +935,14 @@ Returns the sum of depression depths multiplied by cell area. """ function depression_volume( dem::AbstractMatrix; - filled = filldepressions(dem), + filled = filldepressions(dem; flats = false), cellsize = cellsize(dem), ) sum(filter(isfinite, depression_depth(dem; filled))) * prod(abs.(cellsize)) end function depression_volume( dem; - filled = filldepressions(dem), + filled = filldepressions(dem; flats = false), cellsize = cellsize(dem), ) depth = depression_depth(dem; filled) diff --git a/test/hydrology.jl b/test/hydrology.jl new file mode 100644 index 0000000..a8b22b1 --- /dev/null +++ b/test/hydrology.jl @@ -0,0 +1,140 @@ +using Geomorphometry +using Geomorphometry: PriorityFlood, Zhou, resolveflats, flowdirection, maxreceivers +using Test +using Random + +"""Random DEM with broad features, so it contains real depressions and flats.""" +function rough_dem(n, m; seed = 42, levels = 6, quantize = nothing) + rng = Random.MersenneTwister(seed) + dem = zeros(Float32, n, m) + for k in 1:levels + f = 2.0^k + a = randn(rng, Float32, ceil(Int, n / f) + 2, ceil(Int, m / f) + 2) + for i in 1:n, j in 1:m + dem[i, j] += a[ceil(Int, i / f), ceil(Int, j / f)] / k + end + end + isnothing(quantize) || (dem = round.(dem ./ quantize) .* quantize) + dem +end + +"""Directions alone, dropping the weights.""" +directionsof(dem; method, cellsize = (1.0, 1.0)) = + first(flowdirection(dem; method, cellsize)) + +const unitcell = (1.0, 1.0) + +@testset "filldepressions" begin + @testset "Zhou reaches the same surface as PriorityFlood" begin + # Different queue strategies, same definition of a filled surface. Quantizing + # makes large flats, where the two are most likely to diverge. + for dem in (rough_dem(64, 64), rough_dem(80, 80; quantize = 0.25f0)) + @test filldepressions(dem; method = Zhou(), flats = false) == + filldepressions(dem; method = PriorityFlood(), flats = false) + end + end + + @testset "the filled surface has no depressions left" begin + filled = filldepressions(rough_dem(50, 50); method = Zhou(), flats = false) + @test all(filled .>= rough_dem(50, 50)) + interior = CartesianIndices((2:49, 2:49)) + offsets = filter(!iszero, CartesianIndices((-1:1, -1:1))) + @test all(any(filled[I + o] <= filled[I] for o in offsets) for I in interior) + end + + @testset "flats are resolved by default" begin + dem = rough_dem(40, 40; quantize = 0.5f0) + @test filldepressions(dem) != filldepressions(dem; flats = false) + # Resolving only nudges by float steps, so the surface is unchanged to Float32. + @test filldepressions(dem) ≈ filldepressions(dem; flats = false) + end +end + +@testset "resolveflats" begin + @testset "a tilted plane has no flats" begin + dem = Float32[100.0f0 - 2i - 3j for i in 1:10, j in 1:10] + @test resolveflats(dem) == dem + end + + @testset "a flat drains towards its outlet" begin + # A 2x2 flat at 5, draining at (3,3) towards the 4 at (4,4). + dem = Float32[9 9 9 9; 9 5 5 9; 9 5 5 9; 9 9 9 4] + z = resolveflats(dem) + flat = dem .== 5 + @test z[.!flat] == dem[.!flat] # only flat cells move + @test minimum(z[flat]) == z[3, 3] # the outlet cell is the lowest of the flat + end + + @testset "every flat cell drains afterwards" begin + z = filldepressions(rough_dem(60, 60; quantize = 0.5f0)) + interior = CartesianIndices((2:59, 2:59)) + offsets = filter(!iszero, CartesianIndices((-1:1, -1:1))) + @test all(any(z[I + o] < z[I] for o in offsets) for I in interior) + end +end + +@testset "flowdirection" begin + @testset "points downslope on a tilted plane" begin + dem = Float32[100.0f0 - 2i - 3j for i in 1:8, j in 1:8] + dirs = directionsof(dem; method = D8()) + @test all( + CartesianIndex(dirs[I]) == CartesianIndex(1, 1) for + I in CartesianIndices((2:7, 2:7)) + ) + end + + @testset "takes the steepest neighbour, not the lowest" begin + # E drops 1.0 over 1.0; NE drops 1.2 over sqrt(2), so it is lower but shallower. + dem = Float32[10 10 10; 10 10 9; 8.8 10 10] + @test CartesianIndex(directionsof(dem; method = D8())[2, 2]) == + CartesianIndex(0, 1) + end + + @testset "a cell with no lower neighbour is a pit" begin + @test all(Geomorphometry.ispit, directionsof(fill(5.0f0, 5, 5); method = D8())) + end + + @testset "receivers stay within what the method allows" begin + dem = rough_dem(40, 40) + for method in (D8(), DInf(), FD8()) + dirs = directionsof(dem; method) + @test all(d -> Geomorphometry.ndirections(d) <= maxreceivers(method), dirs) + end + # A peak sends flow every way FD8 can, and DInf splits in two somewhere. + peak = Float32[1 1 1; 1 2 1; 1 1 1] + @test Geomorphometry.ndirections(directionsof(peak; method = FD8())[2, 2]) == 8 + @test any(d -> Geomorphometry.ndirections(d) == 2, directionsof(dem; method = DInf())) + end + + @testset "all methods agree on which cells are pits" begin + dem = rough_dem(40, 40; quantize = 0.5f0) + pits(method) = Geomorphometry.ispit.(directionsof(dem; method)) + @test pits(D8()) == pits(DInf()) == pits(FD8()) + end + + @testset "weights sum to one wherever a cell drains" begin + dem = rough_dem(40, 40) + for method in (D8(), DInf(), FD8()) + dirs, weights = flowdirection(dem; method, cellsize = unitcell) + for cell in CartesianIndices(dem) + n = Geomorphometry.ndirections(dirs[cell]) + iszero(n) && continue + @test sum(weights[cell, slot] for slot in 1:n) ≈ 1 + end + end + end + + @testset "weights are finite on a truncated neighbourhood" begin + # At the border a cell is missing neighbours, leaving an angular gap of 180 + # degrees. Taken literally the contour length runs to tand(90) = Inf, and that + # neighbour would take the cell's whole flow. + _, weights = flowdirection(rough_dem(20, 20); method = FD8(), cellsize = unitcell) + @test all(isfinite, weights) + end + + @testset "a narrower weight type still routes the flow" begin + dirs, weights = flowdirection(rough_dem(20, 20); method = FD8(), weighttype = Float16) + @test eltype(weights) == Float16 + @test all(isfinite, weights) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 6f436ac..1313cca 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,6 +37,7 @@ end end # Must run after ExplicitImports include("horizon.jl") + include("hydrology.jl") @testset "pmf" begin # Write your own tests here. A = rand(25, 25) @@ -308,8 +309,11 @@ end @test maximum(acc) <= length(conservation_dem) end + # Filling levels the bowl; the default also gives the resulting flat a gradient + # of one float step per cell, so the surface is level only to within that. fdem = filldepressions(dem) - @test all(==(10), fdem) + @test all(≈(10), fdem) + @test all(==(10), filldepressions(dem; flats = false)) @testset "depression_depth" begin bd = depression_depth(dem) From 762e549c9c03a2b17e7a78cdd8a1b894ca1ba43e Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Mon, 9 Nov 2026 21:52:12 +0100 Subject: [PATCH 2/3] Allocate direction and weight grids through hooks, so wrappers keep their 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) Claude-Session: https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf --- ext/GeomorphometryRastersExt.jl | 6 ++++++ src/direction.jl | 11 +++++++---- src/hydrology.jl | 15 +++++++++++---- test/hydrology.jl | 11 +++++++++++ 4 files changed, 35 insertions(+), 8 deletions(-) diff --git a/ext/GeomorphometryRastersExt.jl b/ext/GeomorphometryRastersExt.jl index ecc477e..526c8c4 100644 --- a/ext/GeomorphometryRastersExt.jl +++ b/ext/GeomorphometryRastersExt.jl @@ -52,4 +52,10 @@ end Geomorphometry.outlets(r::Raster) = Geomorphometry.outlets(parent(r)) Geomorphometry.neighbors(r::Raster, cell::CartesianIndex{2}) = Geomorphometry.neighbors(parent(r), cell) +# A direction eltype cannot represent the source missingval, so allocate without one. +Geomorphometry._directiongrid( + dem::Raster, ::Type{C} +) where {C <: Geomorphometry.FlowDirectionConvention} = + similar(dem, FlowDirection{C, UInt8}; missingval = nothing) + end # module diff --git a/src/direction.jl b/src/direction.jl index e538823..eaefece 100644 --- a/src/direction.jl +++ b/src/direction.jl @@ -20,9 +20,12 @@ Share of each cell's flow per receiver, one dimension larger than the DEM, the t axis in [`D8D`](@ref) bit order. A single-receiver method stores nothing: its weights are a [`Fill`](@extref FillArrays.Fill) and cost no memory. """ -_weights(method::FlowDirectionMethod, dem, T) = - hasweights(method) ? zeros(T, size(dem)..., maxreceivers(method)) : - Fill(one(T), size(dem)..., maxreceivers(method)) +function _weights(method::FlowDirectionMethod, dem, T) + n = maxreceivers(method) + # A single receiver takes everything, so there is nothing to store. + hasweights(method) || return Fill(one(T), size(dem)..., n) + return fill!(_alloc_directions(dem, T, n), zero(T)) +end """Record the share of a cell's `slot`-th receiver. Nothing to record without weights.""" @@ -61,7 +64,7 @@ function flowdirection( weighttype::Type{<:AbstractFloat} = Float32, ) C = _convention(method) - directions = similar(dem, FlowDirection{C, UInt8}) + directions = _directiongrid(dem, C) weights = _weights(method, dem, weighttype) state = _directionstate(method, dem, cellsize) @inbounds for cell in eachcell(dem) diff --git a/src/hydrology.jl b/src/hydrology.jl index 2fb1da0..77f31db 100644 --- a/src/hydrology.jl +++ b/src/hydrology.jl @@ -424,19 +424,26 @@ function flowaccumulation!( return acc, dirs end +""" +Allocate a direction grid shaped like `dem`. A hook, because a direction eltype has no +representation for a source `missingval` and some grid types must drop theirs. +""" +_directiongrid(dem, ::Type{C}) where {C <: FlowDirectionConvention} = + similar(dem, FlowDirection{C, UInt8}) + function _accumulate!(::D8, acc, order, dir, R, dem, cellsize) for i in reverse(order) iszero(dir[i]) && continue acc[R[i] + dir[i]] += acc[i] end - output = similar(dem, FlowDirection{LDD, UInt8}) + output = _directiongrid(dem, LDD) output .= FlowDirection{LDD}.(_orient.(dir, Ref(cellsize))) return output end function _accumulate!(::DInf, acc, order, dir, R, dem, cellsize) aspects = _dinf_aspects(dem, cellsize) visited = fill!(similar(dem, Bool), false) - output = similar(dem, FlowDirection{D8D, UInt8}) + output = _directiongrid(dem, D8D) fill!(output, 0) for i in reverse(order) @@ -509,7 +516,7 @@ function _accumulate!(fd8::FD8, acc, order, dir, R::CartesianIndices{2}, dem, ce visited = falses(size(acc)) nb = vec(collect(CartesianIndices(dists)) .- CartesianIndex(2, 2)) - output = similar(dem, FlowDirection{D8D, UInt8}) + output = _directiongrid(dem, D8D) fill!(output, 0) weights = zeros(size(contour_lengths)) @@ -555,7 +562,7 @@ end function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize) visited = fill!(similar(dem, Bool), false) - output = similar(dem, FlowDirection{D8D, UInt8}) + output = _directiongrid(dem, D8D) fill!(output, 0) scratch = _fd8scratch(dem) diff --git a/test/hydrology.jl b/test/hydrology.jl index a8b22b1..d535495 100644 --- a/test/hydrology.jl +++ b/test/hydrology.jl @@ -138,3 +138,14 @@ end @test all(isfinite, weights) end end + +@testset "direction and weight grids follow the input type" begin + # A direction eltype has no representation for a source missingval, and the weights + # need a direction dimension that `similar` would otherwise drop. Both go through + # hooks the wrapper extensions override. + @test Geomorphometry._directiongrid(zeros(Float32, 3, 3), Geomorphometry.LDD) isa + Matrix{<:Geomorphometry.FlowDirection} + dirs, weights = flowdirection(rough_dem(8, 8); method = FD8(), cellsize = unitcell) + @test size(weights) == (8, 8, maxreceivers(FD8())) + @test all(>=(0), weights) +end From 0dd6635d8b23396fa59f86348ace48762ed9d141 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Tue, 10 Nov 2026 14:12:17 +0100 Subject: [PATCH 3/3] Address review on filling, flats and flow directions 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) Claude-Session: https://claude.ai/code/session_01C4DNRaTBnCQ2KF2aCbq6Nf --- docs/src/usage.md | 27 +++++------ src/direction.jl | 112 +++++++++++++++++++++++----------------------- src/flats.jl | 66 ++++++++++++++++++++------- src/hydrology.jl | 78 ++++++++++++++++++++++++-------- test/hydrology.jl | 50 +++++++++++++++++++++ 5 files changed, 228 insertions(+), 105 deletions(-) diff --git a/docs/src/usage.md b/docs/src/usage.md index 7197c2a..72060b1 100644 --- a/docs/src/usage.md +++ b/docs/src/usage.md @@ -383,26 +383,21 @@ depression_volume(dtm; filled=fdtm) [`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: +faster. Filling leaves flats, which have no flow direction, so by default it also gives +them a gradient with [`resolveflats`](@ref) after +[barnesEfficientAssignmentDrainage2014](@citet); `flats=false` gives the filled surface +alone. [`flowdirection`](@ref) then reads directions off that surface, with the share each +receiver takes beside them. ```@example plots -z = filldepressions(cdtm) # fill depressions, resolve flats -directions, weights = flowdirection(z; method=D8()) # directions off the filled surface -nothing # hide +z = filldepressions(cdtm) +directions, weights = flowdirection(z; method=D8()) +parent(directions[400:800, 300:700]) # arrows, as the drainage direction below ``` -[`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. +[`D8`](@ref) has one receiver that takes everything, so its weights are a `Fill` costing +nothing at any grid size. [`DInf`](@ref) and [`FD8`](@ref) store a share per receiver, two +and eight, which `weighttype=Float16` halves. Filling the depressions in a DEM is not necessary to calculate the flow accumulation. Here we use [`flowaccumulation`](@ref) to do so, and it automatically carves out depressions. Note that the local drainage direction is also returned. By default the FD8 algorithm is used, but you can also use the D∞ or D8 algorithm by setting the `method` keyword argument to `DInf()` or `D8()`. diff --git a/src/direction.jl b/src/direction.jl index eaefece..18589ee 100644 --- a/src/direction.jl +++ b/src/direction.jl @@ -1,8 +1,8 @@ # Single-receiver methods use the same LDD encoding as `flowaccumulation`; methods that can # send flow several ways need the bitwise D8D encoding. -_convention(::D8) = LDD -_convention(::DInf) = D8D -_convention(::FD8) = D8D +convention(::D8) = LDD +convention(::DInf) = D8D +convention(::FD8) = D8D # A method that can send flow several ways needs the share each receiver takes; a single # receiver takes everything, so there is nothing to carry. @@ -20,7 +20,7 @@ Share of each cell's flow per receiver, one dimension larger than the DEM, the t axis in [`D8D`](@ref) bit order. A single-receiver method stores nothing: its weights are a [`Fill`](@extref FillArrays.Fill) and cost no memory. """ -function _weights(method::FlowDirectionMethod, dem, T) +function allocweights(method::FlowDirectionMethod, dem, T) n = maxreceivers(method) # A single receiver takes everything, so there is nothing to store. hasweights(method) || return Fill(one(T), size(dem)..., n) @@ -29,14 +29,18 @@ end """Record the share of a cell's `slot`-th receiver. Nothing to record without weights.""" -@inline _setshare!(::Fill, cell, slot, fraction) = fraction -@inline _setshare!(weights::AbstractArray{T}, cell, slot, fraction) where {T} = +@inline setshare!(::Fill, cell, slot, fraction) = fraction +@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. -_directionstate(::FlowDirectionMethod, dem, cellsize) = nothing -_directionstate(::DInf, dem, cellsize) = _dinf_aspects(dem, cellsize) -_directionstate(::FD8, dem, cellsize) = _fd8scratch(dem) +""" +Whatever a method wants to set up once for the whole DEM rather than per cell: a grid of +aspects for [`DInf`](@ref) off a raster, and for [`FD8`](@ref) the scratch buffers it +reuses from cell to cell. Methods that need neither get `nothing`. +""" +directionstate(::FlowDirectionMethod, dem, cellsize) = nothing +directionstate(::DInf, dem, cellsize) = _dinf_aspects(dem, cellsize) +directionstate(::FD8, dem, cellsize) = _fd8scratch(dem) """ flowdirection(dem; method=D8(), cellsize=cellsize(dem), weighttype=Float32) @@ -63,12 +67,12 @@ function flowdirection( cellsize = cellsize(dem), weighttype::Type{<:AbstractFloat} = Float32, ) - C = _convention(method) + C = convention(method) directions = _directiongrid(dem, C) - weights = _weights(method, dem, weighttype) - state = _directionstate(method, dem, cellsize) + weights = allocweights(method, dem, weighttype) + state = directionstate(method, dem, cellsize) @inbounds for cell in eachcell(dem) - code = _flowdirection!(weights, method, state, dem, cell, cellsize) + code = celldirection!(weights, method, state, dem, cell, cellsize) directions[cell] = FlowDirection{C}(code) end return directions, weights @@ -81,7 +85,7 @@ Offset towards `cell`'s steepest downslope neighbour, or zero when it has none: [Jenson and Domingue (1988)](@cite jensonExtractingTopographicStructure1988), and the fallback wherever a method leaves a cell undecided. """ -function _steepestdescent(dem, cell, cellsize) +function steepestdescent(dem, cell, cellsize) offset = cell - cell z = dem[cell] isfinite(z) || return offset @@ -100,55 +104,51 @@ function _steepestdescent(dem, cell, cellsize) return offset end -_steepestcode(dem, cell, cellsize) = - _directioncode(dem, _steepestdescent(dem, cell, cellsize), cellsize) +steepestcode(dem, cell, cellsize) = + _directioncode(dem, steepestdescent(dem, cell, cellsize), cellsize) -_flowdirection!(weights, ::D8, _, dem, cell, cellsize) = - _orient(_steepestdescent(dem, cell, cellsize), cellsize) +celldirection!(weights, ::D8, _, dem, cell, cellsize) = + _orient(steepestdescent(dem, cell, cellsize), cellsize) # DInf: the two neighbours bounding the steepest facet, using the same `_dinfreceivers` # partition as `flowaccumulation`. -function _flowdirection!(weights, ::DInf, state, dem, cell, cellsize) +function celldirection!(weights, ::DInf, state, dem, cell, cellsize) z = dem[cell] isfinite(z) || return zero(UInt8) - flowneighbors = _dinfreceivers(dem, cell, state, cellsize) - isnothing(flowneighbors) && return _steepestshare!(weights, dem, cell, cellsize) - 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. - takesa = afraction > 0 && isfinite(dem[acell]) && dem[acell] < z - takesb = bfraction > 0 && isfinite(dem[bcell]) && dem[bcell] < z - total = (takesa ? afraction : 0.0) + (takesb ? bfraction : 0.0) - # The partition can name terrain that is not lower, most often beside a flat. Fall back - # rather than leave a cell that can drain without a direction. - total > 0 || return _steepestshare!(weights, dem, cell, cellsize) - - abit = takesa ? _directioncode(dem, acell - cell, cellsize) : zero(UInt8) - bbit = takesb ? _directioncode(dem, bcell - cell, cellsize) : zero(UInt8) - # Shares are stored in the order the directions come back out, which is by bit - if iszero(bbit) - _setshare!(weights, cell, 1, afraction / total) - elseif iszero(abit) - _setshare!(weights, cell, 1, bfraction / total) - elseif abit < bbit - _setshare!(weights, cell, 1, afraction / total) - _setshare!(weights, cell, 2, bfraction / total) - else - _setshare!(weights, cell, 1, bfraction / total) - _setshare!(weights, cell, 2, afraction / total) + receivers = _dinfreceivers(dem, cell, state, cellsize) + isnothing(receivers) && return steepestshare!(weights, dem, cell, cellsize) + acell, bcell, afraction, bfraction = receivers + + # The partition can name a cell that is not actually lower, most often beside a flat. + # Such a receiver takes nothing and the other is renormalised to the whole; if that + # leaves nobody, fall back rather than strand a cell that can drain. + lower(ncell) = isfinite(dem[ncell]) && dem[ncell] < z + a = (_directioncode(dem, acell - cell, cellsize), lower(acell) ? afraction : 0.0) + b = (_directioncode(dem, bcell - cell, cellsize), lower(bcell) ? bfraction : 0.0) + total = last(a) + last(b) + total > 0 || return steepestshare!(weights, dem, cell, cellsize) + + # Shares are read back in the order `decompose` hands the directions out, which is by + # bit, so emit the pair in that order. + code = zero(UInt8) + slot = 0 + @inbounds for (bit, fraction) in (first(a) <= first(b) ? (a, b) : (b, a)) + iszero(fraction) && continue + slot += 1 + setshare!(weights, cell, slot, fraction / total) + code |= bit end - return abit | bbit + return code end # FD8: every downslope neighbour receives flow, weighted by slope and contour length, using # the same `_fd8weights!` as `flowaccumulation`. -function _flowdirection!(weights, fd8::FD8, scratch, dem, cell, cellsize) +function celldirection!(weights, fd8::FD8, scratch, dem, cell, cellsize) z = dem[cell] isfinite(z) || return zero(UInt8) - total = _fd8weights!(scratch, fd8, dem, cell, _Unvisited(), cellsize) - iszero(total) && return _steepestshare!(weights, dem, cell, cellsize) + total = _fd8weights!(scratch, fd8, dem, cell, Unvisited(), cellsize) + iszero(total) && return steepestshare!(weights, dem, cell, cellsize) # Gather the shares against their direction bits, then emit them in bit order, which is # the order the directions come back out in. @@ -165,18 +165,18 @@ function _flowdirection!(weights, fd8::FD8, scratch, dem, cell, cellsize) @inbounds for bit in decompose(D8D, code) iszero(bit) && continue slot += 1 - _setshare!(weights, cell, slot, scratch.bitweights[trailing_zeros(Int(bit)) + 1]) + setshare!(weights, cell, slot, scratch.bitweights[trailing_zeros(Int(bit)) + 1]) end return code end """The steepest descent fallback: one receiver, taking everything.""" -function _steepestshare!(weights, dem, cell, cellsize) - code = _steepestcode(dem, cell, cellsize) - iszero(code) || _setshare!(weights, cell, 1, one(Float32)) +function steepestshare!(weights, dem, cell, cellsize) + code = steepestcode(dem, cell, cellsize) + iszero(code) || setshare!(weights, cell, 1, one(Float32)) return code end """A stencil has no traversal order, so nothing is excluded for being already visited.""" -struct _Unvisited end -Base.getindex(::_Unvisited, _...) = false +struct Unvisited end +Base.getindex(::Unvisited, _...) = false diff --git a/src/flats.jl b/src/flats.jl index e6e8a55..3daf990 100644 --- a/src/flats.jl +++ b/src/flats.jl @@ -1,10 +1,11 @@ """ - resolveflats(dem) - resolveflats!(dem) + resolveflats(dem, excluded=nothing) + resolveflats!(dem, excluded=nothing) Give the flats of `dem` a drainage gradient, following [Barnes et al. (2014)](@cite barnesEfficientAssignmentDrainage2014), and return the -elevation model. `resolveflats!` writes into `dem` instead of a copy. +elevation model. `resolveflats!` writes into `dem` instead of a copy. Cells marked in +`excluded` keep their elevation and take no part in a flat. A flat is a connected region of equal elevation, over which no flow direction is defined. The gradient is the superposition of one away from higher terrain, which gives convergent @@ -13,12 +14,13 @@ twice the weight, so it settles cells where the two disagree without iterating. Cells are raised by whole `nextfloat` steps, so the surface is unchanged at any meaningful precision but every flat cell gains a strictly lower neighbour. A flat with no outlet is a -depression, not a flat, and is left alone. [`filldepressions`](@ref) calls this by default. +depression, not a flat, and is left alone. [`filldepressions`](@ref) calls this by default, +which needs a floating point `dem`: integer elevations cannot express the gradient. """ -resolveflats(dem) = resolveflats!(copy(dem)) +resolveflats(dem, excluded = nothing) = resolveflats!(copy(dem), excluded) -function resolveflats!(dem) - increments = _flatmask(dem) +function resolveflats!(dem, excluded = nothing) + increments = flatmask(dem, excluded) @inbounds for cell in eachcell(dem) increment = increments[cell] iszero(increment) || (dem[cell] = nextfloat(dem[cell], increment)) @@ -26,11 +28,29 @@ function resolveflats!(dem) return dem end +"""Was `cell` excluded by the caller? Nothing is excluded when there is no mask.""" +@inline isexcluded(::Nothing, cell) = false +@inline isexcluded(excluded, cell) = @inbounds excluded[cell] + +""" +The largest part of `increment` that keeps `z` strictly below `limit`. `nextfloat(z, n)` is +a constant-time step over `n` representable values, so the usual case of ample headroom +costs one comparison, and the search that follows is bounded by `increment` itself. +""" +function cappedincrement(z, limit, increment) + nextfloat(z, Int(increment)) < limit && return increment + steps = increment + while steps > zero(steps) && !(nextfloat(z, Int(steps)) < limit) + steps -= one(steps) + end + return steps +end + """ Per-cell increment counts for [`resolveflats`](@ref): zero wherever a cell already drains, positive on cells that have no strictly lower neighbour and belong to a flat that drains. """ -function _flatmask(dem) +function flatmask(dem, excluded = nothing) R = eachcell(dem) Cell = typeof(first(R)) @@ -50,7 +70,7 @@ function _flatmask(dem) # drain but touches higher terrain. @inbounds for cell in R z = dem[cell] - if !isfinite(z) + if !isfinite(z) || isexcluded(excluded, cell) drains[cell] = true continue end @@ -110,7 +130,7 @@ function _flatmask(dem) # edges. `flatheights` is how far each flat reaches, which Step 3 needs to invert it. flatheights = zeros(Int32, nlabels) seen = fill!(similar(dem, Bool), false) - _seed!(seen, highedges) + markseen!(seen, highedges) increment = one(Int32) frontier = copy(highedges) next = Cell[] @@ -119,7 +139,7 @@ function _flatmask(dem) for cell in frontier flatmask[cell] = increment flatheights[labels[cell]] = increment - _growflat!(next, seen, dem, labels, drains, cell) + growflat!(next, seen, dem, labels, drains, cell) end frontier, next = next, frontier increment += one(increment) @@ -132,14 +152,14 @@ function _flatmask(dem) flatmask[cell] = -flatmask[cell] end fill!(seen, false) - _seed!(seen, lowedges) + markseen!(seen, lowedges) increment = one(Int32) frontier = copy(lowedges) empty!(next) @inbounds while !isempty(frontier) empty!(next) for cell in frontier - _growflat!(next, seen, dem, labels, drains, cell) + growflat!(next, seen, dem, labels, drains, cell) flatmask[cell] < 0 && (flatmask[cell] += flatheights[labels[cell]]) flatmask[cell] += 2 * increment end @@ -152,10 +172,26 @@ function _flatmask(dem) @inbounds for cell in R drains[cell] && (flatmask[cell] = 0) end + + # The increments assume the surrounding relief dwarfs a few float steps. Where it does + # not, raising a cell could carry it over a neighbour that was strictly higher and take + # away that neighbour's drainage, so each one is capped below its lowest higher + # neighbour. + @inbounds for cell in R + iszero(flatmask[cell]) && continue + z = dem[cell] + limit = typemax(float(typeof(z))) + for ncell in neighbors(dem, cell) + nz = dem[ncell] + (isfinite(nz) && nz > z) && (limit = min(limit, nz)) + end + isinf(limit) && continue + flatmask[cell] = cappedincrement(z, limit, flatmask[cell]) + end return flatmask end -function _seed!(seen, cells) +function markseen!(seen, cells) for cell in cells seen[cell] = true end @@ -163,7 +199,7 @@ function _seed!(seen, cells) end """Queue the not-yet-visited cells of `cell`'s own flat that still need a gradient.""" -@inline function _growflat!(next, seen, dem, labels, drains, cell) +@inline function growflat!(next, seen, dem, labels, drains, cell) label = labels[cell] @inbounds for ncell in neighbors(dem, cell) if !seen[ncell] && labels[ncell] == label && !drains[ncell] diff --git a/src/hydrology.jl b/src/hydrology.jl index 77f31db..c8c2b2c 100644 --- a/src/hydrology.jl +++ b/src/hydrology.jl @@ -52,19 +52,25 @@ Fill the depressions of a digital elevation model, so that every cell drains to algorithm defines an ordering for `NaN`, so which cells a flood reaches, and in what order, is left to the comparison rules rather than decided. Mask nodata out beforehand if the choice of method must not change the result. -- `flats=true`: Also give the resulting flats a drainage gradient with - [`resolveflats`](@ref), so that every cell of the output has a strictly lower neighbour - and a defined flow direction. Pass `false` for the filled surface on its own. +- `flats`: Also give the resulting flats a drainage gradient with [`resolveflats`](@ref), + so that every cell of the output has a strictly lower neighbour and a defined flow + direction. Pass `false` for the filled surface on its own. Defaults to `true` for a + floating point `dem` and `false` otherwise, as integer elevations cannot express a + gradient of single float steps. Cells excluded by `mask` keep their elevation. """ function filldepressions( dem, mask = fill!(similar(dem, Bool), false); method::DepressionFillingMethod = PriorityFlood(), - flats::Bool = true, + flats::Bool = _canresolveflats(dem), ) filldepressions!(copy(dem), mask; method, flats) end +# Resolving flats nudges elevations by single float steps, which an integer DEM cannot +# represent, so it is only on by default where it can work. +_canresolveflats(dem) = eltype(dem) <: AbstractFloat + abstract type FlowDirectionMethod end """D8 Flow Direction method by [Jenson (1988)](@cite jensonExtractingTopographicStructure1988).""" @@ -88,10 +94,13 @@ function filldepressions!( dem, queued = fill!(similar(dem, Bool), false); method::DepressionFillingMethod = PriorityFlood(), - flats::Bool = true, + flats::Bool = _canresolveflats(dem), ) + # The flood writes into `queued`, so the caller's exclusions are kept to hand back to + # `resolveflats!`, which must leave those cells alone too. + excluded = flats ? copy(queued) : nothing _filldepressions!(method, dem, queued) - flats && resolveflats!(dem) + flats && resolveflats!(dem, excluded) return dem end @@ -345,8 +354,10 @@ function _flowneighbors(dem, cell, aspect, cellsize) return lower_cell, upper_cell, upper_gap / gap, lower_gap / gap end +# The signs matter: `_cellbearing` keeps them, so dropping them here would leave the +# aspect and the receiver bearings in different frames and pick the wrong neighbour. _dinf_aspects(dem::AbstractMatrix, cellsize) = - aspect(dem; method = Horn(), cellsize = abs.(cellsize)) + aspect(dem; method = Horn(), cellsize) _dinf_aspects(dem, cellsize) = nothing _dinf_aspect(aspects, dem, cell, cellsize) = aspects[cell] _dinf_aspect(::Nothing, dem, cell, cellsize) = _localaspect(dem, cell, cellsize) @@ -496,17 +507,11 @@ function _accumulate!(::DInf, acc, order, dir, R, dem, cellsize) end function _accumulate!(fd8::FD8, acc, order, dir, R::CartesianIndices{2}, dem, cellsize) - # Derive contour lengths, which is used to calculate the weights - # Uses the algorithm by Quinn et al. (1991), L1=0.5 L2=0.354 for δx=δy=1 - # TODO Check whether just using the angles is enough. - δx, δy = abs.(cellsize) - δxy = sqrt((δx / 4)^2 + (δy / 4)^2) - contour_lengths = @SMatrix [ - δxy δx δxy - δy 0 δy - δxy δx δxy - ] + # Contour lengths come from `_rastercontour`, which `flowdirection` uses too. + offsets = CartesianIndices((-1:1, -1:1)) + contour_lengths = map(o -> iszero(o) ? 0.0 : _rastercontour(o, cellsize), offsets) + δx, δy = abs.(cellsize) δxy = sqrt(δx^2 + δy^2) dists = @SMatrix [ δxy δx δxy @@ -614,6 +619,24 @@ mutable struct FD8Scratch{C, O, T} bitweights::Vector{T} # a cell's shares, gathered against their direction bit end +""" + _rastercontour(offset, cellsize) + +Contour length a step of `offset` crosses on a square grid, after +[Quinn et al. (1991)](@cite quinnPredictionHillslopeFlow1991). + +The eight Moore neighbours are not a Voronoi neighbourhood — a diagonal shares only a +corner, so there is no face there to measure — and the angular construction in +[`_contourlengths!`](@ref) degenerates on them, giving every direction the same weight. This +is the fix-up for that, and both [`flowdirection`](@ref) and [`flowaccumulation`](@ref) take +their raster contour lengths from here so the two cannot disagree. +""" +function _rastercontour(offset, cellsize) + δx, δy = abs.(cellsize) + (!iszero(offset[1]) && !iszero(offset[2])) && return hypot(δx, δy) / 4 + return iszero(offset[2]) ? δx : δy +end + function _fd8scratch(dem) cell = first(eachcell(dem)) C, O = typeof(cell), typeof(cell - cell) @@ -622,7 +645,9 @@ end """ Recompute the contour length and distance of each neighbour for the neighbourhood shape now -in `scratch.offsets`, after [Quinn et al. (1991)](@cite quinnPredictionHillslopeFlow1991). +in `scratch.offsets`. Each neighbour is given the Voronoi edge it would share with `cell`, +derived from the angular gaps to its neighbours on either side, which is exact wherever the +neighbourhood is a Voronoi one. Rasters are the exception and have their own method. """ function _contourlengths!(scratch::FD8Scratch, dem, cell, cellsize) n = length(scratch.offsets) @@ -660,6 +685,23 @@ function _contourlengths!(scratch::FD8Scratch, dem, cell, cellsize) return scratch end +"""Rasters take their contour lengths from [`_rastercontour`](@ref).""" +function _contourlengths!( + scratch::FD8Scratch{<:CartesianIndex{2}}, dem, cell, cellsize +) + n = length(scratch.offsets) + resize!(scratch.contours, n) + resize!(scratch.distances, n) + @inbounds for (index, offset) in enumerate(scratch.offsets) + scratch.distances[index] = + celldistance(dem, cell, scratch.neighbors[index]; cellsize) + scratch.contours[index] = _rastercontour(offset, cellsize) + end + resize!(scratch.shape, n) + copyto!(scratch.shape, scratch.offsets) + return scratch +end + """ _fd8weights!(scratch, fd8, dem, cell, visited, cellsize) diff --git a/test/hydrology.jl b/test/hydrology.jl index d535495..e8777fd 100644 --- a/test/hydrology.jl +++ b/test/hydrology.jl @@ -149,3 +149,53 @@ end @test size(weights) == (8, 8, maxreceivers(FD8())) @test all(>=(0), weights) end + +@testset "filldepressions leaves alone what it should" begin + @testset "an integer DEM keeps its type and is not nudged" begin + # Resolving flats moves cells by single float steps, which an integer cannot + # represent, so it is off by default rather than throwing. + dem = fill(5, 5, 5) + @test filldepressions(dem) == dem + @test eltype(filldepressions(dem)) == Int + end + + @testset "cells the caller excluded are untouched" begin + dem = fill(5.0, 5, 5) + @test filldepressions(dem, trues(5, 5)) == dem + # Excluding just the interior leaves those cells alone while the rim resolves. + mask = falses(5, 5) + mask[2:4, 2:4] .= true + z = filldepressions(dem, mask) + @test z[2:4, 2:4] == dem[2:4, 2:4] + end + + @testset "a flat cannot be raised over terrain that was higher" begin + # One float step of relief: the centre is above its neighbours by the smallest + # amount representable, so raising the flat around it must not overtake it. + dem = fill(1.0, 5, 5) + dem[3, 3] = nextfloat(1.0) + z = filldepressions(dem) + centre = CartesianIndex(3, 3) + offsets = filter(!iszero, CartesianIndices((-1:1, -1:1))) + @test any(z[centre + o] < z[centre] for o in offsets) + end +end + +@testset "FD8 weights agree between the stencil and accumulation" begin + # Both read their raster contour lengths from `_rastercontour`, so the share a cell + # gives each neighbour must not depend on which route derived it. + peak = Float64[1 1 1; 1 2 1; 1 1 1] + centre = CartesianIndex(2, 2) + method = FD8(; p = 1.0) + + directions, weights = flowdirection(peak; method, cellsize = unitcell) + acc, = flowaccumulation(copy(peak); method, cellsize = unitcell) + + offsets = Geomorphometry.decompose(CartesianIndex{2}, directions[centre]) + for (slot, offset) in enumerate(offsets) + iszero(offset) && continue + # The centre starts with one cell of area, so a neighbour's accumulation above its + # own area is exactly the share the centre sent it. + @test weights[centre, slot] ≈ acc[centre + offset] - 1 rtol = 1e-5 + end +end