slimgroup / imagegather.jl Goto Github PK

View Code? Open in Web Editor NEW

11.0 4.0 0.0 8.08 MB

Image gather tools

License: MIT License

Julia 73.82% Python 26.18%

cig seismic-inversion seismic-imaging devito judi julia julia-language julialang

imagegather.jl's People

Contributors

Stargazers

Watchers

imagegather.jl's Issues

Example does not run

ERROR: LoadError: UndefVarError: offset_map not defined
Stacktrace:
 [1] (::Base.var"#770#772")(::Task) at ./asyncmap.jl:178
 [2] foreach(::Base.var"#770#772", ::Array{Any,1}) at ./abstractarray.jl:2009
 [3] maptwice(::Function, ::Channel{Any}, ::Array{Any,1}, ::UnitRange{Int64}) at ./asyncmap.jl:178
 [4] wrap_n_exec_twice at ./asyncmap.jl:154 [inlined]
 [5] async_usemap(::ImageGather.var"#8#9"{Options,Model,judiVector{Float32,Array{Float32,2}},judiVector{Float32,Array{Float32,2}}}, ::UnitRange{Int64}; ntasks::Int64, batch_size::Nothing) at ./asyncmap.jl:103
 [6] #asyncmap#754 at ./asyncmap.jl:81 [inlined]
 [7] asyncmap at ./asyncmap.jl:81 [inlined]
 [8] judipmap(::ImageGather.var"#8#9"{Options,Model,judiVector{Float32,Array{Float32,2}},judiVector{Float32,Array{Float32,2}}}, ::UnitRange{Int64}) at /Users/francisyin/.julia/dev/JUDI/src/TimeModeling/Modeling/utils.jl:7
 [9] surface_gather(::Model, ::judiVector{Float32,Array{Float32,2}}, ::judiVector{Float32,Array{Float32,2}}; offsets::StepRangeLen{Float32,Float64,Float64}, options::Options) at /Users/francisyin/.julia/dev/ImageGather/src/surface_gather.jl:21
 [10] top-level scope at /Users/francisyin/.julia/dev/ImageGather/examples/layers_cig.jl:72
 [11] include(::String) at ./client.jl:457
 [12] top-level scope at REPL[2]:1
in expression starting at /Users/francisyin/.julia/dev/ImageGather/examples/layers_cig.jl:72

Unable to reproduce the example plots

Hello

I got the following plots when running the layers_cig.jl example, where the curves are not flat for true velocity.

Any idea on what's happened? Thanks!

Foundamental questions about the package

Hi,

I'm getting interested in this package and don't have enough understanding how it works.

In seismic after migration we usually want to get migrated common offset or common angle gathers to make some postprocessing or AVO/AVA analisys.
It seems JUDI's RTM can only result in stacked image, it is unable to get prestack migrated gathers.
I guess this package tries to fix that: it allows to compute RTM while getting prestack migrated gathers right?

If so is it possible to get migrated common angle gathers for further AVA analisys?

Error in example script

