Fix DInf and FD8 methods (#53)

* Fix DInf and FD8 methods.

* Remove DEMON draft.

* Duplicate dep.

* Fix cherry pick.

* Fix docs.

Author: evetion

CHANGELOG.md

@@ -60,7 +60,7 @@ All notable changes to this project will be documented in this file.
 - Fixed bug in TPI
 
 ## [0.3.1] - 2022-06-15
-### 
+### Added
 - Added `skewness` filter
 
 ## [0.3.0] - 2022-02-01

Project.toml

@@ -1,7 +1,7 @@
 name = "Geomorphometry"
 uuid = "714e3e49-7933-471a-9334-4a6a65a92f36"
 authors = ["Maarten Pronk <[email protected]>", "Deltares"]
-version = "0.7.0"
+version = "0.7.1"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -20,8 +20,8 @@ Stencils = "264155e8-78a8-466a-aa59-c9b28c34d21a"
 [weakdeps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Eikonal = "a6aab1ba-8f88-4217-b671-4d0788596809"
-GeoArrays = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
+GeoArrays = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
 
 [extensions]
 GeomorphometryEikonalExt = "Eikonal"

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
@@ -12,4 +13,5 @@ Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 
 [compat]
-ColorSchemes = "3.29"
+Makie = "0.22"
+

docs/src/concepts.md

@@ -20,6 +20,7 @@ dtm = coalesce(r, NaN)
 == Rasters (projected)
 ```@example plots
 r = Raster("saba.tif")
+r = coalesce.(r, NaN)  # hide
 Geomorphometry.cellsize(r)
 ```
 ```@example plots
@@ -29,7 +30,7 @@ heatmap(multihillshade(r))
 == GeoArrays (projected)
 ```@example plots
 r = GeoArrays.read("saba.tif")
-r = coalesce(r, NaN)
+r = coalesce(r, NaN)  # hide
 Geomorphometry.cellsize(r)
 ```
 ```@example plots
@@ -39,6 +40,7 @@ heatmap(multihillshade(r))
 == Rasters (geographic)
 ```@example plots
 r = Raster("Copernicus_DSM_10_N52_00_E004_00_DEM.tif")
+r = coalesce.(r, NaN)  # hide
 Geomorphometry.cellsize(r)
 ```
 ```@example plots
@@ -48,7 +50,7 @@ heatmap(multihillshade(r))
 == GeoArrays (geographic)
 ```@example plots
 r = GeoArrays.read("Copernicus_DSM_10_N52_00_E004_00_DEM.tif")
-r = coalesce(r, NaN)
+r = coalesce(r, NaN)  # hide
 Geomorphometry.cellsize(r)
 ```
 ```@example plots

docs/src/guide.md

@@ -1,3 +1 @@
 # Guides
-
-## 

docs/src/refs.bib

@@ -254,3 +254,19 @@ @article{bartelsThresholdfreeObjectGround2010
   langid   = {english},
   keywords = {_tablet,DSM,DTM,GIS,LIDAR,Remote sensing,Skewness Balancing}
 }
+@article{costa-cabralDigitalElevationModel1994,
+  title = {Digital {{Elevation Model Networks}} ({{DEMON}}): {{A}} Model of Flow over Hillslopes for Computation of Contributing and Dispersal Areas},
+  shorttitle = {Digital {{Elevation Model Networks}} ({{DEMON}})},
+  author = {{Costa-Cabral}, Mariza C. and Burges, Stephen J.},
+  year = {1994},
+  journal = {Water Resources Research},
+  volume = {30},
+  number = {6},
+  pages = {1681--1692},
+  issn = {1944-7973},
+  doi = {10.1029/93WR03512},
+  urldate = {2025-06-11},
+  abstract = {Current algorithms for computing contributing areas from a rectangular grid digital elevation model (DEM) use the flow-routing model of O'Callaghan and Mark (1984), which has two major restrictions: (1) flow which originates over a two-dimensional pixel is treated as a point source (nondimensional) and is projected downslope by a line (one dimensional) (Moore and Grayson, 1991), and (2) the flow direction in each pixel is restricted to eight possibilities. We show that large errors in the computed contributing areas result for any terrain topography: divergent, convergent, or planar. We present a new model, called digital elevation model networks (DEMON), which avoids the above problems by representing flow in two dimensions and directed by aspect. DEMON allows computation of both contributing and dispersal areas. DEMON offers the main advantage of contour-based models (e.g., Moore et al., 1988), the representation of varying flow width over nonplanar topography, while having the convenience of using rectangular grid DEMs.},
+  copyright = {Copyright 1994 by the American Geophysical Union.},
+  langid = {english},
+}

src/hydrology.jl

@@ -68,6 +68,17 @@ function filldepressions!(dem::AbstractMatrix, queued = falses(size(dem)))
     return dem
 end
 
+nbs = CartesianIndex.([
+    (-1, -1),
+    (-1, 1),
+    (1, -1),
+    (1, 1),
+    (-1, 0),
+    (0, -1),
+    (0, 1),
+    (1, 0),
+])
+
 function watersheds(dem::AbstractMatrix, queued = falses(size(dem)))
     open = PriorityQueue{CartesianIndex{2}, eltype(dem)}()
     pit = DataStructures.Queue{CartesianIndex{2}}()
@@ -91,7 +102,10 @@ function watersheds(dem::AbstractMatrix, queued = falses(size(dem)))
             labels[cell] = label
             label += 1
         end
-        for ncell in max(I_first, cell - Δ):min(I_last, cell + Δ)
+        # for ncell in max(I_first, cell - Δ):min(I_last, cell + Δ)
+        for nb in nbs
+            ncell = cell + nb
+            ncell in R || continue
             (queued[ncell] || ncell == cell) && continue
             queued[ncell] = true
             if dem[ncell] <= dem[cell]
@@ -116,31 +130,54 @@ function watersheds(dem::AbstractMatrix, queued = falses(size(dem)))
     return dem, labels
 end
 
-const nb =
+# 1  2  3
+# 4  5  6
+# 7  8  9
+
+nbb2 =
     CartesianIndex.([
-        (0, -1),
-        (-1, -1),
-        (-1, 0),
-        (-1, 1),
-        (0, 1),
-        (1, 1),
-        (1, 0),
-        (1, -1),
-        (0, -1),
-        (-1, -1),
+        (-1, 0),  # N  - 270°
+        (-1, 1),  # NE - 315°
+        (0, 1),   # E  - 0° 
+        (1, 1),   # SE - 45°
+        (1, 0),   # S  - 90°
+        (1, -1),  # SW - 135°
+        (0, -1),  # W  - 180°
+        (-1, -1), # NW - 225°
+        (-1, 0),  # N  - 270°
+        (-1, 1),  # NE - 315°
+        (0, 1),   # E  - 0° 
+        (1, 1),   # SE - 45°
     ])
 
 function infc(aspect)
-    aspect = (aspect + 360) % 360
-    i = round(Int, aspect / 45) + 1
-    return nb[i], nb[i + 1]
+    # Normalize aspect to 0-360 range
+    aspect = mod(aspect+90, 360)
+    
+    # Convert aspect to index (0° = North, clockwise)
+    # Each direction covers 45° sector
+    sector = floor(Int, aspect / 45)
+    
+    # Get the two neighboring directions
+    dir1_idx = sector + 1
+    dir2_idx = (sector + 1) % 8 + 1
+    
+    return nbb2[dir1_idx], nbb2[dir2_idx]
 end
 
 function infa(aspect)
-    aspect = (aspect + 360) % 360
-    a1 = aspect % 45
-    a2 = 45 - a1
-    return a2 / 45, a1 / 45
+    # Normalize aspect to 0-360 range
+    aspect = mod(aspect, 360)
+    
+    # Find position within 45° sector
+    sector_angle = mod(aspect, 45)
+    
+    # Calculate weights for the two directions
+    # Weight decreases linearly from 1 to 0 as we move away from the direction
+    weight2 = sector_angle / 45
+    weight1 = 1 - weight2
+    
+    return weight1, weight2
 end
 
 const directions = centered([1 2 3; 4 5 6; 7 8 9])
@@ -157,20 +194,21 @@ function flowaccumulation(
     method = D8(),
     cellsize = cellsize(dem),
 )
-    flowaccumulation!(dem, copy(closed); method, cellsize)
+    acc = similar(dem, Float32)
+    acc .= abs(cellsize[1] * cellsize[2])
+    flowaccumulation!(dem, acc, copy(closed); method, cellsize)
 end
 
 function flowaccumulation!(
     dem::AbstractMatrix,
+    acc::AbstractMatrix{<:Real},
     closed = falses(size(dem));
     method = D8(),
     cellsize = cellsize(dem),
 )
     dir = fill(CartesianIndex{2}(0, 0), size(dem))
     order = ones(Int64, length(closed) - sum(closed))
 
-    acc = similar(dem, Float32)
-    acc .= abs(cellsize[1] * cellsize[2])
     output = similar(dem, eltype(directions))
 
     open = PriorityQueue{CartesianIndex{2}, eltype(dem)}()
@@ -188,7 +226,10 @@ function flowaccumulation!(
         cell = dequeue!(open)
         order[i] = L[cell]
         i += 1
-        for ncell in max(I_first, cell - Δ):min(I_last, cell + Δ)
+        # for ncell in max(I_first, cell - Δ):min(I_last, cell + Δ)
+        for nb in nbs
+            ncell = cell + nb
+            ncell in R || continue
             # skip visited and center cells
             (closed[ncell] || ncell == cell) && continue
 
@@ -210,30 +251,52 @@ function _accumulate!(::D8, acc, order, dir, R, dem, cellsize)
     end
 end
 function _accumulate!(::DInf, acc, order, dir, R, dem, cellsize)
-    asp = aspect(dem; method = Horn())
+    asp = aspect(dem; method = Horn(), cellsize=abs.(cellsize))
+    visited = falses(size(acc))
+
     for i in reverse(order)
         aspect = asp[i]
 
         if !isfinite(aspect)
             acc[R[i] + dir[i]] += acc[i]
+            visited[i] = true
             continue
         end
 
         a, b = infc(aspect)
-        aa, ab = infa(aspect)
+        aa, bb = infa(aspect)
+
+        # a, b = lookup[dir[i] + CartesianIndex(2, 2)]
 
         # Depression
-        if a != dir[i] && b != dir[i]
+        if (a != dir[i] && b != dir[i])
             acc[R[i] + dir[i]] += acc[i]
+            # acc[R[i] + dir[i]] = 5000
+            visited[i] = true
             continue
         end
 
+        # Scale flows correctly at the edges
+        if !(R[i] + a in R) || visited[R[i] + a]
+            aa = 0
+            bb = 1
+        end
+        if !(R[i] + b in R) || visited[R[i] + b]
+            aa = 1
+            bb = 0
+        end
+
         if R[i] + a in R
             acc[R[i] + a] += acc[i] * aa
         end
         if R[i] + b in R
-            acc[R[i] + b] += acc[i] * ab
+            acc[R[i] + b] += acc[i] * bb
+        end
+        if visited[R[i] + a] && visited[R[i] + b]
+            error()
+            acc[R[i] + dir[i]] += acc[i]
         end
+        visited[i] = true
     end
 end
 
@@ -256,6 +319,9 @@ function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize)
         δxy δx δxy
     ]
 
+    visited = falses(size(acc))
+    nb = vec(collect(CartesianIndices(dists)).-CartesianIndex(2,2))
+
     weights = zeros(size(contour_lengths))
     Σw = 0.0
     for i in reverse(order)
@@ -265,18 +331,19 @@ function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize)
             ri == 5 && continue
             I = R[i] + nb[ri]
             I in R || continue
+            visited[I] && continue
             diff = dem[R[i]] - dem[I]
             if diff < 0 || isnan(diff) # neighbor is higher
-                weights[ri] = 0.0
                 continue
             end
             # TODO Check whether this diff/dist is good enough
             weight = (diff / dist * contour_lengths[ri])^fd8.p
             weights[ri] = weight
             Σw += weight
         end
-        if iszero(Σw)  # depression
+        if iszero(Σw) || iszero(weights[CartesianIndex(2,2)+dir[i]])
             acc[R[i] + dir[i]] += acc[i]
+            visited[i] = true
             continue
         end
         for (ri, weight) in enumerate(weights)
@@ -285,6 +352,8 @@ function _accumulate!(fd8::FD8, acc, order, dir, R, dem, cellsize)
             I in R || continue
             acc[I] += acc[i] * (weight / Σw)
         end
+        visited[i] = true
+        fill!(weights, 0)
     end
 end
 

src/relative.jl

@@ -44,7 +44,7 @@ end
     TRI(dem::AbstractMatrix{<:Real})
 
 TRI stands for Terrain Ruggedness Index, which measures the difference between a central pixel and its surrounding cells.
-This algorithm uses the square root of the sum of the square of the difference between a central pixel and its surrounding cells.
+This algorithm uses the square root of the sum of the square of the absolute difference between a central pixel and its surrounding cells.
 This is recommended for terrestrial use cases.
 """
 function TRI(dem::AbstractMatrix{<:Real}; normalize = false, squared = true)