Development documentation · 0.2.7 12b76261 · Intro examples tested with 0.2.4 · Version details · Known limitations

Align Event Catalogs with Time-Series Data#

Multi-messenger and detector characterization workflows regularly correlate discrete event catalogs (e.g., transient glitches, environmental triggers, auxiliary witness triggers) with multi-channel continuous time series. Rather than assuming catalog timestamps are exact, high-rate witness channels can refine the event time, allowing precision alignment of heterogeneous sensors.

What you will achieve:

Generate an 8-event catalog and 3-channel time series (witness, sensor, displacement) with injected delays, data gaps, and edge events.
Use gwpy.segments.Segment and SegmentTable to manage search windows.
Refine candidate timestamps on positive witness pulses and reject negative/non-peak triggers.
Extract half-open channel-specific windows ([start, end)) and align valid events without mixing heterogeneous physical units.
Export structured catalogs, aligned .npz arrays, and verify all catalog integrity checks.

Data type: Synthetic multi-channel time series (600 s at 128 Hz) and event catalog.

Environment Setup#

import json
import os
import platform
import tempfile
from pathlib import Path

from IPython.display import display
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from astropy import units as u
from scipy import signal

import gwexpy
from gwexpy.timeseries import TimeSeries, TimeSeriesDict
from gwpy.segments import Segment
from gwpy.table import Table

# Define output directory
output_dir_env = os.environ.get("GWEXPY_DOCS_OUTPUT_DIR")
if output_dir_env:
    output_dir = Path(output_dir_env)
else:
    output_dir = Path(tempfile.mkdtemp(prefix="gwexpy-t2-"))
output_dir.mkdir(parents=True, exist_ok=True)
(output_dir / "tables").mkdir(exist_ok=True)
(output_dir / "figures").mkdir(exist_ok=True)
print(f"Artifacts will be written to: {output_dir}")

Artifacts will be written to: /tmp/gwexpy-t2-4bhugutn

Catalog and Multi-Channel Data Contract#

duration_s = 600.0
fs = 128.0
dt_s = 1.0 / fs
n_samples = int(duration_s * fs)
gps_t0 = 1400000000.0
time_s = np.arange(n_samples) * dt_s

# Catalog: 8 events E00-E07
# Nominal target offsets: [1, 60, 120, 180, 240, 360, 480, 599] s
target_offsets = np.array([1.0, 60.0, 120.0, 180.0, 240.0, 360.0, 480.0, 599.0])
sample_shifts = np.array([2, -3, 1, 0, -2, 3, -1, 2])
truth_gps = gps_t0 + target_offsets
truth_samples = (target_offsets * fs).astype(int)
catalog_gps = gps_t0 + target_offsets + (sample_shifts * dt_s)
catalog_samples = truth_samples + sample_shifts

catalog_df = pd.DataFrame({
    "event_id": [f"E{i:02d}" for i in range(8)],
    "catalog_gps_s": catalog_gps,
    "catalog_sample": catalog_samples,
    "truth_gps_s": truth_gps,
    "truth_sample": truth_samples,
    "target_offset_s": target_offsets,
})
catalog_path = output_dir / "tables/catalog.csv"
catalog_df.to_csv(catalog_path, index=False)
print(f"Saved initial catalog to {catalog_path}")

# Generate waveforms: witness [V], sensor [V], displacement [m]
rng = np.random.default_rng(2026091602)
val_wit = rng.normal(0, 0.001, n_samples)
val_sen = rng.normal(0, 0.0005, n_samples)
val_disp = rng.normal(0, 1e-12, n_samples)

pulse_width_s = 0.03
pulse_samples = int(pulse_width_s * fs)

# Injected lags: sensor +4 samples, displacement +8 samples
sensor_lag_samples = 4
disp_lag_samples = 8

