Added progressive slope filter

Project.toml

@@ -1,7 +1,7 @@
 name = "GeoArrayOps"
 uuid = "e6005643-8f00-58df-85c5-7d93a3520d70"
 authors = ["Maarten Pronk <git@evetion.nl>", "Deltares"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -11,19 +11,21 @@ ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 DataStructures = "0.18"
 Distances = "0.10"
 FillArrays = "0.12"
-ImageCore = "^0.9"
-ImageFiltering = "^0.7"
-OffsetArrays = "^1.10"
+ImageCore = "0.9"
+ImageFiltering = "0.7"
+OffsetArrays = "1.10"
+ProgressMeter = "1.7"
 PaddedViews = "0.5"
 StaticArrays = "1"
-julia = "1.5, 1.6"
+julia = "1.5, 1.6, 1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/GeoArrayOps.jl

@@ -1,15 +1,18 @@
 __precompile__()
 module GeoArrayOps
+using ProgressMeter
 
 include("utils.jl")
 include("pmf.jl")
 include("smf.jl")
+include("psf.jl")
 include("plot.jl")
 include("spread.jl")
 include("terrain.jl")
 
-export pmf, smf
+export pmf, smf, psf
 export pssm
+export pitremoval
 export spread, spread2
 export roughness, TRI, TPI
 
@@ -26,6 +29,8 @@ export roughness, TRI, TPI
 # end
 
 precompile(pmf, (Matrix{Float64},))
+precompile(psf, (Matrix{Float64},))
+precompile(pitremoval, (Matrix{Float64},))
 precompile(smf, (Matrix{Float64},))
 precompile(pssm, (Matrix{Float64},))
 precompile(spread, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}))

src/psf.jl

@@ -0,0 +1,74 @@
+
+"""
+```
+B, flags = psf(A; ωₘ, slope, dhₘ, dh₀, cellsize)
+```
+Applies a progressive slope filter to `A`.
+
+# Output
+- `B::Array{T,2}` Maximum allowable values based A + slope, dhₘ, dh₀
+- `flags::Array{Float64,2}` A sized array with window sizes if filtered, zero if not filtered.
+
+Afterwards, one can retrieve the resulting mask for `A` by `A .<= B` or `flags .== 0.`.
+
+# Arguments
+- `A::Array{T,2}` Input Array
+- `ωₘ::Float64=20.` Maximum window size [m]
+- `slope::Float64=0.01` Terrain slope [m/m]
+- `dhₘ::Float64=2.5` Maximum elevation threshold [m]
+- `dh₀::Float64=0.2` Initial elevation threshold [m]
+- `cellsize::Float64=1.` Cellsize in [m]
+"""
+function psf(A::AbstractMatrix{T};
+    ωₘ::Float64 = 20.0,
+    slope::Float64 = 0.01,
+    dhₘ::Float64 = 2.5,
+    dh₀::Float64 = 0.2,
+    cellsize::Float64 = 1.0,
+    circular = false) where {T<:Real}
+
+    # Compute windowsizes and thresholds
+    ωₘ = round(Int, ωₘ / cellsize)
+    κ_max = floor(Int, log2(ωₘ - 1))  # determine # iterations based on exp growth
+    windowsizes = Int.(exp2.(1:κ_max)) .+ 1
+
+    # Compute tresholds
+    dwindows = vcat(windowsizes[1], windowsizes)  # prepend first element so we get 0 as diff
+    window_diffs = [dwindows[i] - dwindows[i-1] for i = 2:length(dwindows)]
+    height_tresholds = [min(dhₘ, slope * window_diff * cellsize + dh₀) for window_diff in window_diffs]
+    @info "Using the following thresholds: $height_tresholds for the following windows: $windowsizes"
+
+    # Set up arrays
+    Af = copy(A)  # array to be morphed
+    nan_mask = isnan.(Af)
+    Af[nan_mask] .= Inf  # Replace NaN with Inf, as to always filter these
+
+    B = copy(A)  # max_elevation raster
+
+    flags = zeros(size(A))  # 0 = ground, other values indicate window size
+    flags[nan_mask] .= NaN
+
+    mask = falses(size(A))
+
+    # Iterate over window sizes and height tresholds
+    p = Progress(sum(windowsizes .^ 2))
+    progress = 0
+    for (ωₖ, dhₜ) in zip(windowsizes, height_tresholds)
+        if circular
+            mapwindowcirc!(minimum_mask, A, ωₖ, Af, Inf)
+        else
+            mapwindow!(minimum, A, ωₖ, Af)
+        end
+        mask .= (A .- Af) .> dhₜ
+        for I in eachindex(flags)
+            if mask[I] && flags[I] == 0
+                flags[I] = ωₖ
+            end
+        end
+        B .= min.(B, Af .+ dhₜ)
+        progress += ωₖ^2
+        ProgressMeter.update!(p, progress)
+    end
+
+    B, flags
+end

