Bugfix 2017-03-17

README.md

@@ -5,6 +5,16 @@ A minimalist, platform-agnostic package for working with univariate geophysical
 http://seisio.readthedocs.org
 
 # CHANGELOG
+## Current
+* Moved pol_sort to the Polarization project.
+* SeisIO data files now include searchable indices at the end of each file.
+  + This change is backwards-compatible and won't affect the ability to read existing files.
+  + A file index contains the following information, written in this order:
+    - (length = ∑\_j∑\_i length(S.id[i])\_i) IDs for each trace in each object
+    - (length = 3∑\_j S.n\_j) start and end times and byte indices of each trace in each object. (time unit = integer μs from Unix epoch)
+    - Byte index to start of IDs.
+    - Byte index to start of Ints.
+
 ## 2017-03-15
 * SeisData merge has been rewritten for greater functionality.
   + Merge speed improved by orders of magnitude for multiple SeisData objects
@@ -28,9 +38,6 @@ http://seisio.readthedocs.org
 * Single-object files can now be written by specifying `sf=true` when calling wseis. By default, single-object file names use IRIS-style naming conventions.
 * Promoted the Polarization submodule to a separate GitHub project.
 * Minor bugfixes: randseischannel, randseisdata, randseisevent, autotap!, IRISget, SeisIO.parserec!, SeisIO.ls, SeisIO,autotuk!
-* *In progress*: SeisIO data files will soon include searchable indices at the end of each file.
-  + This change is backwards-compatible and won't affect the ability to read existing files.
-  + The file index will store trace IDs of each SeisIO object written, start times, end times, and byte indices.
 
 ## 2017-01-31
 * Documentation was completely rewritten.

docs/src/seisdata.rst

@@ -97,18 +97,7 @@ If ``i`` is a string, extract the first channel from ``S`` with ``id=i`` and ret
 Remove all empty channels from ``S``. Empty channels are defined as the set of all channel indices ``i`` s.t. ``isempty(S.x[i]) = true``.
 
 
-Search, Sort
-=============
-.. function:: sort!(S, rev=false)
-
-In-place sort of channels in SeisData struct ``S`` by S.id. Specify ``rev=true`` to reverse the sort order.
-
-.. function:: i = findid(S, C)
-
-Return the index of the first channel in ``S`` with id matching ``C``. If ``C`` is a string, ``findid`` is equivalent to ``findfirst(S.id.==C)``; if ``C`` is a SeisChannel, ``findid`` is equivalent to ``findfirst(S.id.==C.id)``.
-
-
-Save, Load
+Read, Write
 ==========
 .. function:: rsac(fname::String)
 
@@ -130,6 +119,19 @@ Specify ``ts=true`` to write time stamps. Time stamped SAC files created by Seis
 
 Write SeisIO data from S to ``fname``. Supports splat expansion for writing multiple objects, e.g. ``wseis(fname, S, T, U)`` writes ``S``, ``T``, and ``U`` to ``fname``.
 
+To write arrays of SeisIO objects to file, use "splat" notation: for example, for an array ``A`` of type ``Array{SeisEvent,1}``, use syntax ``wseis(fname, A...)``.
+
+
+Search, Sort
+=============
+.. function:: sort!(S, rev=false)
+
+In-place sort of channels in SeisData struct ``S`` by S.id. Specify ``rev=true`` to reverse the sort order.
+
+.. function:: i = findid(S, C)
+
+Return the index of the first channel in ``S`` with id matching ``C``. If ``C`` is a string, ``findid`` is equivalent to ``findfirst(S.id.==C)``; if ``C`` is a SeisChannel, ``findid`` is equivalent to ``findfirst(S.id.==C.id)``.
+
 
 Miscellaneous
 =============

src/Misc/processing.jl

