section plot updated

Gradio_app.ipynb

@@ -61,6 +61,27 @@
      "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Error in callback <function _draw_all_if_interactive at 0x1774d0ea0> (for post_execute):\n"
+     ]
+    },
+    {
+     "ename": "KeyboardInterrupt",
+     "evalue": "",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
+      "File \u001b[0;32m~/miniconda3/envs/phasehunter/lib/python3.11/site-packages/matplotlib/pyplot.py:120\u001b[0m, in \u001b[0;36m_draw_all_if_interactive\u001b[0;34m()\u001b[0m\n\u001b[1;32m    118\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m_draw_all_if_interactive\u001b[39m():\n\u001b[1;32m    119\u001b[0m     \u001b[39mif\u001b[39;00m matplotlib\u001b[39m.\u001b[39mis_interactive():\n\u001b[0;32m--> 120\u001b[0m         draw_all()\n",
+      "File \u001b[0;32m~/miniconda3/envs/phasehunter/lib/python3.11/site-packages/matplotlib/_pylab_helpers.py:132\u001b[0m, in \u001b[0;36mGcf.draw_all\u001b[0;34m(cls, force)\u001b[0m\n\u001b[1;32m    130\u001b[0m \u001b[39mfor\u001b[39;00m manager \u001b[39min\u001b[39;00m \u001b[39mcls\u001b[39m\u001b[39m.\u001b[39mget_all_fig_managers():\n\u001b[1;32m    131\u001b[0m     \u001b[39mif\u001b[39;00m force \u001b[39mor\u001b[39;00m manager\u001b[39m.\u001b[39mcanvas\u001b[39m.\u001b[39mfigure\u001b[39m.\u001b[39mstale:\n\u001b[0;32m--> 132\u001b[0m         manager\u001b[39m.\u001b[39;49mcanvas\u001b[39m.\u001b[39;49mdraw_idle()\n",
+      "File \u001b[0;32m~/miniconda3/envs/phasehunter/lib/python3.11/site-packages/matplotlib/backend_bases.py:2082\u001b[0m, in \u001b[0;36mFigureCanvasBase.draw_idle\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2080\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_is_idle_drawing:\n\u001b[1;32m   2081\u001b[0m     \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_idle_draw_cntx():\n\u001b[0;32m-> 2082\u001b[0m         \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mdraw(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n",
+      "File \u001b[0;32m~/miniconda3/envs/phasehunter/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:397\u001b[0m, in \u001b[0;36mFigureCanvasAgg.draw\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    395\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mrenderer\u001b[39m.\u001b[39mclear()\n\u001b[1;32m    396\u001b[0m \u001b[39m# Acquire a lock on the shared font cache.\u001b[39;00m\n\u001b[0;32m--> 397\u001b[0m \u001b[39mwith\u001b[39;49;00m RendererAgg\u001b[39m.\u001b[39;49mlock, \\\n\u001b[1;32m    398\u001b[0m      (\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtoolbar\u001b[39m.\u001b[39;49m_wait_cursor_for_draw_cm() \u001b[39mif\u001b[39;49;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtoolbar\n\u001b[1;32m    399\u001b[0m       \u001b[39melse\u001b[39;49;00m nullcontext()):\n\u001b[1;32m    400\u001b[0m     \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfigure\u001b[39m.\u001b[39;49mdraw(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrenderer)\n\u001b[1;32m    401\u001b[0m     \u001b[39m# A GUI class may be need to update a window using this draw, so\u001b[39;49;00m\n\u001b[1;32m    402\u001b[0m     \u001b[39m# don't forget to call the superclass.\u001b[39;49;00m\n",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+     ]
     }
    ],
    "source": [

app.py

@@ -21,6 +21,8 @@ from obspy.clients.fdsn.header import URL_MAPPINGS
 import matplotlib.pyplot as plt
 import matplotlib.dates as mdates
 
+from glob import glob
+
 def make_prediction(waveform):
     waveform = np.load(waveform)
     processed_input = prepare_waveform(waveform)
@@ -80,6 +82,15 @@ def mark_phases(waveform):
     plt.close(fig)
     return image
 
+def variance_coefficient(residuals):
+    # calculate the variance of the residuals
+    var = residuals.var()
+    
+    # scale the variance to a coefficient between 0 and 1
+    coeff = 1 - (var / (residuals.max() - residuals.min()))
+    
+    return coeff
+
 def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source_depth_km, velocity_model):
     distances, t0s, st_lats, st_lons, waveforms = [], [], [], [], []
     
@@ -101,6 +112,8 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
                           level='station')
     
     waveforms = []
+    cached_waveforms = glob("data/cached/*.mseed")
+
     for network in inv:
         for station in network:
             try:
@@ -115,8 +128,12 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
                     starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15
                     endtime = starttime + 60
 
