Offshore pickers on ocean-bottom seismometers¶

QuakeScope will process ocean-bottom data, and ocean-bottom data is not land data with more noise. The instrument sits under a water column that reverberates, often on soft sediment, with horizontals that are arbitrarily oriented because nobody was there to point them north. Most stations also carry a hydrophone, a channel that has no land equivalent at all.

SeisBench ships three ocean-bottom pickers, and the differences between them are the interesting part:

Model Backbone Components Hydrophone
PickBlue(base="phasenet") PhaseNet Z12H yes
PickBlue(base="eqtransformer") EQTransformer Z12H yes
OBSTransformer (Niksejel & Zhang, 2024) EQTransformer ZNE no

PickBlue is a constructor rather than a model — it returns the obs weights on whichever backbone you ask for. OBSTransformer is the useful control: trained on ocean-bottom data but taking only three components, which separates trained offshore from uses the fourth channel.

Two land models quantify what running them offshore costs. All five see identical windows on three deployments in contrasting settings.

Deployment Network When Setting
Cascadia Initiative 7D 2012–13 Subduction margin, offshore northern California
AACSE XO 2018–19 Alaska Peninsula subduction, shelf to trench
Blanco X9 2012–13 Oceanic transform, strike-slip in young crust

How detection is scored, and why not against analyst picks¶

Regional catalogs do not pick temporary OBS deployments, so there are no analyst arrivals at these stations to score against. Instead each catalog event gets a predicted P arrival from iasp91, and a model counts as having detected it if it places a P pick within a tolerance of that time.

That prediction is the weak link and the tolerance has to absorb it: iasp91 has no water layer and no sediment column, both of which delay the true arrival relative to the model. Treat the numbers as relative between weight sets on identical data, not as absolute detection rates.

In [1]:
from collections import defaultdict

import matplotlib.pyplot as plt
import numpy as np
import obspy
import pandas as pd
import seisbench.models as sbm
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.geodetics import gps2dist_azimuth, locations2degrees
from obspy.taup import TauPyModel

%matplotlib inline

1. Configuration¶

In [2]:
# Three deployments, chosen for contrasting settings and instrumentation.
# Stations were selected by confirming all four components actually return
# data, not from metadata - and the hydrophone sampling rate differs by an
# order of magnitude between experiments, which section 6 tests the effect of.
EXPERIMENTS = {
    "Cascadia (7D)": dict(
        network="7D", start=UTCDateTime("2012-09-01"), end=UTCDateTime("2013-05-01"),
        lat=40.3, lon=-124.8, min_mag=3.5, max_radius_deg=2.0,
        stations=[("FS07B", "HH", "HDH"), ("FS06B", "BH", "BDH"), ("FS05B", "BH", "BDH")],
        note="Cascadia subduction margin, offshore northern California",
    ),
    "AACSE (XO)": dict(
        network="XO", start=UTCDateTime("2018-08-01"), end=UTCDateTime("2019-06-01"),
        lat=54.8, lon=-155.5, min_mag=3.5, max_radius_deg=2.0,
        stations=[("LA23", "HH", "EDH"), ("LA25", "HH", "EDH"), ("LD36", "HH", "HDH")],
        note="Alaska Peninsula subduction zone, shelf and trench",
    ),
    "Blanco (X9)": dict(
        network="X9", start=UTCDateTime("2012-10-01"), end=UTCDateTime("2013-08-01"),
        lat=43.1, lon=-126.4, min_mag=3.5, max_radius_deg=2.0,
        stations=[("BB060", "HH", "BDH"), ("BB090", "HH", "BDH"), ("BB070", "HH", "BDH")],
        note="Blanco oceanic transform - strike-slip in young oceanic crust",
    ),
}

# SeisBench ships three ocean-bottom pickers and they differ in ways that
# matter here.
#
#   PickBlue is a constructor, not a model: PickBlue(base=...) returns the
#   "obs" weights on either a PhaseNet or an EQTransformer backbone. Both take
#   four components, Z12H, the last being the hydrophone.
#
#   OBSTransformer (Niksejel & Zhang, 2024) is also OBS-trained but takes only
#   three components - no hydrophone. It is the control that separates
#   "trained on ocean-bottom data" from "uses the fourth channel".
#
# The land models quantify what running them offshore costs.
MODEL_SPECS = {
    "pickblue_phasenet": dict(build=lambda: sbm.PickBlue(base="phasenet"),
                              hydrophone=True),
    "pickblue_eqt": dict(build=lambda: sbm.PickBlue(base="eqtransformer"),
                         hydrophone=True),
    "obstransformer": dict(build=lambda: sbm.OBSTransformer.from_pretrained("obst2024"),
                           hydrophone=False),
    "quakescope2026": dict(build=lambda: sbm.PhaseNet.from_pretrained("quakescope2026"),
                           hydrophone=False),
    "original": dict(build=lambda: sbm.PhaseNet.from_pretrained("original"),
                     hydrophone=False),
}