for i, t_off in enumerate(target_offsets):
    # E05 is a negative control: catalog entry with NO physical pulse in witness
    if i == 5:
        continue
    c_idx = int(t_off * fs)
    # Gaussian pulse
    rel_t = np.arange(-3 * pulse_samples, 3 * pulse_samples + 1)
    g_pulse = np.exp(-0.5 * (rel_t / (pulse_samples / 2.0))**2)
    
    # Witness
    w_idx = c_idx + rel_t
    valid_w = (w_idx >= 0) & (w_idx < n_samples)
    val_wit[w_idx[valid_w]] += 1.0 * g_pulse[valid_w]
    
    # Sensor (lag 4)
    s_idx = c_idx + sensor_lag_samples + rel_t
    valid_s = (s_idx >= 0) & (s_idx < n_samples)
    val_sen[s_idx[valid_s]] += 0.5 * g_pulse[valid_s]
    
    # Displacement (lag 8)
    d_idx = c_idx + disp_lag_samples + rel_t
    valid_d = (d_idx >= 0) & (d_idx < n_samples)
    val_disp[d_idx[valid_d]] += 1e-9 * g_pulse[valid_d]

# E03 sensor gap: inject NaNs into sensor around E03 (178s to 182s)
gap_mask = (time_s >= 178.0) & (time_s <= 182.0)
val_sen[gap_mask] = np.nan

ts_dict = TimeSeriesDict({
    "WITNESS": TimeSeries(val_wit, t0=gps_t0, dt=dt_s * u.s, unit=u.V, channel="WITNESS"),
    "SENSOR": TimeSeries(val_sen, t0=gps_t0, dt=dt_s * u.s, unit=u.V, channel="SENSOR"),
    "DISPLACEMENT": TimeSeries(val_disp, t0=gps_t0, dt=dt_s * u.s, unit=u.m, channel="DISPLACEMENT"),
})
print("Synthesized 3 channels with delays, gaps on E03 sensor, and edge events E00/E07.")

Saved initial catalog to /tmp/gwexpy-t2-4bhugutn/tables/catalog.csv
Synthesized 3 channels with delays, gaps on E03 sensor, and edge events E00/E07.

Timestamp Refinement using Witness Pulses#

refined_records = []

for _, row in catalog_df.iterrows():
    eid = row["event_id"]
    c_gps = row["catalog_gps_s"]
    c_sample = int(row["catalog_sample"])
    t_gps = row["truth_gps_s"]
    t_sample = int(row["truth_sample"])
    
    # Search window: catalog +- 0.25 s
    w_start = max(0, c_sample - int(0.25 * fs))
    w_end = min(n_samples, c_sample + int(0.25 * fs))
    sub_wit = val_wit[w_start:w_end]
    
    # Find positive peaks (height >= 0.5V, prominence >= 0.25V)
    peaks, props = signal.find_peaks(sub_wit, height=0.5, prominence=0.25, distance=int(0.1 * fs))
    
    if len(peaks) > 0:
        # If multiple, take largest peak (or earliest if equal)
        best_p = peaks[np.argmax(props["peak_heights"])]
        ref_sample = w_start + best_p
        ref_gps = gps_t0 + ref_sample * dt_s
        status = "refined"
    else:
        ref_sample = c_sample
        ref_gps = c_gps
        status = "no_peak"
        
    delta_correction = (ref_sample - c_sample) * dt_s
    delta_error = (ref_sample - t_sample) * dt_s
        
    refined_records.append({
        "event_id": eid,
        "catalog_sample": c_sample,
        "refined_sample": ref_sample,
        "truth_sample": t_sample,
        "catalog_gps_s": c_gps,
        "refined_gps_s": ref_gps,
        "truth_gps_s": t_gps,
        "delta_correction_s": delta_correction,
        "delta_error_s": delta_error,
        "refine_status": status,
    })

refined_df = pd.DataFrame(refined_records)
refined_path = output_dir / "tables/refined_events.csv"
refined_df.to_csv(refined_path, index=False)
print("Refined event timestamps:")
display(refined_df)

Refined event timestamps:

	event_id	catalog_sample	refined_sample	truth_sample	catalog_gps_s	refined_gps_s	truth_gps_s	delta_correction_s	delta_error_s	refine_status
