fixed autotap! for non-ts data

docs/src/seisdata.rst

@@ -19,13 +19,17 @@ SeisData and SeisObj containers can be created in three ways:
 
 #. ``S = SeisData()`` initializes an empty SeisData container
 
-#. ``S = SeisObj()``  initializes a new SeisObj container
+#. ``S = SeisObj()`` initializes a new SeisObj container
 
-#. ``S = SeisData(s1, s2, s3)`` creates a SeisData container by merging s1, s2, s3. If a merge variable isn't of type SeisData or SeisObj, it's ignored and a warning is given to STDOUT.
+#. ``S = SeisData(s1, s2, s3)`` creates a SeisData container by merging s1, s2, s3. If a variable passed to SeisData in this way isn't of type SeisData or SeisObj, it's ignored with a warning.
 
 Fields can be initialized by name when a new SeisObj container is created; for example,``S = SeisObj(name="DEAD CHANNEL", fs=50)`` initialize a SeisObj with name "DEAD CHANNEL", fs = 50.
 
 
+Example
+--------
+```S = SeisData(SeisObj(name="BRASIL", id="IU.SAML.00.BHZ"), SeisObj(name="UKRAINE", id="IU.KIEV.00.BHE"), SeisObj())`` creates* a new SeisData structure with three channels; the first is named "BRASIL", the second "UKRAINE", the third is blank. This syntax requires unique IDs for each new channel created.
+
 
 Adding Data
 ===========
@@ -40,36 +44,36 @@ Data can be merged directly from the output of any SeisIO command that outputs a
 In a merge operation, pairs of non-NaN data `x`:sub:`i`, `x`:sub:`j` with overlapping time stamps (i.e. `t`:sub:`i` - `t`:sub:`j` < 1/fs) are *averaged*.
 
 
+
 Appending without merging
 -------------------------
-``push!`` adds channels without attempting to merge.
-
-* ``push!(S,T)``
+``push!`` adds a SeisObj instance to a SeisData instance without attempting to merge the channel with existing data. ``append!`` adds a SeisData instance to another SeisData instance without attempting to merge identical channels.
 
-  assign each channel in T to a new channel in S, even if it creates redundant channel IDs.
+``push!(S,T)`` assigns each channel in T to a new channel in S, even if it creates redundant channel IDs.
 
-* ``push!(S,S)``
-
-  append a duplicate of each channel in S to S.
+``push!(S,S)`` appends a duplicate of each channel in S to S.
 
 
 
 Deleting Data
 =============
-``-`` (the subtraction operator) is the standard way to delete unwanted data channels. It's generally safest to use e.g. ``T = S - i``, but in-place deletion (e.g. ``S -= i``) is valid.
-
-* ``S - k``
+``-`` (the subtraction operator) is the standard way to delete unwanted data channels. It's generally safest to use e.g. ``T = S - i``, but in-place deletion (e.g. ``S -= i``) is valid. The exact behavior for a general operation ``S - K`` depends on the type of the addend:
 
-  * If k is an integer, channel k is deleted from S.
+  * If K is an integer, channel k is deleted from S.
 
-  * If k is a string, all channels whose names and ids match k are deleted.
+  * If K is a string, all channels whose names and ids match k are deleted.
 
-  * If k is a SeisObj instance, all channels from S with the same id as k are deleted.
+  * If K is a SeisObj instance, all channels from S with the same id as k are deleted.
 
-* ``deleteat!(S,i)``
+``deleteat!(S,i)`` is identical to ``S-=K`` for an integer K.
 
-  identical to ``S-=i`` for an integer i.
 
+Safe deletion with ``pull``
+---------------------------
+The ``pull`` command extracts a channel from a SeisData instance and returns it as a SeisObj.
+```T = pull(S, name)``, where ``name`` is a string, creates a SeisObj ``T`` from the first channel with name=``name``. The channel is removed from `S`.
++ `T = pull(S, i)`, where `i` is an integer, creates a SeisObj `T` from
+channel `i` and removes channel `i` from `S`.
 
 
 Index, Search, Sort
@@ -106,21 +110,31 @@ The following commands offer tests for equality:
 * ``==`` and ``isequal`` test for equality.
 
 