# Detection is scored against a predicted arrival rather than an analyst pick,
# because regional catalogs do not pick temporary OBS stations. iasp91 has no
# water layer and no sediments, so the prediction is systematically early at
# these sites and the tolerance has to be generous.
ABLATION_MODEL = "pickblue_phasenet"   # the 4-component model used in section 6
GALLERY_N = 10          # windows shown per deployment in section 7

TOLERANCE = 10.0        # seconds around the predicted P
PRE, POST = 60, 120     # window around the predicted arrival

DETECT_FLOOR = 0.02     # run once here, threshold offline
REPORT_THRESHOLD = 0.3

COLORS = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]
C_P, C_S = "#2a78d6", "#eb6834"

2. Data access and arrival prediction¶

These are temporary networks with no public S3 bucket, so everything comes over EarthScope FDSN.

In [3]:
_es = Client("EARTHSCOPE", timeout=300)
_usgs = Client("USGS", timeout=300)
_taup = TauPyModel(model="iasp91")

_station_coords = {}


def station_coords(net, sta, t0, t1):
    """Coordinates for an OBS station, cached."""
    key = (net, sta)
    if key not in _station_coords:
        inv = _es.get_stations(network=net, station=sta, level="station",
                               starttime=t0, endtime=t1)
        s = inv[0][0]
        _station_coords[key] = (s.latitude, s.longitude)
    return _station_coords[key]


def predicted_p(origin, slat, slon):
    """iasp91 P arrival at the station. None if no ray reaches it."""
    deg = locations2degrees(slat, slon, origin.latitude, origin.longitude)
    depth = max((origin.depth or 0) / 1000.0, 0.0)
    arrivals = _taup.get_travel_times(source_depth_in_km=depth,
                                      distance_in_degree=deg, phase_list=["p", "P"])
    if not arrivals:
        return None, deg
    return origin.time + arrivals[0].time, deg


def fetch_obs(net, sta, band, hydro, t0, t1):
    """Four-component OBS stream: Z, 1, 2 and the hydrophone."""
    try:
        st = _es.get_waveforms(net, sta, "*", f"{band}?,{hydro}", t0, t1)
    except Exception:
        return None
    st.merge(fill_value=0)
    comps = {tr.stats.channel[-1] for tr in st}
    if not {"Z", "1", "2"} <= comps:
        return None
    return st


def as_three_component(st):
    """Land models expect ZNE. OBS horizontals are 1 and 2, arbitrarily
    oriented; renaming them is the usual convention and costs nothing here
    because neither model uses absolute orientation."""
    out = obspy.Stream()
    for tr in st:
        if tr.stats.channel[-1] == "H":          # drop the hydrophone
            continue
        tr = tr.copy()
        if tr.stats.channel[-1] == "1":
            tr.stats.channel = tr.stats.channel[:-1] + "N"
        elif tr.stats.channel[-1] == "2":
            tr.stats.channel = tr.stats.channel[:-1] + "E"
        out += tr
    return out

3. Load the weight sets¶

In [4]:
models, needs_h = {}, {}
for name, spec in MODEL_SPECS.items():
    try:
        m = spec['build']()
    except Exception as exc:
        print(f"{name:<20} could not load ({type(exc).__name__}) - skipping")
        continue
    models[name] = m
    needs_h[name] = spec['hydrophone']
    print(f"{name:<20} {type(m).__name__:<15} comps={m.component_order:<6} "
          f"in_samples={m.in_samples:<5} hydrophone={needs_h[name]}")
names = list(models)

pickblue_phasenet    PhaseNet        comps=Z12H   in_samples=3001  hydrophone=True
pickblue_eqt         EQTransformer   comps=Z12H   in_samples=6000  hydrophone=True
obstransformer       OBSTransformer  comps=ZNE    in_samples=6000  hydrophone=False
quakescope2026       PhaseNet        comps=ZNE    in_samples=3001  hydrophone=False
original             PhaseNet        comps=ENZ    in_samples=3001  hydrophone=False

