crimeacs commited on
Commit
05d0cc5
1 Parent(s): 871cd60

WIP on section plot

Browse files
.DS_Store CHANGED
Binary files a/.DS_Store and b/.DS_Store differ
 
Gradio_app.ipynb ADDED
@@ -0,0 +1,506 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ {
2
+ "cells": [
3
+ {
4
+ "cell_type": "code",
5
+ "execution_count": 8,
6
+ "metadata": {},
7
+ "outputs": [
8
+ {
9
+ "name": "stdout",
10
+ "output_type": "stream",
11
+ "text": [
12
+ "Inventory created at 2023-03-31T20:48:41.380800Z\n",
13
+ "\tCreated by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.52\n",
14
+ "\t\t http://service.iris.edu/fdsnws/station/1/query?starttime=2019-07-...\n",
15
+ "\tSending institution: IRIS-DMC (IRIS-DMC)\n",
16
+ "\tContains:\n",
17
+ "\t\tNetworks (7):\n",
18
+ "\t\t\t8P, CI, LB, NN, NP, PB, SY\n",
19
+ "\t\tStations (85):\n",
20
+ "\t\t\t8P.CAU08 (Monache Meadows, CA, USA)\n",
21
+ "\t\t\tCI.APL (Apollo)\n",
22
+ "\t\t\tCI.CCA (California City Airport)\n",
23
+ "\t\t\tCI.CCC (Christmas Canyon China Lake)\n",
24
+ "\t\t\tCI.CGO (Cerro Gordo)\n",
25
+ "\t\t\tCI.CLC (China Lake)\n",
26
+ "\t\t\tCI.CWC (Cottonwood Creek)\n",
27
+ "\t\t\tCI.DAW (Darwin)\n",
28
+ "\t\t\tCI.DTP (Desert Tortoise Park)\n",
29
+ "\t\t\tCI.GSC (Goldstone)\n",
30
+ "\t\t\tCI.HAR (Harper Dry Lake bed)\n",
31
+ "\t\t\tCI.ISA (Isabella)\n",
32
+ "\t\t\tCI.JRC2 (Joshua Ridge: China Lake)\n",
33
+ "\t\t\tCI.LMR2 (Leuhmann Ridge Extension)\n",
34
+ "\t\t\tCI.LRL (Laurel Mtn)\n",
35
+ "\t\t\tCI.MPM (Manuel Prospect Mine)\n",
36
+ "\t\t\tCI.MRS (Mars)\n",
37
+ "\t\t\tCI.Q0068 (Redwood Blvd, California City CA)\n",
38
+ "\t\t\tCI.Q0072 (Lakeland Street, Ridgecrest CA)\n",
39
+ "\t\t\tCI.SLA (Slate Mountain)\n",
40
+ "\t\t\tCI.SRT (Snort)\n",
41
+ "\t\t\tCI.TEH (Cattani Ranch)\n",
42
+ "\t\t\tCI.TOW2 (Tower 2)\n",
43
+ "\t\t\tCI.WAS2 (Alta Sierra 2)\n",
44
+ "\t\t\tCI.WBM (Bowman Road)\n",
45
+ "\t\t\tCI.WBS (Bird Springs)\n",
46
+ "\t\t\tCI.WCS2 (Coso Hot Springs 2)\n",
47
+ "\t\t\tCI.WHF (Hanning Flat)\n",
48
+ "\t\t\tCI.WLH2 (Little Horse 2)\n",
49
+ "\t\t\tCI.WMF (Mccloud Flat)\n",
50
+ "\t\t\tCI.WNM (Nine Mile Canyon)\n",
51
+ "\t\t\tCI.WOR (Onyx Ranch)\n",
52
+ "\t\t\tCI.WRC2 (Renegade Canyon)\n",
53
+ "\t\t\tCI.WRV2 (Rose Valley Canyon 2)\n",
54
+ "\t\t\tCI.WVP2 (Volcano Peak 2)\n",
55
+ "\t\t\tLB.DAC (Inyo County, Darwin, CA, USA)\n",
56
+ "\t\t\tNN.GWY (Greenwater Valley, CA. (GPS 12/06/2000) w84gm)\n",
57
+ "\t\t\tNN.PAN (Panamint Range. (GPS 12/06/2000) w84gm)\n",
58
+ "\t\t\tNN.QSM (Queen of Sheba Mine, CA. (GPS 01/17/2001) w84gm)\n",
59
+ "\t\t\tNP.1035 (CA:Lake Isabella Dam)\n",
60
+ "\t\t\tNP.1809 (CA:Haiwee Rsvr;Pump Pl)\n",
61
+ "\t\t\tNP.5419 (CA:China Lake;Nav Weapon Ctr)\n",
62
+ "\t\t\tPB.B916 (marips916bcs2008, China Lake, CA, USA)\n",
63
+ "\t\t\tPB.B917 (tonyso917bcs2008, China Lake, CA, USA)\n",
64
+ "\t\t\tPB.B918 (mtsprn918bcs2008, China Lake, CA, USA)\n",
65
+ "\t\t\tPB.B921 (randsb921bcs2008, China Lake, CA, USA)\n",
66
+ "\t\t\tSY.CCA (CCA synthetic)\n",
67
+ "\t\t\tSY.CCC (CCC synthetic)\n",
68
+ "\t\t\tSY.CGO (CGO synthetic)\n",
69
+ "\t\t\tSY.CLC (CLC synthetic)\n",
70
+ "\t\t\tSY.CWC (CWC synthetic)\n",
71
+ "\t\t\tSY.DAC (DAC synthetic)\n",
72
+ "\t\t\tSY.DAW (DAW synthetic)\n",
73
+ "\t\t\tSY.DTP (DTP synthetic)\n",
74
+ "\t\t\tSY.FPC (FPC synthetic)\n",
75
+ "\t\t\tSY.FSR (FSR synthetic)\n",
76
+ "\t\t\tSY.GPO (GPO synthetic)\n",
77
+ "\t\t\tSY.GSC (GSC synthetic)\n",
78
+ "\t\t\tSY.HAR (HAR synthetic)\n",
79
+ "\t\t\tSY.ISA (ISA synthetic)\n",
80
+ "\t\t\tSY.JRC (JRC synthetic)\n",
81
+ "\t\t\tSY.JRC2 (JRC2 synthetic)\n",
82
+ "\t\t\tSY.KRV3 (KRV3 synthetic)\n",
83
+ "\t\t\tSY.LMR (LMR synthetic)\n",
84
+ "\t\t\tSY.LMR2 (LMR2 synthetic)\n",
85
+ "\t\t\tSY.LRL (LRL synthetic)\n",
86
+ "\t\t\tSY.MPM (MPM synthetic)\n",
87
+ "\t\t\tSY.OVRO (OVRO synthetic)\n",
88
+ "\t\t\tSY.RRC (RRC synthetic)\n",
89
+ "\t\t\tSY.SEV (SEV synthetic)\n",
90
+ "\t\t\tSY.SLA (SLA synthetic)\n",
91
+ "\t\t\tSY.SRT (SRT synthetic)\n",
92
+ "\t\t\tSY.TEH (TEH synthetic)\n",
93
+ "\t\t\tSY.TOW2 (TOW2 synthetic)\n",
94
+ "\t\t\tSY.WAS2 (WAS2 synthetic)\n",
95
+ "\t\t\tSY.WBM (WBM synthetic)\n",
96
+ "\t\t\tSY.WBP (WBP synthetic)\n",
97
+ "\t\t\tSY.WBS (WBS synthetic)\n",
98
+ "\t\t\tSY.WCS2 (WCS2 synthetic)\n",
99
+ "\t\t\tSY.WHF (WHF synthetic)\n",
100
+ "\t\t\tSY.WLH2 (WLH2 synthetic)\n",
101
+ "\t\t\tSY.WMF (WMF synthetic)\n",
102
+ "\t\t\tSY.WNM (WNM synthetic)\n",
103
+ "\t\t\tSY.WOR (WOR synthetic)\n",
104
+ "\t\t\tSY.WRC2 (WRC2 synthetic)\n",
105
+ "\t\tChannels (0):\n",
106
+ "\n"
107
+ ]
108
+ }
109
+ ],
110
+ "source": [
111
+ "import obspy\n",
112
+ "from obspy.clients.fdsn import Client\n",
113
+ "\n",
114
+ "client_name = 'SCEDC'\n",
115
+ "radius_km = 100\n",
116
+ "timestamp = '2019-07-04 17:33:49'\n",
117
+ "eq_lat = 35.766\n",
118
+ "eq_lon = -117.605\n",
119
+ "\n",
120
+ "origin_time = obspy.UTCDateTime(timestamp)\n",
121
+ "\n",
122
+ "client = Client(\"IRIS\")\n",
123
+ "inventory = client.get_stations(network=\"*\", station=\"*\", channel=\"*\",\n",
124
+ " starttime=origin_time, endtime=origin_time + 120,\n",
125
+ " latitude=eq_lat, longitude=eq_lon, maxradius=radius_km/111.2)\n",
126
+ "print(inventory)\n"
127
+ ]
128
+ },
129
+ {
130
+ "cell_type": "code",
131
+ "execution_count": null,
132
+ "metadata": {},
133
+ "outputs": [],
134
+ "source": [
135
+ "\n",
136
+ "\n",
137
+ "client = Client(client_name)\n",
138
+ "window = radius_km / 111.2\n",
139
+ "\n",
140
+ "assert eq_lat - window > -90 and eq_lat + window < 90, \"Latitude out of bounds\"\n",
141
+ "assert eq_lon - window > -180 and eq_lon + window < 180, \"Longitude out of bounds\"\n",
142
+ "\n",
143
+ "starttime = obspy.UTCDateTime(timestamp)\n",
144
+ "endtime = starttime + 120\n",
145
+ "\n",
146
+ "inv = client.get_stations(network=\"*\", station=\"*\", location=\"*\", channel=\"*H*\", \n",
147
+ " starttime=starttime, endtime=endtime, \n",
148
+ " minlatitude=(eq_lat-window), maxlatitude=(eq_lat+window),\n",
149
+ " minlongitude=(eq_lon-window), maxlongitude=(eq_lon+window), \n",
150
+ " level='station')"
151
+ ]
152
+ },
153
+ {
154
+ "cell_type": "code",
155
+ "execution_count": 64,
156
+ "metadata": {},
157
+ "outputs": [
158
+ {
159
+ "name": "stderr",
160
+ "output_type": "stream",
161
+ "text": [
162
+ "/Users/anovosel/miniconda3/envs/phasehunter/lib/python3.11/site-packages/gradio/outputs.py:43: UserWarning: Usage of gradio.outputs is deprecated, and will not be supported in the future, please import your components from gradio.components\n",
163
+ " warnings.warn(\n"
164
+ ]
165
+ },
166
+ {
167
+ "name": "stdout",
168
+ "output_type": "stream",
169
+ "text": [
170
+ "Running on local URL: http://127.0.0.1:7914\n",
171
+ "\n",
172
+ "To create a public link, set `share=True` in `launch()`.\n"
173
+ ]
174
+ },
175
+ {
176
+ "data": {
177
+ "text/html": [
178
+ "<div><iframe src=\"http://127.0.0.1:7914/\" width=\"100%\" height=\"500\" allow=\"autoplay; camera; microphone; clipboard-read; clipboard-write;\" frameborder=\"0\" allowfullscreen></iframe></div>"
179
+ ],
180
+ "text/plain": [
181
+ "<IPython.core.display.HTML object>"
182
+ ]
183
+ },
184
+ "metadata": {},
185
+ "output_type": "display_data"
186
+ },
187
+ {
188
+ "data": {
189
+ "text/plain": []
190
+ },
191
+ "execution_count": 64,
192
+ "metadata": {},
193
+ "output_type": "execute_result"
194
+ },
195
+ {
196
+ "name": "stdout",
197
+ "output_type": "stream",
198
+ "text": [
199
+ "torch.Size([256])\n"
200
+ ]
201
+ }
202
+ ],
203
+ "source": [
204
+ "# Gradio app that takes seismic waveform as input and marks 2 phases on the waveform as output.\n",
205
+ "\n",
206
+ "import gradio as gr\n",
207
+ "import numpy as np\n",
208
+ "from phasehunter.model import Onset_picker, Updated_onset_picker\n",
209
+ "from phasehunter.data_preparation import prepare_waveform\n",
210
+ "import torch\n",
211
+ "\n",
212
+ "from scipy.stats import gaussian_kde\n",
213
+ "\n",
214
+ "import obspy\n",
215
+ "from obspy.clients.fdsn import Client\n",
216
+ "from obspy.clients.fdsn.header import FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException\n",
217
+ "from obspy.geodetics.base import locations2degrees\n",
218
+ "from obspy.taup import TauPyModel\n",
219
+ "from obspy.taup.helper_classes import SlownessModelError\n",
220
+ "\n",
221
+ "from obspy.clients.fdsn.header import URL_MAPPINGS\n",
222
+ "\n",
223
+ "import matplotlib.pyplot as plt\n",
224
+ "\n",
225
+ "def make_prediction(waveform):\n",
226
+ " waveform = np.load(waveform)\n",
227
+ " processed_input = prepare_waveform(waveform)\n",
228
+ " \n",
229
+ " # Make prediction\n",
230
+ " with torch.no_grad():\n",
231
+ " output = model(processed_input)\n",
232
+ "\n",
233
+ " p_phase = output[:, 0]\n",
234
+ " s_phase = output[:, 1]\n",
235
+ "\n",
236
+ " return processed_input, p_phase, s_phase\n",
237
+ "\n",
238
+ "def mark_phases(waveform):\n",
239
+ " processed_input, p_phase, s_phase = make_prediction(waveform)\n",
240
+ "\n",
241
+ " # Create a plot of the waveform with the phases marked\n",
242
+ " if sum(processed_input[0][2] == 0): #if input is 1C\n",
243
+ " fig, ax = plt.subplots(nrows=2, figsize=(10, 2), sharex=True)\n",
244
+ "\n",
245
+ " ax[0].plot(processed_input[0][0])\n",
246
+ " ax[0].set_ylabel('Norm. Ampl.')\n",
247
+ "\n",
248
+ " else: #if input is 3C\n",
249
+ " fig, ax = plt.subplots(nrows=4, figsize=(10, 6), sharex=True)\n",
250
+ " ax[0].plot(processed_input[0][0])\n",
251
+ " ax[1].plot(processed_input[0][1])\n",
252
+ " ax[2].plot(processed_input[0][2])\n",
253
+ "\n",
254
+ " ax[0].set_ylabel('Z')\n",
255
+ " ax[1].set_ylabel('N')\n",
256
+ " ax[2].set_ylabel('E')\n",
257
+ "\n",
258
+ " p_phase_plot = p_phase*processed_input.shape[-1]\n",
259
+ " p_kde = gaussian_kde(p_phase_plot)\n",
260
+ " p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )\n",
261
+ " ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')\n",
262
+ "\n",
263
+ " s_phase_plot = s_phase*processed_input.shape[-1]\n",
264
+ " s_kde = gaussian_kde(s_phase_plot)\n",
265
+ " s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )\n",
266
+ " ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')\n",
267
+ "\n",
268
+ " for a in ax:\n",
269
+ " a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P')\n",
270
+ " a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S')\n",
271
+ "\n",
272
+ " ax[-1].set_xlabel('Time, samples')\n",
273
+ " ax[-1].set_ylabel('Uncert.')\n",
274
+ " ax[-1].legend()\n",
275
+ "\n",
276
+ " plt.subplots_adjust(hspace=0., wspace=0.)\n",
277
+ "\n",
278
+ " # Convert the plot to an image and return it\n",
279
+ " fig.canvas.draw()\n",
280
+ " image = np.array(fig.canvas.renderer.buffer_rgba())\n",
281
+ " plt.close(fig)\n",
282
+ " return image\n",
283
+ "\n",
284
+ "def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source_depth_km, velocity_model):\n",
285
+ " distances, t0s, st_lats, st_lons, waveforms = [], [], [], [], []\n",
286
+ " \n",
287
+ " taup_model = TauPyModel(model=velocity_model)\n",
288
+ " client = Client(client_name)\n",
289
+ "\n",
290
+ " window = radius_km / 111.2\n",
291
+ "\n",
292
+ " assert eq_lat - window > -90 and eq_lat + window < 90, \"Latitude out of bounds\"\n",
293
+ " assert eq_lon - window > -180 and eq_lon + window < 180, \"Longitude out of bounds\"\n",
294
+ "\n",
295
+ " starttime = obspy.UTCDateTime(timestamp)\n",
296
+ " endtime = starttime + 120\n",
297
+ "\n",
298
+ " inv = client.get_stations(network=\"*\", station=\"*\", location=\"*\", channel=\"*H*\", \n",
299
+ " starttime=starttime, endtime=endtime, \n",
300
+ " minlatitude=(eq_lat-window), maxlatitude=(eq_lat+window),\n",
301
+ " minlongitude=(eq_lon-window), maxlongitude=(eq_lon+window), \n",
302
+ " level='station')\n",
303
+ " \n",
304
+ " waveforms = []\n",
305
+ " for network in inv:\n",
306
+ " for station in network:\n",
307
+ " try:\n",
308
+ " distance = locations2degrees(eq_lat, eq_lon, station.latitude, station.longitude)\n",
309
+ "\n",
310
+ " arrivals = taup_model.get_travel_times(source_depth_in_km=source_depth_km, \n",
311
+ " distance_in_degree=distance, \n",
312
+ " phase_list=[\"P\", \"S\"])\n",
313
+ "\n",
314
+ " if len(arrivals) > 0:\n",
315
+ "\n",
316
+ " starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15\n",
317
+ " endtime = starttime + 60\n",
318
+ "\n",
319
+ " waveform = client.get_waveforms(network=network.code, station=station.code, location=\"*\", channel=\"*\", \n",
320
+ " starttime=starttime, endtime=endtime)\n",
321
+ " \n",
322
+ " waveform = waveform.select(channel=\"H[BH][ZNE]\")\n",
323
+ " waveform = waveform.merge(fill_value=0)\n",
324
+ " waveform = waveform[:3]\n",
325
+ " \n",
326
+ " len_check = [len(x.data) for x in waveform]\n",
327
+ " if len(set(len_check)) > 1:\n",
328
+ " continue\n",
329
+ "\n",
330
+ " if len(waveform) == 3:\n",
331
+ " waveform = prepare_waveform(np.stack([x.data for x in waveform]))\n",
332
+ " \n",
333
+ " distances.append(distance)\n",
334
+ " t0s.append(starttime)\n",
335
+ " st_lats.append(station.latitude)\n",
336
+ " st_lons.append(station.longitude)\n",
337
+ " waveforms.append(waveform)\n",
338
+ "\n",
339
+ " except (IndexError, FDSNNoDataException, FDSNTimeoutException):\n",
340
+ " continue\n",
341
+ "\n",
342
+ " with torch.no_grad():\n",
343
+ " waveforms_torch = torch.vstack(waveforms)\n",
344
+ " output = model(waveforms_torch)\n",
345
+ "\n",
346
+ " p_phases = output[:, 0]\n",
347
+ " s_phases = output[:, 1]\n",
348
+ "\n",
349
+ "\n",
350
+ " print(p_phases.shape)\n",
351
+ " # for i in range(len(waveforms)):\n",
352
+ " # current_P = P_batch[i::len(waveforms)].cpu()\n",
353
+ " # current_S_batch = S_batch[i::len(waveforms)].cpu()\n",
354
+ " # current_Pg_batch = Pg_batch[i::len(waveforms)].cpu()\n",
355
+ " # current_Sg_batch = Sg_batch[i::len(waveforms)].cpu()\n",
356
+ " # current_Pn_batch = Pn_batch[i::len(waveforms)].cpu()\n",
357
+ " # current_Sn_batch = Sn_batch[i::len(waveforms)].cpu()\n",
358
+ " \n",
359
+ " fig,ax = plt.subplots()\n",
360
+ " ax.scatter(st_lats, st_lons)\n",
361
+ " fig.canvas.draw()\n",
362
+ " image = np.array(fig.canvas.renderer.buffer_rgba())\n",
363
+ " plt.close(fig)\n",
364
+ "\n",
365
+ " return image\n",
366
+ "\n",
367
+ "\n",
368
+ "model = Onset_picker.load_from_checkpoint(\"./weights.ckpt\",\n",
369
+ " picker=Updated_onset_picker(),\n",
370
+ " learning_rate=3e-4)\n",
371
+ "model.eval()\n",
372
+ "\n",
373
+ "\n",
374
+ "\n",
375
+ "# # Create the Gradio interface\n",
376
+ "# gr.Interface(mark_phases, inputs, outputs, title='PhaseHunter').launch()\n",
377
+ "\n",
378
+ "\n",
379
+ "with gr.Blocks() as demo:\n",
380
+ " gr.Markdown(\"# PhaseHunter\")\n",
381
+ " gr.Markdown(\"\"\"This app allows one to detect P and S seismic phases along with uncertainty of the detection. \n",
382
+ " The app can be used in three ways: either by selecting one of the sample waveforms;\n",
383
+ " or by selecting an earthquake from the global earthquake catalogue;\n",
384
+ " or by uploading a waveform of interest.\n",
385
+ " \"\"\")\n",
386
+ " with gr.Tab(\"Default example\"):\n",
387
+ " # Define the input and output types for Gradio\n",
388
+ " inputs = gr.Dropdown(\n",
389
+ " [\"data/sample/sample_0.npy\", \n",
390
+ " \"data/sample/sample_1.npy\", \n",
391
+ " \"data/sample/sample_2.npy\"], \n",
392
+ " label=\"Sample waveform\", \n",
393
+ " info=\"Select one of the samples\",\n",
394
+ " value = \"data/sample/sample_0.npy\"\n",
395
+ " )\n",
396
+ "\n",
397
+ " button = gr.Button(\"Predict phases\")\n",
398
+ " outputs = gr.outputs.Image(label='Waveform with Phases Marked', type='numpy')\n",
399
+ " \n",
400
+ " button.click(mark_phases, inputs=inputs, outputs=outputs)\n",
401
+ " \n",
402
+ " with gr.Tab(\"Select earthquake from catalogue\"):\n",
403
+ " gr.Markdown('TEST')\n",
404
+ " \n",
405
+ " client_inputs = gr.Dropdown(\n",
406
+ " choices = list(URL_MAPPINGS.keys()), \n",
407
+ " label=\"FDSN Client\", \n",
408
+ " info=\"Select one of the available FDSN clients\",\n",
409
+ " value = \"IRIS\",\n",
410
+ " interactive=True\n",
411
+ " )\n",
412
+ " with gr.Row(): \n",
413
+ "\n",
414
+ " timestamp_inputs = gr.Textbox(value='2019-07-04 17:33:49',\n",
415
+ " placeholder='YYYY-MM-DD HH:MM:SS',\n",
416
+ " label=\"Timestamp\",\n",
417
+ " info=\"Timestamp of the earthquake\",\n",
418
+ " max_lines=1,\n",
419
+ " interactive=True)\n",
420
+ " \n",
421
+ " eq_lat_inputs = gr.Number(value=35.766, \n",
422
+ " label=\"Latitude\", \n",
423
+ " info=\"Latitude of the earthquake\",\n",
424
+ " interactive=True)\n",
425
+ " \n",
426
+ " eq_lon_inputs = gr.Number(value=-117.605,\n",
427
+ " label=\"Longitude\",\n",
428
+ " info=\"Longitude of the earthquake\",\n",
429
+ " interactive=True)\n",
430
+ " \n",
431
+ " source_depth_inputs = gr.Number(value=10,\n",
432
+ " label=\"Source depth (km)\",\n",
433
+ " info=\"Depth of the earthquake\",\n",
434
+ " interactive=True)\n",
435
+ " \n",
436
+ " radius_inputs = gr.Slider(minimum=1, \n",
437
+ " maximum=150, \n",
438
+ " value=50, label=\"Radius (km)\", \n",
439
+ " info=\"Select the radius around the earthquake to download data from\",\n",
440
+ " interactive=True)\n",
441
+ " \n",
442
+ " velocity_inputs = gr.Dropdown(\n",
443
+ " choices = ['1066a', '1066b', 'ak135', 'ak135f', 'herrin', 'iasp91', 'jb', 'prem', 'pwdk'], \n",
444
+ " label=\"1D velocity model\", \n",
445
+ " info=\"Velocity model for station selection\",\n",
446
+ " value = \"1066a\",\n",
447
+ " interactive=True\n",
448
+ " )\n",
449
+ " \n",
450
+ " \n",
451
+ " button = gr.Button(\"Predict phases\")\n",
452
+ " outputs_section = gr.outputs.Image(label='Waveforms with Phases Marked', type='numpy')\n",
453
+ " \n",
454
+ " button.click(predict_on_section, \n",
455
+ " inputs=[client_inputs, timestamp_inputs, \n",
456
+ " eq_lat_inputs, eq_lon_inputs, \n",
457
+ " radius_inputs, source_depth_inputs, velocity_inputs],\n",
458
+ " outputs=outputs_section)\n",
459
+ "\n",
460
+ " with gr.Tab(\"Predict on your own waveform\"):\n",
461
+ " gr.Markdown(\"\"\"\n",
462
+ " Please upload your waveform in .npy (numpy) format. \n",
463
+ " Your waveform should be sampled at 100 sps and have 3 (Z, N, E) or 1 (Z) channels.\n",
464
+ " \"\"\")\n",
465
+ "\n",
466
+ "\n",
467
+ "\n",
468
+ "demo.launch()"
469
+ ]
470
+ },
471
+ {
472
+ "cell_type": "code",
473
+ "execution_count": null,
474
+ "metadata": {},
475
+ "outputs": [],
476
+ "source": []
477
+ }
478
+ ],
479
+ "metadata": {
480
+ "kernelspec": {
481
+ "display_name": "phasehunter",
482
+ "language": "python",
483
+ "name": "python3"
484
+ },
485
+ "language_info": {
486
+ "codemirror_mode": {
487
+ "name": "ipython",
488
+ "version": 3
489
+ },
490
+ "file_extension": ".py",
491
+ "mimetype": "text/x-python",
492
+ "name": "python",
493
+ "nbconvert_exporter": "python",
494
+ "pygments_lexer": "ipython3",
495
+ "version": "3.11.2"
496
+ },
497
+ "orig_nbformat": 4,
498
+ "vscode": {
499
+ "interpreter": {
500
+ "hash": "6bf57068982d7b420bddaaf1d0614a7795947176033057024cf47d8ca2c1c4cd"
501
+ }
502
+ }
503
+ },
504
+ "nbformat": 4,
505
+ "nbformat_minor": 2
506
+ }
phasehunter/.DS_Store ADDED
Binary file (6.15 kB). View file
 
