Spaces:
Runtime error
Runtime error
Added velocity model builder 0.01
Browse files- .DS_Store +0 -0
- Gradio_app.ipynb +351 -110
- app.py +391 -455
- data/.DS_Store +0 -0
- data/velocity/35.766_-117.605_10.0_2019-07-04 17:33:49_3.csv +4 -0
- data/velocity/35.766_-117.605_2019-07-04 17:33:49_3.csv +3 -3
- data/velocity/current_vel_model.csv +4 -0
- phasehunter/__pycache__/data_preparation.cpython-311.pyc +0 -0
.DS_Store
CHANGED
Binary files a/.DS_Store and b/.DS_Store differ
|
|
Gradio_app.ipynb
CHANGED
@@ -2,14 +2,14 @@
|
|
2 |
"cells": [
|
3 |
{
|
4 |
"cell_type": "code",
|
5 |
-
"execution_count":
|
6 |
"metadata": {},
|
7 |
"outputs": [
|
8 |
{
|
9 |
"name": "stdout",
|
10 |
"output_type": "stream",
|
11 |
"text": [
|
12 |
-
"Running on local URL: http://127.0.0.1:
|
13 |
"\n",
|
14 |
"To create a public link, set `share=True` in `launch()`.\n"
|
15 |
]
|
@@ -17,7 +17,7 @@
|
|
17 |
{
|
18 |
"data": {
|
19 |
"text/html": [
|
20 |
-
"<div><iframe src=\"http://127.0.0.1:
|
21 |
],
|
22 |
"text/plain": [
|
23 |
"<IPython.core.display.HTML object>"
|
@@ -30,7 +30,7 @@
|
|
30 |
"data": {
|
31 |
"text/plain": []
|
32 |
},
|
33 |
-
"execution_count":
|
34 |
"metadata": {},
|
35 |
"output_type": "execute_result"
|
36 |
},
|
@@ -38,16 +38,215 @@
|
|
38 |
"name": "stdout",
|
39 |
"output_type": "stream",
|
40 |
"text": [
|
41 |
-
"
|
42 |
-
"
|
43 |
-
"
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
44 |
]
|
45 |
},
|
46 |
{
|
47 |
"name": "stderr",
|
48 |
"output_type": "stream",
|
49 |
"text": [
|
50 |
-
"
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
51 |
]
|
52 |
}
|
53 |
],
|
@@ -455,9 +654,10 @@
|
|
455 |
" fig.canvas.draw();\n",
|
456 |
" image = np.array(fig.canvas.renderer.buffer_rgba())\n",
|
457 |
" plt.close(fig)\n",
|
458 |
-
" output_picks.to_csv(f'data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv', index=False)\n",
|
459 |
-
" output_csv = f'data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv'\n",
|
460 |
"\n",
|
|
|
|
|
|
|
461 |
" return image, output_picks, output_csv\n",
|
462 |
"\n",
|
463 |
"import numpy as np\n",
|
@@ -468,17 +668,16 @@
|
|
468 |
" return np.argmin(np.abs(array - value))\n",
|
469 |
"\n",
|
470 |
"def compute_velocity_model(azimuth, elevation):\n",
|
471 |
-
" filename = output_csv\n",
|
|
|
472 |
" df = pd.read_csv(filename)\n",
|
473 |
-
"\n",
|
|
|
|
|
474 |
" # Current EQ location\n",
|
475 |
-
" filename = filename.split(\"/\")[-1]\n",
|
476 |
" eq_lat = float(filename.split(\"_\")[0])\n",
|
477 |
" eq_lon = float(filename.split(\"_\")[1])\n",
|
478 |
-
"\n",
|
479 |
-
" ## FIX THIS LATTER\n",
|
480 |
-
" eq_depth = 10 ##FIX THIS LATTER\n",
|
481 |
-
" ## FIX THIS LATTER\n",
|
482 |
"\n",
|
483 |
" # Define the region of interest (latitude, longitude, and depth ranges)\n",
|
484 |
" lat_range = (np.min([df.st_lat.min(), eq_lat]), np.max([df.st_lat.max(), eq_lat]))\n",
|
@@ -595,7 +794,53 @@
|
|
595 |
" </ul>\n",
|
596 |
"</div>\n",
|
597 |
"\"\"\")\n",
|
598 |
-
"
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
599 |
" with gr.Tab(\"Select earthquake from catalogue\"):\n",
|
600 |
"\n",
|
601 |
" gr.HTML(\"\"\"\n",
|
@@ -717,53 +962,6 @@
|
|
717 |
" inputs=inputs_vel_model, \n",
|
718 |
" outputs=outputs_vel_model)\n",
|
719 |
"\n",
|
720 |
-
" with gr.Tab(\"Try on a single station\"):\n",
|
721 |
-
" with gr.Row(): \n",
|
722 |
-
" # Define the input and output types for Gradio\n",
|
723 |
-
" inputs = gr.Dropdown(\n",
|
724 |
-
" [\"data/sample/sample_0.npy\", \n",
|
725 |
-
" \"data/sample/sample_1.npy\", \n",
|
726 |
-
" \"data/sample/sample_2.npy\"], \n",
|
727 |
-
" label=\"Sample waveform\", \n",
|
728 |
-
" info=\"Select one of the samples\",\n",
|
729 |
-
" value = \"data/sample/sample_0.npy\"\n",
|
730 |
-
" )\n",
|
731 |
-
" with gr.Column(scale=1):\n",
|
732 |
-
" P_thres_inputs = gr.Slider(minimum=0.01,\n",
|
733 |
-
" maximum=1,\n",
|
734 |
-
" value=0.1,\n",
|
735 |
-
" label=\"P uncertainty threshold, s\",\n",
|
736 |
-
" step=0.01,\n",
|
737 |
-
" info=\"Acceptable uncertainty for P picks expressed in std() seconds\",\n",
|
738 |
-
" interactive=True,\n",
|
739 |
-
" )\n",
|
740 |
-
" \n",
|
741 |
-
" S_thres_inputs = gr.Slider(minimum=0.01,\n",
|
742 |
-
" maximum=1,\n",
|
743 |
-
" value=0.2,\n",
|
744 |
-
" label=\"S uncertainty threshold, s\",\n",
|
745 |
-
" step=0.01,\n",
|
746 |
-
" info=\"Acceptable uncertainty for S picks expressed in std() seconds\",\n",
|
747 |
-
" interactive=True,\n",
|
748 |
-
" )\n",
|
749 |
-
" with gr.Column(scale=1):\n",
|
750 |
-
" upload = gr.File(label=\"Or upload your own waveform\")\n",
|
751 |
-
" sampling_rate_inputs = gr.Slider(minimum=10,\n",
|
752 |
-
" maximum=1000,\n",
|
753 |
-
" value=100,\n",
|
754 |
-
" label=\"Samlping rate, Hz\",\n",
|
755 |
-
" step=10,\n",
|
756 |
-
" info=\"Sampling rate of the waveform\",\n",
|
757 |
-
" interactive=True,\n",
|
758 |
-
" )\n",
|
759 |
-
"\n",
|
760 |
-
" button = gr.Button(\"Predict phases\")\n",
|
761 |
-
" outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)\n",
|
762 |
-
" \n",
|
763 |
-
" button.click(mark_phases, inputs=[inputs, upload, \n",
|
764 |
-
" P_thres_inputs, S_thres_inputs,\n",
|
765 |
-
" sampling_rate_inputs], \n",
|
766 |
-
" outputs=outputs)\n",
|
767 |
"\n",
|
768 |
" \n",
|
769 |
"\n",
|
@@ -840,54 +1038,102 @@
|
|
840 |
},
|
841 |
{
|
842 |
"cell_type": "code",
|
843 |
-
"execution_count":
|
844 |
"metadata": {},
|
845 |
"outputs": [
|
846 |
-
{
|
847 |
-
"name": "stdout",
|
848 |
-
"output_type": "stream",
|
849 |
-
"text": [
|
850 |
-
"Running on local URL: http://127.0.0.1:7869\n",
|
851 |
-
"\n",
|
852 |
-
"To create a public link, set `share=True` in `launch()`.\n"
|
853 |
-
]
|
854 |
-
},
|
855 |
{
|
856 |
"data": {
|
857 |
-
"text/html": [
|
858 |
-
"<div><iframe src=\"http://127.0.0.1:7869/\" width=\"100%\" height=\"500\" allow=\"autoplay; camera; microphone; clipboard-read; clipboard-write;\" frameborder=\"0\" allowfullscreen></iframe></div>"
|
859 |
-
],
|
860 |
"text/plain": [
|
861 |
-
"
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
862 |
]
|
863 |
},
|
864 |
-
"
|
865 |
-
"output_type": "display_data"
|
866 |
-
},
|
867 |
-
{
|
868 |
-
"data": {
|
869 |
-
"text/plain": []
|
870 |
-
},
|
871 |
-
"execution_count": 24,
|
872 |
"metadata": {},
|
873 |
"output_type": "execute_result"
|
874 |
-
},
|
875 |
-
{
|
876 |
-
"name": "stderr",
|
877 |
-
"output_type": "stream",
|
878 |
-
"text": [
|
879 |
-
"/var/folders/ky/4j6xbvhs5m583jflkhyzxf9h0000gn/T/ipykernel_19576/4006067033.py:95: MatplotlibDeprecationWarning: Unable to determine Axes to steal space for Colorbar. Using gca(), but will raise in the future. Either provide the *cax* argument to use as the Axes for the Colorbar, provide the *ax* argument to steal space from it, or add *mappable* to an Axes.\n",
|
880 |
-
" plt.colorbar(m)\n"
|
881 |
-
]
|
882 |
}
|
883 |
],
|
884 |
"source": [
|
885 |
-
"
|
886 |
-
"import numpy as np\n",
|
887 |
-
"import gradio as gr\n",
|
888 |
-
" \n",
|
889 |
-
"# Define the Gradio interface\n",
|
890 |
-
"\n"
|
891 |
]
|
892 |
},
|
893 |
{
|
@@ -916,14 +1162,9 @@
|
|
916 |
}
|
917 |
],
|
918 |
"source": [
|
919 |
-
"
|
920 |
-
"
|
921 |
-
"
|
922 |
-
"\n",
|
923 |
-
"y_resampled = resample_waveform(y, 1000, 100)\n",
|
924 |
-
"# plot sin\n",
|
925 |
-
"plt.plot(x, y)\n",
|
926 |
-
"plt.plot(x, y_resampled)"
|
927 |
]
|
928 |
}
|
929 |
],
|
@@ -943,7 +1184,7 @@
|
|
943 |
"name": "python",
|
944 |
"nbconvert_exporter": "python",
|
945 |
"pygments_lexer": "ipython3",
|
946 |
-
"version": "3.
|
947 |
},
|
948 |
"orig_nbformat": 4
|
949 |
},
|
|
|
2 |
"cells": [
|
3 |
{
|
4 |
"cell_type": "code",
|
5 |
+
"execution_count": 17,
|
6 |
"metadata": {},
|
7 |
"outputs": [
|
8 |
{
|
9 |
"name": "stdout",
|
10 |
"output_type": "stream",
|
11 |
"text": [
|
12 |
+
"Running on local URL: http://127.0.0.1:7864\n",
|
13 |
"\n",
|
14 |
"To create a public link, set `share=True` in `launch()`.\n"
|
15 |
]
|
|
|
17 |
{
|
18 |
"data": {
|
19 |
"text/html": [
|
20 |
+
"<div><iframe src=\"http://127.0.0.1:7864/\" width=\"100%\" height=\"500\" allow=\"autoplay; camera; microphone; clipboard-read; clipboard-write;\" frameborder=\"0\" allowfullscreen></iframe></div>"
|
21 |
],
|
22 |
"text/plain": [
|
23 |
"<IPython.core.display.HTML object>"
|
|
|
30 |
"data": {
|
31 |
"text/plain": []
|
32 |
},
|
33 |
+
"execution_count": 17,
|
34 |
"metadata": {},
|
35 |
"output_type": "execute_result"
|
36 |
},
|
|
|
38 |
"name": "stdout",
|
39 |
"output_type": "stream",
|
40 |
"text": [
|
41 |
+
"Starting to download inventory\n",
|
42 |
+
"Finished downloading inventory\n",
|
43 |
+
"Processing CI.CCC...\n",
|
44 |
+
"Downloading waveform for CI_CCC_2019-07-04T17:33:40.494920Z\n",
|
45 |
+
"Skipping CI_CCC_2019-07-04T17:33:40.494920Z\n",
|
46 |
+
"Processing CI.CLC...\n",
|
47 |
+
"Processing CI.JRC2...\n",
|
48 |
+
"Reading cached waveform\n",
|
49 |
+
"Added CI.JRC2 to the list of waveforms\n",
|
50 |
+
"Processing CI.LRL...\n",
|
51 |
+
"Reading cached waveform\n",
|
52 |
+
"Added CI.LRL to the list of waveforms\n",
|
53 |
+
"Processing CI.MPM...\n",
|
54 |
+
"Reading cached waveform\n",
|
55 |
+
"Processing CI.Q0072...\n",
|
56 |
+
"Reading cached waveform\n",
|
57 |
+
"Processing CI.SLA...\n",
|
58 |
+
"Reading cached waveform\n",
|
59 |
+
"Added CI.SLA to the list of waveforms\n",
|
60 |
+
"Processing CI.SRT...\n",
|
61 |
+
"Reading cached waveform\n",
|
62 |
+
"Added CI.SRT to the list of waveforms\n",
|
63 |
+
"Processing CI.TOW2...\n",
|
64 |
+
"Reading cached waveform\n",
|
65 |
+
"Added CI.TOW2 to the list of waveforms\n",
|
66 |
+
"Processing CI.WBM...\n",
|
67 |
+
"Downloading waveform for CI_WBM_2019-07-04T17:33:40.063616Z\n",
|
68 |
+
"Skipping CI_WBM_2019-07-04T17:33:40.063616Z\n",
|
69 |
+
"Processing CI.WCS2...\n",
|
70 |
+
"Downloading waveform for CI_WCS2_2019-07-04T17:33:40.200958Z\n",
|
71 |
+
"Skipping CI_WCS2_2019-07-04T17:33:40.200958Z\n",
|
72 |
+
"Processing CI.WMF...\n",
|
73 |
+
"Reading cached waveform\n",
|
74 |
+
"Added CI.WMF to the list of waveforms\n",
|
75 |
+
"Processing CI.WNM...\n",
|
76 |
+
"Reading cached waveform\n",
|
77 |
+
"Processing CI.WRC2...\n",
|
78 |
+
"Downloading waveform for CI_WRC2_2019-07-04T17:33:38.698099Z\n",
|
79 |
+
"Skipping CI_WRC2_2019-07-04T17:33:38.698099Z\n",
|
80 |
+
"Processing CI.WRV2...\n",
|
81 |
+
"Reading cached waveform\n",
|
82 |
+
"Processing CI.WVP2...\n",
|
83 |
+
"Downloading waveform for CI_WVP2_2019-07-04T17:33:39.650402Z\n",
|
84 |
+
"Skipping CI_WVP2_2019-07-04T17:33:39.650402Z\n",
|
85 |
+
"Processing NP.1809...\n",
|
86 |
+
"Reading cached waveform\n",
|
87 |
+
"Processing NP.5419...\n",
|
88 |
+
"Reading cached waveform\n",
|
89 |
+
"Processing PB.B916...\n",
|
90 |
+
"Reading cached waveform\n",
|
91 |
+
"Processing PB.B917...\n",
|
92 |
+
"Reading cached waveform\n",
|
93 |
+
"Processing PB.B918...\n",
|
94 |
+
"Reading cached waveform\n",
|
95 |
+
"Processing PB.B921...\n",
|
96 |
+
"Reading cached waveform\n",
|
97 |
+
"Starting to run predictions\n"
|
98 |
]
|
99 |
},
|
100 |
{
|
101 |
"name": "stderr",
|
102 |
"output_type": "stream",
|
103 |
"text": [
|
104 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:273: FutureWarning: The input object of type 'Tensor' is an array-like implementing one of the corresponding protocols (`__array__`, `__array_interface__` or `__array_struct__`); but not a sequence (or 0-D). In the future, this object will be coerced as if it was first converted using `np.array(obj)`. To retain the old behaviour, you have to either modify the type 'Tensor', or assign to an empty array created with `np.empty(correct_shape, dtype=object)`.\n",
|
105 |
+
" waveforms = np.array(waveforms)[selection_indexes]\n",
|
106 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:273: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
|
107 |
+
" waveforms = np.array(waveforms)[selection_indexes]\n",
|
108 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:280: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
|
109 |
+
" waveforms = [torch.tensor(waveform) for waveform in waveforms]\n"
|
110 |
+
]
|
111 |
+
},
|
112 |
+
{
|
113 |
+
"name": "stdout",
|
114 |
+
"output_type": "stream",
|
115 |
+
"text": [
|
116 |
+
"Starting plotting 3 waveforms\n",
|
117 |
+
"Fetching topography\n",
|
118 |
+
"Plotting topo\n"
|
119 |
+
]
|
120 |
+
},
|
121 |
+
{
|
122 |
+
"name": "stderr",
|
123 |
+
"output_type": "stream",
|
124 |
+
"text": [
|
125 |
+
"/Users/anovosel/miniconda3/envs/phasehunter/lib/python3.11/site-packages/bmi_topography/api_key.py:49: UserWarning: You are using a demo key to fetch data from OpenTopography, functionality will be limited. See https://bmi-topography.readthedocs.io/en/latest/#api-key for more information.\n",
|
126 |
+
" warnings.warn(\n"
|
127 |
+
]
|
128 |
+
},
|
129 |
+
{
|
130 |
+
"name": "stdout",
|
131 |
+
"output_type": "stream",
|
132 |
+
"text": [
|
133 |
+
"Plotting waveform 1/3\n",
|
134 |
+
"Station 35.98249, -117.80885 has P velocity 4.13660431013202 and S velocity 2.2622770044299756\n",
|
135 |
+
"Plotting waveform 2/3\n"
|
136 |
+
]
|
137 |
+
},
|
138 |
+
{
|
139 |
+
"name": "stderr",
|
140 |
+
"output_type": "stream",
|
141 |
+
"text": [
|
142 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:365: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.\n",
|
143 |
+
" output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]],\n",
|
144 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:365: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.\n",
|
145 |
+
" output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]],\n"
|
146 |
+
]
|
147 |
+
},
|
148 |
+
{
|
149 |
+
"name": "stdout",
|
150 |
+
"output_type": "stream",
|
151 |
+
"text": [
|
152 |
+
"Station 35.69235, -117.75051 has P velocity 3.4155476453388767 and S velocity 1.67967367867923\n",
|
153 |
+
"Plotting waveform 3/3\n",
|
154 |
+
"Station 36.11758, -117.85486 has P velocity 4.745724852828504 and S velocity 2.6483289549749593\n",
|
155 |
+
"Plotting stations\n"
|
156 |
+
]
|
157 |
+
},
|
158 |
+
{
|
159 |
+
"name": "stderr",
|
160 |
+
"output_type": "stream",
|
161 |
+
"text": [
|
162 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:365: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.\n",
|
163 |
+
" output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]],\n",
|
164 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:385: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
|
165 |
+
" ax[i].set_xticklabels(ax[i].get_xticks(), rotation = 50)\n"
|
166 |
+
]
|
167 |
+
},
|
168 |
+
{
|
169 |
+
"name": "stdout",
|
170 |
+
"output_type": "stream",
|
171 |
+
"text": [
|
172 |
+
" station_name st_lat st_lon starttime p_phase, s \\\n",
|
173 |
+
"0 CI.JRC2 35.98249 -117.80885 2019-07-04T17:33:39.947494Z 7.320212 \n",
|
174 |
+
"1 CI.SRT 35.69235 -117.75051 2019-07-04T17:33:38.029990Z 4.532020 \n",
|
175 |
+
"2 CI.WMF 36.11758 -117.85486 2019-07-04T17:33:41.867962Z 9.504385 \n",
|
176 |
+
"\n",
|
177 |
+
" p_uncertainty, s s_phase, s s_uncertainty, s velocity_p, km/s \\\n",
|
178 |
+
"0 0.020417 13.385108 0.028439 4.136604 \n",
|
179 |
+
"1 0.017490 9.215676 0.019568 3.415548 \n",
|
180 |
+
"2 0.015920 17.031569 0.046738 4.745725 \n",
|
181 |
+
"\n",
|
182 |
+
" velocity_s, km/s \n",
|
183 |
+
"0 2.262277 \n",
|
184 |
+
"1 1.679674 \n",
|
185 |
+
"2 2.648329 \n"
|
186 |
+
]
|
187 |
+
},
|
188 |
+
{
|
189 |
+
"name": "stderr",
|
190 |
+
"output_type": "stream",
|
191 |
+
"text": [
|
192 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:503: MatplotlibDeprecationWarning: Unable to determine Axes to steal space for Colorbar. Using gca(), but will raise in the future. Either provide the *cax* argument to use as the Axes for the Colorbar, provide the *ax* argument to steal space from it, or add *mappable* to an Axes.\n",
|
193 |
+
" plt.colorbar(m)\n"
|
194 |
+
]
|
195 |
+
},
|
196 |
+
{
|
197 |
+
"name": "stdout",
|
198 |
+
"output_type": "stream",
|
199 |
+
"text": [
|
200 |
+
" station_name st_lat st_lon starttime p_phase, s \\\n",
|
201 |
+
"0 CI.JRC2 35.98249 -117.80885 2019-07-04T17:33:39.947494Z 7.320212 \n",
|
202 |
+
"1 CI.SRT 35.69235 -117.75051 2019-07-04T17:33:38.029990Z 4.532020 \n",
|
203 |
+
"2 CI.WMF 36.11758 -117.85486 2019-07-04T17:33:41.867962Z 9.504385 \n",
|
204 |
+
"\n",
|
205 |
+
" p_uncertainty, s s_phase, s s_uncertainty, s velocity_p, km/s \\\n",
|
206 |
+
"0 0.020417 13.385108 0.028439 4.136604 \n",
|
207 |
+
"1 0.017490 9.215676 0.019568 3.415548 \n",
|
208 |
+
"2 0.015920 17.031569 0.046738 4.745725 \n",
|
209 |
+
"\n",
|
210 |
+
" velocity_s, km/s \n",
|
211 |
+
"0 2.262277 \n",
|
212 |
+
"1 1.679674 \n",
|
213 |
+
"2 2.648329 \n"
|
214 |
+
]
|
215 |
+
},
|
216 |
+
{
|
217 |
+
"name": "stderr",
|
218 |
+
"output_type": "stream",
|
219 |
+
"text": [
|
220 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:503: MatplotlibDeprecationWarning: Unable to determine Axes to steal space for Colorbar. Using gca(), but will raise in the future. Either provide the *cax* argument to use as the Axes for the Colorbar, provide the *ax* argument to steal space from it, or add *mappable* to an Axes.\n",
|
221 |
+
" plt.colorbar(m)\n"
|
222 |
+
]
|
223 |
+
},
|
224 |
+
{
|
225 |
+
"name": "stdout",
|
226 |
+
"output_type": "stream",
|
227 |
+
"text": [
|
228 |
+
" station_name st_lat st_lon starttime p_phase, s \\\n",
|
229 |
+
"0 CI.JRC2 35.98249 -117.80885 2019-07-04T17:33:39.947494Z 7.320212 \n",
|
230 |
+
"1 CI.SRT 35.69235 -117.75051 2019-07-04T17:33:38.029990Z 4.532020 \n",
|
231 |
+
"2 CI.WMF 36.11758 -117.85486 2019-07-04T17:33:41.867962Z 9.504385 \n",
|
232 |
+
"\n",
|
233 |
+
" p_uncertainty, s s_phase, s s_uncertainty, s velocity_p, km/s \\\n",
|
234 |
+
"0 0.020417 13.385108 0.028439 4.136604 \n",
|
235 |
+
"1 0.017490 9.215676 0.019568 3.415548 \n",
|
236 |
+
"2 0.015920 17.031569 0.046738 4.745725 \n",
|
237 |
+
"\n",
|
238 |
+
" velocity_s, km/s \n",
|
239 |
+
"0 2.262277 \n",
|
240 |
+
"1 1.679674 \n",
|
241 |
+
"2 2.648329 \n"
|
242 |
+
]
|
243 |
+
},
|
244 |
+
{
|
245 |
+
"name": "stderr",
|
246 |
+
"output_type": "stream",
|
247 |
+
"text": [
|
248 |
+
"/var/folders/_g/3q5q8_dj0ydcpktxlwxb5vrh0000gq/T/ipykernel_3502/4124724611.py:503: MatplotlibDeprecationWarning: Unable to determine Axes to steal space for Colorbar. Using gca(), but will raise in the future. Either provide the *cax* argument to use as the Axes for the Colorbar, provide the *ax* argument to steal space from it, or add *mappable* to an Axes.\n",
|
249 |
+
" plt.colorbar(m)\n"
|
250 |
]
|
251 |
}
|
252 |
],
|
|
|
654 |
" fig.canvas.draw();\n",
|
655 |
" image = np.array(fig.canvas.renderer.buffer_rgba())\n",
|
656 |
" plt.close(fig)\n",
|
|
|
|
|
657 |
"\n",
|
658 |
+
" output_csv = f'data/velocity/{eq_lat}_{eq_lon}_{source_depth_km}_{timestamp}_{len(waveforms)}.csv'\n",
|
659 |
+
" output_picks.to_csv(output_csv, index=False)\n",
|
660 |
+
" \n",
|
661 |
" return image, output_picks, output_csv\n",
|
662 |
"\n",
|
663 |
"import numpy as np\n",
|
|
|
668 |
" return np.argmin(np.abs(array - value))\n",
|
669 |
"\n",
|
670 |
"def compute_velocity_model(azimuth, elevation):\n",
|
671 |
+
" filename = list(output_csv.temp_files)[0]\n",
|
672 |
+
" \n",
|
673 |
" df = pd.read_csv(filename)\n",
|
674 |
+
" print(df)\n",
|
675 |
+
" filename = filename.split('/')[-1]\n",
|
676 |
+
" \n",
|
677 |
" # Current EQ location\n",
|
|
|
678 |
" eq_lat = float(filename.split(\"_\")[0])\n",
|
679 |
" eq_lon = float(filename.split(\"_\")[1])\n",
|
680 |
+
" eq_depth = float(filename.split(\"_\")[2])\n",
|
|
|
|
|
|
|
681 |
"\n",
|
682 |
" # Define the region of interest (latitude, longitude, and depth ranges)\n",
|
683 |
" lat_range = (np.min([df.st_lat.min(), eq_lat]), np.max([df.st_lat.max(), eq_lat]))\n",
|
|
|
794 |
" </ul>\n",
|
795 |
"</div>\n",
|
796 |
"\"\"\")\n",
|
797 |
+
" with gr.Tab(\"Try on a single station\"):\n",
|
798 |
+
" with gr.Row(): \n",
|
799 |
+
" # Define the input and output types for Gradio\n",
|
800 |
+
" inputs = gr.Dropdown(\n",
|
801 |
+
" [\"data/sample/sample_0.npy\", \n",
|
802 |
+
" \"data/sample/sample_1.npy\", \n",
|
803 |
+
" \"data/sample/sample_2.npy\"], \n",
|
804 |
+
" label=\"Sample waveform\", \n",
|
805 |
+
" info=\"Select one of the samples\",\n",
|
806 |
+
" value = \"data/sample/sample_0.npy\"\n",
|
807 |
+
" )\n",
|
808 |
+
" with gr.Column(scale=1):\n",
|
809 |
+
" P_thres_inputs = gr.Slider(minimum=0.01,\n",
|
810 |
+
" maximum=1,\n",
|
811 |
+
" value=0.1,\n",
|
812 |
+
" label=\"P uncertainty threshold, s\",\n",
|
813 |
+
" step=0.01,\n",
|
814 |
+
" info=\"Acceptable uncertainty for P picks expressed in std() seconds\",\n",
|
815 |
+
" interactive=True,\n",
|
816 |
+
" )\n",
|
817 |
+
" \n",
|
818 |
+
" S_thres_inputs = gr.Slider(minimum=0.01,\n",
|
819 |
+
" maximum=1,\n",
|
820 |
+
" value=0.2,\n",
|
821 |
+
" label=\"S uncertainty threshold, s\",\n",
|
822 |
+
" step=0.01,\n",
|
823 |
+
" info=\"Acceptable uncertainty for S picks expressed in std() seconds\",\n",
|
824 |
+
" interactive=True,\n",
|
825 |
+
" )\n",
|
826 |
+
" with gr.Column(scale=1):\n",
|
827 |
+
" upload = gr.File(label=\"Or upload your own waveform\")\n",
|
828 |
+
" sampling_rate_inputs = gr.Slider(minimum=10,\n",
|
829 |
+
" maximum=1000,\n",
|
830 |
+
" value=100,\n",
|
831 |
+
" label=\"Samlping rate, Hz\",\n",
|
832 |
+
" step=10,\n",
|
833 |
+
" info=\"Sampling rate of the waveform\",\n",
|
834 |
+
" interactive=True,\n",
|
835 |
+
" )\n",
|
836 |
+
"\n",
|
837 |
+
" button = gr.Button(\"Predict phases\")\n",
|
838 |
+
" outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)\n",
|
839 |
+
" \n",
|
840 |
+
" button.click(mark_phases, inputs=[inputs, upload, \n",
|
841 |
+
" P_thres_inputs, S_thres_inputs,\n",
|
842 |
+
" sampling_rate_inputs], \n",
|
843 |
+
" outputs=outputs) \n",
|
844 |
" with gr.Tab(\"Select earthquake from catalogue\"):\n",
|
845 |
"\n",
|
846 |
" gr.HTML(\"\"\"\n",
|
|
|
962 |
" inputs=inputs_vel_model, \n",
|
963 |
" outputs=outputs_vel_model)\n",
|
964 |
"\n",
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
965 |
"\n",
|
966 |
" \n",
|
967 |
"\n",
|
|
|
1038 |
},
|
1039 |
{
|
1040 |
"cell_type": "code",
|
1041 |
+
"execution_count": 2,
|
1042 |
"metadata": {},
|
1043 |
"outputs": [
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1044 |
{
|
1045 |
"data": {
|
|
|
|
|
|
|
1046 |
"text/plain": [
|
1047 |
+
"['DEFAULT_TEMP_DIR',\n",
|
1048 |
+
" '__abstractmethods__',\n",
|
1049 |
+
" '__class__',\n",
|
1050 |
+
" '__delattr__',\n",
|
1051 |
+
" '__dict__',\n",
|
1052 |
+
" '__dir__',\n",
|
1053 |
+
" '__doc__',\n",
|
1054 |
+
" '__eq__',\n",
|
1055 |
+
" '__format__',\n",
|
1056 |
+
" '__ge__',\n",
|
1057 |
+
" '__getattribute__',\n",
|
1058 |
+
" '__getstate__',\n",
|
1059 |
+
" '__gt__',\n",
|
1060 |
+
" '__hash__',\n",
|
1061 |
+
" '__init__',\n",
|
1062 |
+
" '__init_subclass__',\n",
|
1063 |
+
" '__le__',\n",
|
1064 |
+
" '__lt__',\n",
|
1065 |
+
" '__module__',\n",
|
1066 |
+
" '__ne__',\n",
|
1067 |
+
" '__new__',\n",
|
1068 |
+
" '__reduce__',\n",
|
1069 |
+
" '__reduce_ex__',\n",
|
1070 |
+
" '__repr__',\n",
|
1071 |
+
" '__setattr__',\n",
|
1072 |
+
" '__sizeof__',\n",
|
1073 |
+
" '__slots__',\n",
|
1074 |
+
" '__str__',\n",
|
1075 |
+
" '__subclasshook__',\n",
|
1076 |
+
" '__weakref__',\n",
|
1077 |
+
" '_abc_impl',\n",
|
1078 |
+
" '_id',\n",
|
1079 |
+
" '_skip_init_processing',\n",
|
1080 |
+
" '_style',\n",
|
1081 |
+
" 'as_example',\n",
|
1082 |
+
" 'attach_load_event',\n",
|
1083 |
+
" 'base64_to_temp_file_if_needed',\n",
|
1084 |
+
" 'change',\n",
|
1085 |
+
" 'clear',\n",
|
1086 |
+
" 'deserialize',\n",
|
1087 |
+
" 'download_temp_copy_if_needed',\n",
|
1088 |
+
" 'elem_classes',\n",
|
1089 |
+
" 'elem_id',\n",
|
1090 |
+
" 'file_count',\n",
|
1091 |
+
" 'file_types',\n",
|
1092 |
+
" 'get_block_name',\n",
|
1093 |
+
" 'get_config',\n",
|
1094 |
+
" 'get_expected_parent',\n",
|
1095 |
+
" 'get_load_fn_and_initial_value',\n",
|
1096 |
+
" 'get_specific_update',\n",
|
1097 |
+
" 'hash_base64',\n",
|
1098 |
+
" 'hash_file',\n",
|
1099 |
+
" 'hash_url',\n",
|
1100 |
+
" 'info',\n",
|
1101 |
+
" 'interactive',\n",
|
1102 |
+
" 'label',\n",
|
1103 |
+
" 'load_event',\n",
|
1104 |
+
" 'load_event_to_attach',\n",
|
1105 |
+
" 'make_temp_copy_if_needed',\n",
|
1106 |
+
" 'parent',\n",
|
1107 |
+
" 'postprocess',\n",
|
1108 |
+
" 'preprocess',\n",
|
1109 |
+
" 'render',\n",
|
1110 |
+
" 'root',\n",
|
1111 |
+
" 'root_url',\n",
|
1112 |
+
" 'save_uploaded_file',\n",
|
1113 |
+
" 'select',\n",
|
1114 |
+
" 'selectable',\n",
|
1115 |
+
" 'serialize',\n",
|
1116 |
+
" 'set_event_trigger',\n",
|
1117 |
+
" 'share_token',\n",
|
1118 |
+
" 'show_label',\n",
|
1119 |
+
" 'style',\n",
|
1120 |
+
" 'temp_files',\n",
|
1121 |
+
" 'test_input',\n",
|
1122 |
+
" 'type',\n",
|
1123 |
+
" 'unrender',\n",
|
1124 |
+
" 'update',\n",
|
1125 |
+
" 'upload',\n",
|
1126 |
+
" 'value',\n",
|
1127 |
+
" 'visible']"
|
1128 |
]
|
1129 |
},
|
1130 |
+
"execution_count": 2,
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1131 |
"metadata": {},
|
1132 |
"output_type": "execute_result"
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1133 |
}
|
1134 |
],
|
1135 |
"source": [
|
1136 |
+
"dir(output_csv)"
|
|
|
|
|
|
|
|
|
|
|
1137 |
]
|
1138 |
},
|
1139 |
{
|
|
|
1162 |
}
|
1163 |
],
|
1164 |
"source": [
|
1165 |
+
"filename.split(\"/\")[-1]\n",
|
1166 |
+
" eq_lat = float(filename.split(\"_\")[0])\n",
|
1167 |
+
" eq_lon = float(filename.split(\"_\")[1])"
|
|
|
|
|
|
|
|
|
|
|
1168 |
]
|
1169 |
}
|
1170 |
],
|
|
|
1184 |
"name": "python",
|
1185 |
"nbconvert_exporter": "python",
|
1186 |
"pygments_lexer": "ipython3",
|
1187 |
+
"version": "3.11.2"
|
1188 |
},
|
1189 |
"orig_nbformat": 4
|
1190 |
},
|
app.py
CHANGED
@@ -7,18 +7,14 @@ from phasehunter.data_preparation import prepare_waveform
|
|
7 |
import torch
|
8 |
import io
|
9 |
|
10 |
-
from scipy.signal import resample
|
11 |
from scipy.stats import gaussian_kde
|
|
|
12 |
from bmi_topography import Topography
|
13 |
import earthpy.spatial as es
|
14 |
|
15 |
import obspy
|
16 |
from obspy.clients.fdsn import Client
|
17 |
-
from obspy.clients.fdsn.header import
|
18 |
-
FDSNNoDataException,
|
19 |
-
FDSNTimeoutException,
|
20 |
-
FDSNInternalServerException,
|
21 |
-
)
|
22 |
from obspy.geodetics.base import locations2degrees
|
23 |
from obspy.taup import TauPyModel
|
24 |
from obspy.taup.helper_classes import SlownessModelError
|
@@ -35,12 +31,12 @@ from glob import glob
|
|
35 |
def resample_waveform(waveform, original_freq, target_freq):
|
36 |
"""
|
37 |
Resample a waveform from original frequency to target frequency using SciPy's resample function.
|
38 |
-
|
39 |
Args:
|
40 |
waveform (numpy.ndarray): The input waveform as a 1D array.
|
41 |
original_freq (float): The original sampling frequency of the waveform.
|
42 |
target_freq (float): The target sampling frequency of the waveform.
|
43 |
-
|
44 |
Returns:
|
45 |
resampled_waveform (numpy.ndarray): The resampled waveform as a 1D array.
|
46 |
"""
|
@@ -50,20 +46,19 @@ def resample_waveform(waveform, original_freq, target_freq):
|
|
50 |
resampled_length = int(waveform.shape[-1] * resampling_ratio)
|
51 |
# Resample the waveform using SciPy's resample function
|
52 |
resampled_waveform = resample(waveform, resampled_length, axis=-1)
|
53 |
-
|
54 |
return resampled_waveform
|
55 |
|
56 |
-
|
57 |
def make_prediction(waveform, sampling_rate):
|
58 |
waveform = np.load(waveform)
|
59 |
-
print(
|
60 |
|
61 |
if len(waveform.shape) == 1:
|
62 |
waveform = waveform.reshape(1, waveform.shape[0])
|
63 |
-
print(
|
64 |
if sampling_rate != 100:
|
65 |
waveform = resample_waveform(waveform, sampling_rate, 100)
|
66 |
-
print(
|
67 |
|
68 |
orig_waveform = waveform[:, :6000].copy()
|
69 |
processed_input = prepare_waveform(waveform)
|
@@ -83,80 +78,53 @@ def mark_phases(waveform, uploaded_file, p_thres, s_thres, sampling_rate):
|
|
83 |
if uploaded_file is not None:
|
84 |
waveform = uploaded_file.name
|
85 |
|
86 |
-
processed_input, p_phase, s_phase, orig_waveform = make_prediction(
|
87 |
-
waveform, sampling_rate
|
88 |
-
)
|
89 |
|
90 |
# Create a plot of the waveform with the phases marked
|
91 |
-
if sum(processed_input[0][2] == 0):
|
92 |
fig, ax = plt.subplots(nrows=2, figsize=(10, 2), sharex=True)
|
93 |
|
94 |
-
ax[0].plot(orig_waveform[0], color=
|
95 |
-
ax[0].set_ylabel(
|
96 |
|
97 |
-
else:
|
98 |
fig, ax = plt.subplots(nrows=4, figsize=(10, 6), sharex=True)
|
99 |
-
ax[0].plot(orig_waveform[0], color=
|
100 |
-
ax[1].plot(orig_waveform[1], color=
|
101 |
-
ax[2].plot(orig_waveform[2], color=
|
|
|
|
|
|
|
|
|
102 |
|
103 |
-
ax[0].set_ylabel("Z")
|
104 |
-
ax[1].set_ylabel("N")
|
105 |
-
ax[2].set_ylabel("E")
|
106 |
|
107 |
-
do_we_have_p = p_phase.std().item()
|
108 |
if do_we_have_p:
|
109 |
-
p_phase_plot = p_phase
|
110 |
p_kde = gaussian_kde(p_phase_plot)
|
111 |
-
p_dist_space = np.linspace(min(p_phase_plot)
|
112 |
-
ax[-1].plot(p_dist_space, p_kde(p_dist_space), color=
|
113 |
else:
|
114 |
-
ax[-1].text(
|
115 |
-
|
116 |
-
|
117 |
-
"No P phase detected",
|
118 |
-
horizontalalignment="center",
|
119 |
-
verticalalignment="center",
|
120 |
-
transform=ax[-1].transAxes,
|
121 |
-
)
|
122 |
-
|
123 |
-
do_we_have_s = s_phase.std().item() * 60 < s_thres
|
124 |
if do_we_have_s:
|
125 |
-
s_phase_plot = s_phase
|
126 |
s_kde = gaussian_kde(s_phase_plot)
|
127 |
-
s_dist_space = np.linspace(min(s_phase_plot)
|
128 |
-
ax[-1].plot(s_dist_space, s_kde(s_dist_space), color=
|
129 |
|
130 |
for a in ax:
|
131 |
-
a.axvline(
|
132 |
-
|
133 |
-
color="r",
|
134 |
-
linestyle="--",
|
135 |
-
label="P",
|
136 |
-
alpha=do_we_have_p,
|
137 |
-
)
|
138 |
-
a.axvline(
|
139 |
-
s_phase.mean() * processed_input.shape[-1],
|
140 |
-
color="b",
|
141 |
-
linestyle="--",
|
142 |
-
label="S",
|
143 |
-
alpha=do_we_have_s,
|
144 |
-
)
|
145 |
else:
|
146 |
-
ax[-1].text(
|
147 |
-
|
148 |
-
|
149 |
-
|
150 |
-
horizontalalignment="center",
|
151 |
-
verticalalignment="center",
|
152 |
-
transform=ax[-1].transAxes,
|
153 |
-
)
|
154 |
-
|
155 |
-
ax[-1].set_xlabel("Time, samples")
|
156 |
-
ax[-1].set_ylabel("Uncert., samples")
|
157 |
ax[-1].legend()
|
158 |
|
159 |
-
plt.subplots_adjust(hspace=0
|
160 |
|
161 |
# Convert the plot to an image and return it
|
162 |
fig.canvas.draw()
|
@@ -164,7 +132,6 @@ def mark_phases(waveform, uploaded_file, p_thres, s_thres, sampling_rate):
|
|
164 |
plt.close(fig)
|
165 |
return image
|
166 |
|
167 |
-
|
168 |
def bin_distances(distances, bin_size=10):
|
169 |
# Bin the distances into groups of `bin_size` kilometers
|
170 |
binned_distances = {}
|
@@ -180,10 +147,9 @@ def bin_distances(distances, bin_size=10):
|
|
180 |
for bin_index in binned_distances:
|
181 |
first_distance, first_distance_index = binned_distances[bin_index]
|
182 |
first_distances.append(first_distance_index)
|
183 |
-
|
184 |
return first_distances
|
185 |
|
186 |
-
|
187 |
def variance_coefficient(residuals):
|
188 |
# calculate the variance of the residuals
|
189 |
var = residuals.var()
|
@@ -191,21 +157,9 @@ def variance_coefficient(residuals):
|
|
191 |
coeff = 1 - (var / (residuals.max() - residuals.min()))
|
192 |
return coeff
|
193 |
|
194 |
-
|
195 |
-
def predict_on_section(
|
196 |
-
client_name,
|
197 |
-
timestamp,
|
198 |
-
eq_lat,
|
199 |
-
eq_lon,
|
200 |
-
radius_km,
|
201 |
-
source_depth_km,
|
202 |
-
velocity_model,
|
203 |
-
max_waveforms,
|
204 |
-
conf_thres_P,
|
205 |
-
conf_thres_S,
|
206 |
-
):
|
207 |
distances, t0s, st_lats, st_lons, waveforms, names = [], [], [], [], [], []
|
208 |
-
|
209 |
taup_model = TauPyModel(model=velocity_model)
|
210 |
client = Client(client_name)
|
211 |
|
@@ -219,92 +173,60 @@ def predict_on_section(
|
|
219 |
endtime = starttime + 120
|
220 |
|
221 |
try:
|
222 |
-
print(
|
223 |
-
inv = client.get_stations(
|
224 |
-
|
225 |
-
|
226 |
-
|
227 |
-
|
228 |
-
|
229 |
-
|
230 |
-
|
231 |
-
maxlatitude=(eq_lat + window),
|
232 |
-
minlongitude=(eq_lon - window),
|
233 |
-
maxlongitude=(eq_lon + window),
|
234 |
-
level="station",
|
235 |
-
)
|
236 |
-
print("Finished downloading inventory")
|
237 |
-
|
238 |
-
except (
|
239 |
-
IndexError,
|
240 |
-
FDSNNoDataException,
|
241 |
-
FDSNTimeoutException,
|
242 |
-
FDSNInternalServerException,
|
243 |
-
):
|
244 |
fig, ax = plt.subplots()
|
245 |
-
ax.text(0.5,
|
246 |
-
fig.canvas.draw()
|
247 |
image = np.array(fig.canvas.renderer.buffer_rgba())
|
248 |
plt.close(fig)
|
249 |
return image
|
250 |
-
|
251 |
waveforms = []
|
252 |
cached_waveforms = glob("data/cached/*.mseed")
|
253 |
|
254 |
for network in inv:
|
255 |
-
if network.code ==
|
256 |
continue
|
257 |
for station in network:
|
258 |
print(f"Processing {network.code}.{station.code}...")
|
259 |
-
distance = locations2degrees(
|
260 |
-
eq_lat, eq_lon, station.latitude, station.longitude
|
261 |
-
)
|
262 |
|
263 |
-
arrivals = taup_model.get_travel_times(
|
264 |
-
|
265 |
-
|
266 |
-
phase_list=["P", "S"],
|
267 |
-
)
|
268 |
|
269 |
if len(arrivals) > 0:
|
270 |
|
271 |
starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15
|
272 |
endtime = starttime + 60
|
273 |
try:
|
274 |
-
filename
|
275 |
if f"data/cached/{filename}.mseed" not in cached_waveforms:
|
276 |
-
print(f
|
277 |
-
waveform = client.get_waveforms(
|
278 |
-
|
279 |
-
|
280 |
-
|
281 |
-
channel="*",
|
282 |
-
starttime=starttime,
|
283 |
-
endtime=endtime,
|
284 |
-
)
|
285 |
-
waveform.write(
|
286 |
-
f"data/cached/{network.code}_{station.code}_{starttime}.mseed",
|
287 |
-
format="MSEED",
|
288 |
-
)
|
289 |
-
print("Finished downloading and caching waveform")
|
290 |
else:
|
291 |
-
print(
|
292 |
-
waveform = obspy.read(
|
293 |
-
|
294 |
-
)
|
295 |
|
296 |
-
except (
|
297 |
-
|
298 |
-
FDSNNoDataException,
|
299 |
-
FDSNTimeoutException,
|
300 |
-
FDSNInternalServerException,
|
301 |
-
):
|
302 |
-
print(f"Skipping {network.code}_{station.code}_{starttime}")
|
303 |
continue
|
304 |
-
|
305 |
waveform = waveform.select(channel="H[BH][ZNE]")
|
306 |
waveform = waveform.merge(fill_value=0)
|
307 |
-
waveform = waveform[:3].sort(keys=[
|
308 |
|
309 |
len_check = [len(x.data) for x in waveform]
|
310 |
if len(set(len_check)) > 1:
|
@@ -312,9 +234,7 @@ def predict_on_section(
|
|
312 |
|
313 |
if len(waveform) == 3:
|
314 |
try:
|
315 |
-
waveform = prepare_waveform(
|
316 |
-
np.stack([x.data for x in waveform])
|
317 |
-
)
|
318 |
|
319 |
distances.append(distance)
|
320 |
t0s.append(starttime)
|
@@ -323,32 +243,32 @@ def predict_on_section(
|
|
323 |
waveforms.append(waveform)
|
324 |
names.append(f"{network.code}.{station.code}")
|
325 |
|
326 |
-
print(
|
327 |
-
f"Added {network.code}.{station.code} to the list of waveforms"
|
328 |
-
)
|
329 |
|
330 |
except:
|
331 |
continue
|
332 |
-
|
|
|
333 |
# If there are no waveforms, return an empty plot
|
334 |
if len(waveforms) == 0:
|
335 |
-
print(
|
336 |
fig, ax = plt.subplots()
|
337 |
-
ax.text(0.5,
|
338 |
-
fig.canvas.draw()
|
339 |
image = np.array(fig.canvas.renderer.buffer_rgba())
|
340 |
plt.close(fig)
|
341 |
output_picks = pd.DataFrame()
|
342 |
-
output_picks.to_csv(
|
343 |
-
output_csv =
|
344 |
return image, output_picks, output_csv
|
|
|
345 |
|
346 |
-
first_distances = bin_distances(distances, bin_size=10
|
347 |
|
348 |
# Edge case when there are way too many waveforms to process
|
349 |
-
selection_indexes = np.random.choice(
|
350 |
-
|
351 |
-
|
352 |
|
353 |
waveforms = np.array(waveforms)[selection_indexes]
|
354 |
distances = np.array(distances)[selection_indexes]
|
@@ -359,7 +279,7 @@ def predict_on_section(
|
|
359 |
|
360 |
waveforms = [torch.tensor(waveform) for waveform in waveforms]
|
361 |
|
362 |
-
print(
|
363 |
with torch.no_grad():
|
364 |
waveforms_torch = torch.vstack(waveforms)
|
365 |
output = model(waveforms_torch)
|
@@ -367,183 +287,232 @@ def predict_on_section(
|
|
367 |
p_phases = output[:, 0]
|
368 |
s_phases = output[:, 1]
|
369 |
|
370 |
-
p_phases = p_phases.reshape(len(waveforms)
|
371 |
-
s_phases = s_phases.reshape(len(waveforms)
|
372 |
|
373 |
-
# Max confidence - min variance
|
374 |
p_max_confidence = p_phases.std(axis=-1).min()
|
375 |
s_max_confidence = s_phases.std(axis=-1).min()
|
376 |
|
377 |
print(f"Starting plotting {len(waveforms)} waveforms")
|
378 |
fig, ax = plt.subplots(ncols=3, figsize=(10, 3))
|
379 |
-
|
380 |
# Plot topography
|
381 |
-
print(
|
382 |
params = Topography.DEFAULT.copy()
|
383 |
extra_window = 0.5
|
384 |
-
params["south"] = np.min([st_lats.min(), eq_lat])
|
385 |
-
params["north"] = np.max([st_lats.max(), eq_lat])
|
386 |
-
params["west"] = np.min([st_lons.min(), eq_lon])
|
387 |
-
params["east"] = np.max([st_lons.max(), eq_lon])
|
388 |
|
389 |
topo_map = Topography(**params)
|
390 |
topo_map.fetch()
|
391 |
topo_map.load()
|
392 |
|
393 |
-
print(
|
394 |
hillshade = es.hillshade(topo_map.da[0], altitude=10)
|
395 |
-
|
396 |
-
topo_map.da.plot(ax=ax[1], cmap=
|
397 |
-
topo_map.da.plot(ax=ax[2], cmap=
|
398 |
ax[1].imshow(hillshade, cmap="Greys", alpha=0.5)
|
399 |
|
400 |
-
output_picks = pd.DataFrame(
|
401 |
-
|
402 |
-
|
403 |
-
|
404 |
-
|
405 |
-
|
406 |
-
|
407 |
-
"p_uncertainty, s": [],
|
408 |
-
"s_phase, s": [],
|
409 |
-
"s_uncertainty, s": [],
|
410 |
-
"velocity_p, km/s": [],
|
411 |
-
"velocity_s, km/s": [],
|
412 |
-
}
|
413 |
-
)
|
414 |
-
|
415 |
for i in range(len(waveforms)):
|
416 |
print(f"Plotting waveform {i+1}/{len(waveforms)}")
|
417 |
current_P = p_phases[i]
|
418 |
current_S = s_phases[i]
|
419 |
-
|
420 |
-
x = [t0s[i] + pd.Timedelta(seconds=k
|
421 |
x = mdates.date2num(x)
|
422 |
|
423 |
# Normalize confidence for the plot
|
424 |
-
p_conf = 1
|
425 |
-
s_conf = 1
|
426 |
|
427 |
delta_t = t0s[i].timestamp - obspy.UTCDateTime(timestamp).timestamp
|
428 |
|
429 |
-
ax[0].plot(
|
430 |
-
x,
|
431 |
-
waveforms[i][0, 0] * 10 + distances[i] * 111.2,
|
432 |
-
color="black",
|
433 |
-
alpha=0.5,
|
434 |
-
lw=1,
|
435 |
-
)
|
436 |
-
|
437 |
-
if (current_P.std().item() * 60 < conf_thres_P) or (
|
438 |
-
current_S.std().item() * 60 < conf_thres_S
|
439 |
-
):
|
440 |
-
ax[0].scatter(
|
441 |
-
x[int(current_P.mean() * waveforms[i][0].shape[-1])],
|
442 |
-
waveforms[i][0, 0].mean() + distances[i] * 111.2,
|
443 |
-
color="r",
|
444 |
-
alpha=p_conf,
|
445 |
-
marker="|",
|
446 |
-
)
|
447 |
-
ax[0].scatter(
|
448 |
-
x[int(current_S.mean() * waveforms[i][0].shape[-1])],
|
449 |
-
waveforms[i][0, 0].mean() + distances[i] * 111.2,
|
450 |
-
color="b",
|
451 |
-
alpha=s_conf,
|
452 |
-
marker="|",
|
453 |
-
)
|
454 |
|
455 |
-
|
456 |
-
|
457 |
-
).
|
458 |
-
|
459 |
-
|
460 |
-
).item()
|
461 |
|
462 |
# Generate an array from st_lat to eq_lat and from st_lon to eq_lon
|
463 |
x = np.linspace(st_lons[i], eq_lon, 50)
|
464 |
y = np.linspace(st_lats[i], eq_lat, 50)
|
465 |
-
|
466 |
# Plot the array
|
467 |
-
ax[1].scatter(
|
468 |
-
|
469 |
-
)
|
470 |
-
ax[2].scatter(
|
471 |
-
x, y, c=np.zeros_like(x) + velocity_s, alpha=0.1, vmin=0, vmax=8
|
472 |
-
)
|
473 |
|
474 |
else:
|
475 |
velocity_p = np.nan
|
476 |
velocity_s = np.nan
|
477 |
-
|
478 |
-
ax[0].set_ylabel(
|
479 |
-
print(
|
480 |
-
|
481 |
-
|
482 |
-
|
483 |
-
|
484 |
-
|
485 |
-
|
486 |
-
|
487 |
-
|
488 |
-
|
489 |
-
"starttime": [str(t0s[i])],
|
490 |
-
"p_phase, s": [(delta_t + current_P.mean() * 60).item()],
|
491 |
-
"p_uncertainty, s": [current_P.std().item() * 60],
|
492 |
-
"s_phase, s": [(delta_t + current_S.mean() * 60).item()],
|
493 |
-
"s_uncertainty, s": [current_S.std().item() * 60],
|
494 |
-
"velocity_p, km/s": [velocity_p],
|
495 |
-
"velocity_s, km/s": [velocity_s],
|
496 |
-
}
|
497 |
-
)
|
498 |
-
)
|
499 |
-
|
500 |
# Add legend
|
501 |
-
ax[0].scatter(None, None, color=
|
502 |
-
ax[0].scatter(None, None, color=
|
503 |
-
ax[0].xaxis.set_major_formatter(mdates.DateFormatter(
|
504 |
ax[0].xaxis.set_major_locator(mdates.SecondLocator(interval=20))
|
505 |
ax[0].legend()
|
506 |
|
507 |
-
print(
|
508 |
-
for i in range(1,
|
509 |
-
ax[i].scatter(st_lons, st_lats, color=
|
510 |
-
ax[i].scatter(eq_lon, eq_lat, color=
|
511 |
-
ax[i].set_aspect(
|
512 |
-
ax[i].set_xticklabels(ax[i].get_xticks(), rotation=50)
|
513 |
-
|
514 |
-
fig.subplots_adjust(
|
515 |
-
bottom=0.1, top=0.9, left=0.1, right=0.8, wspace=0.02, hspace=0.02
|
516 |
-
)
|
517 |
|
|
|
|
|
|
|
518 |
cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
|
519 |
-
cbar = fig.colorbar(
|
520 |
-
ax[2].scatter(None, None, c=velocity_p, alpha=0.5, vmin=0, vmax=8), cax=cb_ax
|
521 |
-
)
|
522 |
|
523 |
-
cbar.set_label(
|
524 |
-
ax[1].set_title(
|
525 |
-
ax[2].set_title(
|
526 |
|
527 |
for a in ax:
|
528 |
-
a.tick_params(axis=
|
529 |
-
|
530 |
-
plt.subplots_adjust(hspace=0
|
531 |
-
fig.canvas.draw()
|
532 |
image = np.array(fig.canvas.renderer.buffer_rgba())
|
533 |
plt.close(fig)
|
534 |
-
output_picks.to_csv(
|
535 |
-
f"data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv", index=False
|
536 |
-
)
|
537 |
-
output_csv = f"data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv"
|
538 |
|
|
|
|
|
|
|
539 |
return image, output_picks, output_csv
|
540 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
541 |
|
542 |
model = torch.jit.load("model.pt")
|
543 |
|
544 |
with gr.Blocks() as demo:
|
545 |
-
gr.HTML(
|
546 |
-
"""
|
547 |
<div style="padding: 20px; border-radius: 10px;">
|
548 |
<h1 style="font-size: 30px; text-align: center; margin-bottom: 20px;">PhaseHunter <span style="animation: arrow-anim 10s linear infinite; display: inline-block; transform: rotate(45deg) translateX(-20px);">🏹</span>
|
549 |
|
@@ -571,210 +540,177 @@ with gr.Blocks() as demo:
|
|
571 |
<li>Waveforms should be sampled at 100 samples/sec and have 3 (Z, N, E) or 1 (Z) channels. PhaseHunter analyzes the first 6000 samples of your file.</li>
|
572 |
</ul>
|
573 |
</div>
|
574 |
-
"""
|
575 |
-
)
|
576 |
-
|
577 |
with gr.Tab("Try on a single station"):
|
578 |
-
with gr.Row():
|
579 |
# Define the input and output types for Gradio
|
580 |
inputs = gr.Dropdown(
|
581 |
-
[
|
582 |
-
|
583 |
-
|
584 |
-
|
585 |
-
],
|
586 |
-
label="Sample waveform",
|
587 |
info="Select one of the samples",
|
588 |
-
value="data/sample/sample_0.npy"
|
589 |
)
|
590 |
with gr.Column(scale=1):
|
591 |
-
P_thres_inputs = gr.Slider(
|
592 |
-
|
593 |
-
|
594 |
-
|
595 |
-
|
596 |
-
|
597 |
-
|
598 |
-
|
599 |
-
|
600 |
-
|
601 |
-
|
602 |
-
|
603 |
-
|
604 |
-
|
605 |
-
|
606 |
-
|
607 |
-
|
608 |
-
interactive=True,
|
609 |
-
)
|
610 |
with gr.Column(scale=1):
|
611 |
upload = gr.File(label="Or upload your own waveform")
|
612 |
-
sampling_rate_inputs = gr.Slider(
|
613 |
-
|
614 |
-
|
615 |
-
|
616 |
-
|
617 |
-
|
618 |
-
|
619 |
-
|
620 |
-
)
|
621 |
|
622 |
button = gr.Button("Predict phases")
|
623 |
-
outputs = gr.Image(
|
624 |
-
|
625 |
-
|
626 |
-
|
627 |
-
|
628 |
-
|
629 |
-
inputs=[
|
630 |
-
inputs,
|
631 |
-
upload,
|
632 |
-
P_thres_inputs,
|
633 |
-
S_thres_inputs,
|
634 |
-
sampling_rate_inputs,
|
635 |
-
],
|
636 |
-
outputs=outputs,
|
637 |
-
)
|
638 |
-
|
639 |
with gr.Tab("Select earthquake from catalogue"):
|
640 |
|
641 |
-
gr.HTML(
|
642 |
-
"""
|
643 |
<div style="padding: 20px; border-radius: 10px; font-size: 16px;">
|
644 |
<p style="font-weight: bold; font-size: 24px; margin-bottom: 20px;">Using PhaseHunter to Analyze Seismic Waveforms</p>
|
645 |
<p>Select an earthquake from the global earthquake catalogue (e.g. <a href="https://earthquake.usgs.gov/earthquakes/map">USGS</a>) and the app will download the waveform from the FDSN client of your choice. The app will use a velocity model of your choice to select appropriate time windows for each station within a specified radius of the earthquake.</p>
|
646 |
<p>The app will then analyze the waveforms and mark the detected phases on the waveform. Pick data for each waveform is reported in seconds from the start of the waveform.</p>
|
647 |
<p>Velocities are derived from distance and travel time determined by PhaseHunter picks (<span style="font-style: italic;">v = distance/predicted_pick_time</span>). The background of the velocity plot is colored by DEM.</p>
|
648 |
</div>
|
649 |
-
"""
|
650 |
-
)
|
651 |
-
with gr.Row():
|
652 |
with gr.Column(scale=2):
|
653 |
client_inputs = gr.Dropdown(
|
654 |
-
choices=list(URL_MAPPINGS.keys()),
|
655 |
-
label="FDSN Client",
|
656 |
info="Select one of the available FDSN clients",
|
657 |
-
value="IRIS",
|
658 |
-
interactive=True
|
659 |
)
|
660 |
|
661 |
velocity_inputs = gr.Dropdown(
|
662 |
-
choices=[
|
663 |
-
|
664 |
-
|
665 |
-
|
666 |
-
"ak135f",
|
667 |
-
"herrin",
|
668 |
-
"iasp91",
|
669 |
-
"jb",
|
670 |
-
"prem",
|
671 |
-
"pwdk",
|
672 |
-
],
|
673 |
-
label="1D velocity model",
|
674 |
info="Velocity model for station selection",
|
675 |
-
value="1066a",
|
676 |
-
interactive=True
|
677 |
)
|
678 |
|
679 |
with gr.Column(scale=2):
|
680 |
-
timestamp_inputs = gr.Textbox(
|
681 |
-
|
682 |
-
|
683 |
-
|
684 |
-
|
685 |
-
|
686 |
-
|
687 |
-
|
688 |
-
|
689 |
-
source_depth_inputs = gr.Number(
|
690 |
-
value=10,
|
691 |
label="Source depth (km)",
|
692 |
info="Depth of the earthquake",
|
693 |
-
interactive=True
|
694 |
-
|
695 |
-
|
696 |
with gr.Column(scale=2):
|
697 |
-
eq_lat_inputs = gr.Number(
|
698 |
-
|
699 |
-
|
700 |
-
|
701 |
-
|
702 |
-
|
703 |
-
|
704 |
-
|
705 |
-
|
706 |
-
|
707 |
-
info="Longitude of the earthquake",
|
708 |
-
interactive=True,
|
709 |
-
)
|
710 |
-
|
711 |
with gr.Column(scale=2):
|
712 |
-
radius_inputs = gr.Slider(
|
713 |
-
|
714 |
-
|
715 |
-
|
716 |
-
|
717 |
-
|
718 |
-
info="""Select the radius around the earthquake to download data from.\n
|
719 |
Note that the larger the radius, the longer the app will take to run.""",
|
720 |
-
|
721 |
-
|
722 |
-
|
723 |
-
|
724 |
-
|
725 |
-
|
726 |
-
|
727 |
-
|
728 |
-
|
729 |
-
|
730 |
-
interactive=True,
|
731 |
-
)
|
732 |
with gr.Column(scale=2):
|
733 |
-
P_thres_inputs = gr.Slider(
|
734 |
-
|
735 |
-
|
736 |
-
|
737 |
-
|
738 |
-
|
739 |
-
|
740 |
-
|
741 |
-
|
742 |
-
|
743 |
-
|
744 |
-
|
745 |
-
|
746 |
-
|
747 |
-
|
748 |
-
|
749 |
-
|
750 |
-
)
|
751 |
-
|
752 |
button = gr.Button("Predict phases")
|
753 |
-
output_image = gr.Image(
|
754 |
-
label="Waveforms with Phases Marked", type="numpy", interactive=False
|
755 |
-
)
|
756 |
|
757 |
with gr.Row():
|
758 |
-
output_picks = gr.Dataframe(
|
759 |
-
|
760 |
-
|
761 |
output_csv = gr.File(label="Output File", file_types=[".csv"])
|
762 |
|
763 |
-
button.click(
|
764 |
-
|
765 |
-
|
766 |
-
|
767 |
-
|
768 |
-
|
769 |
-
|
770 |
-
|
771 |
-
|
772 |
-
|
773 |
-
|
774 |
-
|
775 |
-
|
776 |
-
|
777 |
-
|
778 |
-
|
779 |
-
|
780 |
-
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
7 |
import torch
|
8 |
import io
|
9 |
|
|
|
10 |
from scipy.stats import gaussian_kde
|
11 |
+
from scipy.signal import resample
|
12 |
from bmi_topography import Topography
|
13 |
import earthpy.spatial as es
|
14 |
|
15 |
import obspy
|
16 |
from obspy.clients.fdsn import Client
|
17 |
+
from obspy.clients.fdsn.header import FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException
|
|
|
|
|
|
|
|
|
18 |
from obspy.geodetics.base import locations2degrees
|
19 |
from obspy.taup import TauPyModel
|
20 |
from obspy.taup.helper_classes import SlownessModelError
|
|
|
31 |
def resample_waveform(waveform, original_freq, target_freq):
|
32 |
"""
|
33 |
Resample a waveform from original frequency to target frequency using SciPy's resample function.
|
34 |
+
|
35 |
Args:
|
36 |
waveform (numpy.ndarray): The input waveform as a 1D array.
|
37 |
original_freq (float): The original sampling frequency of the waveform.
|
38 |
target_freq (float): The target sampling frequency of the waveform.
|
39 |
+
|
40 |
Returns:
|
41 |
resampled_waveform (numpy.ndarray): The resampled waveform as a 1D array.
|
42 |
"""
|
|
|
46 |
resampled_length = int(waveform.shape[-1] * resampling_ratio)
|
47 |
# Resample the waveform using SciPy's resample function
|
48 |
resampled_waveform = resample(waveform, resampled_length, axis=-1)
|
49 |
+
|
50 |
return resampled_waveform
|
51 |
|
|
|
52 |
def make_prediction(waveform, sampling_rate):
|
53 |
waveform = np.load(waveform)
|
54 |
+
print('Loaded', waveform.shape)
|
55 |
|
56 |
if len(waveform.shape) == 1:
|
57 |
waveform = waveform.reshape(1, waveform.shape[0])
|
58 |
+
print('Reshaped', waveform.shape)
|
59 |
if sampling_rate != 100:
|
60 |
waveform = resample_waveform(waveform, sampling_rate, 100)
|
61 |
+
print('Resampled', waveform.shape)
|
62 |
|
63 |
orig_waveform = waveform[:, :6000].copy()
|
64 |
processed_input = prepare_waveform(waveform)
|
|
|
78 |
if uploaded_file is not None:
|
79 |
waveform = uploaded_file.name
|
80 |
|
81 |
+
processed_input, p_phase, s_phase, orig_waveform = make_prediction(waveform, sampling_rate)
|
|
|
|
|
82 |
|
83 |
# Create a plot of the waveform with the phases marked
|
84 |
+
if sum(processed_input[0][2] == 0): #if input is 1C
|
85 |
fig, ax = plt.subplots(nrows=2, figsize=(10, 2), sharex=True)
|
86 |
|
87 |
+
ax[0].plot(orig_waveform[0], color='black', lw=1)
|
88 |
+
ax[0].set_ylabel('Norm. Ampl.')
|
89 |
|
90 |
+
else: #if input is 3C
|
91 |
fig, ax = plt.subplots(nrows=4, figsize=(10, 6), sharex=True)
|
92 |
+
ax[0].plot(orig_waveform[0], color='black', lw=1)
|
93 |
+
ax[1].plot(orig_waveform[1], color='black', lw=1)
|
94 |
+
ax[2].plot(orig_waveform[2], color='black', lw=1)
|
95 |
+
|
96 |
+
ax[0].set_ylabel('Z')
|
97 |
+
ax[1].set_ylabel('N')
|
98 |
+
ax[2].set_ylabel('E')
|
99 |
|
|
|
|
|
|
|
100 |
|
101 |
+
do_we_have_p = (p_phase.std().item()*60 < p_thres)
|
102 |
if do_we_have_p:
|
103 |
+
p_phase_plot = p_phase*processed_input.shape[-1]
|
104 |
p_kde = gaussian_kde(p_phase_plot)
|
105 |
+
p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )
|
106 |
+
ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')
|
107 |
else:
|
108 |
+
ax[-1].text(0.5, 0.75, 'No P phase detected', horizontalalignment='center', verticalalignment='center', transform=ax[-1].transAxes)
|
109 |
+
|
110 |
+
do_we_have_s = (s_phase.std().item()*60 < s_thres)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
111 |
if do_we_have_s:
|
112 |
+
s_phase_plot = s_phase*processed_input.shape[-1]
|
113 |
s_kde = gaussian_kde(s_phase_plot)
|
114 |
+
s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )
|
115 |
+
ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')
|
116 |
|
117 |
for a in ax:
|
118 |
+
a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P', alpha=do_we_have_p)
|
119 |
+
a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S', alpha=do_we_have_s)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
120 |
else:
|
121 |
+
ax[-1].text(0.5, 0.25, 'No S phase detected', horizontalalignment='center', verticalalignment='center', transform=ax[-1].transAxes)
|
122 |
+
|
123 |
+
ax[-1].set_xlabel('Time, samples')
|
124 |
+
ax[-1].set_ylabel('Uncert., samples')
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
125 |
ax[-1].legend()
|
126 |
|
127 |
+
plt.subplots_adjust(hspace=0., wspace=0.)
|
128 |
|
129 |
# Convert the plot to an image and return it
|
130 |
fig.canvas.draw()
|
|
|
132 |
plt.close(fig)
|
133 |
return image
|
134 |
|
|
|
135 |
def bin_distances(distances, bin_size=10):
|
136 |
# Bin the distances into groups of `bin_size` kilometers
|
137 |
binned_distances = {}
|
|
|
147 |
for bin_index in binned_distances:
|
148 |
first_distance, first_distance_index = binned_distances[bin_index]
|
149 |
first_distances.append(first_distance_index)
|
150 |
+
|
151 |
return first_distances
|
152 |
|
|
|
153 |
def variance_coefficient(residuals):
|
154 |
# calculate the variance of the residuals
|
155 |
var = residuals.var()
|
|
|
157 |
coeff = 1 - (var / (residuals.max() - residuals.min()))
|
158 |
return coeff
|
159 |
|
160 |
+
def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source_depth_km, velocity_model, max_waveforms, conf_thres_P, conf_thres_S):
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
161 |
distances, t0s, st_lats, st_lons, waveforms, names = [], [], [], [], [], []
|
162 |
+
|
163 |
taup_model = TauPyModel(model=velocity_model)
|
164 |
client = Client(client_name)
|
165 |
|
|
|
173 |
endtime = starttime + 120
|
174 |
|
175 |
try:
|
176 |
+
print('Starting to download inventory')
|
177 |
+
inv = client.get_stations(network="*", station="*", location="*", channel="*H*",
|
178 |
+
starttime=starttime, endtime=endtime,
|
179 |
+
minlatitude=(eq_lat-window), maxlatitude=(eq_lat+window),
|
180 |
+
minlongitude=(eq_lon-window), maxlongitude=(eq_lon+window),
|
181 |
+
level='station')
|
182 |
+
print('Finished downloading inventory')
|
183 |
+
|
184 |
+
except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
185 |
fig, ax = plt.subplots()
|
186 |
+
ax.text(0.5,0.5,'Something is wrong with the data provider, try another')
|
187 |
+
fig.canvas.draw();
|
188 |
image = np.array(fig.canvas.renderer.buffer_rgba())
|
189 |
plt.close(fig)
|
190 |
return image
|
191 |
+
|
192 |
waveforms = []
|
193 |
cached_waveforms = glob("data/cached/*.mseed")
|
194 |
|
195 |
for network in inv:
|
196 |
+
if network.code == 'SY':
|
197 |
continue
|
198 |
for station in network:
|
199 |
print(f"Processing {network.code}.{station.code}...")
|
200 |
+
distance = locations2degrees(eq_lat, eq_lon, station.latitude, station.longitude)
|
|
|
|
|
201 |
|
202 |
+
arrivals = taup_model.get_travel_times(source_depth_in_km=source_depth_km,
|
203 |
+
distance_in_degree=distance,
|
204 |
+
phase_list=["P", "S"])
|
|
|
|
|
205 |
|
206 |
if len(arrivals) > 0:
|
207 |
|
208 |
starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15
|
209 |
endtime = starttime + 60
|
210 |
try:
|
211 |
+
filename=f'{network.code}_{station.code}_{starttime}'
|
212 |
if f"data/cached/{filename}.mseed" not in cached_waveforms:
|
213 |
+
print(f'Downloading waveform for {filename}')
|
214 |
+
waveform = client.get_waveforms(network=network.code, station=station.code, location="*", channel="*",
|
215 |
+
starttime=starttime, endtime=endtime)
|
216 |
+
waveform.write(f"data/cached/{network.code}_{station.code}_{starttime}.mseed", format="MSEED")
|
217 |
+
print('Finished downloading and caching waveform')
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
218 |
else:
|
219 |
+
print('Reading cached waveform')
|
220 |
+
waveform = obspy.read(f"data/cached/{network.code}_{station.code}_{starttime}.mseed")
|
221 |
+
|
|
|
222 |
|
223 |
+
except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):
|
224 |
+
print(f'Skipping {network.code}_{station.code}_{starttime}')
|
|
|
|
|
|
|
|
|
|
|
225 |
continue
|
226 |
+
|
227 |
waveform = waveform.select(channel="H[BH][ZNE]")
|
228 |
waveform = waveform.merge(fill_value=0)
|
229 |
+
waveform = waveform[:3].sort(keys=['channel'], reverse=True)
|
230 |
|
231 |
len_check = [len(x.data) for x in waveform]
|
232 |
if len(set(len_check)) > 1:
|
|
|
234 |
|
235 |
if len(waveform) == 3:
|
236 |
try:
|
237 |
+
waveform = prepare_waveform(np.stack([x.data for x in waveform]))
|
|
|
|
|
238 |
|
239 |
distances.append(distance)
|
240 |
t0s.append(starttime)
|
|
|
243 |
waveforms.append(waveform)
|
244 |
names.append(f"{network.code}.{station.code}")
|
245 |
|
246 |
+
print(f"Added {network.code}.{station.code} to the list of waveforms")
|
|
|
|
|
247 |
|
248 |
except:
|
249 |
continue
|
250 |
+
|
251 |
+
|
252 |
# If there are no waveforms, return an empty plot
|
253 |
if len(waveforms) == 0:
|
254 |
+
print('No waveforms found')
|
255 |
fig, ax = plt.subplots()
|
256 |
+
ax.text(0.5,0.5,'No waveforms found')
|
257 |
+
fig.canvas.draw();
|
258 |
image = np.array(fig.canvas.renderer.buffer_rgba())
|
259 |
plt.close(fig)
|
260 |
output_picks = pd.DataFrame()
|
261 |
+
output_picks.to_csv('data/picks.csv', index=False)
|
262 |
+
output_csv = 'data/picks.csv'
|
263 |
return image, output_picks, output_csv
|
264 |
+
|
265 |
|
266 |
+
first_distances = bin_distances(distances, bin_size=10/111.2)
|
267 |
|
268 |
# Edge case when there are way too many waveforms to process
|
269 |
+
selection_indexes = np.random.choice(first_distances,
|
270 |
+
np.min([len(first_distances), max_waveforms]),
|
271 |
+
replace=False)
|
272 |
|
273 |
waveforms = np.array(waveforms)[selection_indexes]
|
274 |
distances = np.array(distances)[selection_indexes]
|
|
|
279 |
|
280 |
waveforms = [torch.tensor(waveform) for waveform in waveforms]
|
281 |
|
282 |
+
print('Starting to run predictions')
|
283 |
with torch.no_grad():
|
284 |
waveforms_torch = torch.vstack(waveforms)
|
285 |
output = model(waveforms_torch)
|
|
|
287 |
p_phases = output[:, 0]
|
288 |
s_phases = output[:, 1]
|
289 |
|
290 |
+
p_phases = p_phases.reshape(len(waveforms),-1)
|
291 |
+
s_phases = s_phases.reshape(len(waveforms),-1)
|
292 |
|
293 |
+
# Max confidence - min variance
|
294 |
p_max_confidence = p_phases.std(axis=-1).min()
|
295 |
s_max_confidence = s_phases.std(axis=-1).min()
|
296 |
|
297 |
print(f"Starting plotting {len(waveforms)} waveforms")
|
298 |
fig, ax = plt.subplots(ncols=3, figsize=(10, 3))
|
299 |
+
|
300 |
# Plot topography
|
301 |
+
print('Fetching topography')
|
302 |
params = Topography.DEFAULT.copy()
|
303 |
extra_window = 0.5
|
304 |
+
params["south"] = np.min([st_lats.min(), eq_lat])-extra_window
|
305 |
+
params["north"] = np.max([st_lats.max(), eq_lat])+extra_window
|
306 |
+
params["west"] = np.min([st_lons.min(), eq_lon])-extra_window
|
307 |
+
params["east"] = np.max([st_lons.max(), eq_lon])+extra_window
|
308 |
|
309 |
topo_map = Topography(**params)
|
310 |
topo_map.fetch()
|
311 |
topo_map.load()
|
312 |
|
313 |
+
print('Plotting topo')
|
314 |
hillshade = es.hillshade(topo_map.da[0], altitude=10)
|
315 |
+
|
316 |
+
topo_map.da.plot(ax = ax[1], cmap='Greys', add_colorbar=False, add_labels=False)
|
317 |
+
topo_map.da.plot(ax = ax[2], cmap='Greys', add_colorbar=False, add_labels=False)
|
318 |
ax[1].imshow(hillshade, cmap="Greys", alpha=0.5)
|
319 |
|
320 |
+
output_picks = pd.DataFrame({'station_name' : [],
|
321 |
+
'st_lat' : [], 'st_lon' : [],
|
322 |
+
'starttime' : [],
|
323 |
+
'p_phase, s' : [], 'p_uncertainty, s' : [],
|
324 |
+
's_phase, s' : [], 's_uncertainty, s' : [],
|
325 |
+
'velocity_p, km/s' : [], 'velocity_s, km/s' : []})
|
326 |
+
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
327 |
for i in range(len(waveforms)):
|
328 |
print(f"Plotting waveform {i+1}/{len(waveforms)}")
|
329 |
current_P = p_phases[i]
|
330 |
current_S = s_phases[i]
|
331 |
+
|
332 |
+
x = [t0s[i] + pd.Timedelta(seconds=k/100) for k in np.linspace(0,6000,6000)]
|
333 |
x = mdates.date2num(x)
|
334 |
|
335 |
# Normalize confidence for the plot
|
336 |
+
p_conf = 1/(current_P.std()/p_max_confidence).item()
|
337 |
+
s_conf = 1/(current_S.std()/s_max_confidence).item()
|
338 |
|
339 |
delta_t = t0s[i].timestamp - obspy.UTCDateTime(timestamp).timestamp
|
340 |
|
341 |
+
ax[0].plot(x, waveforms[i][0, 0]*10+distances[i]*111.2, color='black', alpha=0.5, lw=1)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
342 |
|
343 |
+
if (current_P.std().item()*60 < conf_thres_P) or (current_S.std().item()*60 < conf_thres_S):
|
344 |
+
ax[0].scatter(x[int(current_P.mean()*waveforms[i][0].shape[-1])], waveforms[i][0, 0].mean()+distances[i]*111.2, color='r', alpha=p_conf, marker='|')
|
345 |
+
ax[0].scatter(x[int(current_S.mean()*waveforms[i][0].shape[-1])], waveforms[i][0, 0].mean()+distances[i]*111.2, color='b', alpha=s_conf, marker='|')
|
346 |
+
|
347 |
+
velocity_p = (distances[i]*111.2)/(delta_t+current_P.mean()*60).item()
|
348 |
+
velocity_s = (distances[i]*111.2)/(delta_t+current_S.mean()*60).item()
|
349 |
|
350 |
# Generate an array from st_lat to eq_lat and from st_lon to eq_lon
|
351 |
x = np.linspace(st_lons[i], eq_lon, 50)
|
352 |
y = np.linspace(st_lats[i], eq_lat, 50)
|
353 |
+
|
354 |
# Plot the array
|
355 |
+
ax[1].scatter(x, y, c=np.zeros_like(x)+velocity_p, alpha=0.1, vmin=0, vmax=8)
|
356 |
+
ax[2].scatter(x, y, c=np.zeros_like(x)+velocity_s, alpha=0.1, vmin=0, vmax=8)
|
|
|
|
|
|
|
|
|
357 |
|
358 |
else:
|
359 |
velocity_p = np.nan
|
360 |
velocity_s = np.nan
|
361 |
+
|
362 |
+
ax[0].set_ylabel('Z')
|
363 |
+
print(f"Station {st_lats[i]}, {st_lons[i]} has P velocity {velocity_p} and S velocity {velocity_s}")
|
364 |
+
|
365 |
+
output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]],
|
366 |
+
'st_lat' : [st_lats[i]], 'st_lon' : [st_lons[i]],
|
367 |
+
'starttime' : [str(t0s[i])],
|
368 |
+
'p_phase, s' : [(delta_t+current_P.mean()*60).item()], 'p_uncertainty, s' : [current_P.std().item()*60],
|
369 |
+
's_phase, s' : [(delta_t+current_S.mean()*60).item()], 's_uncertainty, s' : [current_S.std().item()*60],
|
370 |
+
'velocity_p, km/s' : [velocity_p], 'velocity_s, km/s' : [velocity_s]}))
|
371 |
+
|
372 |
+
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
373 |
# Add legend
|
374 |
+
ax[0].scatter(None, None, color='r', marker='|', label='P')
|
375 |
+
ax[0].scatter(None, None, color='b', marker='|', label='S')
|
376 |
+
ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))
|
377 |
ax[0].xaxis.set_major_locator(mdates.SecondLocator(interval=20))
|
378 |
ax[0].legend()
|
379 |
|
380 |
+
print('Plotting stations')
|
381 |
+
for i in range(1,3):
|
382 |
+
ax[i].scatter(st_lons, st_lats, color='b', label='Stations')
|
383 |
+
ax[i].scatter(eq_lon, eq_lat, color='r', marker='*', label='Earthquake')
|
384 |
+
ax[i].set_aspect('equal')
|
385 |
+
ax[i].set_xticklabels(ax[i].get_xticks(), rotation = 50)
|
|
|
|
|
|
|
|
|
386 |
|
387 |
+
fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
|
388 |
+
wspace=0.02, hspace=0.02)
|
389 |
+
|
390 |
cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
|
391 |
+
cbar = fig.colorbar(ax[2].scatter(None, None, c=velocity_p, alpha=0.5, vmin=0, vmax=8), cax=cb_ax)
|
|
|
|
|
392 |
|
393 |
+
cbar.set_label('Velocity (km/s)')
|
394 |
+
ax[1].set_title('P Velocity')
|
395 |
+
ax[2].set_title('S Velocity')
|
396 |
|
397 |
for a in ax:
|
398 |
+
a.tick_params(axis='both', which='major', labelsize=8)
|
399 |
+
|
400 |
+
plt.subplots_adjust(hspace=0., wspace=0.5)
|
401 |
+
fig.canvas.draw();
|
402 |
image = np.array(fig.canvas.renderer.buffer_rgba())
|
403 |
plt.close(fig)
|
|
|
|
|
|
|
|
|
404 |
|
405 |
+
output_csv = f'data/velocity/{eq_lat}_{eq_lon}_{source_depth_km}_{timestamp}_{len(waveforms)}.csv'
|
406 |
+
output_picks.to_csv(output_csv, index=False)
|
407 |
+
|
408 |
return image, output_picks, output_csv
|
409 |
|
410 |
+
import numpy as np
|
411 |
+
from matplotlib import colors, cm
|
412 |
+
|
413 |
+
# Function to find the closest index for a given value in an array
|
414 |
+
def find_closest_index(array, value):
|
415 |
+
return np.argmin(np.abs(array - value))
|
416 |
+
|
417 |
+
def compute_velocity_model(azimuth, elevation):
|
418 |
+
filename = list(output_csv.temp_files)[0]
|
419 |
+
|
420 |
+
df = pd.read_csv(filename)
|
421 |
+
print(df)
|
422 |
+
filename = filename.split('/')[-1]
|
423 |
+
|
424 |
+
# Current EQ location
|
425 |
+
eq_lat = float(filename.split("_")[0])
|
426 |
+
eq_lon = float(filename.split("_")[1])
|
427 |
+
eq_depth = float(filename.split("_")[2])
|
428 |
+
|
429 |
+
# Define the region of interest (latitude, longitude, and depth ranges)
|
430 |
+
lat_range = (np.min([df.st_lat.min(), eq_lat]), np.max([df.st_lat.max(), eq_lat]))
|
431 |
+
lon_range = (np.min([df.st_lon.min(), eq_lon]), np.max([df.st_lon.max(), eq_lon]))
|
432 |
+
depth_range = (0, 50)
|
433 |
+
|
434 |
+
# Define the number of nodes in each dimension
|
435 |
+
n_lat = 10
|
436 |
+
n_lon = 10
|
437 |
+
n_depth = 10
|
438 |
+
num_points = 100
|
439 |
+
|
440 |
+
# Create the grid
|
441 |
+
lat_values = np.linspace(lat_range[0], lat_range[1], n_lat)
|
442 |
+
lon_values = np.linspace(lon_range[0], lon_range[1], n_lon)
|
443 |
+
depth_values = np.linspace(depth_range[0], depth_range[1], n_depth)
|
444 |
+
|
445 |
+
# Initialize the velocity model with constant values
|
446 |
+
initial_velocity = 0 # km/s, this can be P-wave or S-wave velocity
|
447 |
+
velocity_model = np.full((n_lat, n_lon, n_depth), initial_velocity, dtype=float)
|
448 |
+
|
449 |
+
# Loop through the stations and update the velocity model
|
450 |
+
for i in range(len(df)):
|
451 |
+
if ~np.isnan(df['velocity_p, km/s'].iloc[i]):
|
452 |
+
# Interpolate coordinates along the great circle path between the earthquake and the station
|
453 |
+
lon_deg = np.linspace(df.st_lon.iloc[i], eq_lon, num_points)
|
454 |
+
lat_deg = np.linspace(df.st_lat.iloc[i], eq_lat, num_points)
|
455 |
+
depth_interpolated = np.interp(np.linspace(0, 1, num_points), [0, 1], [eq_depth, 0])
|
456 |
+
|
457 |
+
# Loop through the interpolated coordinates and update the grid cells with the average P-wave velocity
|
458 |
+
for lat, lon, depth in zip(lat_deg, lon_deg, depth_interpolated):
|
459 |
+
lat_index = find_closest_index(lat_values, lat)
|
460 |
+
lon_index = find_closest_index(lon_values, lon)
|
461 |
+
depth_index = find_closest_index(depth_values, depth)
|
462 |
+
|
463 |
+
if velocity_model[lat_index, lon_index, depth_index] == initial_velocity:
|
464 |
+
velocity_model[lat_index, lon_index, depth_index] = df['velocity_p, km/s'].iloc[i]
|
465 |
+
else:
|
466 |
+
velocity_model[lat_index, lon_index, depth_index] = (velocity_model[lat_index, lon_index, depth_index] +
|
467 |
+
df['velocity_p, km/s'].iloc[i]) / 2
|
468 |
+
|
469 |
+
# Create the figure and axis
|
470 |
+
fig = plt.figure(figsize=(8, 8))
|
471 |
+
ax = fig.add_subplot(111, projection='3d')
|
472 |
+
|
473 |
+
# Set the plot limits
|
474 |
+
ax.set_xlim3d(lat_range[0], lat_range[1])
|
475 |
+
ax.set_ylim3d(lon_range[0], lon_range[1])
|
476 |
+
ax.set_zlim3d(depth_range[1], depth_range[0])
|
477 |
+
|
478 |
+
ax.set_xlabel('Latitude')
|
479 |
+
ax.set_ylabel('Longitude')
|
480 |
+
ax.set_zlabel('Depth (km)')
|
481 |
+
ax.set_title('Velocity Model')
|
482 |
+
|
483 |
+
# Create the meshgrid
|
484 |
+
x, y, z = np.meshgrid(
|
485 |
+
np.linspace(lat_range[0], lat_range[1], velocity_model.shape[0]+1),
|
486 |
+
np.linspace(lon_range[0], lon_range[1], velocity_model.shape[1]+1),
|
487 |
+
np.linspace(depth_range[0], depth_range[1], velocity_model.shape[2]+1),
|
488 |
+
indexing='ij'
|
489 |
+
)
|
490 |
+
|
491 |
+
# Create the color array
|
492 |
+
norm = plt.Normalize(vmin=velocity_model.min(), vmax=velocity_model.max())
|
493 |
+
colors = plt.cm.plasma(norm(velocity_model))
|
494 |
+
|
495 |
+
# Plot the voxels
|
496 |
+
ax.voxels(x, y, z, velocity_model > 0, facecolors=colors, alpha=0.5, edgecolor='k')
|
497 |
+
|
498 |
+
# Set the view angle
|
499 |
+
ax.view_init(elev=elevation, azim=azimuth)
|
500 |
+
|
501 |
+
m = cm.ScalarMappable(cmap=plt.cm.plasma, norm=norm)
|
502 |
+
m.set_array([])
|
503 |
+
plt.colorbar(m)
|
504 |
+
|
505 |
+
# Show the plot
|
506 |
+
fig.canvas.draw();
|
507 |
+
image = np.array(fig.canvas.renderer.buffer_rgba())
|
508 |
+
plt.close(fig)
|
509 |
+
|
510 |
+
return image
|
511 |
|
512 |
model = torch.jit.load("model.pt")
|
513 |
|
514 |
with gr.Blocks() as demo:
|
515 |
+
gr.HTML("""
|
|
|
516 |
<div style="padding: 20px; border-radius: 10px;">
|
517 |
<h1 style="font-size: 30px; text-align: center; margin-bottom: 20px;">PhaseHunter <span style="animation: arrow-anim 10s linear infinite; display: inline-block; transform: rotate(45deg) translateX(-20px);">🏹</span>
|
518 |
|
|
|
540 |
<li>Waveforms should be sampled at 100 samples/sec and have 3 (Z, N, E) or 1 (Z) channels. PhaseHunter analyzes the first 6000 samples of your file.</li>
|
541 |
</ul>
|
542 |
</div>
|
543 |
+
""")
|
|
|
|
|
544 |
with gr.Tab("Try on a single station"):
|
545 |
+
with gr.Row():
|
546 |
# Define the input and output types for Gradio
|
547 |
inputs = gr.Dropdown(
|
548 |
+
["data/sample/sample_0.npy",
|
549 |
+
"data/sample/sample_1.npy",
|
550 |
+
"data/sample/sample_2.npy"],
|
551 |
+
label="Sample waveform",
|
|
|
|
|
552 |
info="Select one of the samples",
|
553 |
+
value = "data/sample/sample_0.npy"
|
554 |
)
|
555 |
with gr.Column(scale=1):
|
556 |
+
P_thres_inputs = gr.Slider(minimum=0.01,
|
557 |
+
maximum=1,
|
558 |
+
value=0.1,
|
559 |
+
label="P uncertainty threshold, s",
|
560 |
+
step=0.01,
|
561 |
+
info="Acceptable uncertainty for P picks expressed in std() seconds",
|
562 |
+
interactive=True,
|
563 |
+
)
|
564 |
+
|
565 |
+
S_thres_inputs = gr.Slider(minimum=0.01,
|
566 |
+
maximum=1,
|
567 |
+
value=0.2,
|
568 |
+
label="S uncertainty threshold, s",
|
569 |
+
step=0.01,
|
570 |
+
info="Acceptable uncertainty for S picks expressed in std() seconds",
|
571 |
+
interactive=True,
|
572 |
+
)
|
|
|
|
|
573 |
with gr.Column(scale=1):
|
574 |
upload = gr.File(label="Or upload your own waveform")
|
575 |
+
sampling_rate_inputs = gr.Slider(minimum=10,
|
576 |
+
maximum=1000,
|
577 |
+
value=100,
|
578 |
+
label="Samlping rate, Hz",
|
579 |
+
step=10,
|
580 |
+
info="Sampling rate of the waveform",
|
581 |
+
interactive=True,
|
582 |
+
)
|
|
|
583 |
|
584 |
button = gr.Button("Predict phases")
|
585 |
+
outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)
|
586 |
+
|
587 |
+
button.click(mark_phases, inputs=[inputs, upload,
|
588 |
+
P_thres_inputs, S_thres_inputs,
|
589 |
+
sampling_rate_inputs],
|
590 |
+
outputs=outputs)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
591 |
with gr.Tab("Select earthquake from catalogue"):
|
592 |
|
593 |
+
gr.HTML("""
|
|
|
594 |
<div style="padding: 20px; border-radius: 10px; font-size: 16px;">
|
595 |
<p style="font-weight: bold; font-size: 24px; margin-bottom: 20px;">Using PhaseHunter to Analyze Seismic Waveforms</p>
|
596 |
<p>Select an earthquake from the global earthquake catalogue (e.g. <a href="https://earthquake.usgs.gov/earthquakes/map">USGS</a>) and the app will download the waveform from the FDSN client of your choice. The app will use a velocity model of your choice to select appropriate time windows for each station within a specified radius of the earthquake.</p>
|
597 |
<p>The app will then analyze the waveforms and mark the detected phases on the waveform. Pick data for each waveform is reported in seconds from the start of the waveform.</p>
|
598 |
<p>Velocities are derived from distance and travel time determined by PhaseHunter picks (<span style="font-style: italic;">v = distance/predicted_pick_time</span>). The background of the velocity plot is colored by DEM.</p>
|
599 |
</div>
|
600 |
+
""")
|
601 |
+
with gr.Row():
|
|
|
602 |
with gr.Column(scale=2):
|
603 |
client_inputs = gr.Dropdown(
|
604 |
+
choices = list(URL_MAPPINGS.keys()),
|
605 |
+
label="FDSN Client",
|
606 |
info="Select one of the available FDSN clients",
|
607 |
+
value = "IRIS",
|
608 |
+
interactive=True
|
609 |
)
|
610 |
|
611 |
velocity_inputs = gr.Dropdown(
|
612 |
+
choices = ['1066a', '1066b', 'ak135',
|
613 |
+
'ak135f', 'herrin', 'iasp91',
|
614 |
+
'jb', 'prem', 'pwdk'],
|
615 |
+
label="1D velocity model",
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
616 |
info="Velocity model for station selection",
|
617 |
+
value = "1066a",
|
618 |
+
interactive=True
|
619 |
)
|
620 |
|
621 |
with gr.Column(scale=2):
|
622 |
+
timestamp_inputs = gr.Textbox(value='2019-07-04 17:33:49',
|
623 |
+
placeholder='YYYY-MM-DD HH:MM:SS',
|
624 |
+
label="Timestamp",
|
625 |
+
info="Timestamp of the earthquake",
|
626 |
+
max_lines=1,
|
627 |
+
interactive=True)
|
628 |
+
|
629 |
+
source_depth_inputs = gr.Number(value=10,
|
|
|
|
|
|
|
630 |
label="Source depth (km)",
|
631 |
info="Depth of the earthquake",
|
632 |
+
interactive=True)
|
633 |
+
|
|
|
634 |
with gr.Column(scale=2):
|
635 |
+
eq_lat_inputs = gr.Number(value=35.766,
|
636 |
+
label="Latitude",
|
637 |
+
info="Latitude of the earthquake",
|
638 |
+
interactive=True)
|
639 |
+
|
640 |
+
eq_lon_inputs = gr.Number(value=-117.605,
|
641 |
+
label="Longitude",
|
642 |
+
info="Longitude of the earthquake",
|
643 |
+
interactive=True)
|
644 |
+
|
|
|
|
|
|
|
|
|
645 |
with gr.Column(scale=2):
|
646 |
+
radius_inputs = gr.Slider(minimum=1,
|
647 |
+
maximum=200,
|
648 |
+
value=50,
|
649 |
+
label="Radius (km)",
|
650 |
+
step=10,
|
651 |
+
info="""Select the radius around the earthquake to download data from.\n
|
|
|
652 |
Note that the larger the radius, the longer the app will take to run.""",
|
653 |
+
interactive=True)
|
654 |
+
|
655 |
+
max_waveforms_inputs = gr.Slider(minimum=1,
|
656 |
+
maximum=100,
|
657 |
+
value=10,
|
658 |
+
label="Max waveforms per section",
|
659 |
+
step=1,
|
660 |
+
info="Maximum number of waveforms to show per section\n (to avoid long prediction times)",
|
661 |
+
interactive=True,
|
662 |
+
)
|
|
|
|
|
663 |
with gr.Column(scale=2):
|
664 |
+
P_thres_inputs = gr.Slider(minimum=0.01,
|
665 |
+
maximum=1,
|
666 |
+
value=0.1,
|
667 |
+
label="P uncertainty threshold, s",
|
668 |
+
step=0.01,
|
669 |
+
info="Acceptable uncertainty for P picks expressed in std() seconds",
|
670 |
+
interactive=True,
|
671 |
+
)
|
672 |
+
S_thres_inputs = gr.Slider(minimum=0.01,
|
673 |
+
maximum=1,
|
674 |
+
value=0.2,
|
675 |
+
label="S uncertainty threshold, s",
|
676 |
+
step=0.01,
|
677 |
+
info="Acceptable uncertainty for S picks expressed in std() seconds",
|
678 |
+
interactive=True,
|
679 |
+
)
|
680 |
+
|
|
|
|
|
681 |
button = gr.Button("Predict phases")
|
682 |
+
output_image = gr.Image(label='Waveforms with Phases Marked', type='numpy', interactive=False)
|
|
|
|
|
683 |
|
684 |
with gr.Row():
|
685 |
+
output_picks = gr.Dataframe(label='Pick data',
|
686 |
+
type='pandas',
|
687 |
+
interactive=False)
|
688 |
output_csv = gr.File(label="Output File", file_types=[".csv"])
|
689 |
|
690 |
+
button.click(predict_on_section,
|
691 |
+
inputs=[client_inputs, timestamp_inputs,
|
692 |
+
eq_lat_inputs, eq_lon_inputs,
|
693 |
+
radius_inputs, source_depth_inputs,
|
694 |
+
velocity_inputs, max_waveforms_inputs,
|
695 |
+
P_thres_inputs, S_thres_inputs],
|
696 |
+
outputs=[output_image, output_picks, output_csv])
|
697 |
+
|
698 |
+
with gr.Row():
|
699 |
+
with gr.Column(scale=2):
|
700 |
+
inputs_vel_model = [
|
701 |
+
## FIX FILE NAME ISSUE
|
702 |
+
gr.Slider(minimum=-180, maximum=180, value=0, step=5, label="Azimuth", interactive=True),
|
703 |
+
gr.Slider(minimum=-90, maximum=90, value=30, step=5, label="Elevation", interactive=True)
|
704 |
+
]
|
705 |
+
button = gr.Button("Look at 3D Velocities")
|
706 |
+
outputs_vel_model = gr.Image(label="3D Velocity Model")
|
707 |
+
|
708 |
+
button.click(compute_velocity_model,
|
709 |
+
inputs=inputs_vel_model,
|
710 |
+
outputs=outputs_vel_model)
|
711 |
+
|
712 |
+
|
713 |
+
|
714 |
+
|
715 |
+
|
716 |
+
demo.launch()
|
data/.DS_Store
CHANGED
Binary files a/data/.DS_Store and b/data/.DS_Store differ
|
|
data/velocity/35.766_-117.605_10.0_2019-07-04 17:33:49_3.csv
ADDED
@@ -0,0 +1,4 @@
|
|
|
|
|
|
|
|
|
|
|
1 |
+
station_name,st_lat,st_lon,starttime,"p_phase, s","p_uncertainty, s","s_phase, s","s_uncertainty, s","velocity_p, km/s","velocity_s, km/s"
|
2 |
+
CI.JRC2,35.98249,-117.80885,2019-07-04T17:33:39.947494Z,7.320212364196777,0.020417090272530913,13.38510799407959,0.028438671142794192,4.13660431013202,2.2622770044299756
|
3 |
+
CI.SRT,35.69235,-117.75051,2019-07-04T17:33:38.029990Z,4.53201961517334,0.01748959010001272,9.215676307678223,0.019567650742828846,3.4155476453388767,1.67967367867923
|
4 |
+
CI.WMF,36.11758,-117.85486,2019-07-04T17:33:41.867962Z,9.504384994506836,0.015920251607894897,17.03156852722168,0.04673808580264449,4.745724852828504,2.6483289549749593
|
data/velocity/35.766_-117.605_2019-07-04 17:33:49_3.csv
CHANGED
@@ -1,4 +1,4 @@
|
|
1 |
station_name,st_lat,st_lon,starttime,"p_phase, s","p_uncertainty, s","s_phase, s","s_uncertainty, s","velocity_p, km/s","velocity_s, km/s"
|
2 |
-
CI.
|
3 |
-
CI.
|
4 |
-
CI.
|
|
|
1 |
station_name,st_lat,st_lon,starttime,"p_phase, s","p_uncertainty, s","s_phase, s","s_uncertainty, s","velocity_p, km/s","velocity_s, km/s"
|
2 |
+
CI.WMF,36.11758,-117.85486,2019-07-04T17:33:41.867962Z,9.503650665283203,0.016685163136571646,17.022592544555664,0.04997979383915663,4.746091546067712,2.649725414106055
|
3 |
+
CI.SRT,35.69235,-117.75051,2019-07-04T17:33:38.029990Z,4.53201961517334,0.01748959010001272,9.215676307678223,0.019567650742828846,3.4155476453388767,1.67967367867923
|
4 |
+
CI.JRC2,35.98249,-117.80885,2019-07-04T17:33:39.947494Z,7.3213396072387695,0.014792646397836506,13.395279884338379,0.02523316943552345,4.135967410510336,2.2605591132307814
|
data/velocity/current_vel_model.csv
ADDED
@@ -0,0 +1,4 @@
|
|
|
|
|
|
|
|
|
|
|
1 |
+
station_name,st_lat,st_lon,starttime,"p_phase, s","p_uncertainty, s","s_phase, s","s_uncertainty, s","velocity_p, km/s","velocity_s, km/s"
|
2 |
+
CI.WMF,36.11758,-117.85486,2019-07-04T17:33:41.867962Z,9.503650665283203,0.016685163136571646,17.022592544555664,0.04997979383915663,4.746091546067712,2.649725414106055
|
3 |
+
CI.SRT,35.69235,-117.75051,2019-07-04T17:33:38.029990Z,4.53201961517334,0.01748959010001272,9.215676307678223,0.019567650742828846,3.4155476453388767,1.67967367867923
|
4 |
+
CI.JRC2,35.98249,-117.80885,2019-07-04T17:33:39.947494Z,7.3213396072387695,0.014792646397836506,13.395279884338379,0.02523316943552345,4.135967410510336,2.2605591132307814
|
phasehunter/__pycache__/data_preparation.cpython-311.pyc
ADDED
Binary file (9.14 kB). View file
|
|