4. Run every model over every deployment¶

One short window per catalogued event per station. The slow part is the network, not the inference.

In [5]:
records = []
examples, gallery = {}, {}
for label, exp in EXPERIMENTS.items():
    print(f"{label} - {exp['note']}")
    try:
        cat = _usgs.get_events(starttime=exp['start'], endtime=exp['end'],
                               latitude=exp['lat'], longitude=exp['lon'],
                               maxradius=exp['max_radius_deg'],
                               minmagnitude=exp['min_mag'])
    except Exception as exc:
        print(f"    catalog failed: {type(exc).__name__}")
        continue
    print(f"    {len(cat)} catalogued events M>={exp['min_mag']}")

    for sta, band, hydro in exp['stations']:
        try:
            slat, slon = station_coords(exp['network'], sta, exp['start'], exp['end'])
        except Exception as exc:
            print(f"    {sta}: metadata failed ({type(exc).__name__})")
            continue
        n_win = 0
        for ev in cat:
            origin = ev.preferred_origin() or (ev.origins[0] if ev.origins else None)
            if origin is None:
                continue
            tp, deg = predicted_p(origin, slat, slon)
            if tp is None:
                continue
            st = fetch_obs(exp['network'], sta, band, hydro, tp - PRE, tp + POST)
            if st is None:
                continue
            n_win += 1
            st3 = as_three_component(st)
            mag = ev.preferred_magnitude() or ev.magnitudes[0]
            has_h = any(tr.stats.channel[-1] == 'H' for tr in st)
            win_picks = {}
            for name, model in models.items():
                use = st if needs_h[name] else st3
                try:
                    out = model.classify(use, P_threshold=DETECT_FLOOR,
                                         S_threshold=DETECT_FLOOR)
                except Exception:
                    continue
                best = None
                for p in out.picks:
                    if p.phase != 'P':
                        continue
                    dt = p.peak_time - tp
                    if abs(dt) <= TOLERANCE and (best is None or
                                                float(p.peak_value) > best[1]):
                        best = (dt, float(p.peak_value))
                if best:
                    win_picks[name] = best
                records.append(dict(
                    experiment=label, station=sta, weights=name,
                    mag=round(mag.mag, 1), dist_deg=round(deg, 3),
                    has_hydrophone=has_h,
                    dt=round(best[0], 2) if best else np.nan,
                    conf=round(best[1], 3) if best else 0.0,
                    n_picks=len(out.picks),
                ))
                if name == ABLATION_MODEL and best and best[1] >= REPORT_THRESHOLD:
                    examples.setdefault(label, []).append(
                        (st, tp, sta, mag.mag, best))
            per_sta = sum(1 for g in gallery.get(label, []) if g['sta'] == sta)
            room = int(np.ceil(GALLERY_N / max(len(exp['stations']), 1)))
            if per_sta < room and len(gallery.get(label, [])) < GALLERY_N + room:
                keep = st.slice(tp - 30, tp + 70).copy()
                gallery.setdefault(label, []).append(dict(
                    stream=keep, tp=tp, sta=sta, mag=mag.mag,
                    deg=deg, picks=dict(win_picks)))
        print(f"    {sta:<6} {n_win} windows with data")
    print()

det = pd.DataFrame(records)
print(f"{len(det)} (event, station, model) rows")

Cascadia (7D) - Cascadia subduction margin, offshore northern California

    19 catalogued events M>=3.5

    FS07B  19 windows with data

    FS06B  19 windows with data

    FS05B  19 windows with data

AACSE (XO) - Alaska Peninsula subduction zone, shelf and trench

    12 catalogued events M>=3.5

    LA23   12 windows with data

    LA25   12 windows with data

    LD36   12 windows with data

Blanco (X9) - Blanco oceanic transform - strike-slip in young oceanic crust

    15 catalogued events M>=3.5

    BB060  15 windows with data

    BB090  15 windows with data

    BB070  15 windows with data

690 (event, station, model) rows

5. Detection rate¶

A detection is a P pick within the tolerance of the predicted arrival, counted at the reporting threshold.

In [6]:
hit = det[det.conf >= REPORT_THRESHOLD]
rows = []
for label in EXPERIMENTS:
    sub_all = det[det.experiment == label]
    if not len(sub_all):
        continue
    for name in names:
        total = len(sub_all[sub_all.weights == name])
        found = len(hit[(hit.experiment == label) & (hit.weights == name)])
        if not total:
            continue
        res = hit[(hit.experiment == label) & (hit.weights == name)]['dt']
        rows.append(dict(experiment=label, weights=name,
                         windows=total, detected=found,
                         rate=round(found / total, 3),
                         median_dt=round(float(res.median()), 2) if len(res) else np.nan))
