Compatibility with Raster matrix

[pierrethiriet]

I am struggling to replicate the minimal example, including Rasters from the Concept chapter of the documentation.

In the case of the slope function, for example, the Raster doesn't seem recognized with the following output:

MethodError: no method matching slope(::Raster{…}; method::Horn)

But, the Geomorphometry.cellsize(raster) command worked as expected.
(Julia 1.11.5, Rasters v0.14.4, Geomorphometry v0.7.0)

Minimal example :
file downloaded from "https://github.com/Deltares/Geomorphometry.jl/releases/download/v0.6.0/saba.tif" and located in the script folder.

using Rasters
using Geomorphometry
import ArchGDAL

path = @__DIR__
# file from url : "https://github.com/Deltares/Geomorphometry.jl/releases/download/v0.6.0/saba.tif"

r = Raster(joinpath(path, "saba.tif"))
Geomorphometry.cellsize(r)
slope(r, method=Horn())

Thanks a lot
Pierre

[evetion]

Thanks for making an issue! Can you post the whole error and stacktrace?

[pierrethiriet]

Here is the full error stacktrace:

LoadError: MethodError: no method matching slope(::Raster{Union{Missing, Float32}, 2, Tuple{X{Projected{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Intervals{DimensionalData.Dimensions.Lookups.Start}, DimensionalData.Dimensions.Lookups.NoMetadata, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, DimensionalData.Dimensions.Lookups.ReverseOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Intervals{DimensionalData.Dimensions.Lookups.Start}, DimensionalData.Dimensions.Lookups.NoMetadata, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}}, Tuple{Band{DimensionalData.Dimensions.Lookups.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Base.ReshapedArray{Union{Missing, Float32}, 2, Array{Union{Missing, Float32}, 3}, Tuple{}}, Symbol, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GDALsource, Dict{String, Any}}, Missing}; method::Horn)
The function `slope` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  slope(::AbstractMatrix{<:Real}; cellsize, method, exaggeration, direction)
   @ Geomorphometry C:\Users\pthiriet\.julia\packages\Geomorphometry\qqdTK\src\terrain.jl:29

Stacktrace: 
  [1] top-level scope
    @ xx\Typology\test.jl:10
  [2] eval
    @ .\boot.jl:430 [inlined]
  [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:2734
  [4] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Base .\essentials.jl:1055
  [5] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base .\essentials.jl:1052
  [6] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\eval.jl:271
  [7] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\eval.jl:181
  [8] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)  
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\repl.jl:276
  [9] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [10] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [11] #67
    @ c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\eval.jl:150 [inlined]
 [12] with_logstate(f::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, logstate::Base.CoreLogging.LogState)
    @ Base.CoreLogging .\logging\logging.jl:522
 [13] with_logger
    @ .\logging\logging.jl:632 [inlined]
 [14] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [15] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [16] invokelatest(::Any)
    @ Base .\essentials.jl:1052
 [17] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\pthiriet\.vscode\extensions\julialang.language-julia-1.140.2\scripts\packages\VSCodeServer\src\eval.jl:34
in expression starting at xx\Typology\test.jl:10

[rafaqz]

Looks like you having missing values and the alg can't handle that?

You can maybe replace them with zeros if they are e.g. the sea.

Raster(filename; missingval=0) should do that.

Otherwise use missingval=NaN but be aware the NaNs may smear into the neighbourhood cells. @evetion will know better what the alg needs.

[pierrethiriet]

Thanks !
Indeed, missingval=0 solved the issue.
But, theoretically speaking, if the target pixel has a value and at least one of the surrounding ones also, the slope should be calculated, even if others are missing, no?

[evetion]

Yes, all operations in Geomorphometry currently need a Matrix{<:Real}, so no missing values. That said, it's on the backlog to add it.

But, theoretically speaking, if the target pixel has a value and at least one of the surrounding ones also, the slope should be calculated, even if others are missing, no?

I fear not, for most calculations all neighbor pixels are taken into account, and any non-finite value will propagate. Thus, it's really best to replace/interpolate any missing/nan values.

[rafaqz]

It would be good to add an error at least explaining this for Union{<:Real,Mising}, mostly because its the default in Rasters.jl now and I imagine this might happen a bit. We could suggest specifying missingval or interpolating.