Skip to content

Commit

Permalink
time stamps switched to microsec.
Browse files Browse the repository at this point in the history
  • Loading branch information
jpjones76 committed May 25, 2016
1 parent 768042c commit 7ed72c6
Show file tree
Hide file tree
Showing 18 changed files with 173 additions and 160 deletions.
2 changes: 1 addition & 1 deletion docs/src/seisdata.rst
Original file line number Diff line number Diff line change
Expand Up @@ -149,4 +149,4 @@ Native File I/O
===============
Use ``rseis`` and ``wseis`` to read and save in native format. Use ``writesac(S)`` to save trace data in ``S`` to single-channel :ref:`SAC <sac1>` files.

SAC is widely used and well-supported, but writes in single-precision format. Rudimentary time stamping is enabled by default. Time stamped SAC data are treated by SAC itself as unevenly spaced, generic ``x-y`` data (``LEVEN=0, IFTYPE=4``). However, third-party SAC readers interpret such files less predictably: timestamped data *might* be loaded as the real part of a complex time series, with the time values themselves as the imaginary part...or the other way around...or not at all.
The SAC file format is widely used and well-supported, but writes in single-precision format. Rudimentary time stamping is enabled by default. Time stamped SAC files from SeisIO are treated by the SAC program as unevenly spaced, generic ``x-y`` data (``LEVEN=0, IFTYPE=4``). However, third-party readers might interpret timestamped files less predictably: timestamped data *might* be loaded as the real part of a complex time series, with the time values stored as the imaginary part...or they might load the other way around...or they might not load at all.
12 changes: 7 additions & 5 deletions docs/src/seisdata_fields.rst
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ Field Names and Meanings

* ``notes``: array of notes

* ``t``: pseudo-sparse two-column array of times
* ``t``: pseudo-sparse two-column array of times :ref:`(see below) <seisdata_t>`

* ``x``: vector of sample data

Expand Down Expand Up @@ -82,14 +82,16 @@ Scalars and single-type arrays of the follwoing types are OK: ``Char, Complex64,
---------
The ``notes`` array grows automatically as data are modified. Any ``merge`` or ``sync`` operation is automatically noted. Additional notes can be written manually using the ``note`` command.

.. _seisdata_t:

``t``
-----
For time series data (``fs > 0``), ``t`` is a sparse delta-encoded representation of *time gaps* in the correspoding data ``x``. The first column stores *indices*; the second column stores gap lengths in *seconds*.
For *time series data* ``x``, ``t`` is a sparse delta-compressed representation of *time gaps in microseconds* in ``x``. The first column stores *indices*; the second column stores *gap lengths*.

* The second column of the first row, i.e. ``t[1,2]``, always stores the time of the first sample in ``x`` relative to the Unix epoch (0.0 = 1970-01-01T00:00:00).
* The second column of the first row, i.e. ``t[1,2]``, always stores the *time of the first sample* in ``x`` relative to the Unix epoch (e.g. ``t[1,2] = 10`` means the time series begins at 1970-01-01T00:00:00.000010 UTC).

* The last row of ``t`` should always take the form ``[length(x) 0.0]``.

* Other rows should take the form ``[(start of gap in x) (length of gap)]``.
* Other rows should take the form ``[(starting index of gap in x) (length of gap)]``.

For irregularly sampled data (``fs = 0``), ``t`` is a sparse delta-encoded representation of time stamps, one per sample. In this case, only the second column of ``t`` is used. The first column is usually a column of zeros.
For *irregularly sampled data* (``fs = 0``), ``t`` is a dense delta-compressed representation of *time stamps in microseconds for each sample*. ``t`` for irregularly sampled data has a single column, but is treated by Julia as a two-dimensiona array (not a one-dimensional vector).
6 changes: 3 additions & 3 deletions docs/src/seisdata_fileformat.rst
Original file line number Diff line number Diff line change
Expand Up @@ -181,10 +181,10 @@ As with string arrays in ``misc``, split notes using the value of ``a`` as a cha
:widths: 5, 32, 5, 5

``nt``, number of gaps, i64, 1
``ti``, time gap indices, f64, ``nt``
``tv``, time gap values, f64, ``nt``
``ti``, time gap indices, i64, ``nt``
``tv``, time gap values, i64, ``nt``