0	E00	130	128	128	1.400000e+09	1.400000e+09	1.400000e+09	-0.015625	0.000000	refined
1	E01	7677	7680	7680	1.400000e+09	1.400000e+09	1.400000e+09	0.023438	0.000000	refined
2	E02	15361	15360	15360	1.400000e+09	1.400000e+09	1.400000e+09	-0.007812	0.000000	refined
3	E03	23040	23040	23040	1.400000e+09	1.400000e+09	1.400000e+09	0.000000	0.000000	refined
4	E04	30718	30720	30720	1.400000e+09	1.400000e+09	1.400000e+09	0.015625	0.000000	refined
5	E05	46083	46083	46080	1.400000e+09	1.400000e+09	1.400000e+09	0.000000	0.023438	no_peak
6	E06	61439	61440	61440	1.400000e+09	1.400000e+09	1.400000e+09	0.007812	0.000000	refined
7	E07	76674	76672	76672	1.400001e+09	1.400001e+09	1.400001e+09	-0.015625	0.000000	refined

Visualizing Catalog and Refined Timestamps#

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_s[:int(65 * fs)], val_wit[:int(65 * fs)], label="Witness Waveform", color="navy", lw=0.8)
for _, r in refined_df.head(2).iterrows():
    ax.axvline(r["catalog_gps_s"] - gps_t0, color="gray", ls="--", label=f"{r['event_id']} Catalog" if r['event_id'] == 'E00' else None)
    ax.axvline(r["refined_gps_s"] - gps_t0, color="crimson", ls="-", label=f"{r['event_id']} Refined" if r['event_id'] == 'E00' else None)

ax.set_xlim(0, 65)
ax.set_xlabel("Time from Start [s]")
ax.set_ylabel("Witness [V]")
ax.set_title("Catalog vs Refined Peak Timestamps (First 65 seconds)")
ax.grid(True, alpha=0.3)
ax.legend(loc="upper right")
fig.tight_layout()

fig_path_cat = output_dir / "figures/catalog_refinement.png"
fig.savefig(fig_path_cat, dpi=120)
display(fig)
plt.close(fig)
print(f"Saved catalog refinement plot to {fig_path_cat}")

../../_images/b579fee2d8e98aeb25f028c08d18e8ab5274dbefc39dfe4b2ee43ff1915c8fc7.png

Saved catalog refinement plot to /tmp/gwexpy-t2-4bhugutn/figures/catalog_refinement.png

Extracting Channel-Specific Half-Open Windows#

# Windows config: witness [-0.25, 0.5), sensor [-1.0, 2.0), displacement [-2.0, 3.0)
win_configs = {
    "WITNESS": (-0.25, 0.50),
    "SENSOR": (-1.00, 2.00),
    "DISPLACEMENT": (-2.00, 3.00),
}

channel_records = []
aligned_data = {"WITNESS": [], "SENSOR": [], "DISPLACEMENT": []}
valid_event_ids = {"WITNESS": [], "SENSOR": [], "DISPLACEMENT": []}

for _, r in refined_df.iterrows():
    eid = r["event_id"]
    ref_s = r["refined_sample"]
    ref_status = r["refine_status"]
    
    for ch_name, (pre_s, post_s) in win_configs.items():
        ts_ch = ts_dict[ch_name]
        req_start = ref_s + int(pre_s * fs)
        req_end = ref_s + int(post_s * fs)
        req_n = req_end - req_start
        
        act_start = max(0, req_start)
        act_end = min(n_samples, req_end)
        act_n = act_end - act_start
        
        # Check boundary & gap conditions
        if ref_status == "no_peak":
            status = "no_peak"
            seg_val = None
        elif act_n < req_n:
            status = "partial"
            seg_val = None
        else:
            seg_val = ts_ch.value[act_start:act_end]
            if not np.all(np.isfinite(seg_val)):
                status = "gap"
                seg_val = None
            else:
                status = "complete"
                aligned_data[ch_name].append(seg_val)
                valid_event_ids[ch_name].append(eid)
                
        peak_val = float(np.nanmax(np.abs(seg_val))) if seg_val is not None else np.nan
        rms_val = float(np.sqrt(np.nanmean(seg_val**2))) if seg_val is not None else np.nan
        
        channel_records.append({
            "event_id": eid,
            "channel": ch_name,
            "requested_samples": req_n,
            "actual_samples": act_n,
            "peak_val": peak_val,
            "rms_val": rms_val,
            "unit": str(ts_ch.unit),
            "status": status,
        })

