added missing folder

Gradio_app.ipynb

@@ -152,7 +152,7 @@
     "    coeff = 1 - (var / (residuals.max() - residuals.min()))\n",
     "    return coeff\n",
     "\n",
-    "def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source_depth_km, velocity_model, max_waveforms):\n",
+    "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):\n",
     "    distances, t0s, st_lats, st_lons, waveforms, names = [], [], [], [], [], []\n",
     "    \n",
     "    taup_model = TauPyModel(model=velocity_model)\n",
@@ -175,6 +175,7 @@
     "                            minlongitude=(eq_lon-window), maxlongitude=(eq_lon+window), \n",
     "                            level='station')\n",
     "        print('Finished downloading inventory')\n",
+    "        \n",
     "    except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):\n",
     "        fig, ax = plt.subplots()\n",
     "        ax.text(0.5,0.5,'Something is wrong with the data provider, try another')\n",
@@ -187,7 +188,6 @@
     "    cached_waveforms = glob(\"data/cached/*.mseed\")\n",
     "\n",
     "    for network in inv:\n",
-    "        # Skip the SYntetic networks\n",
     "        if network.code == 'SY':\n",
     "            continue\n",
     "        for station in network:\n",
@@ -203,8 +203,9 @@
     "                starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15\n",
     "                endtime = starttime + 60\n",
     "                try:\n",
-    "                    if f\"data/cached/{network.code}_{station.code}_{starttime}.mseed\" not in cached_waveforms:\n",
-    "                        print('Downloading waveform')\n",
+    "                    filename=f'{network.code}_{station.code}_{starttime}'\n",
+    "                    if f\"data/cached/{filename}.mseed\" not in cached_waveforms:\n",
+    "                        print(f'Downloading waveform for {filename}')\n",
     "                        waveform = client.get_waveforms(network=network.code, station=station.code, location=\"*\", channel=\"*\", \n",
     "                                                    starttime=starttime, endtime=endtime)\n",
     "                        waveform.write(f\"data/cached/{network.code}_{station.code}_{starttime}.mseed\", format=\"MSEED\")\n",
@@ -245,12 +246,16 @@
     "    \n",
     "    # If there are no waveforms, return an empty plot\n",
     "    if len(waveforms) == 0:\n",
+    "        print('No waveforms found')\n",
     "        fig, ax = plt.subplots()\n",
     "        ax.text(0.5,0.5,'No waveforms found')\n",
     "        fig.canvas.draw();\n",
     "        image = np.array(fig.canvas.renderer.buffer_rgba())\n",
     "        plt.close(fig)\n",
-    "        return image\n",
+    "        output_picks = pd.DataFrame()\n",
+    "        output_picks.to_csv('data/picks.csv', index=False)\n",
+    "        output_csv = 'data/picks.csv'\n",
+    "        return image, output_picks, output_csv\n",
     "    \n",
     "\n",
     "    first_distances = bin_distances(distances, bin_size=10/111.2)\n",
@@ -277,9 +282,12 @@
     "    p_phases = output[:, 0]\n",
     "    s_phases = output[:, 1]\n",
     "\n",
-    "    # Max confidence - min variance     \n",
-    "    p_max_confidence = np.min([p_phases[i::len(waveforms)].std() for i in range(len(waveforms))]) \n",
-    "    s_max_confidence = np.min([s_phases[i::len(waveforms)].std() for i in range(len(waveforms))])\n",
+    "    p_phases = p_phases.reshape(len(waveforms),-1)\n",
+    "    s_phases = s_phases.reshape(len(waveforms),-1)\n",
+    "\n",
+    "    # Max confidence - min variance    \n",
+    "    p_max_confidence = p_phases.std(axis=-1).min()\n",
+    "    s_max_confidence = s_phases.std(axis=-1).min()\n",
     "\n",
     "    print(f\"Starting plotting {len(waveforms)} waveforms\")\n",
     "    fig, ax = plt.subplots(ncols=3, figsize=(10, 3))\n",
@@ -304,17 +312,18 @@
     "    topo_map.da.plot(ax = ax[2], cmap='Greys', add_colorbar=False, add_labels=False)\n",
     "    ax[1].imshow(hillshade, cmap=\"Greys\", alpha=0.5)\n",
     "\n",
