search

jpjones76 / SeisIO.jl

Julia language support for geophysical time series data
http://seisio.readthedocs.org
Other
47 stars 20 forks source link

Causal instrument response correction #46

Open dylanmikesell opened 4 years ago

dylanmikesell commented 4 years ago

I would like to implement the response deconvolution algorithm in this paper from Matt Haney et al. (2012) in BSSA. I was reading through the developer's guide on processing algorithms. This is a published algorithm and I will be translating directly from Matt's MATLAB code. This is an algorithm that ensures causality and has worked great for my purposes of signal correlations across a heterogeneous network of seismic stations. I am curious if this is of use for the SeisIO.jl package and if so, how should I go about adding this (keyword for "type" in remove_resp!). Should I add it to SeisIO.jl/src/Processing/Resp/0_resp_aux.jl?

jpjones76 commented 4 years ago

It's definitely useful for SeisIO and I encourage it.

Recommendations:

If you have a Julia implementation of Matt's algorithm, we can work out insertion points for the new function.

dylanmikesell commented 4 years ago

Okay. I am working on the Julia version now. Will probably be late tomorrow or Tuesday before I have something worth sharing.

Working through his Cleveland Volcano example has brought up another issue with SeisData.resp.f0, but I will post that in another issue....if it turns out I can't figure it out.

dylanmikesell commented 4 years ago

I am back to translating this MATLAB code. Any opinions on interpolation packages?

Interpolations.jl Dierckx.jl GridInterpolations.jl ApproXD.jl

I need to do a 1D nearest neighbor interpolation (to follow Matt's code). I don't see an interpolation dep in SeisIO yet. It will be an interpolation on evenly gridded data. Happy to use whatever people suggest. I am playing with Interpolations.jl right now.

jpjones76 commented 4 years ago

I'm checking each package. Nothing currently in SeisIO requires interpolation, but a robust interpolation package will eventually be needed for many processing operations. I want to be certain that whatever we use is robust and well-maintained.

One thing that keeps me away from R, and which I feel Julia must actively prevent, is a huge ecosystem of packages with tremendous redundancy. By extension I feel like it falls on us to avoid that in SeisIO. (Julia is much better about this than R right now, but that's not setting the bar high enough; R's package ecosystem is a superfund site. Python machine learning might be even worse...)

dylanmikesell commented 4 years ago

Any thoughts on which library? Interpolations.jl seemed pretty good to me after looking through the documentation. Last commit was 15 days ago so it looks like they are still maintaining this package.

tclements commented 4 years ago

Just to chime in - I've used Interpolations.jl and it's had everything I've needed. Dierckx.jl wraps the same Fortran library that Scipy uses for interpolation. Interpolations.jl is probably easier to debug and the author is one of the most active contributors to Julia.

dylanmikesell commented 4 years ago

I have a working version of the interpolation I need using Interpolations.jl. It is actually very slick in my opinion. I am going to keep using this package unless I hear otherwise.

jpjones76 commented 4 years ago

Aha, yes, I think Interpolations.jl is best. Sorry for not responding sooner. I swear I replied to this yesterday, but there's no sign of my comment...??

dylanmikesell commented 4 years ago

No worries. I got the basics finished today. See image that compares new Julia code to Matt's MATLAB output. This is an example that compares colocated broadband and short period sensors in the AV network. Station name is CLES, and this is what Matt gives an example data set when he shares the MATLAB code. You can see that in Julia we have exactly the same results. The BB and SP signals match from MATLAB or Julia, and the small discrepancies between the BHZ and SHZ channels are consistent in MATLAB or Julia. decon_compare_zoom

I need to clean up a few more things and then I will share with @jpjones76 so you can help me figure out where to integrate into the SeisIO.jl base. I will start with a new file in src/Processing/Resp/.

jpjones76 commented 3 years ago

Hi, I haven't heard back from you about this since my email of August 12. Haney's methods are elegant but prohibitively memory-intensive; can you please provide an example to demonstrate the need? I'm looking specifically for something where your code works, but ordinary instrument response translation fails: for example, if you can show me a case where detrend! ==> taper! ==> translate_resp! adds acausal artifacts that bias subsequent analysis, that's a pretty compelling argument.

dylanmikesell commented 3 years ago

Hi Josh. Sorry the semester started and I had field work in the latter half of August. I will work on an example, but it won't be for another 2-3 weeks once the semester is over.

On Wed, Dec 2, 2020 at 3:00 PM Joshua Jones notifications@github.com wrote:

Hi, I haven't heard back from you about this since my email of August 12. Haney's methods are elegant but prohibitively memory-intensive; can you please provide an example to demonstrate the need? I'm looking specifically for something where your code works, but ordinary instrument response translation fails: for example, if you can show me a case where detrend! ==> taper! ==> translate_resp! adds acausal artifacts that bias subsequent analysis, that's a pretty compelling argument.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/jpjones76/SeisIO.jl/issues/46#issuecomment-737521680, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABVIE4TX5UPIEQDRPSYWQPDSS22IFANCNFSM4NPKUENQ .