diff --git a/Project.toml b/Project.toml index fd56951..f54e005 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoArrayOps" uuid = "e6005643-8f00-58df-85c5-7d93a3520d70" authors = ["Maarten Pronk ", "Deltares"] -version = "0.2.0" +version = "0.3.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" @@ -11,6 +11,7 @@ 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" @@ -18,12 +19,13 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" 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" diff --git a/src/GeoArrayOps.jl b/src/GeoArrayOps.jl index 64ebb1e..9050794 100644 --- a/src/GeoArrayOps.jl +++ b/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})) diff --git a/src/psf.jl b/src/psf.jl new file mode 100644 index 0000000..9f5fab1 --- /dev/null +++ b/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 diff --git a/src/terrain.jl b/src/terrain.jl index 568f254..06b47ae 100644 --- a/src/terrain.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index 3e37d25..2c80d83 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 93c7ca7..744cec9 100644 --- a/test/runtests.jl +++ b/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]