User manual
This section describes how to use Seis to perform various basic processing operations on seismic data.
Preamble
The examples in this manual all assume that you have first used the module like so:
julia> using SeisSample data
Seis comes with a small selection of sample data which will be used for the examples in this manual. See the help for sample_data for details.
Types
Working with Seis requires understanding a little bit about the types it uses. Fundamentally, there are three important types:
Each of these types have fields which either must contain data, or allow data to be missing.
Traces
A Trace holds a single continuous record of evenly-sampled data. As such, it cannot contain gaps or overlaps. It may be referenced to an absolute start time, and by convention all times of day in Seis are UTC. Equally, however, a trace may have no absolute reference time.
Traces have a small number of fields which are directly accessible by users:
deltais the samping interval (1 / sampling rate) in s;bis the start time of the trace, relative to any absolute time if any;stais aStationholding information about the site where this recording was made;evtis anEventholding information about any event related to this recording, such as an earthquake;picksholds a set of times of importance such as arrival picks; andmetaholds any other information.
None of these fields can be missing.
To access the data stored in a trace, use the trace function to return the underlying array.
The user-accessible fields can be retrieved or modified as usual for composite types:
julia> # A trace starting a 0 s with sampling interval 1 s made of 100 random samples
t = Trace(0, 1, rand(100))
Seis.Trace{Float64,Vector{Float64},Seis.Geographic{Float64}}:
b: 0.0
delta: 1.0
GeogStation{Float64}:
sta.meta: Seis.SeisDict{Symbol, Any}()
GeogEvent{Float64}:
evt.meta: Seis.SeisDict{Symbol, Any}()
Trace:
picks: 0
meta:
julia> # Update the sampling interval to 0.1 s
t.delta = 0.1
0.1Mutability
Traces are mutable. This means that the values they hold can be changed by the user at will. It is often useful to copy a trace after you have performed an operation on it. Thoughout Seis, as with Julia generally, functions which modify a Trace have an exclamation mark (!) at the end of their name.
Traces are 'parameterised' on the type of floating point numbers they use, the type of data storage, and the geometry in which they are defined (i.e., either in geographic or Cartesian coordinates). A Trace is technically a Trace{T, V, P} where {T, V, P}, where T is the number type used, V is the type of vector holding the trace data, and P is the geometry type. This is why a trace has a scary-looking type like Trace{Float32, Array{Float32, 1}, Seis.Geographic{Float32}}. You only need to deal with these type parameters if you want to change these defaults. See Geometry for more information.
Collections of Traces
In Seis, methods are usually written to accept Traces. There is no special container for several traces; instead you can use Arrays of traces just as you would any other collection of types.
For example, to find the epicentral distance between the earthquake and a number of recordings of it at different stations, you could use Julia's broadcasting feature to call the distance_deg function on each Trace. (Note that the sample_data function returns a simple Vector{Trace} of recordings.)
julia> t = sample_data(:array)
60-element Vector{Trace{Float32, Vector{Float32}, Seis.Geographic{Float32}}}:
Seis.Trace(.ABA..SHZ: delta=0.05, b=996.73, nsamples=6002)
Seis.Trace(.APA..SHZ: delta=0.05, b=996.87, nsamples=6002)
Seis.Trace(.AWI..SHZ: delta=0.05, b=996.86, nsamples=6002)
Seis.Trace(.BBH..SHZ: delta=0.05, b=996.34, nsamples=6002)
Seis.Trace(.BBO..SHZ: delta=0.05, b=996.51, nsamples=6002)
Seis.Trace(.BDL..SHZ: delta=0.05, b=995.81, nsamples=6002)
Seis.Trace(.BTA..SHZ: delta=0.05, b=995.8, nsamples=6002)
Seis.Trace(.BWH..SHZ: delta=0.05, b=996.1, nsamples=6002)
Seis.Trace(.CRA..SHZ: delta=0.05, b=999.85, nsamples=6002)
Seis.Trace(.CSF..SHZ: delta=0.05, b=995.22, nsamples=6002)
⋮
Seis.Trace(.WCB..SHZ: delta=0.05, b=987.78, nsamples=6002)
Seis.Trace(.WME..SHZ: delta=0.05, b=987.89, nsamples=6002)
Seis.Trace(.WPM..SHZ: delta=0.05, b=987.8, nsamples=6002)
Seis.Trace(.XAL..SHZ: delta=0.05, b=994.08, nsamples=6002)
Seis.Trace(.XDE..SHZ: delta=0.05, b=994.94, nsamples=6002)
Seis.Trace(.YEL..SHZ: delta=0.05, b=996.79, nsamples=6002)
Seis.Trace(.YLL..SHZ: delta=0.05, b=988.17, nsamples=6002)
Seis.Trace(.YRC..SHZ: delta=0.05, b=988.03, nsamples=6002)
Seis.Trace(.YRE..SHZ: delta=0.05, b=988.29, nsamples=6002)
julia> typeof(t)
Vector{Trace{Float32, Vector{Float32}, Seis.Geographic{Float32}}} (alias for Array{Trace{Float32, Array{Float32, 1}, Seis.Geographic{Float32}}, 1})
julia> eltype(t)
Trace{Float32, Vector{Float32}, Seis.Geographic{Float32}}
julia> distance_deg.(t)
60-element Vector{Float64}:
150.77944621926457
151.35719900309945
150.82652363302728
148.47569756225323
148.85352129112246
148.80415824159704
148.7142987502766
148.39243651154584
153.2137411026592
149.1411942978774
⋮
150.1093711892705
150.11145961826915
150.27951531203723
148.7781191276318
149.06873665863068
143.1208080990685
150.3761745994614
150.2332219600003
150.51277958154284To find the nearest station, you could write the following:
julia> t[argmin(distance_deg.(t))].sta
Seis.Station{Float32,Seis.Geographic{Float32}}:
lon: -1.083
lat: 60.5509
dep: missing
net: missing
sta: YEL
loc: missing
cha: SHZ
elev: missing
azi: 0.0
inc: 0.0
meta:Stations
Stations define a stationary point in space where a recording was made. They have the following fields which you can access and modify:
net,sta,locandchaare used respectively for the network, station, channel and location names of a single recording channel. They correspond to the SEED channel naming convention, but all can bemissing.lon,latandelevare the longitude and latitude (in °, positive eastwards and northwards) of the station, whilstelevis the elevation above sea level in m.aziandinctogether define the orientation of the channel; that is, in which spatial direction positive values point.aziis the local azimuth in ° measured clockwise from north, whilstincis the local inclination in ° measured from the vertical downwards. For seismic data, vertical channels will haveinc == 0andazi == 0, whilst horizontal channels will haveinc == 90andazi == 0for north andazi == 90for east.metaholds all other information.
See Geometry for more details on other coordinate systems in which objects can be placed in Seis.
Events
An Event marks the source of a Trace. Typically, it represents the source of energy for a recording, or it may simply indicate the start time of a section of data with no energy source implied.
Events have the following accessible fields, all of which can be missing:
timeis aDates.DateTimegiving an absolute date and time in UTC against which aTrace'sbfield (beginning time) is relative. For an earthquake or other source of energy, it should be the origin time of that event.lon,latanddepare the longitude and latitude (in °) and depth (in km) of any identified source for the data, such as an earthquake.idis aStringgiving some identifying information about an event, and could be a catalogue ID or otherwise.metaholds all other information. By convention, the following fields inmetamight be used:magholds an event magnitude;quakemlholds information about an event in QuakeML form, if any;cataloggives the name of the catalogue for this event.
Getting data in
Data can be loaded into Julia via a number of means:
Reading data from disk in SAC or miniSEED format.
For this use either
read_sacorread_mseed.To read a single file, use the one-argument form:
julia> t1 = read_sac("single_file.sac") Seis.Trace{Float32,Array{Float32,1},Seis.Geographic{Float32}}: b: 52.66 delta: 0.01 Station{Float32,Seis.Geographic{Float32}}: sta.lon: -120.0 sta.lat: 48.0 sta.sta: CDV sta.azi: 0.0 sta.inc: 0.0 sta.meta: Seis.SeisDict{Symbol,Any}() Event{Float32,Seis.Geographic{Float32}}: evt.lon: -125.0 evt.lat: 48.0 evt.dep: 0.0 evt.time: 1981-03-29T10:38:14 evt.id: K8108838 evt.meta: Seis.SeisDict{Symbol,Any}() Trace: picks: 2 meta: SAC_lpspol => true SAC_nevid => 0 SAC_iftype => 1 file => "../test/test_data/seis.sac" SAC_idep => 50 SAC_iztype => 9 SAC_lcalda => true SAC_unused18 => false SAC_lovrok => true SAC_norid => 0 SAC_ievtyp => 42 julia> t2 = read_mseed("file.mseed") 2-element Array{Trace{Float32,Array{Float32,1},Seis.Geographic{Float32}},1}: Seis.Trace(GB.CWF..BHZ: delta=0.02, b=0.0, nsamples=3000) Seis.Trace(GB.CWF..HHZ: delta=0.01, b=0.0, nsamples=6000)Note SAC files can contain only one single continuous data channel, whilst miniSEED files can contain more than one, and so an array of
Traces are returned.Reading several files using a globbing pattern:
julia> read_sac("Event_*/*Z.sac", "DATA")Downloading data from a remote server.
For this, install and use SeisRequests.
julia> using SeisRequests julia> t = get_data(code="IU.ANMO.00.BH?", starttime="2018-02-02", endtime="2018-02-02T01:00:00") # an hour of data [ Info: Request status: Successful request, results follow 3-element Array{Trace{Float64,Array{Float64,1},Seis.Geographic{Float64}},1}: Seis.Trace(IU.ANMO.00.BH1: delta=0.05, b=0.0, nsamples=72000) Seis.Trace(IU.ANMO.00.BH2: delta=0.05, b=0.0, nsamples=72000) Seis.Trace(IU.ANMO.00.BHZ: delta=0.05, b=0.0, nsamples=72000)
Creating data from scratch
Traces can be constructed using the constructors. For example:
julia> b = 0
0
julia> delta = 0.01
0.01
julia> data = randn(1000)
1000-element Vector{Float64}:
0.5401856622505432
-0.3492438171703137
0.2883887799457867
-1.9347726382482167
-0.23171919888298964
0.19176078111150419
-0.36921079828000164
0.7669061174009048
-1.2594193802034581
-0.1101340008722498
⋮
2.7257358191654424
0.14469796969865903
-1.6505614584001924
0.2057503705160407
1.5435005909794588
0.022089478606524257
-0.9553612148485425
-1.1827823968452387
-0.6147305351467764
julia> t = Trace(b, delta, data) # Fill trace with data already available
Seis.Trace{Float64,Vector{Float64},Seis.Geographic{Float64}}:
b: 0.0
delta: 0.01
GeogStation{Float64}:
sta.meta: Seis.SeisDict{Symbol, Any}()
GeogEvent{Float64}:
evt.meta: Seis.SeisDict{Symbol, Any}()
Trace:
picks: 0
meta:Writing data out
Seis currently supports writing in SAC format, using the write_sac function.
julia> t = sample_data()
Seis.Trace{Float32,Vector{Float32},Seis.Geographic{Float32}}:
b: 52.66
delta: 0.01
GeogStation{Float32}:
sta.lon: -120.0
sta.lat: 48.0
sta.sta: CDV
sta.azi: 0.0
sta.inc: 0.0
sta.meta: Seis.SeisDict{Symbol, Any}()
GeogEvent{Float32}:
evt.lon: -125.0
evt.lat: 48.0
evt.dep: 0.0
evt.time: 1981-03-29T10:38:14
evt.id: K8108838
evt.meta: Seis.SeisDict{Symbol, Any}()
Trace:
picks: 2
meta: SAC_lpspol => true
SAC_nevid => 0
SAC_iftype => 1
file => "/home/runner/work/Seis.jl/Seis.jl/src/../data/seis.sac"
SAC_idep => 50
SAC_iztype => 9
SAC_lcalda => true
SAC_unused18 => false
SAC_lovrok => true
SAC_norid => 0
SAC_ievtyp => 42
julia> write_sac(t, "outfile.sac")
4000Basic processing
Once read in or obtained some other way, data can be processed in a number of ways.
Note that, as mentioned above (Mutability), traces can be modified in place, or copies taken when processing steps are applied. Modifying in place is usually faster as it does not involve copying all the information, but it is often more convenient to copy traces as well. For this reason, there is always a modifying version of a function (ending in !), and a copying version, without.
- Cut the data to a certain window, in either relative time, absolute time, or relative to a time pick, with
cutorcut!. - Decimate the data, with or without an antialiasing filter (
decimate,decimate!). - Remove the mean value of the trace with
remove_meanorremove_mean!. - Remove a linear trend in the data with
remove_trendorremove_trend!. - Taper the signal with
taperortaper!. - Filter the data, performing a low-pass (
lowpass,lowpass!), high-pass (highpass,highpass!), band-pass (bandpass,bandpass!) or band-stop (bandstop,bandstop!) filter. For each, an acausal version of the filter can be obtained by passing thetwopass=trueoption. - Differentiate (
differentiate,differentiate!) or integrate (integrate,integrate!) the trace. - Normalise the data to a certain value with
normalizeornormalize!(ornormalisefor us Brits). - Take the envelope of a trace with
envelopeorenvelope!.
In addition, pairs of traces can be rotated:
rotate_throughandrotate_through!rotate a pair of horizontal traces by a certain angle.rotate_to_gcpandrotate_to_gcp!rotate pairs to radial (pointing away from the source at the station) and transverse (90° clockwise from this) components.
Accessing raw trace data
Beyond these basic functions, Traces are designed to be used to build your own processing workflows on top. For example, shear wave splitting analysis can be done with SeisSplit and array processing with Beamforming.
To access the raw trace, use the trace function, which returns the data:
julia> t = Trace(0, 0.1, [1, 2, 3, 4, 5])
Seis.Trace{Float64,Vector{Float64},Seis.Geographic{Float64}}:
b: 0.0
delta: 0.1
GeogStation{Float64}:
sta.meta: Seis.SeisDict{Symbol, Any}()
GeogEvent{Float64}:
evt.meta: Seis.SeisDict{Symbol, Any}()
Trace:
picks: 0
meta:
julia> data = trace(t)
5-element Vector{Float64}:
1.0
2.0
3.0
4.0
5.0As is usual with Julia, data is a variable bound to the same object as the underlying data of the trace t. Modifying data will also modify the values in t:
julia> data[1] = 0 # Set the first point to be zero
0
julia> trace(t) # Note the first point is now 0
5-element Vector{Float64}:
0.0
2.0
3.0
4.0
5.0data can be empty!d, or push!ed or append!ed to, and so on—it is simply a vector of data points.
Accessing other properties
Times
The first sample of a trace t occurs at time t.b, and each sample is t.delta seconds after the previous one. Therefore, each sample of the trace has a time, and this array (actually a subtype of AbstractRange) can be obtained with times:
julia> times(t)
0.0:0.1:0.4For consistency, the start and end times of a trace are given by starttime and endtime:
julia> starttime(t), endtime(t)
(0.0, 0.4)All these times are relative to the event date if that is set. If it is, then startdate and enddate give the start and end date of the trace in UTC:
julia> using Dates
julia> t.evt.time = DateTime(2013, 5, 24, 5, 44, 49, 880) # Set the origin time
2013-05-24T05:44:49.880
julia> startdate(t), enddate(t)
(Dates.DateTime("2013-05-24T05:44:49.880"), Dates.DateTime("2013-05-24T05:44:50.280"))Like before, we can get a set of dates for every sample with dates:
julia> dates(t)
Dates.DateTime("2013-05-24T05:44:49.880"):Dates.Millisecond(100):Dates.DateTime("2013-05-24T05:44:50.280")There should be 5 values here, since we gave the trace a set of 5 data points. nsamples tells us this:
julia> nsamples(t)
5The closest sample to a certain time or date can be found with nearest_sample:
julia> nearest_sample(t, 0.11)
2
julia> nearest_sample(t, DateTime(2013, 5, 24, 5, 44, 50))
2Picks
Seis allows you to assign arbitrary time picks to a Trace and use these in flexible ways.
Picks are simply markers which label a certain point in time with a label. They are useful for marking particular phase arrivals and for other purposes. All picks in Seis are relative to the origin time for the trace if any. So if a trace's starttime is 3 s and a pick is at 4 s, the pick points to a time 1 s after the trace starts.
Picks have two fields: time and name. time is the time in seconds after zero time, which is the event time if set. name may be a String such as "PcP", or missing if no name is needed for this pick.
Getting picks
Picks are accessible in the picks field of a trace:
julia> t = sample_data();
julia> t.picks
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 2 entries:
:F => Seis.Pick{Float32}(time=60.980003, name=missing)
:A => Seis.Pick{Float32}(time=53.670002, name=missing)You can see that there are two picks defined: one at 53.67 s with key :A, and one at 60.98 s with key :F. Each pick is associated with a key, which can be a Symbol, or an Int. Picks with a Symbol key are named picks and can be accessed in two ways.
You can access a named pick with dot notation if you know the literal key:
julia> t.picks.A
Seis.Pick{Float32}(time=53.670002, name=missing)You can also use getindex or [] notation to access a named pick using a variable bound to a Symbol or a literal value:
julia> t.picks[:A]
Seis.Pick{Float32}(time=53.670002, name=missing)
julia> pickkey = :A; t.picks[pickkey]
Seis.Pick{Float32}(time=53.670002, name=missing)For any pick, access its time and name like so:
julia> t.picks.A.time, t.picks.A.name
(53.670002f0, missing)Setting picks
Setting picks is done similarly. To add a named pick with name :S, time 2 s and name "Sg", you can do
julia> t.picks.S = 2, "Sg"
(2, "Sg")
julia> t.picks
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 3 entries:
:F => Seis.Pick{Float32}(time=60.980003, name=missing)
:A => Seis.Pick{Float32}(time=53.670002, name=missing)
:S => Seis.Pick{Float32}(time=2.0, name="Sg")Or, you can use setindex! ([]):
julia> t.picks[:S] = 2.5
2.5
julia> t.picks
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 3 entries:
:F => Seis.Pick{Float32}(time=60.980003, name=missing)
:A => Seis.Pick{Float32}(time=53.670002, name=missing)
:S => Seis.Pick{Float32}(time=2.5, name=missing)Note that the :S pick is overwritten. Note also here that by not providing a name for the pick, its name defaults to missing.
Removing picks
Remove picks as usual for Dict-like collections of things, with delete!:
julia> delete!(t.picks, :S)
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 2 entries:
:F => Seis.Pick{Float32}(time=60.980003, name=missing)
:A => Seis.Pick{Float32}(time=53.670002, name=missing)
julia> t.picks
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 2 entries:
:F => Seis.Pick{Float32}(time=60.980003, name=missing)
:A => Seis.Pick{Float32}(time=53.670002, name=missing)You can also use the special method for the picks field of setting a pick to missing:
julia> t.picks.F = missing
missing
julia> t.picks
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 1 entry:
:A => Seis.Pick{Float32}(time=53.670002, name=missing)Remove all picks with the clear_picks! function:
julia> clear_picks!(t)
Dict{Union{Int64, Symbol}, Seis.Pick{Float32}}()
julia> t.picks
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}}()Arrays of picks
Because t.picks is a special type of Dict (a Seis.SeisDict), we can also access the time and name fields of sets of picks easily.
For example, if dealing with a set of stations where we have a pick for each, we can extract the vector of all pick times simply:
julia> t = sample_data(:array)
60-element Vector{Trace{Float32, Vector{Float32}, Seis.Geographic{Float32}}}:
Seis.Trace(.ABA..SHZ: delta=0.05, b=996.73, nsamples=6002)
Seis.Trace(.APA..SHZ: delta=0.05, b=996.87, nsamples=6002)
Seis.Trace(.AWI..SHZ: delta=0.05, b=996.86, nsamples=6002)
Seis.Trace(.BBH..SHZ: delta=0.05, b=996.34, nsamples=6002)
Seis.Trace(.BBO..SHZ: delta=0.05, b=996.51, nsamples=6002)
Seis.Trace(.BDL..SHZ: delta=0.05, b=995.81, nsamples=6002)
Seis.Trace(.BTA..SHZ: delta=0.05, b=995.8, nsamples=6002)
Seis.Trace(.BWH..SHZ: delta=0.05, b=996.1, nsamples=6002)
Seis.Trace(.CRA..SHZ: delta=0.05, b=999.85, nsamples=6002)
Seis.Trace(.CSF..SHZ: delta=0.05, b=995.22, nsamples=6002)
⋮
Seis.Trace(.WCB..SHZ: delta=0.05, b=987.78, nsamples=6002)
Seis.Trace(.WME..SHZ: delta=0.05, b=987.89, nsamples=6002)
Seis.Trace(.WPM..SHZ: delta=0.05, b=987.8, nsamples=6002)
Seis.Trace(.XAL..SHZ: delta=0.05, b=994.08, nsamples=6002)
Seis.Trace(.XDE..SHZ: delta=0.05, b=994.94, nsamples=6002)
Seis.Trace(.YEL..SHZ: delta=0.05, b=996.79, nsamples=6002)
Seis.Trace(.YLL..SHZ: delta=0.05, b=988.17, nsamples=6002)
Seis.Trace(.YRC..SHZ: delta=0.05, b=988.03, nsamples=6002)
Seis.Trace(.YRE..SHZ: delta=0.05, b=988.29, nsamples=6002)
julia> t[1].picks # Picks of the first trace
Seis.SeisDict{Union{Int64, Symbol}, Seis.Pick{Float32}} with 2 entries:
:A => Seis.Pick{Float32}(time=1134.08, name=missing)
:T1 => Seis.Pick{Float32}(time=137.35, name=missing)
julia> t.picks.A.time # Times of all A picks for all traces
60-element Vector{Float32}:
1134.08
1135.62
1134.61
1128.9
1130.0
1129.5
1129.3
1128.4
1141.9
1130.67
⋮
1132.43
1132.69
1133.04
1129.33
1130.14
1114.24
1133.36
1132.82
1133.64Geometry
Seis allows you to have traces defined in any coordinate system you like. By default, it provides two main types:
Geographic, where coordinates are given as longitudes and latitudes in °, and depths in km below sea level, andCartesian, where coordinates are $x$, $y$ and $z$ in a right-handed system and each is given in m. Conventionally in Seis, $x$ is a local easting and $y$ is a local northing, meaning $z$ is the upwards direction and usually $z$ is $0$ at sea level or the local reference level.
By default in Seis, all Events, Stations and Traces are geographic, and use the Geographic type. Therefore you do not need to do anything special to work in geographic sense with longitudes, latitudes, and so on.
Cartesian geometry
However, if you are dealing with data where it makes more sense to consider positions in a Cartesian system, that is straightforward. Use the special constructors for Cartesian objects to create them:
CartTracecreates a newTracein an $x, y, z$ system.CartStationcreates a Cartesian-referencedStation.CartEventmakes a Cartesian-referenced event.
Where it makes sense, certain accessor functions like distance_km are defined for these types, and these are necessarily different to those for geographic objects.
For example, distance_km will calculate the straight-line distance between the event and station when in Cartesian coordinates:
julia> t = CartTrace(0, 1e-4, 1000)
Seis.Trace{Float64,Vector{Float64},Seis.Cartesian{Float64}}:
b: 0.0
delta: 0.0001
CartStation{Float64}:
sta.meta: Seis.SeisDict{Symbol, Any}()
CartEvent{Float64}:
evt.meta: Seis.SeisDict{Symbol, Any}()
Trace:
picks: 0
meta:
julia> t.sta.x, t.sta.y, t.sta.z = 101, -35, 12
(101, -35, 12)
julia> t.evt.x, t.evt.y, t.evt.z = 12, 14, -100
(12, 14, -100)
julia> distance_km(t)
0.10159724405711013Geometry type hierarchy and custom geometry types
Geographic and Cartesian are subtypes of Seis.Position. This means that if you want to define your own coordinate system, you are able to do so by subtyping Seis.Position and you can then write methods taking a Trace{T, V, P} where {T, V, P}, where P is your new type. These will replace the generic methods allowing special processing to be done in cases where the generic calculations need to be changed, without requiring writing duplicate methods where they do not. An example of a different coordinate system might be when stations have positions which vary in time.