@@ -22,69 +22,6 @@ namestrip!(S::String) = (strip(S, ['\,', '\\', '!', '\@', '\#', '\$',
 namestrip!(S::Array{String,1}) = [namestrip!(i) for i in S]
 namestrip!(S::SeisData) = [namestrip!(i) for i in S.name]
 
-"""
-    T = pol_sort(S::SeisData)
-
-Sort data in `S` by channel ID; return a SeisData structure with seismic data channels arranged {Z, {N,1}, {E,2}, ... }.
-
-    T = pol_sort(S::SeisData, sort_ord = CS::Array{String,1})
-
-As above, with custom sort order `CS`.
-
-    T = pol_sort(S::SeisData, sort_ord = CS::Array{String,1}, inst_codes = IC::Array{Char,1})
-
-As above for channels with instrument codes in `IC`. The instrument code is the second character of the channel field of each S.id string.
-
-Note that `inst_codes` is a Char array, but `sort_ord` is a String array.
-
-In-place channel sorting is not possible as pol_sort discards non-seismic data channels.
-"""
-
-function pol_sort(S::SeisData;
-  sort_ord = ["Z","N","E","0","1","2","3","4","5","6","7","8","9"]::Array{String,1},
-  inst_codes = ['G', 'H', 'L', 'M', 'N', 'P']::Array{Char,1})
-  id = Array{String,1}(S.n)
-  j = Array{Int64,1}(S.n)
-  for i = 1:1:S.n
-    id[i] = S.id[i][1:end-1]
-  end
-
-  ids = sort(unique(id))
-  k = 0
-  for i in ids
-    c = find(ids.==i)
-    d = find(id.==i)
-    L = length(d)
-    if i[end] in inst_codes
-      for m in sort_ord
-        cc = findfirst(S.id .== i*m)
-        if cc > 0
-           k += 1; j[k] = cc; deleteat!(d, findfirst(d.==cc))
-        end
-      end
-
-      # assign other components
-      L = length(d)
-      if L == 1
-        k += 1
-        j[k] = d[1]
-      elseif L > 0
-        j[k+1:k+L] = sortperm(S.id[d])
-        k += L
-      end
-    else
-      j[k+1:k+L] = sortperm(S.id[d])
-      k += L
-    end
-  end
-
-  T = SeisData(length(j))
-  for i = 1:length(j)
-    T[i] = S[j[i]]
-  end
-  return T
-end
-
 """
     autotap!(S::SeisData)
 

src/SeisIO.jl

@@ -26,7 +26,7 @@ gcdist, getbandcode, ls,                                      # Utils/
 fctopz, translate_resp, equalize_resp!,                       # Utils/resp.jl
 d2u, j2md, md2j, parsetimewin, timestamp, u2d,                # Utils/time.jl
 distaz!,                                                      # Misc/event_misc.jl
-autotap!, namestrip!, pol_sort, purge!,                       # Misc/processing.jl
+autotap!, namestrip!, purge!,                                 # Misc/processing.jl
 randseischannel, randseisdata, randseisevent, randseishdr     # Misc/randseis.jl
 
 # Everything depends on time.jl, which calls SeisIO.ls

src/Types/write.jl

@@ -312,8 +312,12 @@ function wseis(S...; sf=false::Bool, path="./"::String, pref=""::String, name=""
         f0 = isempty(pref) ? autoname(getfield(seis,:t)) : pref
         C = UInt8['D']
         ID = join(seis.id,'\0').data
-        TS = vcat([seis.t[j][1,2] for j=1:1:seis.n]...)
-        TE = TS .+ vcat([seis.t[j][2:end,2] for j=1:1:seis.n]...) .+ map(Int64, round(1.0e6.*[length(seis.x[j]) for j=1:1:seis.n]./seis.fs))
+        TS = Array{Int64,1}(seis.n)
+        TE = Array{Int64,1seis.n}()
+        for j = 1:1:seis.n
+          TS[j] = seis.t[j][1,2]
+          TE[j] = TS[j] + sum(seis.t[j][2:end,2]) + round(Int64, 1.0e6*length(seis.x[j])/seis.fs[j])
+        end
 
       elseif typeof(seis) == SeisHdr
         f0 = isempty(pref) ? autoname(getfield(seis),:ot) : pref
@@ -326,9 +330,12 @@ function wseis(S...; sf=false::Bool, path="./"::String, pref=""::String, name=""
         f0 = isempty(pref) ? (d2u(seis.hdr.ot) == 0.0 ? autoname(getfield(getfield(seis,:data),:t)) : autoname(getfield(getfield(seis,:hdr),:ot))) : pref
         C = UInt8['E']
         ID = join(seis.data.id,'\0').data
-        TS = vcat([seis.data.t[j][1,2] for j=1:1:seis.data.n]...)
-        TE = TS .+ vcat([seis.data.t[j][2:end,2] for j=1:1:seis.data.n]...) .+ map(Int64, round(1.0e6.*[length(seis.data.x[j]) for j=1:1:seis.data.n]./seis.data.fs))
-
+        TS = Array{Int64,1}(seis.n)
+        TE = Array{Int64,1seis.n}()
+        for j = 1:1:seis.n
+          TS[j] = seis.data.t[j][1,2]
+          TE[j] = TS[j] + sum(seis.data.t[j][2:end,2]) + round(Int64, 1.0e6*length(seis.data.x[j])/seis.data.fs[j])
+        end
       end
       io = open(path*"/"*join([f0, n0, suff],'.'), "w")
       write(io, "SEISIO".data)