Hello everyone,
I have 20 seismic nodes with channels DP1 (North), DP2 (East), and DPZ (Vertical) recording continuous data. I want to rotate the horizontal components (DP1, DP2) from NE to RT frame relative to a reference station, so that DP1→DPR and DP2→DPT.
The problem: after rotation, the comparison plot shows DPR matching DP2 and DPT matching DP1 — as if the components are swapped.
I verified in SAC that DP1 is indeed oriented North and the rotation there gives the correct result (DP1→DPR). So I suspect something is wrong in my ObsPy workflow. Has anyone encountered this? Is there something in ObsPy’s NE->RT convention I am missing?
Here are the steps I follow:
- Load and merge DP1 and DP2 for each station
- Sort by channel (gives
[DP1, DP2]order, confirmed) - Rename:
DP1→DPN,DP2→DPE - Call
stream.rotate(method="NE->RT", back_azimuth=baz) - Plot
data_before['DP1']vsDPR.data— they don’t match; insteaddata_before['DP2']matchesDPR.data
import obspy
import numpy as np
import matplotlib.pyplot as plt
from obspy.geodetics import gps2dist_azimuth
from collections import defaultdict
from pathlib import Path
base_path = Path("path/to/sac/files")
# --- load and group horizontal components by station ---
st_h_group = defaultdict(obspy.Stream)
for h_comp in ['1', '2']:
for f in base_path.glob(f"XO.*DP{h_comp}.*.sac"):
tr = obspy.read(str(f))[0]
st_h_group[tr.stats.station] += tr
# --- merge ---
for station, stream in st_h_group.items():
stream.merge(method=1, fill_value=0)
stream.sort(keys=['channel']) # ensures [DP1, DP2] order, confirmed by printing
# --- coordinates and BAZ ---
coords = {sta: (st[0].stats.sac.stla, st[0].stats.sac.stlo) for sta, st in st_h_group.items()}
ref_station = '101' # example
targ_station = '101' # self-pair: BAZ=0, so DPR should equal DP1 exactly
_, _, baz = gps2dist_azimuth(coords[ref_station][0], coords[ref_station][1],
coords[targ_station][0], coords[targ_station][1])
print(f"BAZ = {baz}") # confirmed 0.0
# --- rotate one window ---
h_stream = st_h_group[targ_station].copy()
h_stream.trim(starttime=h_stream[0].stats.starttime,
endtime=h_stream[0].stats.starttime + 20)
h_stream.detrend('demean')
h_stream.detrend('linear')
# capture data before rotation by index (not channel name) to avoid rename-across-windows bug
data_before_N = h_stream[0].data.copy() # DP1, North
data_before_E = h_stream[1].data.copy() # DP2, East
# rename DP1->DPN, DP2->DPE
for tr in h_stream:
if tr.stats.channel.endswith('1'):
tr.stats.channel = tr.stats.channel[:-1] + 'N'
elif tr.stats.channel.endswith('2'):
tr.stats.channel = tr.stats.channel[:-1] + 'E'
h_stream.rotate(method="NE->RT", back_azimuth=baz)
r_tr = h_stream.select(channel="*DPR")[0]
t_tr = h_stream.select(channel="*DPT")[0]
# --- plot: for BAZ=0, DPR should equal DP1 (North) ---
t = np.arange(h_stream[0].stats.npts) * h_stream[0].stats.delta
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
axes[0].plot(t, data_before_N, 'b', label='DP1/North (before)')
axes[0].plot(t, r_tr.data, 'r', label='DPR (after)')
axes[0].set_title(f'Expected: DPR == DP1 for BAZ=0')
axes[0].legend()
axes[1].plot(t, data_before_E, 'b', label='DP2/East (before)')
axes[1].plot(t, t_tr.data, 'r', label='DPT (after)')
axes[1].set_title(f'Expected: DPT == DP2 for BAZ=0')
axes[1].legend()
plt.tight_layout()
plt.show()
