New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Convolutional model for seismic data (using a depth model) #190

Merged
merged 65 commits into from Feb 2, 2016

Conversation

Projects
None yet
3 participants
@victortxa
Contributor

victortxa commented Mar 30, 2015

Generate a convolutional model from a depth geological model (Vp; density optional). It's useful to put direct in geology instead of a model in time, avoiding lose the main information. It can be helpful in situations when it's boring to convert from depth to time.

Checklist:

  • Make tests for new code
  • Create/update docstrings
  • Include relevant equations and citations in docstrings
  • Code follows PEP8 style conventions
  • Code and docs have been spellchecked
  • Include new dependencies in docs, requirements.txt, README, and .travis.yml
  • Documentation builds properly
  • All tests pass
  • Can be merged
  • Changelog entry (leave for last)
  • First-time contributor? Add yourself to doc/contributors.rst (leave for last)
@leouieda

This comment has been minimized.

leouieda commented on fatiando/seismic/conv.py in 81d9dbf Mar 30, 2015

Please leave 2 black lines between the imports and the first function.

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Mar 30, 2015

Hi @victortxa. I am really curious about what does it really do. Does it create a synthetic seismic ?
I know that documentation will be provided soon, just wanted to know the general ideia because I might have done something similar in the past.

Cheers and wellcome

@victortxa

This comment has been minimized.

Contributor

victortxa commented Mar 30, 2015