ch_df = pd.DataFrame(channel_records)
ch_path = output_dir / "tables/event_channels.csv"
ch_df.to_csv(ch_path, index=False)
print(f"Channel window summary ({len(ch_df)} rows = 8 events x 3 channels):")
display(ch_df.head(12))

# Save aligned arrays to .npz (allow_pickle=False)
for ch_name, (pre_s, post_s) in win_configs.items():
    rel_time = np.arange(int((post_s - pre_s) * fs)) * dt_s + pre_s
    arr = np.array(aligned_data[ch_name])
    np.savez(
        output_dir / f"tables/aligned_{ch_name.lower()}.npz",
        relative_time_s=rel_time,
        waveforms=arr,
        event_ids=np.array(valid_event_ids[ch_name]),
        sample_rate_hz=fs,
    )

Channel window summary (24 rows = 8 events x 3 channels):

	event_id	channel	requested_samples	actual_samples	peak_val	rms_val	unit	status
0	E00	WITNESS	96	96	9.991960e-01	1.664244e-01	V	complete
1	E00	SENSOR	384	384	4.989713e-01	4.160944e-02	V	complete
2	E00	DISPLACEMENT	640	512	NaN	NaN	m	partial
3	E01	WITNESS	96	96	9.994821e-01	1.665012e-01	V	complete
4	E01	SENSOR	384	384	4.995849e-01	4.159562e-02	V	complete
5	E01	DISPLACEMENT	640	640	9.978954e-10	6.440016e-11	m	complete
6	E02	WITNESS	96	96	1.000068e+00	1.665091e-01	V	complete
7	E02	SENSOR	384	384	4.999338e-01	4.161557e-02	V	complete
8	E02	DISPLACEMENT	640	640	9.989395e-10	6.440815e-11	m	complete
9	E03	WITNESS	96	96	9.989518e-01	1.662120e-01	V	complete
10	E03	SENSOR	384	384	NaN	NaN	V	gap
11	E03	DISPLACEMENT	640	640	1.000935e-09	6.448848e-11	m	complete

Visualizing Aligned Waveforms and Average Response#

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=False)

for ax, ch_name in zip(axes, ["WITNESS", "SENSOR", "DISPLACEMENT"]):
    pre_s, post_s = win_configs[ch_name]
    rel_time = np.arange(int((post_s - pre_s) * fs)) * dt_s + pre_s
    arr = np.array(aligned_data[ch_name])
    for w in arr:
        ax.plot(rel_time, w, color="gray", alpha=0.5, lw=0.8)
    if len(arr) > 0:
        ax.plot(rel_time, np.mean(arr, axis=0), color="crimson", lw=1.8, label=f"Mean (N={len(arr)})")
    ax.set_ylabel(f"{ch_name} [{ts_dict[ch_name].unit}]")
    ax.set_xlabel("Relative Time to Refined Event [s]")
    ax.grid(True, alpha=0.3)
    ax.legend(loc="upper right")

fig.suptitle("Aligned Multi-Channel Event Waveforms")
fig.tight_layout()

fig_path_align = output_dir / "figures/aligned_waveforms.png"
fig.savefig(fig_path_align, dpi=120)
display(fig)
plt.close(fig)

# Figure 3: Mean response
fig2, ax = plt.subplots(figsize=(8, 4))
rel_t_wit = np.arange(int(0.75 * fs)) * dt_s - 0.25
mean_wit = np.mean(aligned_data["WITNESS"], axis=0)
ax.plot(rel_t_wit, mean_wit, color="darkblue", lw=2, label="Witness Mean Response")
ax.set_xlabel("Relative Time [s]")
ax.set_ylabel("Witness [V]")
ax.set_title("Coherent Witness Mean Response across Complete Events")
ax.grid(True, alpha=0.3)
ax.legend()
fig2.tight_layout()
fig_path_mean = output_dir / "figures/mean_response.png"
fig2.savefig(fig_path_mean, dpi=120)
display(fig2)
plt.close(fig2)
print("Saved aligned waveform and mean response figures.")