src/terrain.jl

@@ -1,4 +1,5 @@
 const neib_8 = @SMatrix [1.0 1 1; 1 0 1; 1 1 1]
+const neib_8_inf = @SMatrix [1.0 1 1; 1 Inf 1; 1 1 1]
 
 function buffer_array(A::AbstractMatrix{<:Real})
     oA = OffsetMatrix(fill(zero(eltype(A)), size(A) .+ 2), UnitRange.(0, size(A) .+ 1))
@@ -48,6 +49,7 @@ function TPI(dem::AbstractMatrix{<:Real})
     ex_dem = buffer_array(dem)
     tpi = similar(dem)
     x = @MMatrix zeros(3, 3)
+    Δ = CartesianIndex(1, 1)
 
     @inbounds for I ∈ CartesianIndices(size(tpi))
         patch = I-Δ:I+Δ
@@ -71,6 +73,7 @@ function TRI(dem::AbstractMatrix{<:Real})
     ex_dem = buffer_array(dem)
     tri = similar(dem)
     x = @MMatrix zeros(3, 3)
+    Δ = CartesianIndex(1, 1)
 
     @inbounds for I ∈ CartesianIndices(size(tri))
         patch = I-Δ:I+Δ
@@ -81,3 +84,34 @@ function TRI(dem::AbstractMatrix{<:Real})
     end
     return tri
 end
+
+
+"""
+pitremoval(dem::Matrix{<:Real})
+
+Remove pits from a DEM Array if the center cell of a 3x3 patch is `limit` lower or than the surrounding cells.
+"""
+function pitremoval(dem::AbstractMatrix{<:Real}, limit = 2.0)
+
+    ex_dem = buffer_array(dem)
+    tri = similar(dem)
+    x = @MMatrix zeros(3, 3)
+    Δ = CartesianIndex(1, 1)
+
+    @inbounds for I ∈ CartesianIndices(size(tri))
+        patch = I-Δ:I+Δ
+        x .= view(ex_dem, patch)
+        tri[I] = _pitremoval(x, limit)
+    end
+    return tri
+end
+
+@inline function _pitremoval(x, limit)
+    A = vec(x)
+    order = sortperm(A)
+    ifelse(
+        (order[1] == 5) &&
+            ((A[order[2]] - A[order[1]]) >= limit),
+        Inf,
+        x[2, 2])
+end

src/utils.jl

@@ -104,6 +104,7 @@ function mapwindowcirc!(f, img, window, out, fill = Inf)
     out
 end
 
+
 function mapwindowcirc_approx!(f, img, window, out, fill = Inf)
     R = CartesianIndices(img)
     Δ = CartesianIndex(1, 1)
@@ -124,22 +125,18 @@ function mapwindowcirc_approx!(f, img, window, out, fill = Inf)
     out
 end
 
-function maximum_mask(x, m)
+@inline function maximum_mask(x, m)
     o = -Inf
-    for I in eachindex(x)
-        if m[I]
-            o = max(o, x[I])
-        end
+    @inbounds for I in eachindex(x)
+        o = ifelse(m[I], max(o, x[I]), o)
     end
     o
 end
 
-function minimum_mask(x, m)
+@inline function minimum_mask(x, m)
     o = Inf
-    for I in eachindex(x)
-        if m[I]
-            o = min(o, x[I])
-        end
+    @inbounds for I in eachindex(x)
+        o = ifelse(m[I], min(o, x[I]), o)
     end
     o
 end

test/runtests.jl

@@ -19,6 +19,12 @@ end
 @testset "pssm" begin
     B = pssm(rand(25, 25))
 end
+@testset "psf" begin
+    B = psf(rand(25, 25))
+end
+@testset "pitremoval" begin
+    B = pitremoval(rand(25, 25))
+end
 @testset "spread" begin
     points = [0.0 0 0 0 2; 0 0 0 0 0; 0 0 0 0 0; 0 1 0 0 0; 0 0 0 0 0]
     initial = [8.0 8 8 8 4; 8 8 8 8 8; 8 8 8 8 8; 0 0 8 8 8; 0 0 8 8 8]