zyin62@eas-coda-fherr08 examples]$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.1 (2021-12-22)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> include("layers_sscig.jl")
Building born operator
Operator `born` ran in 0.14 s
Building forward operator
Operator `forward` ran in 0.08 s
Building adjoint born operator
Operator `gradient` ran in 0.10 s
Building forward operator
Operator `forward` ran in 0.07 s
ERROR: LoadError: PyError ($(Expr(:escape, :(ccall(#= /data/home/zyin62/.julia/packages/PyCall/7a7w0/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'TypeError'>
TypeError("__new__() got an unexpected keyword argument 'evaluate'")
  File "/data/home/zyin62/.julia/dev/ImageGather/src/implementation.py", line 56, in cig_grad
    op = Operator(pde + go_expr + g_expr,subs=subs, name="cig_sso", opt=opt_op(model))
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/operator/operator.py", line 159, in __new__
    op = cls._build(expressions, **kwargs)
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/operator/operator.py", line 183, in _build
    expressions = cls._lower_exprs(expressions, **kwargs)
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/tools/timing.py", line 76, in __call__
    retval = self.func(*args, **kwargs)
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/operator/operator.py", line 292, in _lower_exprs
    processed = [LoweredEq(i) for i in expressions]
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/operator/operator.py", line 292, in <listcomp>
    processed = [LoweredEq(i) for i in expressions]
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/ir/equations/equation.py", line 148, in __new__
    rhs = diff2sympy(expr.rhs)
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/finite_differences/differentiable.py", line 591, in diff2sympy
    return _diff2sympy(expr)[0]
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/finite_differences/differentiable.py", line 575, in _diff2sympy
    ax, af = _diff2sympy(a)
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/finite_differences/differentiable.py", line 575, in _diff2sympy
    ax, af = _diff2sympy(a)
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/finite_differences/differentiable.py", line 575, in _diff2sympy
    ax, af = _diff2sympy(a)
  [Previous line repeated 1 more time]
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/devito/finite_differences/differentiable.py", line 587, in _diff2sympy
    return obj.func(*args, evaluate=False), True
  File "/data/home/zyin62/.julia/adcme/lib/python3.7/site-packages/sympy/tensor/indexed.py", line 168, in __new__
    obj = Expr.__new__(cls, base, *args, **kw_args)

Stacktrace:
  [1] pyerr_check
    @ ~/.julia/packages/PyCall/7a7w0/src/exception.jl:62 [inlined]
  [2] pyerr_check
    @ ~/.julia/packages/PyCall/7a7w0/src/exception.jl:66 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/.julia/packages/PyCall/7a7w0/src/exception.jl:83
  [4] macro expansion
    @ ~/.julia/packages/PyCall/7a7w0/src/exception.jl:97 [inlined]
  [5] #107
    @ ~/.julia/packages/PyCall/7a7w0/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/.julia/packages/PyCall/7a7w0/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Vector{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/.julia/packages/PyCall/7a7w0/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Vector{Float32}}, kwargs::Base.Pairs{Symbol, Integer, Tuple{Symbol, Symbol}, NamedTuple{(:isic, :space_order), Tuple{Bool, Int64}}})
    @ PyCall ~/.julia/packages/PyCall/7a7w0/src/pyfncall.jl:11
 [10] pycall(::PyCall.PyObject, ::Type{PyCall.PyArray}, ::PyCall.PyObject, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Integer, Tuple{Symbol, Symbol}, NamedTuple{(:isic, :space_order), Tuple{Bool, Int64}}})
    @ PyCall ~/.julia/packages/PyCall/7a7w0/src/pyfncall.jl:80
 [11] propagate(J::judiExtendedJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ ImageGather ~/.julia/dev/ImageGather/src/subsurface_gather.jl:94
 [12] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#202#203"{judiExtendedJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, judiVector{Float32, Matrix{Float32}}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:37
 [13] multi_src_propagate(F::judiExtendedJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:66
 [14] *(F::judiExtendedJacobian{Float32, :adjoint_born, judiDataSourceModeling{Float32, :forward}}, q::judiVector{Float32, Matrix{Float32}})
    @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/LinearOperators/operators.jl:172
 [15] top-level scope
    @ ~/.julia/dev/ImageGather/examples/layers_sscig.jl:65
 [16] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [17] top-level scope
    @ REPL[1]:1
in expression starting at /data/home/zyin62/.julia/dev/ImageGather/examples/layers_sscig.jl:65

[zyin62@eas-coda-fherr08 examples]$ pip show devito
Name: devito
Version: 4.6.2
Summary: Finite Difference DSL for symbolic computation.
Home-page: http://www.devitoproject.org
Author: Imperial College London
Author-email: [email protected]
License: MIT
Location: /data/home/zyin62/.local/lib/python3.6/site-packages
Requires: anytree, cached-property, cgen, click, codecov, codepy, distributed, flake8, multidict, nbval, numpy, pip, psutil, py-cpuinfo, pyrevolve, pytest, pytest-cov, pytest-runner, scipy, sympy
Required-by: 
[zyin62@eas-coda-fherr08 examples]$ pip show sympy
Name: sympy
Version: 1.9
Summary: Computer algebra system (CAS) in Python
Home-page: https://sympy.org
Author: SymPy development team
Author-email: [email protected]
License: BSD
Location: /data/home/zyin62/.local/lib/python3.6/site-packages
Requires: mpmath
Required-by: devito

any thought?

Error when using example `layers_cig.jl`

Hi,

Started using this package.
When running layers_cig.jl example I get error:

ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/twYvK/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'>
ValueError('too many values to unpack (expected 3)')
  File "/home/kerim/Documents/Colada/r/julia-1.6/.julia/dev/ImageGather/src/implementation.py", line 19, in double_rtm
    _, u, _ = forward(model, src_coords, None, wavelet, space_order=space_order,

receiver data in layers_sscig.jl is NaN when space_order=4

Hi @mloubout,

When I set space_order=4 in opt = Options(space_order=4, sum_padding=true), I am getting

This doesn't happen with other values of space_order (e.g. 6, 8, 10 and 12). Any idea of what is going on? This is for layers_sscig.jl.

Thanks.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Trying to run on real data (Viking Graben)

Hi,

I'm trying to compute prestack RTM gathers for Viking Graben line 12 (the data is taken from slimgroup/JUDI.jl#181).
Field data contains 1001 shot with offsets from -262:-25:-3267m.
Shot/Rec step = 25m.
120 rec per shot.

My attempt is:

container = segy_scan(prestk_dir, prestk_file, ["SourceX", "SourceY", "GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(container; segy_depth_key = segy_depth_key_rec)

# JUDI options
jopt = JUDI.Options(
    space_order=32,
    limit_m = true,
    buffer_size = buffer_size,
    optimal_checkpointing=false)

# Left-hand side preconditioners
Ml = judiDataMute(q.geometry, d_obs.geometry, vp=1100f0, t0=0.001f0, mode=:reflection) # keep reflections

# Setup operators
Pr = judiProjection(d_obs.geometry)
F = judiModeling(model0; options=jopt)
Ps = judiProjection(q.geometry)
J = judiJacobian(Pr*F*adjoint(Ps), q)

shot_from = 1
shot_to = length(d_obs)
shot_step = 1

indsrc = rand(shot_from:shot_from+shot_step-1):shot_step:shot_to

# Topmute
d_obs = Ml*d_obs

# PRESTACK RTM
# Common surface offset image gather
offsets = 262f0:25f0:3237f0
CIG = surface_gather(model0, q[indsrc], d_obs[indsrc]; offsets=offsets, options=jopt)

figure()
imshow(CIG[50,:,:], vmin=-1, vmax=1, cmap="PuOr")
gcf()

and the picture I get is something like:

What do I do wrong? Because the result is not something that I expect (the model is computed using FWI).
Does the result should be CDP sorted arrays? If so how to get X coordinate of each CDP?
If my data contains negative offsets -262:-25:-3267m what offsets should I better use?

additional info about size:

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.