``t`` is a two-column array, ``t = [ti tv]``.
If ``fs>0, t`` is a two-column array, ``t = [ti tv]``; otherwise, Julia reshapes ``t`` to a one-column, two-dimensional array.

``x``
^^^^^
Expand Down
4 changes: 2 additions & 2 deletions src/FileFormats/LennartzAsc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,15 @@ function rlennasc(fname::ASCIIString)

sta = replace(h[3],"\'", "")
S.fs = 1000/parse(h[5])
ts = Dates.datetime2unix(DateTime(join([h[8],"T",h[9]])))
ts = round(Int, Dates.datetime2unix(DateTime(join([h[8],"T",h[9]])))/μs)

cmp = split(fname,'.')[end]
x = readdlm(fid)
close(fid)

S.name = join([sta, cmp], '.')
S.id = join(["", sta, "", cmp], '.')
S.t = map(Float64, [1.0 ts; length(S.x) 0.0])
S.t = [1 ts; length(S.x) 0]
S.x = x[:,1]
S.src = "lennartz ascii"
return S
Expand Down
2 changes: 1 addition & 1 deletion src/FileFormats/SAC.jl
Original file line number Diff line number Diff line change
Expand Up @@ -280,7 +280,7 @@ function sactoseis(D::Dict{ASCIIString,Any})
loc = [D["stla"], D["stlo"], D["stel"], 0, 0]
end
x = map(Float64, D["data"])
t = map(Float64, [1.0 sac2epoch(D); D["npts"] 0.0])
t = map(Float64, [1 sac2epoch(D)/μs; D["npts"] 0])