summary = pd.DataFrame(rows)
print(summary.to_string(index=False))
print(f'\ndetection = a P pick within {TOLERANCE:g} s of the iasp91 prediction, '
      f'at confidence >= {REPORT_THRESHOLD}')
print('median_dt = median offset from the prediction. What matters is that it '
      'is small and similar across models, which says they lock onto the same '
      'arrival and the offset belongs to the prediction rather than the picker.')

   experiment           weights  windows  detected  rate  median_dt
Cascadia (7D) pickblue_phasenet       57        42 0.737      -1.22
Cascadia (7D)      pickblue_eqt       57        45 0.789      -1.26
Cascadia (7D)    obstransformer       57        35 0.614      -1.10
Cascadia (7D)    quakescope2026       57        41 0.719      -0.93
Cascadia (7D)          original       57        37 0.649      -0.93
   AACSE (XO) pickblue_phasenet       36        28 0.778      -1.49
   AACSE (XO)      pickblue_eqt       36        26 0.722      -1.59
   AACSE (XO)    obstransformer       36        32 0.889      -1.69
   AACSE (XO)    quakescope2026       36        28 0.778      -1.43
   AACSE (XO)          original       36        24 0.667      -1.32
  Blanco (X9) pickblue_phasenet       45        24 0.533      -3.42
  Blanco (X9)      pickblue_eqt       45        24 0.533      -3.39
  Blanco (X9)    obstransformer       45        23 0.511      -3.11
  Blanco (X9)    quakescope2026       45        24 0.533      -2.83
  Blanco (X9)          original       45        22 0.489      -3.31

detection = a P pick within 10 s of the iasp91 prediction, at confidence >= 0.3
median_dt = median offset from the prediction. What matters is that it is small and similar across models, which says they lock onto the same arrival and the offset belongs to the prediction rather than the picker.
In [7]:
if len(summary):
    labels = [l for l in EXPERIMENTS if l in set(summary.experiment)]
    fig, ax = plt.subplots(figsize=(10.5, 4.3))
    width = 0.8 / max(len(names), 1)
    for i, name in enumerate(names):
        xs, vals = [], []
        for j, lab in enumerate(labels):
            r = summary[(summary.experiment == lab) & (summary.weights == name)]
            if len(r):
                xs.append(j + (i - (len(names) - 1) / 2) * width)
                vals.append(float(r['rate'].iloc[0]))
        ax.bar(xs, vals, width=width * 0.9, color=COLORS[i % len(COLORS)], label=name)
        for x, v in zip(xs, vals):
            ax.text(x, v + 0.015, f'{v:.2f}', ha='center', fontsize=8)
    ax.set_xticks(range(len(labels)))
    ax.set_xticklabels(labels, fontsize=9)
    ax.set_ylabel('fraction of catalogued events detected')
    ax.set_ylim(0, 1.05)
    ax.grid(alpha=0.25, lw=0.5, axis='y')
    ax.set_axisbelow(True)
    ax.legend(frameon=False, fontsize=9, ncol=len(names))
    ax.set_title('Detection of catalogued events on ocean-bottom stations',
                 fontsize=11, loc='left')
    fig.tight_layout()
    plt.show()
No description has been provided for this image

Where the picks land relative to the prediction¶

What matters is that the cluster is tight and in the same place for every model — that indicates they are all locking onto the same physical arrival, and that the offset is a property of the prediction rather than of any weight set.

The offset direction is not a clean diagnostic here. A water column and sediments delay the true arrival relative to iasp91, but offshore catalog locations and depths are themselves poorly constrained — many are fixed rather than solved — and that error enters the prediction directly and in either direction. A consistent offset of a second or two says the prediction is biased, not that the picks are wrong.

In [8]:
fig, ax = plt.subplots(figsize=(10.5, 3.8))
bins = np.linspace(-TOLERANCE, TOLERANCE, 33)
for name, color in zip(names, COLORS):
    d = hit[hit.weights == name]['dt'].dropna()
    if not len(d):
        continue
    ax.hist(d, bins=bins, histtype='step', lw=1.8, color=color,
            label=f'{name} (n={len(d)})')
