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 Seis

Sample 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:

  • delta is the samping interval (1 / sampling rate) in s;
  • b is the start time of the trace, relative to any absolute time if any;
  • sta is a Station holding information about the site where this recording was made;
  • evt is an Event holding information about any event related to this recording, such as an earthquake;
  • picks holds a set of times of importance such as arrival picks; and
  • meta holds 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.1

Mutability

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.

Note

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.51277958154284

To 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, loc and cha are 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 be missing.
  • lon, lat and elev are the longitude and latitude (in °, positive eastwards and northwards) of the station, whilst elev is the elevation above sea level in m.
  • azi and inc together define the orientation of the channel; that is, in which spatial direction positive values point. azi is the local azimuth in ° measured clockwise from north, whilst inc is the local inclination in ° measured from the vertical downwards. For seismic data, vertical channels will have inc == 0 and azi == 0, whilst horizontal channels will have inc == 90 and azi == 0 for north and azi == 90 for east.
  • meta holds 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:

  • time is a Dates.DateTime giving an absolute date and time in UTC against which a Trace's b field (beginning time) is relative. For an earthquake or other source of energy, it should be the origin time of that event.
  • lon, lat and dep are the longitude and latitude (in °) and depth (in km) of any identified source for the data, such as an earthquake.
  • id is a String giving some identifying information about an event, and could be a catalogue ID or otherwise.
  • meta holds all other information. By convention, the following fields in meta might be used:
    • mag holds an event magnitude;
    • quakeml holds information about an event in QuakeML form, if any;
    • catalog gives 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_sac or read_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")
4000

Basic 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.

In addition, pairs of traces can be rotated:

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.0

As 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.0

data 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.4

For 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)
5

The 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))
2

Picks

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.64

Geometry

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, and
  • Cartesian, 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:

  • CartTrace creates a new Trace in an $x, y, z$ system.
  • CartStation creates a Cartesian-referenced Station.
  • CartEvent makes 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.10159724405711013

Geometry 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.

Seis.Geographic

Seis.Cartesian