-    "    output_picks = pd.DataFrame({'station_name' : [], 'starttime' : [], \n",
+    "    output_picks = pd.DataFrame({'station_name' : [], \n",
+    "                                'st_lat' : [], 'st_lon' : [],\n",
+    "                                 'starttime' : [], \n",
     "                                 'p_phase, s' : [], 'p_uncertainty, s' : [], \n",
     "                                 's_phase, s' : [], 's_uncertainty, s' : [],\n",
     "                                 'velocity_p, km/s' : [], 'velocity_s, km/s' : []})\n",
-    "                                  \n",
-    "                                  \n",
+    "                        \n",
     "    for i in range(len(waveforms)):\n",
     "        print(f\"Plotting waveform {i+1}/{len(waveforms)}\")\n",
-    "        current_P = p_phases[i::len(waveforms)]\n",
-    "        current_S = s_phases[i::len(waveforms)]\n",
-    "\n",
+    "        current_P = p_phases[i]\n",
+    "        current_S = s_phases[i]\n",
+    "        \n",
     "        x = [t0s[i] + pd.Timedelta(seconds=k/100) for k in np.linspace(0,6000,6000)]\n",
     "        x = mdates.date2num(x)\n",
     "\n",
@@ -322,32 +331,40 @@
     "        p_conf = 1/(current_P.std()/p_max_confidence).item()\n",
     "        s_conf = 1/(current_S.std()/s_max_confidence).item()\n",
     "\n",
-    "        ax[0].plot(x, waveforms[i][0, 0]*10+distances[i]*111.2, color='black', alpha=0.5, lw=1)\n",
+    "        delta_t = t0s[i].timestamp - obspy.UTCDateTime(timestamp).timestamp\n",
     "\n",
-    "        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='|')\n",
-    "        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='|')\n",
-    "        ax[0].set_ylabel('Z')\n",
+    "        ax[0].plot(x, waveforms[i][0, 0]*10+distances[i]*111.2, color='black', alpha=0.5, lw=1)\n",
     "\n",
-    "        delta_t = t0s[i].timestamp - obspy.UTCDateTime(timestamp).timestamp\n",
+    "        if (current_P.std().item()*60 < conf_thres_P) or (current_S.std().item()*60 < conf_thres_S):\n",
+    "            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='|')\n",
+    "            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='|')\n",
+    "        \n",
+    "            velocity_p = (distances[i]*111.2)/(delta_t+current_P.mean()*60).item()\n",
+    "            velocity_s = (distances[i]*111.2)/(delta_t+current_S.mean()*60).item()\n",
     "\n",
-    "        velocity_p = (distances[i]*111.2)/(delta_t+current_P.mean()*60).item()\n",
-    "        velocity_s = (distances[i]*111.2)/(delta_t+current_S.mean()*60).item()\n",
+    "            # Generate an array from st_lat to eq_lat and from st_lon to eq_lon\n",
+    "            x = np.linspace(st_lons[i], eq_lon, 50)\n",
+    "            y = np.linspace(st_lats[i], eq_lat, 50)\n",
+    "            \n",
+    "            # Plot the array\n",
+    "            ax[1].scatter(x, y, c=np.zeros_like(x)+velocity_p, alpha=0.1, vmin=0, vmax=8)\n",
+    "            ax[2].scatter(x, y, c=np.zeros_like(x)+velocity_s, alpha=0.1, vmin=0, vmax=8)\n",
     "\n",
+    "        else:\n",
+    "            velocity_p = np.nan\n",
+    "            velocity_s = np.nan\n",
+    "        \n",
+    "        ax[0].set_ylabel('Z')\n",
     "        print(f\"Station {st_lats[i]}, {st_lons[i]} has P velocity {velocity_p} and S velocity {velocity_s}\")\n",
     "\n",
-    "        output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]], 'starttime' : [str(t0s[i])], \n",
+    "        output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]], \n",
+    "                                                        'st_lat' : [st_lats[i]], 'st_lon' : [st_lons[i]],\n",
+    "                                                        'starttime' : [str(t0s[i])], \n",
     "                                                        'p_phase, s' : [(delta_t+current_P.mean()*60).item()], 'p_uncertainty, s' : [current_P.std().item()*60], \n",
     "                                                        's_phase, s' : [(delta_t+current_S.mean()*60).item()], 's_uncertainty, s' : [current_S.std().item()*60],\n",
     "                                                        'velocity_p, km/s' : [velocity_p], 'velocity_s, km/s' : [velocity_s]}))\n",
     "        \n",
-    "        # Generate an array from st_lat to eq_lat and from st_lon to eq_lon\n",
-    "        x = np.linspace(st_lons[i], eq_lon, 50)\n",
-    "        y = np.linspace(st_lats[i], eq_lat, 50)\n",
     "        \n",