ax.axvline(0, color='#8a8a8a', lw=1)
ax.set_xlabel('pick minus iasp91 prediction (s)')
ax.set_ylabel('detections')
ax.set_title('Offset from the predicted arrival', fontsize=11, loc='left')
ax.grid(alpha=0.25, lw=0.5)
ax.legend(frameon=False, fontsize=9)
fig.tight_layout()
plt.show()
No description has been provided for this image

6. Does the hydrophone earn its place?¶

obs is the only model here that can use the fourth channel. The hydrophone is also the channel whose sampling rate varies most between these deployments — 100 Hz at AACSE, 40 Hz at Cascadia, 10 Hz at Blanco — while the model resamples everything to 100 Hz. An upsampled 10 Hz trace carries very little of the band a P onset lives in.

Re-running obs on the same windows with the hydrophone withheld separates the value of the channel from the value of the training.

In [9]:
if ABLATION_MODEL in models:
    rows = []
    for label, entries in examples.items():
      for st, tp, sta, mag, best in entries:
        with_h = st
        without_h = obspy.Stream([tr for tr in st if tr.stats.channel[-1] != 'H'])
        entry = {'experiment': label, 'station': sta, 'M': mag}
        for tag, stream in (('with_H', with_h), ('without_H', without_h)):
            out = models[ABLATION_MODEL].classify(stream, P_threshold=DETECT_FLOOR,
                                         S_threshold=DETECT_FLOOR)
            cand = [float(p.peak_value) for p in out.picks
                    if p.phase == 'P' and abs(p.peak_time - tp) <= TOLERANCE]
            entry[tag] = round(max(cand), 3) if cand else 0.0
        entry['hydrophone_rate'] = next(
            (tr.stats.sampling_rate for tr in st if tr.stats.channel[-1] == 'H'), np.nan)
        rows.append(entry)
    abl = pd.DataFrame(rows)
    if len(abl):
        abl['delta'] = (abl.with_H - abl.without_H).round(3)
        print(f'{len(abl)} detected windows re-run with the hydrophone withheld\n')
        print(abl.groupby(['experiment', 'hydrophone_rate'])[['with_H', 'without_H', 'delta']]
              .agg(['mean', 'count']).round(3).to_string())
        print(f'\nmean change in P confidence when the hydrophone is included: '
              f'{abl.delta.mean():+.4f}')
        worse = int((abl.delta < -0.01).sum()); better = int((abl.delta > 0.01).sum())
        print(f'windows where it helped by >0.01: {better}; hurt by >0.01: {worse}; '
              f'unchanged: {len(abl) - better - worse}')
    else:
        print('no detections captured for the ablation')

94 detected windows re-run with the hydrophone withheld

                              with_H       without_H        delta      
                                mean count      mean count   mean count
experiment    hydrophone_rate                                          
AACSE (XO)    100.0            0.807    28     0.866    28 -0.059    28
Blanco (X9)   10.0             0.855    24     0.876    24 -0.022    24
Cascadia (7D) 40.0             0.772    30     0.693    30  0.079    30
              125.0            0.699    12     0.714    12 -0.015    12

mean change in P confidence when the hydrophone is included: +0.0002
windows where it helped by >0.01: 32; hurt by >0.01: 35; unchanged: 27

7. Scanning the records¶

Ten windows per deployment, spread across its stations rather than drawn from one, vertical component, with every model's pick overlaid so disagreements are visible at a glance. The dotted line is the iasp91 prediction, which per section 5 carries the catalog's location error and should be read as an approximate marker rather than truth.

What to look for: picks clustering on a visible onset is the model working; picks spread across seconds of an emergent arrival is the hard case that ocean-bottom noise creates; and a pick with nothing visible under it is worth following up.