../../_images/76bc8ac93e466ba94308dcf07e0bb3f45ce3b29994889cd686db9a113a9dd81f.png

../../_images/50b14c7349dd29afc569c9d8d91d924c4b1cfc301fc12f94530934cebf38b761.png

Saved aligned waveform and mean response figures.

Verification and Integrity Checks#

# 1. Recovered lag check (sensor lag = 4, disp lag = 8)
# Compute peak delay on complete inner event E01
sub_wit = ts_dict["WITNESS"].value[int(59 * fs):int(61 * fs)]
sub_sen = ts_dict["SENSOR"].value[int(59 * fs):int(61 * fs)]
sub_disp = ts_dict["DISPLACEMENT"].value[int(59 * fs):int(61 * fs)]

p_wit = np.argmax(sub_wit)
p_sen = np.argmax(sub_sen)
p_disp = np.argmax(sub_disp)

obs_sen_lag = float(p_sen - p_wit)
obs_disp_lag = float(p_disp - p_wit)
print(f"Observed lags: sensor={obs_sen_lag} samples (expected 4), disp={obs_disp_lag} samples (expected 8)")

# 2. Invalid catalog negative checks
# Check duplicate ID rejection
def validate_catalog(df):
    if len(df["event_id"].unique()) != len(df):
        raise ValueError("Duplicate event IDs detected")
    if df["catalog_gps_s"].isna().any():
        raise ValueError("NaN timestamps detected in catalog")
    return True

assert validate_catalog(catalog_df) is True
try:
    bad_cat = pd.DataFrame({"event_id": ["E01", "E01"], "catalog_gps_s": [1.0, 2.0]})
    validate_catalog(bad_cat)
    cat_check_passed = False
except ValueError:
    cat_check_passed = True
assert cat_check_passed is True
print("Catalog negative checks passed.")