-    "        # Plot the array\n",
-    "        ax[1].scatter(x, y, c=np.zeros_like(x)+velocity_p, alpha=0.1, vmin=0, vmax=8)\n",
-    "        ax[2].scatter(x, y, c=np.zeros_like(x)+velocity_s, alpha=0.1, vmin=0, vmax=8)\n",
-    "\n",
     "    # Add legend\n",
     "    ax[0].scatter(None, None, color='r', marker='|', label='P')\n",
     "    ax[0].scatter(None, None, color='b', marker='|', label='S')\n",
@@ -379,18 +396,11 @@
     "    fig.canvas.draw();\n",
     "    image = np.array(fig.canvas.renderer.buffer_rgba())\n",
     "    plt.close(fig)\n",
-    "\n",
-    "    output_picks.to_csv('data/picks.csv', index=False)\n",
-    "    output_csv = 'data/picks.csv'\n",
+    "    output_picks.to_csv(f'data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv', index=False)\n",
+    "    output_csv = f'data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv'\n",
     "\n",
     "    return image, output_picks, output_csv\n",
     "\n",
-    "def download_picks(output_picks):\n",
-    "    output_csv = io.BytesIO()\n",
-    "    output_picks.to_csv(output_csv, index=False)\n",
-    "    output_csv.seek(0)\n",
-    "    return output_csv\n",
-    "\n",
     "model = torch.jit.load(\"model.pt\")\n",
     "\n",
     "with gr.Blocks() as demo:\n",
@@ -448,7 +458,7 @@
     "        gr.HTML(\"\"\"\n",
     "        <div style=\"padding: 20px; border-radius: 10px; font-size: 16px;\">\n",
     "        <p style=\"font-weight: bold; font-size: 24px; margin-bottom: 20px;\">Using PhaseHunter to Analyze Seismic Waveforms</p>\n",
-    "        <p>Select an earthquake from the global earthquake catalogue 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>\n",
+    "        <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>\n",
     "        <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>\n",
     "        <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>\n",
     "        </div>\n",
@@ -500,7 +510,8 @@
     "            with gr.Column(scale=2):\n",
     "                radius_inputs = gr.Slider(minimum=1, \n",
     "                                        maximum=200, \n",
-    "                                        value=50, label=\"Radius (km)\", \n",
+    "                                        value=50, \n",
+    "                                        label=\"Radius (km)\", \n",
     "                                        step=10,\n",
     "                                        info=\"\"\"Select the radius around the earthquake to download data from.\\n \n",
     "                                        Note that the larger the radius, the longer the app will take to run.\"\"\",\n",
@@ -514,6 +525,23 @@
     "                                info=\"Maximum number of waveforms to show per section\\n (to avoid long prediction times)\",\n",
     "                                interactive=True,\n",
     "                                )\n",
+    "            with gr.Column(scale=2):\n",
+    "                P_thres_inputs = gr.Slider(minimum=0.01,\n",
+    "                                maximum=1,\n",
+    "                                value=0.1,\n",
+    "                                label=\"P uncertainty threshold, s\",\n",
+    "                                step=0.01,\n",
+    "                                info=\"Acceptable uncertainty for P picks expressed in std() seconds\",\n",
+    "                                interactive=True,\n",
+    "                                )\n",
+    "                S_thres_inputs = gr.Slider(minimum=0.01,\n",
+    "                                maximum=1,\n",
+    "                                value=0.2,\n",
+    "                                label=\"S uncertainty threshold, s\",\n",
+    "                                step=0.01,\n",
+    "                                info=\"Acceptable uncertainty for S picks expressed in std() seconds\",\n",
+    "                                interactive=True,\n",
+    "                                )\n",
     "            \n",
     "        button = gr.Button(\"Predict phases\")\n",
     "        output_image = gr.Image(label='Waveforms with Phases Marked', type='numpy', interactive=False)\n",
@@ -528,7 +556,8 @@
     "                 inputs=[client_inputs, timestamp_inputs, \n",
     "                         eq_lat_inputs, eq_lon_inputs, \n",
     "                         radius_inputs, source_depth_inputs, \n",
-    "                         velocity_inputs, max_waveforms_inputs],\n",
+    "                         velocity_inputs, max_waveforms_inputs,\n",
+    "                         P_thres_inputs, S_thres_inputs],\n",
     "                 outputs=[output_image, output_picks, output_csv])\n",
     "\n",
     "demo.launch()"