In [10]:
def gallery_figure(label, entries, window=(-25, 60)):
    n = len(entries)
    fig, axes = plt.subplots(n, 1, figsize=(11.5, 1.35 * n + 1.0), sharex=True)
    axes = np.atleast_1d(axes)
    for ax, e in zip(axes, entries):
        tr = e['stream'].select(component='Z')
        if not tr:
            continue
        tr = tr[0]
        t = tr.times(reftime=e['tp'])
        x = tr.data.astype(float)
        peak = np.abs(x).max()
        if peak > 0:
            x = x / peak
        ax.plot(t, x, color='#3d3d3d', lw=0.5)
        ax.axvline(0, color='#8a8a8a', lw=1.1, ls=':')
        for i, name in enumerate(names):
            if name not in e['picks']:
                continue
            dt, conf = e['picks'][name]
            ax.axvline(dt, color=COLORS[i % len(COLORS)], lw=1.5,
                       alpha=0.85 if conf >= REPORT_THRESHOLD else 0.35)
        found = len(e['picks'])
        ax.set_ylabel(f"{e['sta']}\nM{e['mag']:.1f} {e['deg']*111:.0f}km",
                      fontsize=7.5, rotation=0, ha='right', va='center', labelpad=32)
        ax.set_yticks([])
        ax.grid(alpha=0.18, lw=0.4, axis='x')
        ax.text(0.995, 0.82, f'{found}/{len(names)} models',
                transform=ax.transAxes, ha='right', fontsize=7, color='#7a7973')
    handles = [plt.Line2D([], [], color=COLORS[i % len(COLORS)], lw=1.6, label=n)
               for i, n in enumerate(names)]
    handles.append(plt.Line2D([], [], color='#8a8a8a', lw=1.1, ls=':',
                              label='iasp91 prediction'))
    axes[0].legend(handles=handles, fontsize=7.5, frameon=False, ncol=3,
                   loc='lower left', bbox_to_anchor=(0, 1.05))
    axes[-1].set_xlim(*window)
    axes[-1].set_xlabel('seconds from the predicted P')
    fig.suptitle(label, fontsize=11, x=0.01, ha='left', y=0.998)
    fig.tight_layout(rect=[0, 0, 1, 0.985])
    return fig


for label in EXPERIMENTS:
    entries = gallery.get(label) or []
    if not entries:
        continue
    gallery_figure(label, entries[:GALLERY_N])
    plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

One record in full¶

All four components for a single detected event per deployment, including the hydrophone, at its native sampling rate.

In [11]:
def plot_obs(label, entry, window=(-30, 90)):
    st, tp, sta, mag, best = entry
    order = ['Z', '1', '2', 'H']
    traces = [next((tr for tr in st if tr.stats.channel[-1] == c), None) for c in order]
    traces = [t for t in traces if t is not None]
    fig, axes = plt.subplots(len(traces), 1, figsize=(11, 1.5 * len(traces) + 1.1),
                             sharex=True)
    axes = np.atleast_1d(axes)
    for ax, tr in zip(axes, traces):
        ax.plot(tr.times(reftime=tp), tr.data, color='#3d3d3d', lw=0.5)
        ax.axvline(0, color='#8a8a8a', lw=1.2, ls=':')
        ax.axvline(best[0], color=C_P, lw=1.6)
        ax.set_ylabel(f'{tr.stats.channel}\n{tr.stats.sampling_rate:g} Hz', fontsize=8)
        ax.grid(alpha=0.22, lw=0.5)
        ax.tick_params(labelsize=8)
    axes[-1].set_xlim(*window)
    axes[-1].set_xlabel('seconds from the iasp91 predicted P '
                        '(dotted); solid line is the pick')
    axes[0].set_title(f'{label}  -  {sta}  -  M{mag}  -  '
                      f'pick {best[0]:+.1f} s at confidence {best[1]:.2f}',
                      fontsize=10, loc='left')
    fig.tight_layout()
    return fig


for label in EXPERIMENTS:
    entries = examples.get(label) or []
    if entries:
        plot_obs(label, max(entries, key=lambda e: e[4][1]))
        plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Reading the result¶

The comparison is relative, not absolute. Detection is scored against an iasp91 prediction with no water layer and no sediments, on a tolerance wide enough to absorb that. What the numbers support is which weight set does better on identical windows; what they do not support is a claim about how complete an OBS catalog would be.

The land models are being used out of domain on purpose. That is the measurement: QuakeScope will encounter ocean-bottom data, and the question is what running a land-trained picker over it costs. Renaming the horizontals from 1/2 to N/E is the usual convention and is harmless here, since none of these models uses absolute orientation.

Instrumentation varies more than the settings do. Hydrophone sampling runs from 100 Hz to 10 Hz across three deployments, and seismometer bands differ too. A weight set that depends on the fourth channel will behave differently between experiments for reasons that have nothing to do with the seismicity, which is what section 6 is for.

Thresholds still belong to the weight set. Everything above is reported at a shared cutoff, and the shared-threshold comparison was already shown to mislead in the land benchmark — see phasenet_sequence_comparison.ipynb. Note that obs ships different defaults from the rest, P 0.2 and S 0.1, which is itself a signal that its probabilities are not on the same scale.