fixed sync bug

src/Utils/event_utils.jl

@@ -201,6 +201,11 @@ gc_unwrap!(t::Array{Float64,1}) = (t[t .< 0] .+= 2.0*π; return t)
 
 function mkevthdr(evt_line::String)
   evt = split(evt_line,'|')
+  CONTRIB_ID = try
+      parse(Int64, evt[9])
+    catch
+      -1
+    end
   return SeisHdr( id = parse(Int64, evt[1]),
                   time = Dates.DateTime(evt[2]),
                   lat = parse(Float64, evt[3]),
@@ -209,7 +214,7 @@ function mkevthdr(evt_line::String)
                   auth = evt[6],
                   cat = evt[7],
                   contrib = evt[8],
-                  contrib_id = parse(Int64, evt[9]),
+                  contrib_id = CONTRIB_ID,
                   mag_typ = evt[10],
                   mag = parse(Float32, evt[11]),
                   mag_auth = evt[12],
@@ -290,6 +295,7 @@ function getevt(evt::String, cc::String;
   # If the phase string supplied is "all", request window is spad s before P to twice the last phase arrival
   # If a phase name is supplied, request window is spad s before that phase to epad s after next phase
   pstr = Array{String,1}(S.data.n)
+  bads = falses(S.data.n)
   for i = 1:1:S.data.n
     pdat = getpha(S.data.misc[i]["dist"], S.hdr.dep, to=to, v=v, vv=vv)
     if pha == "all"
@@ -298,10 +304,22 @@ function getevt(evt::String, cc::String;
       t = 2.0*parse(Float64,pdat[getPhaEn(pdat),4])
       S.data.misc[i]["PhaseWindow"] = string(pdat[j,3], " : Coda")
     else
-      s = parse(Float64,pdat[find(pdat[:,3].==pha)[1],4]) - spad
-      (p2,t) = getNextPhase(pha, pdat)
+      # Note: at Δ > ~90, we must use Pdiff; we can't use P
+      p1 = pha
+      j = findfirst(pdat[:,3].==p1)
+      if j == 0
+        p1 = pha*"diff"
+        j = findfirst(pdat[:,3].==p1)
+        if j == 0
+          error(string("Neither ", pha, " nor ", pha, "diff found!"))
+        else
+          warn(string(pha, "diff substituted for ", pha, " at ", S.data.id[i]))
+        end
+      end
+      s = parse(Float64,pdat[j,4]) - spad
+      (p2,t) = getNextPhase(p1, pdat)
       t += epad
-      S.data.misc[i]["PhaseWindow"] = string(pha, " : ", p2)
+      S.data.misc[i]["PhaseWindow"] = string(p1, " : ", p2)
     end
     s = string(u2d(d2u(S.hdr.time) + s))
     t = string(u2d(d2u(S.hdr.time) + t))
@@ -312,14 +330,24 @@ function getevt(evt::String, cc::String;
     C = FDSNget(net = NET, sta = STA, loc = LOC, cha = CHA,
                 s = s, t = t, si = false, y = false, v=v, vv=vv)
     vv && println("FDSNget output:\n", C)
-    S.data.t[i] = C.t[1]
-    S.data.x[i] = C.x[1]
-    S.data.notes[i] = C.notes[1]
-    S.data.src[i] = C.src[1]
-    if (v | vv)
-         println(now(), ": data acquired for ", S.data.id[i])
+    if C.n == 0
+      bads[i] = true
+    else
+      S.data.t[i] = C.t[1]
+      S.data.x[i] = C.x[1]
+      S.data.notes[i] = C.notes[1]
+      S.data.src[i] = C.src[1]
+      if (v | vv)
+        println(now(), ": data acquired for ", S.data.id[i])
+      end
     end
   end
+  bad = find(bads.==true)
+  if !isempty(bad)
+    ids = join(S.data.id[bad],',')
+    warn(string("Channels ", ids, " removed (no data were found)."))
+    deleteat!(S.data, bad)
+  end
   return S
 end
 

src/Utils/sync.jl

@@ -95,57 +95,72 @@ Start time can be synchronized to a variety of values:
 """
 function sync!(S::SeisData; resample=false::Bool, fs=0::Real,
   s="max"::Union{String,Real,DateTime},
-  t="max"::Union{String,Real,DateTime})
+  t="max"::Union{String,Real,DateTime},
+  v=false::Bool)
 
-  # Do not edit order of operations
+  for i = 1:S.n
+    S.x[i] -= mean(S.x[i])
+  end
   ungap!(S, m=false, w=false)
   #autotap!(S)
-  start_times = zeros(S.n)
-  end_times = zeros(S.n)
+
+  # Do not edit order of operations -------------------------------------------
+  start_times = zeros(Int64, S.n)
+  end_times = zeros(Int64, S.n)
+
+  # non-timeseries data
   c = find(S.fs .== 0)
   for i in c
     start_times[i] = S.t[i][1,2]
     end_times[i] = sum(S.t[i][:,2])
   end
+
+  # non-null timeseries data
   k = find((S.fs .> 0).*[!isempty(S.x[i]) for i=1:S.n])
   for i in k
-    S.x[i] -= mean(S.x[i])
     start_times[i] = round(Int, S.t[i][1,2]*S.fs[i]) / S.fs[i]
+    end_times[i] = start_times[i] + round(Int, length(S.x[i])/(μs*S.fs[i]))
   end
-  [end_times[i] = length(S.x[i])/S.fs[i] for i in k]
-  # Do not edit order of operations
+  # Do not edit order of operations -------------------------------------------
+  v && println("start_times = ", start_times)
+  v && println("end_times = ", end_times)
 
   # Possible options for s
+  v && println("s = ", s)
   if isa(s,Real)
     t_start = round(Int, s/μs)
   elseif isa(s,DateTime)
     t_start = round(Int, Dates.datetime2unix(s)/μs)
   else
     starts = start_times[k]
+    v && println("starts = ", starts)
     if s == "max"
-      t_start = minimum(starts)
+      t_start = floor(Int, minimum(starts))
     elseif s == "min"
-      t_start = maximum(starts)
+      t_start = ceil(Int, maximum(starts))
     else
       t_start = round(Int, Dates.datetime2unix(DateTime(s))/μs)
     end
   end
+  v && println("t_start =", t_start)
 
   # Possible options for t
+  v && println("t = ", t)
   if isa(t,Real)
     t_end = t_start + round(Int, t/μs)
   elseif isa(t,DateTime)
     t_end = round(Int, Dates.datetime2unix(t)/μs)
   else
-    ends = end_times[k] + start_times[k]
     if t == "max"
-      t_end = maximum(ends)
+      t_end = maximum(end_times)
     elseif t == "min"
-      t_end = minimum(ends)
+      t_end = minimum(end_times)
     else
       t_end = round(Int, Dates.datetime2unix(DateTime(t))/μs)
     end
   end
+  v && println("t_end =", t_end)
+
   t_end <= t_start && error("No time overlap with given start \& end times!")
   @printf(STDOUT, "Synching %.2f seconds of data\n", (t_end - t_start)*μs)
 
@@ -161,13 +176,11 @@ function sync!(S::SeisData; resample=false::Bool, fs=0::Real,
       note(S, i, @sprintf("Resampled to %.1f", f0))
     end
   end
-
+  v && println("t_start = ", t_start)
+  v && println("t_end = ", t_end)
   sstr = string(Dates.unix2datetime(t_start*μs))
   tstr = string(Dates.unix2datetime(t_end*μs))
 
-  # Pad and truncate
-  end_times += start_times
-
   # Loop over non-timeseries data
   for i in c
     if S.fs[i] == 0