+Annotation
+==========
+SeisData and SeisObj are intended to be annotated as data are analyzed; the command ``note`` exists for this purpose. Calling ``note`` appends a timestamped note to the ``notes`` field of the target container.
+
+* ``note(S, i, NOTE)`` adds string ``NOTE`` to channel ``i`` of ```S``.
+
+* ``S += "NAME: NOTE"``, adds a note with a short form date. NAME should be the intended channel name; NOTE should be the intended note content.
+
+* The command ``S += NOTE`` adds the text of string ``NOTE`` to all channels.
+
 
 Utility Functions
 =================
-* ``note``: Add a timestamped note.
+* ``plotseis``: Plot time-aligned data from a SeisData object. Time series data are represented by straight lines; irregularly sampled data (``fs=0``) use normalized stem plots.
 
-* ``plotseis``: Plot time-aligned data. Time series data are represented by straight lines; irregularly sampled data (``fs=0``) use normalized stem plots.
+* ``prune!, prune``: Merge channels with redundant header fields.
 
-* ``prune!, prune``: Merges channels with redundant header fields.
+* ``purge!, purge``: Delete all channels with no data (defined for any channel ``i`` by ``isempty(S.x[i]) == true``).
 
-* ``purge!, purge``: Deletes all channels with no data (defined for any channel i as isempty(S.x[i]) == true).
-
-* ``sync!, sync``: Synchronize time windows for all channels and fill time gaps. Calls ungap at invocation.
+* ``sync!, sync``: Synchronize time windows for all channels and fill time gaps. Calls ``ungap!`` at invocation.
 
 * ``ungap!, ungap``: Fill all time gaps in each channel of regularly sampled data.
 
+* ``autotap!``: Cosine taper time series data around time gaps.
+
 
 
 Native File I/O

src/SeisData/SeisData.jl

@@ -11,48 +11,8 @@ filter, filt!, sort!, sort #, isempty
 """
     S = SeisObj()
 
-Create a single channel instance of univariate geophysical data. A SeisObj has
-the following fields, which can be set by name at creation; for example,
-T = SeisObj(gain=1.0e6) creates a SeisObj with T.gain = 1.0e6.
-
-* name: ASCIIString. Freeform; strings > 26 characters long cannot be saved to
-SeisData files.
-* id: ASCIIString. ID must follow the convention `net`.`sta`.`loc`.`chan`;
-If unsure of a value, leave it blank.
-+ `network` is a 2-character network code.
-+ `station` is a 5-character (maximum) station code. For strict compliance with
-SEED standards, station codes shouldn't contain punctuation or non-standard
-characters.
-+ `loc` is a 2-character SEED-style location code. Not widely used.
-+ `chan` is a 3-character channel code. If you only know the orientation, try
-`getbandcode(fs)` for the first character; if you also know the critical
-frequency FC, `getbandcode(fs, fc=FC)`.
-* fs: Float64. Sampling frequency in Hz.
-+ For irregularly sampled data, such as "campaign" style gas flux, set fs = 0.
-* gain: Float64. Stage 0 scalar gain.
-* loc: 5-element vector. The first three elements should be [latitude [°N],
-longitude [°E], elevation [m]]. For data channels from a three-component
-seismometer with orthogonal channels, loc[4] should be the azimuth of the
-horizontal components [N°E]; loc [5] should be the incidence angle of the
-vertical component [° from vertical].
-* misc: Dictionary with ASCII keys. Can hold any type of value, but only
-characters, strings, numbers, and arrays will be saved to/read from file.
-* notes: An array of strings. A log of channel information as data is processed.
-The command `note` appends a timestamped note.
-* resp: Instrument frequency response, given as a matrix with complex zeros in
-the first column and complex poles in the second.
-* src: ASCII string. The source of the data. Normally this is filled in
-automatically when data are loaded in.
-* t: Sparse two-column array of delta-encoded times.
-+ For regularly sampled data (S.fs > 0), the first column gives the index
-in the data (S.x) of the first value after each time gap. The second column
-gives the gap length in seconds.
-+ For irregularly sampled data (S.fs == 0), the first column is meaningless;
-the second column gives delta-encoded sample times.
-+ For all data, S.t[1,2] is the start of the time series, measured from Unix
-Epoch time (1 January 1970).
-* x: Array of data points.
-* units: Freeform string with the unit type.
+Create a single channel instance of univariate geophysical data. See
+documentation for information on fields.
 """
 type SeisObj
   name::ASCIIString
@@ -90,37 +50,7 @@ end
 
 Create a multichannel structure for univariate geophysical data. A SeisData
 structure has the same fields as a SeisObj, but they cannot be set at creation.
-See the help for `SeisObj` for details of field names and meanings.
-
-### Creating and modifying SeisData structures
-You can create a SeisData structure piecewise from individual channels by
-concatenating SeisObj structures:
-* `S = SeisData(SeisObj(name="BRASIL"), SeisObj(name="UKRAINE"), SeisObj())``
-**creates** a new SeisData structure with three channels; the first is named
-"BRASIL", the second "UKRAINE", the third is blank.
-* You can **merge** SeisObj and SeisData structures into existing SeisData
-structures with the addition operator. `S = SeisData(); T = SeisObj(); S += T`
-merges T into S.
-+ If a merged SeisObj structure has the same ID field as a channel in an existing
-SeisData structure, the data are meged using the `.t` fields. Data pairs separated
-by less than half the sampling interval are *averaged*.
-+ To append a SeisObj `T` as a new channel in `S` without merging, use `append!(S,T)`.
-+ You can create a new SeisData structure by adding two or more SeisObj instances
-together, e.g. `S = SeisObj(); T = SeisObj(); U = SeisObj(); R = S + T + U`
-* You can **remove** a channel from a SeisData structure in three ways:
-+ `T = S[i]` creates a SeisObj out of all data in channel `i` of `S`;
-`T = S[i1:i2]` creates a SeisData object from multiple channels. Both methods
-leave S intact.
-+ `S -= i`, where `i` is an integer, deletes channel `i` from S.
-+ `T = pull(S, name)`, where `name` is a string, creates a SeisObj `T` from
-the first channel with name=`name`. The channel is removed from `S`.
-+ `T = pull(S, i)`, where `i` is an integer, creates a SeisObj `T` from
-channel `i` and removes channel `i` from `S`.
-
-### Differences between SeisData and SeisObj interactions
-* The command `note(S, i, [note])` adds a timestamped note to channel `i` of `S`.
-* The command `S += "[name]: [note]"`, adds a note with a short form date.
-* The command `S += "[note]"` adds a note with short-form date to all channels.
+
 """
 type SeisData
   n::Int64

src/SeisData/misc.jl

@@ -67,7 +67,7 @@ Fill gaps in x, as specified in t, assuming sampling rate fs
 """
 function gapfill!(x::Array{Float64,1}, t::Array{Float64,2}, fs::Float64; m=true::Bool, w=true::Bool)
   (fs == 0 || fs == Inf) && (return x)