misc = prunesac(D)
for k in ["nvhdr", "knetwk", "kstnm", "kcmpnm", "scale",
Expand Down
2 changes: 1 addition & 1 deletion src/FileFormats/SEGY.jl
Original file line number Diff line number Diff line change
Expand Up @@ -337,7 +337,7 @@ function segytoseis(S::Dict{ASCIIString,Any})
id = join(["",sta,"",cmp],'.')
end
x = S["data"]
t = map(Float64, [1.0 sac2epoch(S); length(x) 0])
t = [1 round(Int,sac2epoch(S)/μs); length(x) 0]
gain = 1.0
units = "unknown"
loc = zeros(5)
Expand Down
2 changes: 1 addition & 1 deletion src/FileFormats/UW.jl
Original file line number Diff line number Diff line change
Expand Up @@ -457,7 +457,7 @@ function uwtoseis(S::Dict{ASCIIString,Any})
src = S["src"]
x = map(Float64, S["data"][i])
notes = S["comment"]
t = map(Float64, [1.0 S["start"][i]+S["ot"]; length(x) 0])
t = [1 round(Int,(S["start"][i]+S["ot"])/μs); length(x) 0]
misc = Dict{ASCIIString,Any}()
for k in misc_keys
misc[k] = S[k][i]
Expand Down
2 changes: 1 addition & 1 deletion src/FileFormats/Win32.jl
Original file line number Diff line number Diff line change
Expand Up @@ -221,7 +221,7 @@ function win32toseis(D = Dict{ASCIIString,Any}())
# for an example of the (lack of) correspondence

x = map(Float64, D[k]["data"].*D[k]["scale"])
t = [1.0 D[k]["startTime"]; length(D[k]["data"]) 0.0]
t = [1 round(Int,D[k]["startTime"]/μs); length(D[k]["data"]) 0]
src = "win32"
notes = [string(now, " Record source: ", D[k]["netName"]);
string(now, " Location comment: ", D[k]["comment"]);
Expand Down
8 changes: 4 additions & 4 deletions src/FileFormats/mSEED.jl
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ function parserec(S::SeisData, sid; v=false::Bool, vv=false::Bool, fmt=10::Int)
# Lazy coding; I assume fs doesn't change
channel = channel[1]
dt = 1/S.fs[channel]
te = sum(S.t[channel][:,2]) + length(S.x[channel])*dt
te = sum(S.t[channel][:,2]) + length(S.x[channel])*round(Int,dt/μs)
end

# =========================================================================
Expand Down Expand Up @@ -213,13 +213,13 @@ function parserec(S::SeisData, sid; v=false::Bool, vv=false::Bool, fmt=10::Int)
#te = ts + nsamp*dt
vv && println("ts = ", ts, " te = ", te)
if te == 0
S.t[channel] = [1.0 ts; Float64(nsamp) 0.0]
S.t[channel] = [1 round(Int,ts/μs); nsamp 0]
else
S.t[channel] = S.t[channel][1:end-1,:]
if ts-te > dt
S.t[channel] = [S.t[channel]; [length(S.x[channel])+1.0 ts-te]]
S.t[channel] = [S.t[channel]; [length(S.x[channel])+1 ts-te]]
end
S.t[channel] = [S.t[channel]; [Float64(length(S.x[channel])+nsamp) 0.0]]
S.t[channel] = [S.t[channel]; [length(S.x[channel])+nsamp 0]]
end

# Data Read
Expand Down
75 changes: 6 additions & 69 deletions src/SeisData/SeisData.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ type SeisObj
notes::Array{ASCIIString,1}
resp::Array{Complex{Float64},2}
src::ASCIIString
t::Array{Float64,2}
t::Array{Int64,2}
x::Array{Float64,1}
units::ASCIIString

Expand All @@ -37,7 +37,7 @@ type SeisObj
notes=Array{ASCIIString,1}()::Array{ASCIIString,1},
resp=Array{Complex{Float64},2}(0,2)::Array{Complex{Float64},2},
src=""::ASCIIString,
t=Array{Float64,2}(0,2)::Array{Float64,2},
t=Array{Int64,2}(0,2)::Array{Int64,2},
x=Array{Float64,1}()::Array{Float64,1},
units=""::ASCIIString) = begin
return new(name, id, fs, gain, loc, misc, notes, resp, src, t, x, units)
Expand All @@ -63,7 +63,7 @@ type SeisData
notes::Array{Array{ASCIIString,1},1} # notes
resp::Array{Array{Complex{Float64},2},1} # resp
src::Array{ASCIIString,1} # src
t::Array{Array{Float64,2},1} # time
t::Array{Array{Int64,2},1} # time
x::Array{Array{Float64,1},1} # data
units::Array{ASCIIString,1} # units

Expand All @@ -78,7 +78,7 @@ type SeisData
notes=Array{Array{ASCIIString,1},1}(), # notes
resp=Array{Array{Complex{Float64},2},1}(), # resp
src=Array{ASCIIString,1}(), # src
t=Array{Array{Float64,2},1}(), # time
t=Array{Array{Int64,2},1}(), # time
x=Array{Array{Float64,1},1}(), # data
units=Array{ASCIIString,1}() # units
) = begin
Expand Down Expand Up @@ -275,70 +275,7 @@ delete!(S::SeisData, J::Array{Int,1}) = deleteat!(S, J)
# ============================================================================
# Merge and extract
# Dealing with sparse time difference representations
"""
t = t_expand(T::Array{Float64,2}, dt::Real)
Expand sparse-difference (SeisData-style) time stamp representation t to full
time stamps.
"""
function t_expand(t::Array{Float64,2}, dt::Real)
dt == Inf && return cumsum(t[:,2])
i = round(Int, t[:,1])
tt = dt.*ones(i[end])
tt[i] += t[:,2]
return cumsum(tt)
end

"""
t = t_collapse(T::Array{Float64,1}, dt::Real)
Collapse full time stamp representation t to SeisData sparse-difference
representation t.
"""
function t_collapse(t::Array{Float64,1}, dt::Real)
ts = map(Float32, [dt; diff(t)])
L = length(t)
i = find([!isapprox(ts[i],Float32(dt)) for i = 1:1:length(t)])
tt = cat(1, [1.0 t[1]], [map(Float64, i) ts[i]-dt])
(isempty(i) || i[length(i)] != Float64(L)) && (tt = cat(1, tt,
[Float64(L) 0.0]))
return tt
end

function xtmerge(t1::Array{Float64,2}, x1::Array{Float64,1},
t2::Array{Float64,2}, x2::Array{Float64,1}, dt::Float64)
t = [t_expand(t1, dt); t_expand(t2, dt)]
if dt == Inf
dt = 0
end
x = [x1; x2]