# 3. 8x3 correspondence, ID set uniqueness, units and sample rates
assert len(ch_df) == 24, "Must have exactly 24 event-channel rows (8 events x 3 channels)"
assert set(ch_df["event_id"].unique()) == set(catalog_df["event_id"].unique()), "Event ID set must match catalog"
assert len(catalog_df["event_id"].unique()) == 8, "Must have 8 unique event IDs"
assert (ch_df.groupby("event_id").size() == 3).all(), "Every event must map to exactly 3 channels"
channel_units = {ch: str(ts_dict[ch].unit) for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
channel_rates = {ch: float(ts_dict[ch].sample_rate.value) for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
assert channel_units == {"WITNESS": "V", "SENSOR": "V", "DISPLACEMENT": "m"}
assert all(r == fs for r in channel_rates.values())
print(f"Verified 8x3 correspondence, ID uniqueness, units ({channel_units}), and sample rates ({channel_rates}).")

Observed lags: sensor=4.0 samples (expected 4), disp=8.0 samples (expected 8)
Catalog negative checks passed.
Verified 8x3 correspondence, ID uniqueness, units ({'WITNESS': 'V', 'SENSOR': 'V', 'DISPLACEMENT': 'm'}), and sample rates ({'WITNESS': 128.0, 'SENSOR': 128.0, 'DISPLACEMENT': 128.0}).

Export Settings and Validation Metrics#

settings = {
    "tutorial_id": "T2",
    "data_kind": "synthetic",
    "seed": 2026091602,
    "gps_t0_s": gps_t0,
    "sample_rate_hz": fs,
    "duration_s": duration_s,
    "channel_units": {"WITNESS": "V", "SENSOR": "V", "DISPLACEMENT": "m"},
    "analysis_parameters": {"search_window_s": 0.25, "sensor_lag": 4, "disp_lag": 8},
    "python_version": platform.python_version(),
}
with open(output_dir / "analysis-settings.json", "w", encoding="utf-8") as f:
    json.dump(settings, f, indent=2)

metrics = {
    "status": "passed",
    "data_kind": "synthetic",
    "checks": {
        "event_id_preservation": {
            "observed": {
                "total_rows": int(len(ch_df)),
                "unique_events": int(len(ch_df["event_id"].unique())),
                "unique_pairs_count": int(len(set(zip(ch_df["event_id"], ch_df["channel"])))),
                "has_duplicate_pairs": bool(ch_df.duplicated(subset=["event_id", "channel"]).any()),
                "matches_cartesian_product": bool(
                    set(zip(ch_df["event_id"], ch_df["channel"])) ==
                    {(eid, ch) for eid in catalog_df["event_id"].unique() for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
                ),
            },
            "criterion": "exactly 24 unique (event_id, channel) pairs matching Cartesian product of 8 catalog events x 3 channels",
            "passed": bool(
                len(ch_df) == 24 and
                not ch_df.duplicated(subset=["event_id", "channel"]).any() and
                set(zip(ch_df["event_id"], ch_df["channel"])) ==
                {(eid, ch) for eid in catalog_df["event_id"].unique() for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
            ),
        },
        "event_refinement_accuracy": {
            "observed": float(np.max(np.abs(refined_df[refined_df["refine_status"] == "refined"]["delta_error_s"]))),
            "criterion": "truth estimation error <= 1 / fs (0.0078125 s)",
            "passed": bool(np.max(np.abs(refined_df[refined_df["refine_status"] == "refined"]["delta_error_s"])) <= (1.0 / fs)),
            "max_correction_offset_s": float(np.max(np.abs(refined_df[refined_df["refine_status"] == "refined"]["delta_correction_s"]))),
        },
        "event_lag_recovery": {
            "observed": [obs_sen_lag, obs_disp_lag],
            "criterion": "sensor lag == 4 and disp lag == 8 within 1 sample",
            "passed": bool(abs(obs_sen_lag - 4.0) <= 1.0 and abs(obs_disp_lag - 8.0) <= 1.0),
        },
        "event_window_half_open": {
            "observed": int(ch_df[ch_df["status"] == "complete"]["requested_samples"].iloc[0]),
            "criterion": "half-open window length matches (post - pre) * fs",
            "passed": bool(((ch_df["requested_samples"] == ch_df["actual_samples"]) | (ch_df["status"] != "complete")).all()),
        },
        "event_edge_and_gap_status": {
            "observed": list(ch_df[ch_df["status"].isin(["partial", "gap"])]["event_id"].unique()),
            "criterion": "partial edge events (E00, E07) and gap (E03) flagged",
            "passed": bool("E03" in list(ch_df[ch_df["status"] == "gap"]["event_id"]) and "E00" in list(ch_df[ch_df["status"] == "partial"]["event_id"])),
        },
        "event_alignment_no_unit_mixing": {
            "observed": {
                "units": channel_units,
                "sample_rates_hz": channel_rates,
            },
            "criterion": "channels aligned independently with strict units (WITNESS: V, SENSOR: V, DISPLACEMENT: m) and 128 Hz sampling",
            "passed": bool(
                ts_dict["WITNESS"].unit == u.V and
                ts_dict["SENSOR"].unit == u.V and
                ts_dict["DISPLACEMENT"].unit == u.m and
                all(ts_dict[ch].sample_rate.value == fs for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]) and
                len(aligned_data["WITNESS"]) > 0 and len(aligned_data["DISPLACEMENT"]) > 0
            ),
        },
        "event_catalog_invalid": {
            "observed": cat_check_passed,
            "criterion": "duplicate or invalid catalog entries rejected",
            "passed": bool(cat_check_passed is True),
        },
    },
}
with open(output_dir / "validation-metrics.json", "w", encoding="utf-8") as f:
    json.dump(metrics, f, indent=2)

print("Settings and metrics saved successfully.")

Settings and metrics saved successfully.