phasehunter/__pycache__/data_preparation.cpython-311.pyc ADDED
Binary file (9.14 kB). View file
 
phasehunter/__pycache__/model.cpython-311.pyc ADDED
Binary file (16.4 kB). View file
 
phasehunter/app.py DELETED
@@ -1,188 +0,0 @@
1
- # Gradio app that takes seismic waveform as input and marks 2 phases on the waveform as output.
2
-
3
- import gradio as gr
4
- import numpy as np
5
- from phasehunter.model import Onset_picker, Updated_onset_picker
6
- from phasehunter.data_preparation import prepare_waveform
7
- import torch
8
-
9
- from scipy.stats import gaussian_kde
10
-
11
- import obspy
12
- from obspy.clients.fdsn import Client
13
- from obspy.clients.fdsn.header import FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException
14
- from obspy.geodetics.base import locations2degrees
15
- from obspy.taup import TauPyModel
16
- from obspy.taup.helper_classes import SlownessModelError
17
-
18
- from obspy.clients.fdsn.header import URL_MAPPINGS
19
-
20
- import matplotlib.pyplot as plt
21
-
22
- def make_prediction(waveform):
23
- waveform = np.load(waveform)
24
- processed_input = prepare_waveform(waveform)
25
-
26
- # Make prediction
27
- with torch.no_grad():
28
- output = model(processed_input)
29
-
30
- p_phase = output[:, 0]
31
- s_phase = output[:, 1]
32
-
33
- return processed_input, p_phase, s_phase
34
-
35
- def mark_phases(waveform):
36
- processed_input, p_phase, s_phase = make_prediction(waveform)
37
-
38
- # Create a plot of the waveform with the phases marked
39
- if sum(processed_input[0][2] == 0): #if input is 1C
40
- fig, ax = plt.subplots(nrows=2, figsize=(10, 2), sharex=True)
41
-
42
- ax[0].plot(processed_input[0][0])
43
- ax[0].set_ylabel('Norm. Ampl.')
44
-
45
- else: #if input is 3C
46
- fig, ax = plt.subplots(nrows=4, figsize=(10, 6), sharex=True)
47
- ax[0].plot(processed_input[0][0])
48
- ax[1].plot(processed_input[0][1])
49
- ax[2].plot(processed_input[0][2])
50
-
51
- ax[0].set_ylabel('Z')
52
- ax[1].set_ylabel('N')
53
- ax[2].set_ylabel('E')
54
-
55
- p_phase_plot = p_phase*processed_input.shape[-1]
56
- p_kde = gaussian_kde(p_phase_plot)
57
- p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )
58
- ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')
59
-
60
- s_phase_plot = s_phase*processed_input.shape[-1]
61
- s_kde = gaussian_kde(s_phase_plot)
62
- s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )
63
- ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')
64
-
65
- for a in ax:
66
- a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P')
67
- a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S')
68
-
69
- ax[-1].set_xlabel('Time, samples')
70
- ax[-1].set_ylabel('Uncert.')
71
- ax[-1].legend()
72
-
73
- plt.subplots_adjust(hspace=0., wspace=0.)
74
-
75
- # Convert the plot to an image and return it
76
- fig.canvas.draw()
77
- image = np.array(fig.canvas.renderer.buffer_rgba())
78
- plt.close(fig)
79
- return image
80
-
81
- #??
82
-
83
- def download_data(timestamp, eq_lat, eq_lon, client_name, radius_km):
84
- client = Client(client_name)
85
- window = radius_km / 111.2
86
-
87
- assert eq_lat - window > -90 and eq_lat + window < 90, "Latitude out of bounds"
88
- assert eq_lon - window > -180 and eq_lon + window < 180, "Longitude out of bounds"
89
-
90
- starttime = obspy.UTCDateTime(timestamp)
91
- endtime = startime + 120
92
-
93
- inv = client.get_stations(network="*", station="*", location="*", channel="*H*",
94
- starttime=obspy.UTCDateTime(starttime), endtime=endtime,
95
- minlatitude=eq_lat-window, maxlatitude=eq_lat+window,
96
- minlongitude=eq_lon-window, maxlongitude=eq_lon+window,
97
- level='channel')
98
-
99
- for network in inv:
100
- for station in network:
101
- print(station)
102
-
103
- # waveform = client.get_waveforms(network=network.code, station=station.code, location="*", channel="*",
104
- # starttime=obspy.UTCDateTime(start_date), endtime=obspy.UTCDateTime(end_date))
105
-
106
- return 0
107
-
108
- model = Onset_picker.load_from_checkpoint("./weights.ckpt",
109
- picker=Updated_onset_picker(),
110
- learning_rate=3e-4)
111
- model.eval()
112
-
113
-
114
-
115
- # # Create the Gradio interface
116
- # gr.Interface(mark_phases, inputs, outputs, title='PhaseHunter').launch()
117
-
118
-
119
- with gr.Blocks() as demo:
120
- gr.Markdown("# PhaseHunter")
121
- gr.Markdown("""This app allows one to detect P and S seismic phases along with uncertainty of the detection.
122
- The app can be used in three ways: either by selecting one of the sample waveforms;
123
- or by selecting an earthquake from the global earthquake catalogue;
124
- or by uploading a waveform of interest.
125
- """)
126
- with gr.Tab("Default example"):
127
- # Define the input and output types for Gradio
128
- inputs = gr.Dropdown(
129
- ["data/sample/sample_0.npy",
130
- "data/sample/sample_1.npy",
131
- "data/sample/sample_2.npy"],
132
- label="Sample waveform",
133
- info="Select one of the samples",
134
- value = "data/sample/sample_0.npy"
135
- )
136
-
137
- button = gr.Button("Predict phases")
138
- outputs = gr.outputs.Image(label='Waveform with Phases Marked', type='numpy')
139
-
140
- button.click(mark_phases, inputs=inputs, outputs=outputs)
141
-
142
- with gr.Tab("Select earthquake from catalogue"):
143
- gr.Markdown('TEST')
144
-
145
- client_inputs = gr.Dropdown(
146
- choices = list(URL_MAPPINGS.keys()),
147
- label="FDSN Client",
148
- info="Select one of the available FDSN clients",
149
- value = "IRIS",
150
- interactive=True
151
- )
152
- with gr.Row():
153
-
154
- timestamp_inputs = gr.Textbox(value='2019-07-04 17:33:49',
155
- placeholder='YYYY-MM-DD HH:MM:SS',
156
- label="Timestamp",
157
- info="Timestamp of the earthquake",
158
- max_lines=1,
159
- interactive=True)
160
-
161
- eq_lat_inputs = gr.Number(value=35.766,
162
- label="Latitude",
163
- info="Latitude of the earthquake",
164
- interactive=True)
165
-
166
- eq_lo_inputs = gr.Number(value=117.605,
167
- label="Longitude",
168
- info="Longitude of the earthquake",
169
- interactive=True)
170
-
171
- radius_inputs = gr.Slider(minimum=1,
172
- maximum=150,
173
- value=50, label="Radius (km)",
174
- info="Select the radius around the earthquake to download data from",
175
- interactive=True)
176
-
177
- button = gr.Button("Predict phases")
178
- button.click(mark_phases, inputs=inputs, outputs=outputs)
179
-
180
- with gr.Tab("Predict on your own waveform"):
181
- gr.Markdown("""
182
- Please upload your waveform in .npy (numpy) format.
183
- Your waveform should be sampled at 100 sps and have 3 (Z, N, E) or 1 (Z) channels.
184
- """)
185
-
186
- button.click(download_data, inputs=[timestamp_inputs, eq_lat_inputs,eq_lo_inputs, radius_inputs], outputs=outputs)
187
-
188
- demo.launch()