search

jpjones76 / SeisIO.jl

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

begin time "b" in SAC is not handled correctly #60

Closed taotaokai closed 3 years ago

taotaokai commented 4 years ago

Hi,

On line 101 of SAC.jl instead of adding "b" time as milliseconds

ts += round(Int64, b*1.0f3),

it should be added as seconds

ts += round(Int64, b*1.0f6).

taotaokai commented 4 years ago

Also on line 48 of SAC.jl it should be BUF.sac_fv[6] = rem(ts, 1000)*1.0f-6.

jpjones76 commented 4 years ago

Really? Is the SAC manual wrong? It clearly states "GMT millisecond", not "GMT microsecond"...

http://ds.iris.edu/files/sac-manual/manual/file_format.html

taotaokai commented 4 years ago

Thank you! I do not know if it worths your time, but below is how I found the problem:

  1. Create a sample seismogram "seis.sac" in SAC: funcgen seis write seis.sac quit

  2. In Julia read in this sample sac file and write out a new file as "seis2.sac": using SeisIO s = read_data("sac", "seis.sac") writesac(s[1], "seis2.sac")

  3. compare the SAC headers using saclst (a utility progam from SAC package) $ saclst npts delta b e kztime f seis.sac seis2.sac seis.sac 1000 0.01 9.46 19.45 10:38:14.000 seis2.sac 1000 0.01 0.46 10.46 10:38:14.009

Issues: 1) the seismogram has a duration of 0.01*(1000-1)=9.99 seconds. Since seis.sac is generated by SAC, "b" and "e" should mean seconds. However the "e" time is not set correclty in seis2.sac (should be 10.45); 2) in seis2.sac SeisIO moves 9 seconds from "b", but adds to kztime as 9 milliseconds.

I did the following modifications: 1) SAC.jl line 101: ts is in microsecond, so adding "b" to ts should be ts += round(Int64, b1.0f6) 2) SAC.jl line 48: for the same reason, the sub-millisecond residual time should be put in "b" as BUF.sac_fv[6] = rem(ts, 1000)1.0f-6 3) SAC.jl line 49: the "e" time should be calculated as BUF.sac_fv[7] = BUF.sac_fv[6] + BUF.sac_fv[1]*(nx-1)

Then I ran the above test again and got the correct output: $ saclst npts delta b e kztime f seis.sac seis2.sac seis.sac 1000 0.01 9.46 19.45 10:38:14.000 seis2.sac 1000 0.01 0.000999 9.991 10:38:23.460

jpjones76 commented 3 years ago

Hi again. Thanks for the explanation. I'm currently looking for additional test data to investigate this issue and your proposed fix.

jpjones76 commented 3 years ago

Thank you! I do not know if it worths your time, but below is how I found the problem:

Hi again. I'm going to push a fix live tonight, based on this comment, with my thanks.

I think you are correct that my SAC writer incorrectly handled a non-null begin time. The solution should lead to perfectly re-readable data. Writing the SAC test seismogram no longer produces any time shifts, even subsample ones. I verified that the issue was fixed using SeisIO.jl master branch and SAC v101.6a; it's now part of my automated tests.

One reason this took so long to resolve is that I could find no real test data with b != 0.0f0. I'm not 100% confident that the test file generated by SAC funcgen reflects how "b" and "e" are set in real data. If other users report that these changes break their code, I might need to make adjustments.

Thank you very much!

jpjones76 commented 3 years ago

Fixed in SeisIO v1.2.0.