jpjones76 / SeisIO.jl

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

Can't read multiple SAC files using `read_data` #88

Open adigitoleo opened 2 years ago

adigitoleo commented 2 years ago

I'm unable to read multiple SAC channels into a single SeisData instance using read_data. Tried on v1.2.1 and also latest master. The files are the same as uploaded for #87 .

julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4* (2021-07-14 15:36 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8705G CPU @ 3.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

shell> ls
2019.262.07.06.33.3000.CN.FRB..BHE.M.SAC  2019.262.07.06.33.3000.CN.FRB..BHN.M.SAC  2019.262.07.06.33.3000.CN.FRB..BHZ.M.SAC

julia> s = SeisData()
SeisData with 0 channels (0 shown)
    ID:
  NAME:
   LOC:
    FS:
  GAIN:
  RESP:
 UNITS:
   SRC:
  MISC:
 NOTES:
     T:
     X:

     C: 0 open, 0 total

julia> read_data!(s, "sac", "2019.262.07.06.33.3000.CN.FRB..BH*.SAC")

julia> s
SeisData with 1 channels (1 shown)
    ID: C.FR..BH
  NAME:
   LOC: 63.7469 N, -68.5451 E, 25.0 m
    FS: 40.0
  GAIN: 1.0
  RESP: a0 1.0, f0 1.0, 0z, 0p
 UNITS:
   SRC: /home/admin/hazcode/julia/julia-s…
  MISC: 0 entries
 NOTES: 1 entries
     T: 2019-09-19T07:06:33 (2 gaps)
     X: -6.722e+04
        -6.701e+04
            ...
        -2.850e+04
        (nx = 432000)
     C: 0 open, 0 total
tclements commented 2 years ago

This is being caused by #87 but read_data is working as expected. When I read all three channels at once

s = read_data("sac", "./*.SAC")
SeisData with 1 channels (1 shown)
    ID: C.FR..BH                           
  NAME:                                    
   LOC: 63.7469 N, -68.5451 E, 25.0 m      
    FS: 40.0                               
  GAIN: 1.0                                
  RESP: a0 1.0, f0 1.0, 0z, 0p             
 UNITS:                                    
   SRC: /home/timclements/Downloads/SACTE… 
  MISC: 0 entries                          
 NOTES: 1 entries                          
     T: 2019-09-19T07:06:33 (2 gaps)       
     X: -6.722e+04                         
        -6.701e+04                         
            ...                            
        -2.850e+04                         
        (nx = 432000)                      
     C: 0 open, 0 total

there are two gaps. We can check out the data gaps with timing information S.t:

 s.t
1-element Array{Array{Int64,2},1}:
 [1 1568876793300000; 144001 -3600000000; 288001 -3600000000; 432000 0]

each data segment has 144,000 samples or 1 hour of data. We can check that with some helper functions:

SeisIO.u2d.(SeisIO.t_win(s.t[1],s.fs[1]) .* SeisIO.μs)
3×2 Array{Dates.DateTime,2}:
 2019-09-19T07:06:33.3  2019-09-19T08:06:33.275
 2019-09-19T07:06:33.3  2019-09-19T08:06:33.275
 2019-09-19T07:06:33.3  2019-09-19T08:06:33.275

Since all three channels have the same name and start/end times, calling ungap should leave us with a single channel:

ungap(s)
SeisData with 1 channels (1 shown)
    ID: C.FR..BH                           
  NAME:                                    
   LOC: 63.7469 N, -68.5451 E, 25.0 m      
    FS: 40.0                               
  GAIN: 1.0                                
  RESP: a0 1.0, f0 1.0, 0z, 0p             
 UNITS:                                    
   SRC: /home/timclements/Downloads/SACTE… 
  MISC: 0 entries                          
 NOTES: 2 entries                          
     T: 2019-09-19T07:06:33 (0 gaps)       
     X: -3.430e+04                         
        -3.494e+04                         
            ...                            
        -2.850e+04                         
        (nx = 144000)                      
     C: 0 open, 0 total

which it does. Solving #87 will solve this issue.

jpjones76 commented 2 years ago

Duplicate of issue #87; leaving open while I investigate. If my guess is right, the fix is easy, but I need to do some further checking before I push code live (and I'm still looking for a suitable Travis-CI replacement, whose greed is ultimately the cause of these delays)...