Modeling/Inversion with ground relief above sea level

[kerim371]

Hi,

Lets suppose we work with land data.
Of course we have to take into account the relief (elevations).
If the land elevation is above sea level (say +100 meters) then probably model's vertical origin should start with -100 right? and vertical spacing is positive (say +25 m).

The problem is that I don't know what sign should have RecGroupElevation and SourceSurfaceElevation.
I think they should be -100 m but anyway JUDI reads them by absolute value: link1, link2.
In this case most likely the JUDI will treat them as they located below sea level.

Also there is more difficult case when elevation on the profile vary from positive to negative values.
Taking absolute elevations will lead to incorrect geometry.

If all my elevations are above sea level may I set vertical origin to positive (+100 m) and vertical spacing to negative (-25 m)?

[kerim371]

Oh, it seems negative spacing doesn't work:

ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'>
ValueError('`nt` must be > 0')
  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/interface.py", line 46, in forward_rec
    rec, _, I, _ = forward(model, src_coords, rec_coords, wavelet, save=False,
  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/propagators.py", line 30, in forward
    src, rcv = src_rec(model, u, src_coords, rcv_coords, wavelet, nt)
  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/pysource/geom_utils.py", line 14, in src_rec
    src = PointSource(name="src%s" % namef, grid=model.grid, ntime=nt,
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/types/basic.py", line 873, in __new__
    newobj._shape = cls.__shape_setup__(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/types/sparse.py", line 336, in __shape_setup__
    raise ValueError('`nt` must be > 0')

Stacktrace:
  [1] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:62 [inlined]
  [2] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:66 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:83
  [4] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/exception.jl:97 [inlined]
  [5] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:11
 [10] pycall(::PyCall.PyObject, ::Type{Tuple{PyCall.PyArray, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:80
 [11] (::JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:20
 [12] (::JUDI.var"#1#2"{JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:68
 [13] lock(f::JUDI.var"#1#2"{JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [14] pylock(f::JUDI.var"#189#190"{Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, DataType})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:65
 [15] wrapcall_data(::PyCall.PyObject, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:fw, :space_order, :f0, :illum), Tuple{Bool, Int64, Float32, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:19
 [16] devito_interface(modelPy::PyCall.PyObject, srcGeometry::GeometryIC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryIC{Float32}, recData::Nothing, dm::Nothing, options::JUDIOptions, illum::Bool, fw::Bool)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/python_interface.jl:72
 [17] time_modeling(model_full::Model, srcGeometry::GeometryOOC{Float32}, srcData::Matrix{Float32}, recGeometry::GeometryOOC{Float32}, recData::Nothing, dm::Nothing, op::Symbol, options::JUDIOptions, fw::Bool)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/time_modeling_serial.jl:39
 [18] propagate(F::judiDataSourceModeling{Float32, :forward}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:9
 [19] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#228#229"{judiDataSourceModeling{Float32, :forward}, judiVector{Float32, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:32
 [20] multi_src_propagate(F::judiDataSourceModeling{Float32, :forward}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:67
 [21] *
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:173 [inlined]
 [22] afoldl
    @ ./operators.jl:533 [inlined]
 [23] *(a::judiProjection{Float32}, b::judiModeling{Float32, :forward}, c::JUDI.jAdjoint{judiProjection{Float32}}, xs::judiVector{Float32, Matrix{Float32}})
    @ Base ./operators.jl:560
 [24] top-level scope
    @ /mnt/HDD_2TB/130406/inversion/scripts/tests_modeling.jl:129

[mloubout]

I'm actually not quite sure what value people usually put in the SegY header for these type of geometry, is there some open data you are using that shows what these values usually are?

Negative spacing will be very problematic yeah, negative origin would be the way to go as long as the segy is read correctly and put negative values for the rec depth.

[kerim371]

@mloubout I'm not sure about open source data but in SEGY header positive elevations usually mean above sea level. Even though seismic software may treat this differently and I always look for the documentation whether above sea level is positive or negative.

With JUDI it seems the simplest solution is to treat RecGroupElevation and SourceSurfaceElevation without making hem absolute.
In this way if we have some registered data above sea level then we can set vertical origin to negative value, vertical spacing is positive and RecGroupElevation SourceSurfaceElevation in SEGY header are negative.

I propose to add an argument to JUDI::Geometry to multiply elevations.
For example make something like:

function Geometry(data::SegyIO.SeisBlock; key="source", segy_depth_key="", dt=nothing, t=nothing, elev_mult=1f0)

and use elev_mult as multiplier to elevations.

It is worth to add information to the documentation that JUDI treats negative elevations above sea level and positive elevations below sea level to make it consistent with model origin/spacings.

[mloubout]

I don't think adding this keyword to Geometry is a good idea, I tend to prefer to have a robust reader rather than having a lot of arguments. So if the abs is the issue I would just remove it

[kerim371]

Yes removing abs would help.
But adding this information to the documentation is necessary as it is not easy to guess.