-                    waveform = client.get_waveforms(network=network.code, station=station.code, location="*", channel="*", 
-                                                starttime=starttime, endtime=endtime)
+                    if f"data/cached/{network.code}_{station.code}_{starttime}.mseed" not in cached_waveforms:
+                        waveform = client.get_waveforms(network=network.code, station=station.code, location="*", channel="*", 
+                                                    starttime=starttime, endtime=endtime)
+                        waveform.write(f"data/cached/{network.code}_{station.code}_{starttime}.mseed", format="MSEED")
+                    else:
+                        waveform = obspy.read(f"data/cached/{network.code}_{station.code}_{starttime}.mseed")
                     
                     waveform = waveform.select(channel="H[BH][ZNE]")
                     waveform = waveform.merge(fill_value=0)
@@ -148,29 +165,39 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
     p_phases = output[:, 0]
     s_phases = output[:, 1]
 
-    fig, ax = plt.subplots(nrows=1, figsize=(10, 3), sharex=True)
+    # Max confidence - min variance     
+    p_max_confidence = np.min([p_phases[i::len(waveforms)].std() for i in range(len(waveforms))]) 
+    s_max_confidence = np.min([s_phases[i::len(waveforms)].std() for i in range(len(waveforms))])
+
+    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 3), sharex=True)
     for i in range(len(waveforms)):
         current_P = p_phases[i::len(waveforms)]
         current_S = s_phases[i::len(waveforms)]
+
         x = [t0s[i] + pd.Timedelta(seconds=k/100) for k in np.linspace(0,6000,6000)]
         x = mdates.date2num(x)
-        ax.plot(x, waveforms[i][0, 0]+distances[i]*111.2, color='black', alpha=0.5)
-        ax.scatter(x[int(current_P.mean()*waveforms[i][0].shape[-1])], waveforms[i][0, 0].mean()+distances[i]*111.2, color='r')
-        ax.scatter(x[int(current_S.mean()*waveforms[i][0].shape[-1])], waveforms[i][0, 0].mean()+distances[i]*111.2, color='b')
-        ax.set_ylabel('Z')
 
-        ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))
-        ax.xaxis.set_major_locator(mdates.SecondLocator(interval=10))
+        # Normalize confidence for the plot
+        p_conf = 1/(current_P.std()/p_max_confidence).item()
+        s_conf = 1/(current_S.std()/s_max_confidence).item()
 
-        # for a in ax:
-        #     a.axvline(current_P.mean()*waveforms[i][0].shape[-1], color='r', linestyle='--', label='P')
-        #     a.axvline(current_S.mean()*waveforms[i][0].shape[-1], color='b', linestyle='--', label='S')
+        ax[0].plot(x, waveforms[i][0, 0]*10+distances[i]*111.2, color='black', alpha=0.5, lw=1)
 
-        # ax[-1].set_xlabel('Time, samples')
-        # ax[-1].set_ylabel('Uncert.')
-        # ax[-1].legend()
+        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='|')
+        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='|')
+        ax[0].set_ylabel('Z')
 
-        plt.subplots_adjust(hspace=0., wspace=0.)
+        ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))
+        ax[0].xaxis.set_major_locator(mdates.SecondLocator(interval=5))
+        
+    ax[0].scatter(None, None, color='r', marker='|', label='P')
+    ax[0].scatter(None, None, color='b', marker='|', label='S')
+    ax[0].legend()
+
+    ax[1].scatter(st_lats, st_lons, color='b', marker='d', label='Stations')
+    ax[1].scatter(eq_lat, eq_lon, color='r', marker='*', label='Earthquake')
+    ax[1].legend()
+    plt.subplots_adjust(hspace=0., wspace=0.)
         
     fig.canvas.draw();
     image = np.array(fig.canvas.renderer.buffer_rgba())
@@ -184,12 +211,6 @@ model = Onset_picker.load_from_checkpoint("./weights.ckpt",
                                     learning_rate=3e-4)
 model.eval()
 
-
-
-# # Create the Gradio interface
-# gr.Interface(mark_phases, inputs, outputs, title='PhaseHunter').launch()
-
-
 with gr.Blocks() as demo:
     gr.Markdown("# PhaseHunter")
     gr.Markdown("""This app allows one to detect P and S seismic phases along with uncertainty of the detection. 
@@ -250,7 +271,9 @@ with gr.Blocks() as demo:
             radius_inputs = gr.Slider(minimum=1, 
                                     maximum=150, 
                                     value=50, label="Radius (km)", 
-                                    info="Select the radius around the earthquake to download data from",
+                                    step=10,
+                                    info="""Select the radius around the earthquake to download data from.\n 
+                                    Note that the larger the radius, the longer the app will take to run.""",
                                     interactive=True)
             
             velocity_inputs = gr.Dropdown(
@@ -263,7 +286,7 @@ with gr.Blocks() as demo:
             
             
         button = gr.Button("Predict phases")
-        outputs_section = gr.outputs.Image(label='Waveforms with Phases Marked', type='numpy')
+        outputs_section = gr.Image(label='Waveforms with Phases Marked', type='numpy', interactive=False)
         
         button.click(predict_on_section, 
                  inputs=[client_inputs, timestamp_inputs, 
@@ -277,6 +300,4 @@ with gr.Blocks() as demo:
         Your waveform should be sampled at 100 sps and have 3 (Z, N, E) or 1 (Z) channels.
         """)
 
-
-
 demo.launch()