-  mx = m ? mean(x) : NaN
+  mx = m ? mean(x[!isnan[x]]) : NaN
   u = round(Int, max(20,0.2*fs))
   for i = size(t,1):-1:2
     g = t[i,2]
@@ -296,7 +296,7 @@ Synchronize S and downsample data to the lowest non-null interval in S.fs.
 sync(S::SeisData; r=false::Bool) = (T = deepcopy(S); sync!(T, resample=r);
   return T)
 
-function autotuk!(x, v, u, m)
+function autotuk!(x, v, u)
   g = find(diff(v) .> 1)
   L = length(g)
   if L > 0
@@ -311,50 +311,60 @@ function autotuk!(x, v, u, m)
         x[j+1:k] .*= tukey(N, u/N)
       else
         warn(string(@sprintf("Channel %i: Time window too small, ",i),
-          @sprintf("x[%i:%i]; replaced with mean.", j+1, k)))
-        x[j+1:k] = m
+          @sprintf("x[%i:%i]; replaced with zeros.", j+1, k)))
+        x[j+1:k] = 0
       end
     end
   end
   return x
 end
 
 """
-autotap!(U)
+    !autotap(U)
 
-Automatically cosine taper (Tukey window) any
+Automatically cosine taper (Tukey window) all data in U
 """
 function autotap!(U::SeisObj)
-  mx = mean(U.x[find(!isnan(U.x))])
+  U.fs == 0 && return
+  j = find(!isnan(U.x))
+  mx = mean(U.x[j])
   u = round(Int, max(20,0.2*U.fs))
 
-  # First check for auto-fill values (i.e. values that don't change)
-  autotuk!(U.x, find(diff(U.x).!=0), u, mx)
+  # First remove the mean
+  U.x -= mx
+
+  # Then check for auto-fill values (i.e. values that don't change) and NaNs
+  # autotuk!(U.x, find(diff(U.x).!=0), u)
 
   # Then check for NaNs
-  autotuk!(U.x, find(!isnan(U.x)), u, mx)
+  autotuk!(U.x, find(!isnan(U.x)), u)
 
-  # Then replace NaNs with the mean
-  U.x[find(isnan(U.x))] = mx
+  # Then replace NaNs with zeros
+  U.x[find(isnan(U.x))] = 0
 
   # And note it
-  note(U, "Auto-tapered data.")
+  note(U, "Auto-tapered data; replaced NaNs with zeros.")
   return U
 end
 function autotap!(U::SeisData)
   for i = 1:U.n
-    mx = mean(U.x[i][find(!isnan(U.x[i]))])
+    U.fs[i] == 0 && continue
+    j = find(!isnan(U.x[i]))
+    mx = mean(U.x[i][j])
     u = round(Int, max(20,0.2*U.fs[i]))
 
-    # First check for auto-fill values (i.e. values that don't change)
-    autotuk!(U.x[i], find(diff(U.x[i]).!=0), u, mx)
+    # Remove mean
+    U.x[i] -= mx
+
+    # Check for auto-fill values (i.e. values that don't change); taper around them
+    # autotuk!(U.x[i], find(diff(U.x[i]).!=0), u)
 
-    # Then check for NaNs
-    autotuk!(U.x[i], find(!isnan(U.x[i])), u, mx)
+    # Check for NaNs
+    autotuk!(U.x[i], find(!isnan(U.x[i])), u)
 
-    # Then replace NaNs with the mean
-    U.x[i][find(isnan(U.x[i]))] = mx
-    note(U, i, "Auto-tapered data.")
+    # Then replace NaNs with zeros
+    U.x[i][find(isnan(U.x[i]))] = 0
+    note(U, i, "Auto-tapered data; replaced NaNs with zeros.")
   end
   return U
 end