Hi @eusoubrasileiro thank you.
The goal is to take a Vp (density is optional) model in depth and generates the convolutional seismic response in time. In general, the Vp should be given in time (like here: http://www.agilegeoscience.com/blog/2013/12/12/to-make-a-wedge.html ), what makes 'harder' to think geologically. The source is a ricker. @leouieda is introducing me to Git, so I make this pull request to begin, but I'll detail the function asap.
Cheers.

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Mar 30, 2015

@victortxa nice! some time ago I did this I didn't know about Evan Bianco and agile geoscience. I did as a cookbook for testing obspy... Well but It was too much for just testing obpsy io in a cookbook.

Now It will be really nice to have such a tool, really cool!
Again good initiative keep the good job!

Btw I guess that with so many seismic modules we are heading to need to restruture the seismic package to new inner packages? What do you think @leouieda?

@leouieda leouieda changed the title from generate convolutional model from depth model to Convolutional model for seismic data (using a depth model) Apr 8, 2015

@leouieda leouieda added this to the 0.5 milestone Apr 8, 2015

@leouieda

This comment has been minimized.

leouieda commented on fatiando/seismic/conv.py in 30812ba Apr 8, 2015

seismogramgm

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 24, 2015

@victortxa any progress on this?

Could you write an introductory paragraph at the start of the PR description? It should contain a summary of the changes you're making (maybe an example?). Code from #192 will be useful for this (should be merged soon).

@victortxa

This comment has been minimized.

Contributor

victortxa commented Apr 24, 2015

@leouieda Should I write a description of the code or only about the changes I made? About changes do you mean what the code adds to Fatiando? #192 is going to be merged to Fatiando, right? Sorry, I'm still a little confused here.

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 24, 2015

Just a general description of what this PR adds. Nothing too big (2-3 sentences is fine).

#192 is practically done so it won't take too long to be merged. Once it is, you can updated your branch and use the plotting functions proposed there.

@@ -0,0 +1,126 @@
"""

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

@victortxa the first line of the module docstring should describe in a single sentence what this module does. In this case, a good examples would be "Seismic forward modeling by convolution" or "Zero-offset convolutional seismic modeling".

Next, put a small paragraph describing what the functions here do. Very general.

Finally, put a list of the functions present here with a line saying what each does (usually this is the first line of the function docstring).

"""
"""
import numpy as np

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

First thing you should import is from __future__ import division. This will enable the new division handling of Python 3. In Python 2 1/2 == 0. In Python 3 1/2 == 0.5 and 1//2 == 0. So / denotes division and // denotes integer division. This makes things better because you don't have to call float on everything.

This comment has been minimized.

@victortxa

victortxa May 12, 2015

Contributor

Should I use Python 3 already?

for j in range(0, n_traces):
kk = np.ceil(TWT[0, j]/dt_dwn)
lim = np.ceil(TWT[-1, j]/dt_dwn)-1
# linear interpolation

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

This is nitpicking but you should avoid this kind of comment. I can see from the line below that you are doing linear interpolation. So the comment doesn't contribute much. Instead, use the comment to explain why you are doing a linear interpolation.

# extension of the model repeats the last value of the true model
vel[lim:, j] = vel[lim-1, j]
# first values equal to the first value of the true model
vel[0:kk, j] = model[0, j]

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

Same thing for these two comments above as well.

# first values equal to the first value of the true model
vel[0:kk, j] = model[0, j]
# resampling from dt_dwn to dt
vel_l = np.zeros((np.ceil(TMAX/dt), n_traces))

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

This one is great!

# resampling from dt_dwn to dt
vel_l = np.zeros((np.ceil(TMAX/dt), n_traces))
TWT_ts = np.zeros((np.ceil(TMAX/dt), n_traces))
resmpl = int(dt/dt_dwn)

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

As per the comment in the beginning of the file, if you want the integer division of dt and dt_dwn, use dt//dt_dwn. Use int if this needs to be an integer type (to use for indexing, for example).

This comment has been minimized.

@victortxa

victortxa May 12, 2015

Contributor

I use this variable four lines below as array indexing, so I should keep the calculation like this, right?

This comment has been minimized.

@leouieda

leouieda May 13, 2015

Member

OK, I hadn't seen it there. You can keep it like that or use integer division to make it more explicit.

for jj in range(1, len(TWT_ts)):
rho_l[jj, j] = rho[resmpl*jj, j]
# calculate RC
if ~isinstance(rho, np.ndarray):

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

Is this really a ~? Usually ~ negates numpy arrays. I've never seen it used like this before. Could you use not isinstance(... instead?

This comment has been minimized.

@victortxa

victortxa May 12, 2015

Contributor

I think I looked how to negate something in Python and found this. I thought it was valid for any variable, and it fit my expectations (I prefer a symbol than writing not).

rc[1:, :] = (vel_l[1:, :]-vel_l[:-1, :])/(vel_l[1:, :]+vel_l[:-1, :])
else:
rc[1:, :] = ((vel_l[1:, :]*rho_l[1:, :]-vel_l[:-1, :]*rho_l[:-1, :]) /
(vel_l[1:, :]*rho_l[1:, :]+vel_l[:-1, :]*rho_l[:-1, :]))

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

So all of the above is time model building, right? It might be better to separate it in a different function then, like you did with the wavelet code. This will make this function easier to test. Specially since most of the code is model building, hence that's where it's most likely to have bugs. Fortunately, model building is easier to test that the actual results of convolution.

else:
aux = np.floor(len(w)/2.)
synth_l[:, j] = np.convolve(rc[:, j], w, mode='full')[aux:-aux]
return synth_l, TWT_ts

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

Could we separate the convolution part in another function as well? Imagine someone wants to use a time domain model. Then, this function would simply chain the call to depth->time conversion and time->convolution.

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

This would be good to have so people can interface with existing code to model stuff (you mentioned most people model in the time domain).

return synth_l, TWT_ts
def rickerwave(f, dt):

This comment has been minimized.

@leouieda

leouieda May 8, 2015

Member

In Python, there is no concept of "private" functions (things that can only access from within a module/class). But there is a convention to let people know if a certain function/class should not be used in other places. For "private" functions, we append a _ to the beginning of the function name, e.g. _rickerwave. This marks that this function is for internal use only. Such functions don't need to have very detailed docstrings, just a single line is enough.

So consider if this will be a private function for this module. If not, than you should make a complete docstring for it. It is probably better to have it be public so that people can plot the wavelet if they want.

@victortxa

This comment has been minimized.

Contributor

victortxa commented Aug 31, 2015

@leouieda I think/hope it's everything fine, finally...

@leouieda

This comment has been minimized.

Member

leouieda commented Sep 1, 2015

@victortxa I'll have a look at this as soon as I can. Seems that all tests pass and there are no merge conflicts 👍 I just need to take a quick look at the code and see if there is anything missing that I can spot (probably not).

from scipy import interpolate # linear interpolation of velocity/density
def convolutional_model(rc, f, wavelet, dt=2.e-3):

This comment has been minimized.

@leouieda

leouieda Dec 9, 2015

Member

@victortxa now that I'm thinking about this better, I don't think dt should have a default value. This is the sampling interval of rc and the results would be totally wrong if this value is wrong. So it's better to require the user to think about.

return synth_l
def reflectivity(model_t, rho=1.):

This comment has been minimized.

@leouieda

leouieda Dec 9, 2015

Member

As with dt, I don't think rho should have a default value. Defaults are better when the user doesn't really need to think about this, only for advanced use. This is not the case with rho. Forgetting to pass rho would lead to wrong reflectivity.

This comment has been minimized.

@victortxa

victortxa Jan 5, 2016

Contributor

@leouieda I put a default value in rho because in general only velocity is provided by the user, giving rho is more optional. Do you think it's better to put it as a needed input and just warn that put 1.0 to disconsider this parameter?

(model_t[1:, :]*rho[1:, :]+model_t[:-1, :]*rho[:-1, :]))
except TypeError:
rc[1:, :] = ((model_t[1:, :]-model_t[:-1, :]) /
(model_t[1:, :]+model_t[:-1, :]))

This comment has been minimized.

@leouieda

leouieda Dec 9, 2015

Member

Since rho shouldn't have a default, you don't need the try ... except blocks anymore. What you should have is some asserts in the beginning to check if the shapes of model_t and rho match.

You should also check if both have 2 dimensions (you're assuming this in the calculations). Something like assert rho.ndim == 2.

This comment has been minimized.

@victortxa

victortxa Jan 7, 2016

Contributor

@leouieda Rewriting the code I noticed one thing. Maybe the user is interested in calculate the reflectivity just for one trace.
So, how do you think is the best way to keep this option (if you think it a good idea)? I realized to keep the tryblock, but in the except I would write considering model_t and rho as one dimensional. In this way, assert model_t.shape == rho.shape would replace assert rho.ndim == 2.

This comment has been minimized.

@leouieda

leouieda Jan 8, 2016

Member

@victortxa since you're only indexing in the first dimension, you can omit the : for the column index:

rc[1:] = ((model_t[1:]*rho[1:]-model_t[:-1]*rho[:-1]) /
             (model_t[1:]*rho[1:]+model_t[:-1]*rho[:-1]))

This should work the same in 2D or 1D. No need to checking the dimensions then.

def reflectivity(model_t, rho=1.):
"""
Calculate reflectivity series

This comment has been minimized.

@leouieda

leouieda Dec 9, 2015

Member

Could you add a few sentences here to say that this function expects things in the time-domain and the user should use depth_2_time to convert?

Parameters:
* model : 2D-array
Vp values

This comment has been minimized.

@leouieda

leouieda Dec 9, 2015

Member

Please add something like "in the depth domain" to make it clear.

return rc
def depth_2_time(model, dt=2.e-3, dz=1., rho=1.):

This comment has been minimized.

@leouieda

leouieda Dec 9, 2015

Member

Again, dt and dz should not have defaults because these will vary depending on the model.

Also, model doesn't have to be Vp. It could be any physical property given in depth that you want to convert to time. So instead of passing model as Vp and rho as well, the user can call this function twice to get vp and rho in time:

vp_t = depth_2_time(vp, dt, dz)
rho_t = depth_2_time(rho, dt, dz)

This will simplify the function code and make it easier to reuse.

This comment has been minimized.

@victortxa

victortxa Jan 5, 2016

Contributor

I agree, but I only did like this because I kept rho as an optional parameter. If we modify this (above), this adjust here will also be necessary. The negative side I see will be the necessity to call the depth_2_time twice, but if the user doesn't want (or have) any rho then it will not required. Maybe we could keep the default value in the reflectivityfunction and do this modification here in the depth_2_time.
What do you think?

This comment has been minimized.

@leouieda

leouieda Jan 6, 2016

Member

@victortxa I still prefer a more explicit interface. Having rho as default feels strange because it's a crucial parameter to calculate the impedance. It's usually better to be explicit so that users are required to think about that they are doing. In this case, that would mean calling depth_2_time twice or creating a np.ones array for reflectivity. Both are not too complicated.

@leouieda

This comment has been minimized.

Member

leouieda commented Dec 9, 2015

@victortxa sorry I took so long to look at this. I made a few comments about the function parameters. Not too much to do, just a few tweaks to make the code and the API cleaner. You'll need to adapt the cookbook and tests to the required changes.

After this is done, I think it should be ready to merge. You'll need to merge the fatiando/fatiando master branch into this one before we can merge the PR.

@victortxa

This comment has been minimized.

Contributor

victortxa commented Jan 5, 2016

@leouieda Now I took so long to answer, sorry. Maybe it is easier if I show you the modifications I mentioned above, it will be easier to understand. 😬

@victortxa

This comment has been minimized.

Contributor

victortxa commented Jan 8, 2016

@leouieda Now I guess it's all done.

@leouieda

This comment has been minimized.

Member

leouieda commented Feb 2, 2016

@victortxa all seems good to go! Sorry for the long wait. Merging!

leouieda added a commit that referenced this pull request Feb 2, 2016

Merge pull request #190 from victortxa/convolutional_model
Convolutional model for seismic data (using a depth model)

@leouieda leouieda merged commit 9601c0a into fatiando:master Feb 2, 2016

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details

@leouieda leouieda modified the milestones: 0.4, 0.5 Feb 2, 2016

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Feb 2, 2016

Cool!

On Tue, Feb 2, 2016 at 3:02 PM, Leonardo Uieda notifications@github.com
wrote:

Merged #190 #190.


Reply to this email directly or view it on GitHub
#190 (comment).

Um homem que fosse só homem e dissesse as coisas que Jesus disse não
seria um grande mestre de moral. Seria um lunático no mesmo nível de um
homem que diz ser um ovo cozido ou então seria o próprio diabo. Cada um de
nós precisa tomar a sua decisão. Ou este homem era, e é, o Filho de Deus,
ou então um louco, ou algo pior... Mas não venhamos com nenhum argumento
complacente que diga que ele foi apenas um grande mestre humano. Ele não
nos deu esta escolha. Nunca pretendeu fazê-lo. C. S. Lewis - Cristianismo
Puro e Simples

@leouieda leouieda referenced this pull request Feb 2, 2016

Merged

Pad arrays #239

11 of 11 tasks complete
@victortxa

This comment has been minimized.

Contributor

victortxa commented Feb 2, 2016

👍 Fine! 😄

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment