event utils update

jpjones76 · Nov 18, 2016 · e0ab131 · e0ab131
1 parent 825dbec
commit e0ab131
Show file tree

Hide file tree

Showing 2 changed files with 78 additions and 81 deletions.
diff --git a/src/SeisIO.jl b/src/SeisIO.jl
@@ -17,7 +17,7 @@ FDSNget,                                                       # Web/FDSN.jl
 IRISget, irisws,                                               # Web/IRIS.jl
 SeedLink, SeedLink!, has_sta, SLinfo,                          # Web/SeedLink.jl
 GetSta, chparse,                                               # Web/WebMisc.jl
-evq, gcdist, distaz!, getpha, getevt,                          # Utils/event_utils.jl
+evq, gcdist, distaz!, getpha, getevt, getsta, getNextPhase,    # Utils/event_utils.jl
 pull, getbandcode, prune!, purge, purge!, chan_sort, note,     # Utils/misc.jl
 gapfill, ungap!, ungap, sync!, sync, gapfill!,                 # Utils/sync.jl
 randseischa, randseisdata, randseisevent, randseishdr,         # Utils/randseis.jl

diff --git a/src/Utils/event_utils.jl b/src/Utils/event_utils.jl
@@ -201,11 +201,8 @@ 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
+  cid = evt[9]
+  CONTRIB_ID = isnumber(cid) ? parse(cid) : -1
   return SeisHdr( id = parse(Int64, evt[1]),
                   time = Dates.DateTime(evt[2]),
                   lat = parse(Float64, evt[3]),
@@ -242,6 +239,68 @@ function distaz!(S::SeisEvent)
   end
 end
 
+"""
+    T = getpha(pha::String, Δ::Float64, z::Float64)
+
+Get onset times of phases **pha** relative to origin time for an event at
+distance **Δ** (degrees), depth **z** (km). Returns a matrix of strings.
+
+Detail: getpha is a command-line interface to the IRIS travel time calculator,
+which calls TauP (1,2,3). Specify **pha** as a comma-separated string, e.g.
+"P,S,PKiKP". **pha** also accepts special keywords (e.g. \"ttall\") as described
+on the IRIS web pages.
+
+References:
+(1) IRIS travel time calculator: https://service.iris.edu/irisws/traveltime/1/
+(2) TauP manual: http://www.seis.sc.edu/downloads/TauP/taup.pdf
+(3) Crotwell, H. P., Owens, T. J., & Ritsema, J. (1999). The TauP Toolkit:
+Flexible seismic travel-time and ray-path utilities, SRL 70(2), 154-160.
+"""
+function getpha(Δ::Float64, z::Float64;
+  phases=""::String,
+  model="iasp91"::String,
+  to=10.0::Real, v=false::Bool, vv=false::Bool)
+
+  # Generate URL and do web query
+  if isempty(phases)
+    pq = ""
+  else
+    pq = string("&phases=", phases)
+  end
+
+  url = string("http://service.iris.edu/irisws/traveltime/1/query?",
+  "distdeg=", Δ, "&evdepth=", z, pq, "&model=", model,
+  "&mintimeonly=true&noheader=true")
+  (v | vv) && println("url = ", url)
+  req = readall(get(url, timeout=to, headers=webhdr()))
+  (v | vv) && println("Request result:\n", req)
+
+  # Parse results
+  phase_data = split(req, '\n')
+  sa_prune!(phase_data)
+  Nf = length(split(phase_data[1]))
+  Np = length(phase_data)
+  Pha = Array{String,2}(Np, Nf)
+  for p = 1:Np
+    Pha[p,1:Nf] = split(phase_data[p])
+  end
+  return Pha
+end
+
+getphaseTime(pha::String, Pha::Array{String,2}) = parse(Float64, Pha[find(Pha[:,3].==pha)[1],4])
+getPhaSt(pha::String, Pha::Array{String,2}) = findmin([parse(Float64,i) for i in Pha[:,4]])
+getPhaEn(pha::String, Pha::Array{String,2}) = findmax([parse(Float64,i) for i in Pha[:,4]])
+function getNextPhase(pha::String, Pha::Array{String,2})
+  s = Pha[:,3]
+  t = [parse(Float64,i) for i in Pha[:,4]]
+  j = find(s.==pha)[1]
+  i = t.-t[j].>0
+  tt = t[i]
+  ss = s[i]
+  k = sortperm(tt.-t[j])[1]
+  return(ss[k],tt[k])
+end
+
 """
     getevt(evt::String, cc::String)
 
@@ -303,19 +362,19 @@ function getevt(evt::String, cc::String;
       s = parse(Float64,pdat[j,4]) - spad
       t = 2.0*parse(Float64,pdat[getPhaEn(pdat),4])
       S.data.misc[i]["PhaseWindow"] = string(pdat[j,3], " : Coda")
-    else
-      # 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
+    else
+      # 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
@@ -350,65 +409,3 @@ function getevt(evt::String, cc::String;
   end
   return S
 end
-
-"""
-    T = getpha(pha::String, Δ::Float64, z::Float64)
-
-Get onset times of phases **pha** relative to origin time for an event at
-distance **Δ** (degrees), depth **z** (km). Returns a matrix of strings.
-
-Detail: getpha is a command-line interface to the IRIS travel time calculator,
-which calls TauP (1,2,3). Specify **pha** as a comma-separated string, e.g.
-"P,S,PKiKP". **pha** also accepts special keywords (e.g. \"ttall\") as described
-on the IRIS web pages.
-
-References:
-(1) IRIS travel time calculator: https://service.iris.edu/irisws/traveltime/1/
-(2) TauP manual: http://www.seis.sc.edu/downloads/TauP/taup.pdf
-(3) Crotwell, H. P., Owens, T. J., & Ritsema, J. (1999). The TauP Toolkit:
-Flexible seismic travel-time and ray-path utilities, SRL 70(2), 154-160.
-"""
-function getpha(Δ::Float64, z::Float64;
-  phases=""::String,
-  model="iasp91"::String,
-  to=10.0::Real, v=false::Bool, vv=false::Bool)
-
-  # Generate URL and do web query
-  if isempty(phases)
-    pq = ""
-  else
-    pq = string("&phases=", phases)
-  end
-
-  url = string("http://service.iris.edu/irisws/traveltime/1/query?",
-  "distdeg=", Δ, "&evdepth=", z, pq, "&model=", model,
-  "&mintimeonly=true&noheader=true")
-  (v | vv) && println("url = ", url)
-  req = readall(get(url, timeout=to, headers=webhdr()))
-  (v | vv) && println("Request result:\n", req)
-
-  # Parse results
-  phase_data = split(req, '\n')
-  sa_prune!(phase_data)
-  Nf = length(split(phase_data[1]))
-  Np = length(phase_data)
-  Pha = Array{String,2}(Np, Nf)
-  for p = 1:Np
-    Pha[p,1:Nf] = split(phase_data[p])
-  end
-  return Pha
-end
-
-getphaseTime(pha::String, Pha::Array{String,2}) = parse(Float64, Pha[find(Pha[:,3].==pha)[1],4])
-getPhaSt(pha::String, Pha::Array{String,2}) = findmin([parse(Float64,i) for i in Pha[:,4]])
-getPhaEn(pha::String, Pha::Array{String,2}) = findmax([parse(Float64,i) for i in Pha[:,4]])
-function getNextPhase(pha::String, Pha::Array{String,2})
-  s = Pha[:,3]
-  t = [parse(Float64,i) for i in Pha[:,4]]
-  j = find(s.==pha)[1]
-  i = t.-t[j].>0
-  tt = t[i]
-  ss = s[i]
-  k = sortperm(tt.-t[j])[1]
-  return(ss[k],tt[k])
-end