# Sort
i = sortperm(t)
t1 = t[i]
x1 = x[i]

# Resolve conflicts
if minimum(diff(t1)) < 0.5*dt
J = flipdim(find(diff(t1) .< 0.5*dt), 1)
for j in J
t1[j] = 0.5*(t1[j]+t1[j+1])
if isnan(x1[j])
x1[j] = x1[j+1]
elseif !isnan(x1[j+1])
x1[j] = 0.5*(x1[j]+x1[j+1])
end
deleteat!(t1, j+1)
deleteat!(x1, j+1)
end
end
if 0 < dt < Inf
t1 = t_collapse(t1, dt)
else
t1 = [zeros(length(t1)) t1]
end
return (t1, x1)
end

"""
merge!(S::SeisObj, T::SeisObj)
Expand All @@ -361,7 +298,7 @@ function merge!(S::SeisObj, U::SeisObj)
if !isapprox(S.gain,T.gain)
(T.x .*= (S.gain/T.gain); T.gain = copy(S.gain)) # rescale T.x to match S.x
end
(S.t, S.x) = xtmerge(S.t, S.x, T.t, T.x, 1/T.fs) # merge time and data
(S.t, S.x) = xtmerge(S.t, S.x, T.t, T.x, T.fs) # merge time and data
merge!(S.misc, T.misc) # merge misc
S.notes = cat(1, S.notes, T.notes) # merge notes
note(S, @sprintf("Merged %i samples", length(T.x)))
Expand All @@ -384,7 +321,7 @@ function merge!(S::SeisData, U::SeisObj)
if !isapprox(S.gain[i],T.gain)
(T.x .*= (S.gain[i]/T.gain); T.gain = copy(S.gain[i])) # rescale T.x to match S.x
end
(S.t[i], S.x[i]) = xtmerge(S.t[i], S.x[i], T.t, T.x, 1/T.fs) # time, data
(S.t[i], S.x[i]) = xtmerge(S.t[i], S.x[i], T.t, T.x, T.fs) # time, data
merge!(S.misc[i], T.misc) # misc
S.notes[i] = cat(1, S.notes[i], T.notes) # notes
note(S, i, @sprintf("Merged %i samples", length(T.x)))
Expand Down
16 changes: 10 additions & 6 deletions src/SeisData/fileio.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ function writesac(S::SeisData; ts=true, v=true)
sacFloatVals[1] = Float32(dt)
sacFloatVals[4] = Float32(S.gain[i])
sacFloatVals[6] = Float32(0)
sacFloatVals[7] = Float32(dt*length(S.x[i]) + sum(S.t[i][2:end,2]))
sacFloatVals[7] = Float32(dt*length(S.x[i]) + sum(S.t[i][2:end,2])*μs)
if !isempty(S.loc[i])
sacFloatVals[32] = S.loc[i][1]
sacFloatVals[33] = S.loc[i][2]
Expand Down Expand Up @@ -78,7 +78,7 @@ function writesac(S::SeisData; ts=true, v=true)

# Data
x = map(Float32, S.x[i])
ts && (tdata = map(Float32, t_expand(S.t[i], dt)-S.t[i][1,2]))
ts && (tdata = map(Float32, μs*(t_expand(S.t[i], dt)-S.t[i][1,2])))

# Write to file
sacwrite(fname, sacFloatVals, sacIntVals, sacCharVals, x, t=tdata, ts=ts)
Expand Down Expand Up @@ -388,11 +388,15 @@ function r_seisobj(io::IOStream)
notes = read_string_array(io)
T = read(io, Int64)
if T > 0
ti = read(io, Float64, T)
tv = read(io, Float64, T)
t = [ti tv]
if fs > 0.0
ti = read(io, Int64, T)
tv = read(io, Int64, T)
t = [ti tv]
else
t = reshape(read(io, Int64, T), T, 1)
end
else
t = Array{Float64,2}()
t = Array{Int64,2}()
end
X = read(io, Int64)
if X > 0
Expand Down

0 comments on commit 7ed72c6

Please sign in to comment.