jpjones76 / SeisIO.jl

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

SAC PZ difference: MATLAB IrisFetch.m vs. get_data() #47

Open dylanmikesell opened 4 years ago

dylanmikesell commented 4 years ago

As I am doing this instrument deconvolution code, I have ran into a small difference in behaviors and I would like to know if this is expected. Here is a minimum working example.

MATLAB: (using the latest IrisFetch.m and IRIS-WS.jar files)

javaaddpath("./IRIS-WS-2.0.19.jar");

tc1 = [2017, 7, 4, 11, 18, 0];
tc2 = [2017, 7, 4, 11, 20, 0];

dum = irisFetch.Traces('AV','CLES','--','BHZ',datestr(tc1,31),datestr(tc2,31),'includePZ');

dum.sacpz.zeros = 
   1.0e+03 *

  0.000000000000000 + 0.000000000000000i
  0.000000000000000 + 0.000000000000000i
  0.000000000000000 + 0.000000000000000i
 -0.392000000000000 + 0.000000000000000i
 -1.960000000000000 + 0.000000000000000i
 -1.490000000000000 + 1.740000000000000i
 -1.490000000000000 - 1.740000000000000i

dum.sacpz.poles =
   1.0e+04 *

 -0.000003691000000 + 0.000003702000000i
 -0.000003691000000 - 0.000003702000000i
 -0.034300000000000 + 0.000000000000000i
 -0.037000000000000 + 0.046700000000000i
 -0.037000000000000 - 0.046700000000000i
 -0.083600000000000 + 0.152200000000000i
 -0.083600000000000 - 0.152200000000000i
 -0.490000000000000 + 0.470000000000000i
 -0.490000000000000 - 0.470000000000000i
 -0.690000000000000 + 0.000000000000000i
 -1.500000000000000 + 0.000000000000000i

dum.sacpz.units =

    'M'

The unit indicate meters.

SeisIO.jl

S  = get_data("FDSN", "AV.CLES..BHZ", src="IRIS", s="2017-07-04T11:18:00", t="2017-07-04T11:20:00")

S[1].resp.z = 
6-element Array{Complex{Float32},1}:
     0.0f0 + 0.0f0im   
     0.0f0 + 0.0f0im   
  -392.0f0 + 0.0f0im   
 -1960.0f0 + 0.0f0im   
 -1490.0f0 + 1740.0f0im
 -1490.0f0 - 1740.0f0im

S[1].resp.p = 
11-element Array{Complex{Float32},1}:
 -0.03691f0 + 0.03702f0im
 -0.03691f0 - 0.03702f0im
   -343.0f0 + 0.0f0im    
   -370.0f0 + 467.0f0im  
   -370.0f0 - 467.0f0im  
   -836.0f0 + 1522.0f0im 
   -836.0f0 - 1522.0f0im 
  -4900.0f0 + 4700.0f0im 
  -4900.0f0 - 4700.0f0im 
  -6900.0f0 + 0.0f0im    
 -15000.0f0 + 0.0f0im   

S[1].units = 
"m/s"

So these two calls result in a different number of zeros, which can be related to the different units.

1) I am curious if this is expected behavior for these two different functions? I thought they both used the same IRIS SACPZ web-request tool. I guess they don't. There is no keyword to change the IrisFetch.m behavior. Is there a keyword to change the get_data() behavior so it returns the zeros for units in meter? 2) Do we want to add some functionality to this instrument deconvolution code for SeisIO.jl? Would it be useful to let the user request which units the data should be returned in (e.g. displacement [m], velocity [m/s], acceleration [m/s/s])? I see there is already the function convert_seis!. Should I just use the default units in the SACPZ for the instrument decon and then let the user use convert_seis!?

jpjones76 commented 4 years ago
  1. Same tool. My intent is for SeisIO to match :resp to the units of the channel's dependent variable -- hence, an extra complex zero in IrisFetch.m. I probably need that to document that explicitly, though.

  2. I'd recommend get_data!(S, ...); convert_seis!(S, ...). I can add a keyword to get_data if needed, but I prefer not to, simply because having many keywords grows cumbersome. Much as I love POSIX-like options (probably obvious, eh?), I suspect <100 people on Earth have them memorized for every Linux function.

Would it break your work flow if I added a utility function to change pole-zero representation to user-specified units? Is that something you'd want?

dylanmikesell commented 4 years ago

You can you do whatever you like to the function :). I would just like to be able to get disp or vel when I want it. Either during the instrument decon step, or after with convert_seis!. I will leave it to you to decide. I don't think we should add another keyword. I think convert_seis! is sufficient for my purposes. It's one extra line/function call in my scripts. I am still trying to get the translation of MATLAB's freqz() and DSP.jl's freqz() sorted out and then this decon code will be done. I plan to get the guts of this into a few function, but let you decide how to integrate it into SeisIO. I have now idea what would be the best way, so not going to submit a pull request if that is okay. I would be happy to chat though about how to integrate it if you like. Seems like you have some ideas though.