-
Notifications
You must be signed in to change notification settings - Fork 0
/
observer_waveform.py
executable file
·120 lines (106 loc) · 5.45 KB
/
observer_waveform.py
1
2
3
4
5
6
7
8
9
10
11
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
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# -*- coding: utf-8 -*-
import os
import json
import time
import obspy
import numpy
import asyncio
import threading
import websockets
import matplotlib.pyplot
import matplotlib.animation
station_net = "AS" # Station network code
station_stn = "SHAKE" # Station station code
station_loc = "00" # Station location code
channel_prefix = "E" # Station channel prefix code E.g. EHx == E, BHx == B
station_address = "127.0.0.1:8073" # Observer address
time_span = 120 # Time span in seconds
refresh_time = 1000 # Refresh time in milliseconds
window_size = 2 # Spectrogram window size in seconds
overlap_percent = 86 # Spectrogram overlap in percent
spectrogram_power_range = [20, 120] # Spectrogram power range in dB
fig, axs = matplotlib.pyplot.subplots(6, 1, num="Observer Waveform", figsize=(9.6, 7.0))
matplotlib.pyplot.subplots_adjust(left=0, right=1, top=1, bottom=0, hspace=0, wspace=0)
def resample_trace(trace, target_sampling_rate):
if trace.stats.sampling_rate != target_sampling_rate:
trace.interpolate(target_sampling_rate)
return trace
def make_trace(channel, sps, counts_list, timestamp):
trace = obspy.core.Trace(data=numpy.ma.MaskedArray(counts_list, dtype=numpy.float64))
trace.stats.network = station_net
trace.stats.station = station_stn
trace.stats.location = station_loc
trace.stats.channel = channel
trace.stats.sampling_rate = sps
trace.stats.starttime = obspy.UTCDateTime(timestamp)
return trace
async def get_data():
global bhe_data, bhn_data, bhz_data
while True:
try:
async with websockets.connect(f"ws://{station_address}/api/v1/socket") as websocket:
while True:
data = json.loads(await websocket.recv())
timestamp = data["ts"] / 1000
bhe_data = make_trace(f"{channel_prefix}HE", len(data["ehe"]), data["ehe"], timestamp)
bhn_data = make_trace(f"{channel_prefix}HN", len(data["ehn"]), data["ehn"], timestamp)
bhz_data = make_trace(f"{channel_prefix}HZ", len(data["ehz"]), data["ehz"], timestamp)
except Exception as e:
print(e)
time.sleep(0.5)
continue
def update(frame):
try:
# Resample new data to match the stream sampling rate
bhe_resampled = resample_trace(bhe_data, bhe_stream.stats.sampling_rate)
bhn_resampled = resample_trace(bhn_data, bhn_stream.stats.sampling_rate)
bhz_resampled = resample_trace(bhz_data, bhz_stream.stats.sampling_rate)
# Update streams with fixed length
for stream, new_data in zip([bhe_stream, bhn_stream, bhz_stream], [bhe_resampled, bhn_resampled, bhz_resampled]):
new_samples = int(new_data.stats.npts)
stream_length = int(stream.stats.sampling_rate * time_span)
if len(stream.data) >= stream_length:
stream.data = numpy.roll(stream.data, -new_samples)
stream.data[-new_samples:] = new_data.data
else:
stream.data = numpy.concatenate((stream.data, new_data.data))
if len(stream.data) > stream_length:
stream.data = stream.data[-stream_length:]
stream.stats.starttime = stream.stats.starttime + 1.0
# Plot data
for i, (stream, component) in enumerate(zip([bhe_stream, bhn_stream, bhz_stream], [f"{channel_prefix}HE", f"{channel_prefix}HN", f"{channel_prefix}HZ"])):
axs[i*2].clear()
axs[i*2+1].clear()
times = numpy.arange(stream.stats.npts) / stream.stats.sampling_rate
waveform_data = stream.copy().filter("bandpass", freqmin=0.1, freqmax=10.0, zerophase=True).data
if not numpy.any(numpy.isnan(waveform_data)) and not numpy.any(numpy.isinf(waveform_data)):
axs[i*2].plot(times, waveform_data, label=component, color="blue")
axs[i*2].legend(loc="upper left")
axs[i*2].xaxis.set_visible(False)
axs[i*2].yaxis.set_visible(False)
axs[i*2].set_xlim([times[0], times[-1]])
axs[i*2].set_ylim([numpy.min(waveform_data), numpy.max(waveform_data)])
NFFT = int(stream.stats.sampling_rate * window_size)
noverlap = int(NFFT * (overlap_percent / 100))
spec_data = stream.copy().filter("highpass", freq=0.1, zerophase=True).data
if not numpy.any(numpy.isnan(spec_data)) and not numpy.any(numpy.isinf(spec_data)):
axs[i*2+1].specgram(spec_data, NFFT=NFFT, Fs=stream.stats.sampling_rate, noverlap=noverlap, cmap="jet", vmin=spectrogram_power_range[0], vmax=spectrogram_power_range[1])
axs[i*2+1].set_ylim(0, 15)
axs[i*2+1].yaxis.set_visible(False)
axs[i*2+1].xaxis.set_visible(False)
except Exception as e:
print(f"Error plotting data: {e}")
if __name__ == "__main__":
thread1 = threading.Thread(target=asyncio.run, args=(get_data(),))
thread1.start()
time.sleep(3)
bhe_stream = bhe_data.copy()
bhn_stream = bhn_data.copy()
bhz_stream = bhz_data.copy()
stream_length = int(bhe_stream.stats.sampling_rate * time_span)
bhe_stream.data = numpy.zeros(stream_length)
bhn_stream.data = numpy.zeros(stream_length)
bhz_stream.data = numpy.zeros(stream_length)
ani = matplotlib.animation.FuncAnimation(fig, update, interval=refresh_time, cache_frame_data=False)
matplotlib.pyplot.show()
os._exit(0)