crawtools.clockdrift.driftcorrect

Apply a linear drift correction to one station's arrival picks

Correction is zero at starttime, offset at endtime and linear in between

driftcorrect 

driftcorrect(catalog_file, station, starttime, endtime, offset, out_filename='out_catalog.qml')

correct time residuals on one station

:param catalog_file: name of QuakeML file with events :param station: station :param starttime: sychronization start time (assumed perfectly synchronized) :param endtime: synchronization end time :param offset: seconds to add to pick times at endtime :param out_filename: output catalog name

Source code in crawtools/clockdrift/driftcorrect.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def driftcorrect(catalog_file, station, starttime, endtime, offset,
                 out_filename='out_catalog.qml'):
    """
    correct time residuals on one station

    :param catalog_file: name of QuakeML file with events
    :param station: station
    :param starttime: sychronization start time (assumed perfectly synchronized)
    :param endtime: synchronization end time
    :param offset: seconds to add to pick times at endtime
    :param out_filename: output catalog name
    """
    catalog = read_events(catalog_file)
    starttime = UTCDateTime(starttime)
    endtime = UTCDateTime(endtime)
    assert endtime > starttime
    sync_span = endtime - starttime
    print(f'{offset=}, {sync_span/86400=}')

    firstpick = endtime
    lastpick = starttime
    for event in catalog:
        for pick in event.picks:
            if pick.waveform_id.station_code == station:
                firstpick = min(pick.time, firstpick)
                lastpick = max(pick.time, lastpick)
                pick.time += offset * (pick.time - starttime) / sync_span
    print(f'first, last pick times were {firstpick}, {lastpick}')
    assert firstpick > starttime
    assert lastpick < endtime
    catalog.write(out_filename, format='QUAKEML') 
    print(f'Wrote to {out_filename}')

crawtools.clockdrift.driftplot

Plot EQ catalog time residuals as a function of time

driftplot 

driftplot(catalog_file, stations=None, phases=['[pP].*', '[sS].*'])

Plots time residuals as a function of time

:param catalog_file: name of QuakeML file with events :param stations: list of stations to limit to :param phases: regexps of phases to plot (up to 2)

Source code in crawtools/clockdrift/driftplot.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def driftplot(catalog_file, stations=None, phases=[r'[pP].*', r'[sS].*']):
    """
    Plots time residuals as a function of time

    :param catalog_file: name of QuakeML file with events
    :param stations: list of stations to limit to
    :param phases: regexps of phases to plot (up to 2)
    """
    catalog = read_events(catalog_file)
    if stations is None:
        stations = _get_stations(catalog)
        print(f'found {stations=}')

    fig, axs = plt.subplots(len(stations), 1, sharex=True)
    colors='brgmcyk'
    for station, ax in zip(stations, axs):
        times, tress = [], []
        for phase in phases:
            times.append([])
            tress.append([])
        for event in catalog:
            for pick in event.picks:
                if pick.waveform_id.station_code == station:
                    for i in range(len(phases)):
                        time, tres = _get_pick(pick, event.origins[0], phases[i])
                        if time is not None:
                            times[i].append(time)
                            tress[i].append(tres)
        for i in range(len(phases)):
            ax.scatter(times[i], tress[i], s = 2.5, edgecolor = colors[i],
                       label=phases[i])
        ax.axhline(y=0, color='k', linestyle='-')
        ax.set_ylabel(station)
        ax.set_ylim(-1, 1)
    axs[0].set_title("Time residuals")
    axs[0].legend()
    axs[-1].set_xlabel("Origin Time")
    plt.savefig(f"{catalog_file}_drift.png", format='png')
    plt.show()