fixed uploading bug (1,...) shape

Gradio_app.ipynb

@@ -2,14 +2,14 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 55,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Running on local URL:  http://127.0.0.1:7901\n",
+      "Running on local URL:  http://127.0.0.1:7866\n",
       "\n",
       "To create a public link, set `share=True` in `launch()`.\n"
      ]
@@ -17,7 +17,7 @@
     {
      "data": {
       "text/html": [
-       "<div><iframe src=\"http://127.0.0.1:7901/\" width=\"100%\" height=\"500\" allow=\"autoplay; camera; microphone; clipboard-read; clipboard-write;\" frameborder=\"0\" allowfullscreen></iframe></div>"
+       "<div><iframe src=\"http://127.0.0.1:7866/\" width=\"100%\" height=\"500\" allow=\"autoplay; camera; microphone; clipboard-read; clipboard-write;\" frameborder=\"0\" allowfullscreen></iframe></div>"
       ],
       "text/plain": [
        "<IPython.core.display.HTML object>"
@@ -30,9 +30,23 @@
      "data": {
       "text/plain": []
      },
-     "execution_count": 55,
+     "execution_count": 13,
      "metadata": {},
      "output_type": "execute_result"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0.13119522482156754\n"
+     ]
     }
    ],
    "source": [
@@ -66,6 +80,9 @@
     "\n",
     "def make_prediction(waveform):\n",
     "    waveform = np.load(waveform)\n",
+    "    if len(waveform.shape) == 1:\n",
+    "        waveform = waveform.reshape(1, waveform.shape[0])\n",
+    "\n",
     "    processed_input = prepare_waveform(waveform)\n",
     "    \n",
     "    # Make prediction\n",
@@ -77,7 +94,7 @@
     "\n",
     "    return processed_input, p_phase, s_phase\n",
     "\n",
-    "def mark_phases(waveform, uploaded_file):\n",
+    "def mark_phases(waveform, uploaded_file, p_thres, s_thres):\n",
     "\n",
     "    if uploaded_file is not None:\n",
     "        waveform = uploaded_file.name\n",
@@ -101,19 +118,28 @@
     "        ax[1].set_ylabel('N')\n",
     "        ax[2].set_ylabel('E')\n",
     "\n",
-    "    p_phase_plot = p_phase*processed_input.shape[-1]\n",
-    "    p_kde = gaussian_kde(p_phase_plot)\n",
-    "    p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )\n",
-    "    ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')\n",
-    "\n",
-    "    s_phase_plot = s_phase*processed_input.shape[-1]\n",
-    "    s_kde = gaussian_kde(s_phase_plot)\n",
-    "    s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )\n",
-    "    ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')\n",
-    "\n",
-    "    for a in ax:\n",
-    "        a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P')\n",
-    "        a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S')\n",
+    "    print(p_phase.std().item()*60)\n",
+    "    do_we_have_p = (p_phase.std().item()*60 < p_thres)\n",
+    "    if do_we_have_p:\n",
+    "        p_phase_plot = p_phase*processed_input.shape[-1]\n",
+    "        p_kde = gaussian_kde(p_phase_plot)\n",
+    "        p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )\n",
+    "        ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')\n",
+    "    else:\n",
+    "        ax[-1].text(0.5, 0.75, 'No P phase detected', horizontalalignment='center', verticalalignment='center', transform=ax[-1].transAxes)\n",
+    "\n",
+    "    do_we_have_s = (s_phase.std().item()*60 < s_thres)\n",
+    "    if do_we_have_s:\n",
+    "        s_phase_plot = s_phase*processed_input.shape[-1]\n",
+    "        s_kde = gaussian_kde(s_phase_plot)\n",
+    "        s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )\n",
+    "        ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')\n",
+    "\n",
+    "        for a in ax:\n",
+    "            a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P', alpha=do_we_have_p)\n",
+    "            a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S', alpha=do_we_have_s)\n",
+    "    else:\n",
+    "        ax[-1].text(0.5, 0.25, 'No S phase detected', horizontalalignment='center', verticalalignment='center', transform=ax[-1].transAxes)\n",
     "\n",
     "    ax[-1].set_xlabel('Time, samples')\n",
     "    ax[-1].set_ylabel('Uncert., samples')\n",
@@ -401,6 +427,109 @@
     "\n",
     "    return image, output_picks, output_csv\n",
     "\n",
+    "import numpy as np\n",
+    "from matplotlib import colors, cm\n",
+    "\n",
+    "# Function to find the closest index for a given value in an array\n",
+    "def find_closest_index(array, value):\n",
+    "    return np.argmin(np.abs(array - value))\n",
+    "\n",
+    "def compute_velocity_model(azimuth, elevation):\n",
+    "    filename = output_csv\n",
+    "    df = pd.read_csv(filename)\n",
+    "\n",
+    "    # Current EQ location\n",
+    "    filename = filename.split(\"/\")[-1]\n",
+    "    eq_lat = float(filename.split(\"_\")[0])\n",
+    "    eq_lon = float(filename.split(\"_\")[1])\n",
+    "\n",
+    "    ## FIX THIS LATTER\n",
+    "    eq_depth = 10 ##FIX THIS LATTER\n",
+    "    ## FIX THIS LATTER\n",
+    "\n",
+    "    # Define the region of interest (latitude, longitude, and depth ranges)\n",
+    "    lat_range = (np.min([df.st_lat.min(), eq_lat]), np.max([df.st_lat.max(), eq_lat]))\n",
+    "    lon_range = (np.min([df.st_lon.min(), eq_lon]),  np.max([df.st_lon.max(), eq_lon]))\n",
+    "    depth_range = (0, 50)\n",
+    "\n",
+    "    # Define the number of nodes in each dimension\n",
+    "    n_lat = 10\n",
+    "    n_lon = 10\n",
+    "    n_depth = 10\n",
+    "    num_points = 100\n",
+    "\n",
+    "    # Create the grid\n",
+    "    lat_values = np.linspace(lat_range[0], lat_range[1], n_lat)\n",
+    "    lon_values = np.linspace(lon_range[0], lon_range[1], n_lon)\n",
+    "    depth_values = np.linspace(depth_range[0], depth_range[1], n_depth)\n",
+    "\n",
+    "    # Initialize the velocity model with constant values\n",
+    "    initial_velocity = 0  # km/s, this can be P-wave or S-wave velocity\n",
+    "    velocity_model = np.full((n_lat, n_lon, n_depth), initial_velocity, dtype=float)\n",
+    "\n",
+    "    # Loop through the stations and update the velocity model\n",
+    "    for i in range(len(df)):\n",
+    "        if ~np.isnan(df['velocity_p, km/s'].iloc[i]):\n",
+    "            # Interpolate coordinates along the great circle path between the earthquake and the station\n",
+    "            lon_deg = np.linspace(df.st_lon.iloc[i], eq_lon, num_points)\n",
+    "            lat_deg = np.linspace(df.st_lat.iloc[i], eq_lat, num_points)\n",
+    "            depth_interpolated = np.interp(np.linspace(0, 1, num_points), [0, 1], [eq_depth, 0])\n",
+    "\n",
+    "            # Loop through the interpolated coordinates and update the grid cells with the average P-wave velocity\n",
+    "            for lat, lon, depth in zip(lat_deg, lon_deg, depth_interpolated):\n",
+    "                lat_index = find_closest_index(lat_values, lat)\n",
+    "                lon_index = find_closest_index(lon_values, lon)\n",
+    "                depth_index = find_closest_index(depth_values, depth)\n",
+    "                \n",
+    "                if velocity_model[lat_index, lon_index, depth_index] == initial_velocity:\n",
+    "                    velocity_model[lat_index, lon_index, depth_index] = df['velocity_p, km/s'].iloc[i]\n",
+    "                else:\n",
+    "                    velocity_model[lat_index, lon_index, depth_index] = (velocity_model[lat_index, lon_index, depth_index] +\n",
+    "                                                                        df['velocity_p, km/s'].iloc[i]) / 2\n",
+    "                    \n",
+    "    # Create the figure and axis\n",
+    "    fig = plt.figure(figsize=(8, 8))\n",
+    "    ax = fig.add_subplot(111, projection='3d')\n",
+    "\n",
+    "    # Set the plot limits\n",
+    "    ax.set_xlim3d(lat_range[0], lat_range[1])\n",
+    "    ax.set_ylim3d(lon_range[0], lon_range[1])\n",
+    "    ax.set_zlim3d(depth_range[1], depth_range[0])\n",
+    "\n",
+    "    ax.set_xlabel('Latitude')\n",
+    "    ax.set_ylabel('Longitude')\n",
+    "    ax.set_zlabel('Depth (km)')\n",
+    "    ax.set_title('Velocity Model')\n",
+    "    \n",
+    "    # Create the meshgrid\n",
+    "    x, y, z = np.meshgrid(\n",
+    "        np.linspace(lat_range[0], lat_range[1], velocity_model.shape[0]+1),\n",
+    "        np.linspace(lon_range[0], lon_range[1], velocity_model.shape[1]+1),\n",
+    "        np.linspace(depth_range[0], depth_range[1], velocity_model.shape[2]+1),\n",
+    "        indexing='ij'\n",
+    "    )\n",
+    "\n",
+    "    # Create the color array\n",
+    "    norm = plt.Normalize(vmin=velocity_model.min(), vmax=velocity_model.max())\n",
+    "    colors = plt.cm.plasma(norm(velocity_model))\n",
+    "\n",
+    "    # Plot the voxels\n",
+    "    ax.voxels(x, y, z, velocity_model > 0, facecolors=colors, alpha=0.5, edgecolor='k')\n",
+    "\n",
+    "    # Set the view angle\n",
+    "    ax.view_init(elev=elevation, azim=azimuth)\n",
+    "\n",
+    "    m = cm.ScalarMappable(cmap=plt.cm.plasma, norm=norm)\n",
+    "    m.set_array([])\n",
+    "    plt.colorbar(m)\n",
+    "\n",
+    "    # Show the plot\n",
+    "    fig.canvas.draw();\n",
+    "    image = np.array(fig.canvas.renderer.buffer_rgba())\n",
+    "    plt.close(fig)\n",
+    "\n",
+    "    return image\n",
+    "\n",
     "model = torch.jit.load(\"model.pt\")\n",
     "\n",
     "with gr.Blocks() as demo:\n",
@@ -433,25 +562,6 @@
     "    </ul>\n",
     "</div>\n",
     "\"\"\")\n",
-    "\n",
-    "    with gr.Tab(\"Try on a single station\"):\n",
-    "        with gr.Row(): \n",
-    "            # Define the input and output types for Gradio\n",
-    "            inputs = gr.Dropdown(\n",
-    "                [\"data/sample/sample_0.npy\", \n",
-    "                \"data/sample/sample_1.npy\", \n",
-    "                \"data/sample/sample_2.npy\"], \n",
-    "                label=\"Sample waveform\", \n",
-    "                info=\"Select one of the samples\",\n",
-    "                value = \"data/sample/sample_0.npy\"\n",
-    "            )\n",
-    "\n",
-    "            upload = gr.File(label=\"Or upload your own waveform\")\n",
-    "\n",
-    "        button = gr.Button(\"Predict phases\")\n",
-    "        outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)\n",
-    "    \n",
-    "        button.click(mark_phases, inputs=[inputs, upload], outputs=outputs)\n",
     "        \n",
     "    with gr.Tab(\"Select earthquake from catalogue\"):\n",
     "\n",
@@ -559,10 +669,134 @@
     "                         velocity_inputs, max_waveforms_inputs,\n",
     "                         P_thres_inputs, S_thres_inputs],\n",
     "                 outputs=[output_image, output_picks, output_csv])\n",
+    "        \n",
+    "        with gr.Row():\n",
+    "            with gr.Column(scale=2):\n",
+    "                inputs_vel_model = [\n",
+    "                    ## FIX FILE NAME ISSUE\n",
+    "                    gr.Slider(minimum=-180, maximum=180, value=0, step=5, label=\"Azimuth\", interactive=True),\n",
+    "                    gr.Slider(minimum=-90, maximum=90, value=30, step=5, label=\"Elevation\", interactive=True)\n",
+    "                ]\n",
+    "                button = gr.Button(\"Look at 3D Velocities\")\n",
+    "            outputs_vel_model = gr.Image(label=\"3D Velocity Model\")\n",
+    "\n",
+    "            button.click(compute_velocity_model, \n",
+    "                         inputs=inputs_vel_model, \n",
+    "                         outputs=outputs_vel_model)\n",
+    "\n",
+    "    with gr.Tab(\"Try on a single station\"):\n",
+    "        with gr.Row(): \n",
+    "            # Define the input and output types for Gradio\n",
+    "            inputs = gr.Dropdown(\n",
+    "                [\"data/sample/sample_0.npy\", \n",
+    "                \"data/sample/sample_1.npy\", \n",
+    "                \"data/sample/sample_2.npy\"], \n",
+    "                label=\"Sample waveform\", \n",
+    "                info=\"Select one of the samples\",\n",
+    "                value = \"data/sample/sample_0.npy\"\n",
+    "            )\n",
+    "            with gr.Column(scale=1):\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",
+    "                \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",
+    "            upload = gr.File(label=\"Or upload your own waveform\")\n",
+    "\n",
+    "        button = gr.Button(\"Predict phases\")\n",
+    "        outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)\n",
+    "    \n",
+    "        button.click(mark_phases, inputs=[inputs, upload, P_thres_inputs, S_thres_inputs], outputs=outputs)\n",
+    "\n",
+    "            \n",
+    "\n",
     "\n",
     "demo.launch()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 33,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "output_csv.value"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.save(\"test.npy\", np.random.randint(0,10, size=(6000)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Running on local URL:  http://127.0.0.1:7869\n",
+      "\n",
+      "To create a public link, set `share=True` in `launch()`.\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<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>"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": []
+     },
+     "execution_count": 24,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/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",
+      "  plt.colorbar(m)\n"
+     ]
+    }
+   ],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import gradio as gr\n",
+    "    \n",
+    "# Define the Gradio interface\n",
+    "\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -573,7 +807,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "phasehunter",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -587,14 +821,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.2"
+   "version": "3.9.8"
   },
-  "orig_nbformat": 4,
-  "vscode": {
-   "interpreter": {
-    "hash": "6bf57068982d7b420bddaaf1d0614a7795947176033057024cf47d8ca2c1c4cd"
-   }
-  }
+  "orig_nbformat": 4
  },
  "nbformat": 4,
  "nbformat_minor": 2

app.py

@@ -13,7 +13,11 @@ import earthpy.spatial as es
 
 import obspy
 from obspy.clients.fdsn import Client
-from obspy.clients.fdsn.header import FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException
+from obspy.clients.fdsn.header import (
+    FDSNNoDataException,
+    FDSNTimeoutException,
+    FDSNInternalServerException,
+)
 from obspy.geodetics.base import locations2degrees
 from obspy.taup import TauPyModel
 from obspy.taup.helper_classes import SlownessModelError
@@ -26,10 +30,14 @@ from mpl_toolkits.axes_grid1 import ImageGrid
 
 from glob import glob
 
+
 def make_prediction(waveform):
     waveform = np.load(waveform)
+    if len(waveform.shape) == 1:
+        waveform = waveform.reshape(1, waveform.shape[0])
+
     processed_input = prepare_waveform(waveform)
-    
+
     # Make prediction
     with torch.inference_mode():
         output = model(processed_input)
@@ -39,7 +47,8 @@ def make_prediction(waveform):
 
     return processed_input, p_phase, s_phase
 
-def mark_phases(waveform, uploaded_file):
+
+def mark_phases(waveform, uploaded_file, p_thres, s_thres):
 
     if uploaded_file is not None:
         waveform = uploaded_file.name
@@ -47,41 +56,76 @@ def mark_phases(waveform, uploaded_file):
     processed_input, p_phase, s_phase = make_prediction(waveform)
 
     # Create a plot of the waveform with the phases marked
-    if sum(processed_input[0][2] == 0): #if input is 1C
+    if sum(processed_input[0][2] == 0):  # if input is 1C
         fig, ax = plt.subplots(nrows=2, figsize=(10, 2), sharex=True)
 
-        ax[0].plot(processed_input[0][0], color='black', lw=1)
-        ax[0].set_ylabel('Norm. Ampl.')
+        ax[0].plot(processed_input[0][0], color="black", lw=1)
+        ax[0].set_ylabel("Norm. Ampl.")
 
-    else: #if input is 3C
+    else:  # if input is 3C
         fig, ax = plt.subplots(nrows=4, figsize=(10, 6), sharex=True)
-        ax[0].plot(processed_input[0][0], color='black', lw=1)
-        ax[1].plot(processed_input[0][1], color='black', lw=1)
-        ax[2].plot(processed_input[0][2], color='black', lw=1)
-
-        ax[0].set_ylabel('Z')
-        ax[1].set_ylabel('N')
-        ax[2].set_ylabel('E')
-
-    p_phase_plot = p_phase*processed_input.shape[-1]
-    p_kde = gaussian_kde(p_phase_plot)
-    p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )
-    ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')
-
-    s_phase_plot = s_phase*processed_input.shape[-1]
-    s_kde = gaussian_kde(s_phase_plot)
-    s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )
-    ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')
-
-    for a in ax:
-        a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P')
-        a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S')
-
-    ax[-1].set_xlabel('Time, samples')
-    ax[-1].set_ylabel('Uncert., samples')
+        ax[0].plot(processed_input[0][0], color="black", lw=1)
+        ax[1].plot(processed_input[0][1], color="black", lw=1)
+        ax[2].plot(processed_input[0][2], color="black", lw=1)
+
+        ax[0].set_ylabel("Z")
+        ax[1].set_ylabel("N")
+        ax[2].set_ylabel("E")
+
+    print(p_phase.std().item() * 60)
+    do_we_have_p = p_phase.std().item() * 60 < p_thres
+    if do_we_have_p:
+        p_phase_plot = p_phase * processed_input.shape[-1]
+        p_kde = gaussian_kde(p_phase_plot)
+        p_dist_space = np.linspace(min(p_phase_plot) - 10, max(p_phase_plot) + 10, 500)
+        ax[-1].plot(p_dist_space, p_kde(p_dist_space), color="r")
+    else:
+        ax[-1].text(
+            0.5,
+            0.75,
+            "No P phase detected",
+            horizontalalignment="center",
+            verticalalignment="center",
+            transform=ax[-1].transAxes,
+        )
+
+    do_we_have_s = s_phase.std().item() * 60 < s_thres
+    if do_we_have_s:
+        s_phase_plot = s_phase * processed_input.shape[-1]
+        s_kde = gaussian_kde(s_phase_plot)
+        s_dist_space = np.linspace(min(s_phase_plot) - 10, max(s_phase_plot) + 10, 500)
+        ax[-1].plot(s_dist_space, s_kde(s_dist_space), color="b")
+
+        for a in ax:
+            a.axvline(
+                p_phase.mean() * processed_input.shape[-1],
+                color="r",
+                linestyle="--",
+                label="P",
+                alpha=do_we_have_p,
+            )
+            a.axvline(
+                s_phase.mean() * processed_input.shape[-1],
+                color="b",
+                linestyle="--",
+                label="S",
+                alpha=do_we_have_s,
+            )
+    else:
+        ax[-1].text(
+            0.5,
+            0.25,
+            "No S phase detected",
+            horizontalalignment="center",
+            verticalalignment="center",
+            transform=ax[-1].transAxes,
+        )
+
+    ax[-1].set_xlabel("Time, samples")
+    ax[-1].set_ylabel("Uncert., samples")
     ax[-1].legend()
 
-    plt.subplots_adjust(hspace=0., wspace=0.)
+    plt.subplots_adjust(hspace=0.0, wspace=0.0)
 
     # Convert the plot to an image and return it
     fig.canvas.draw()
@@ -89,6 +133,7 @@ def mark_phases(waveform, uploaded_file):
     plt.close(fig)
     return image
 
+
 def bin_distances(distances, bin_size=10):
     # Bin the distances into groups of `bin_size` kilometers
     binned_distances = {}
@@ -104,9 +149,10 @@ def bin_distances(distances, bin_size=10):
     for bin_index in binned_distances:
         first_distance, first_distance_index = binned_distances[bin_index]
         first_distances.append(first_distance_index)
-    
+
     return first_distances
 
+
 def variance_coefficient(residuals):
     # calculate the variance of the residuals
     var = residuals.var()
@@ -114,9 +160,21 @@ def variance_coefficient(residuals):
     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, max_waveforms, conf_thres_P, conf_thres_S):
+
+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,
+):
     distances, t0s, st_lats, st_lons, waveforms, names = [], [], [], [], [], []
-    
+
     taup_model = TauPyModel(model=velocity_model)
     client = Client(client_name)
 
@@ -130,60 +188,92 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
     endtime = starttime + 120
 
     try:
-        print('Starting to download inventory')
-        inv = client.get_stations(network="*", station="*", location="*", channel="*H*", 
-                            starttime=starttime, endtime=endtime, 
-                            minlatitude=(eq_lat-window), maxlatitude=(eq_lat+window),
-                            minlongitude=(eq_lon-window), maxlongitude=(eq_lon+window), 
-                            level='station')
-        print('Finished downloading inventory')
-        
-    except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):
+        print("Starting to download inventory")
+        inv = client.get_stations(
+            network="*",
+            station="*",
+            location="*",
+            channel="*H*",
+            starttime=starttime,
+            endtime=endtime,
+            minlatitude=(eq_lat - window),
+            maxlatitude=(eq_lat + window),
+            minlongitude=(eq_lon - window),
+            maxlongitude=(eq_lon + window),
+            level="station",
+        )
+        print("Finished downloading inventory")
+
+    except (
+        IndexError,
+        FDSNNoDataException,
+        FDSNTimeoutException,
+        FDSNInternalServerException,
+    ):
         fig, ax = plt.subplots()
-        ax.text(0.5,0.5,'Something is wrong with the data provider, try another')
-        fig.canvas.draw();
+        ax.text(0.5, 0.5, "Something is wrong with the data provider, try another")
+        fig.canvas.draw()
         image = np.array(fig.canvas.renderer.buffer_rgba())
         plt.close(fig)
         return image
-    
+
     waveforms = []
     cached_waveforms = glob("data/cached/*.mseed")
 
     for network in inv:
-        if network.code == 'SY':
+        if network.code == "SY":
             continue
         for station in network:
             print(f"Processing {network.code}.{station.code}...")
-            distance = locations2degrees(eq_lat, eq_lon, station.latitude, station.longitude)
+            distance = locations2degrees(
+                eq_lat, eq_lon, station.latitude, station.longitude
+            )
 
-            arrivals = taup_model.get_travel_times(source_depth_in_km=source_depth_km, 
-                                                    distance_in_degree=distance, 
-                                                    phase_list=["P", "S"])
+            arrivals = taup_model.get_travel_times(
+                source_depth_in_km=source_depth_km,
+                distance_in_degree=distance,
+                phase_list=["P", "S"],
+            )
 
             if len(arrivals) > 0:
 
                 starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15
                 endtime = starttime + 60
                 try:
-                    filename=f'{network.code}_{station.code}_{starttime}'
+                    filename = f"{network.code}_{station.code}_{starttime}"
                     if f"data/cached/{filename}.mseed" not in cached_waveforms:
-                        print(f'Downloading waveform for {filename}')
-                        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")
-                        print('Finished downloading and caching waveform')
+                        print(f"Downloading waveform for {filename}")
+                        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",
+                        )
+                        print("Finished downloading and caching waveform")
                     else:
-                        print('Reading cached waveform')
-                        waveform = obspy.read(f"data/cached/{network.code}_{station.code}_{starttime}.mseed")
-                        
-
-                except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):
-                    print(f'Skipping {network.code}_{station.code}_{starttime}')
+                        print("Reading cached waveform")
+                        waveform = obspy.read(
+                            f"data/cached/{network.code}_{station.code}_{starttime}.mseed"
+                        )
+
+                except (
+                    IndexError,
+                    FDSNNoDataException,
+                    FDSNTimeoutException,
+                    FDSNInternalServerException,
+                ):
+                    print(f"Skipping {network.code}_{station.code}_{starttime}")
                     continue
-            
+
                 waveform = waveform.select(channel="H[BH][ZNE]")
                 waveform = waveform.merge(fill_value=0)
-                waveform = waveform[:3].sort(keys=['channel'], reverse=True)
+                waveform = waveform[:3].sort(keys=["channel"], reverse=True)
 
                 len_check = [len(x.data) for x in waveform]
                 if len(set(len_check)) > 1:
@@ -191,7 +281,9 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
 
                 if len(waveform) == 3:
                     try:
-                        waveform = prepare_waveform(np.stack([x.data for x in waveform]))
+                        waveform = prepare_waveform(
+                            np.stack([x.data for x in waveform])
+                        )
 
                         distances.append(distance)
                         t0s.append(starttime)
@@ -200,32 +292,32 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
                         waveforms.append(waveform)
                         names.append(f"{network.code}.{station.code}")
 
-                        print(f"Added {network.code}.{station.code} to the list of waveforms")
+                        print(
+                            f"Added {network.code}.{station.code} to the list of waveforms"
+                        )
 
                     except:
                         continue
-                
-    
+
     # If there are no waveforms, return an empty plot
     if len(waveforms) == 0:
-        print('No waveforms found')
+        print("No waveforms found")
         fig, ax = plt.subplots()
-        ax.text(0.5,0.5,'No waveforms found')
-        fig.canvas.draw();
+        ax.text(0.5, 0.5, "No waveforms found")
+        fig.canvas.draw()
         image = np.array(fig.canvas.renderer.buffer_rgba())
         plt.close(fig)
         output_picks = pd.DataFrame()
-        output_picks.to_csv('data/picks.csv', index=False)
-        output_csv = 'data/picks.csv'
+        output_picks.to_csv("data/picks.csv", index=False)
+        output_csv = "data/picks.csv"
         return image, output_picks, output_csv
-    
 
-    first_distances = bin_distances(distances, bin_size=10/111.2)
+    first_distances = bin_distances(distances, bin_size=10 / 111.2)
 
     # Edge case when there are way too many waveforms to process
-    selection_indexes = np.random.choice(first_distances, 
-                                         np.min([len(first_distances), max_waveforms]),
-                                         replace=False)
+    selection_indexes = np.random.choice(
+        first_distances, np.min([len(first_distances), max_waveforms]), replace=False
+    )
 
     waveforms = np.array(waveforms)[selection_indexes]
     distances = np.array(distances)[selection_indexes]
@@ -236,7 +328,7 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
 
     waveforms = [torch.tensor(waveform) for waveform in waveforms]
 
-    print('Starting to run predictions')
+    print("Starting to run predictions")
     with torch.no_grad():
         waveforms_torch = torch.vstack(waveforms)
         output = model(waveforms_torch)
@@ -244,129 +336,183 @@ def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source
     p_phases = output[:, 0]
     s_phases = output[:, 1]
 
-    p_phases = p_phases.reshape(len(waveforms),-1)
-    s_phases = s_phases.reshape(len(waveforms),-1)
+    p_phases = p_phases.reshape(len(waveforms), -1)
+    s_phases = s_phases.reshape(len(waveforms), -1)
 
-    # Max confidence - min variance    
+    # Max confidence - min variance
     p_max_confidence = p_phases.std(axis=-1).min()
     s_max_confidence = s_phases.std(axis=-1).min()
 
     print(f"Starting plotting {len(waveforms)} waveforms")
     fig, ax = plt.subplots(ncols=3, figsize=(10, 3))
-    
+
     # Plot topography
-    print('Fetching topography')
+    print("Fetching topography")
     params = Topography.DEFAULT.copy()
     extra_window = 0.5
-    params["south"] = np.min([st_lats.min(), eq_lat])-extra_window
-    params["north"] = np.max([st_lats.max(), eq_lat])+extra_window
-    params["west"] = np.min([st_lons.min(), eq_lon])-extra_window
-    params["east"] = np.max([st_lons.max(), eq_lon])+extra_window
+    params["south"] = np.min([st_lats.min(), eq_lat]) - extra_window
+    params["north"] = np.max([st_lats.max(), eq_lat]) + extra_window
+    params["west"] = np.min([st_lons.min(), eq_lon]) - extra_window
+    params["east"] = np.max([st_lons.max(), eq_lon]) + extra_window
 
     topo_map = Topography(**params)
     topo_map.fetch()
     topo_map.load()
 
-    print('Plotting topo')
+    print("Plotting topo")
     hillshade = es.hillshade(topo_map.da[0], altitude=10)
-    
-    topo_map.da.plot(ax = ax[1], cmap='Greys', add_colorbar=False, add_labels=False)
-    topo_map.da.plot(ax = ax[2], cmap='Greys', add_colorbar=False, add_labels=False)
+
+    topo_map.da.plot(ax=ax[1], cmap="Greys", add_colorbar=False, add_labels=False)
+    topo_map.da.plot(ax=ax[2], cmap="Greys", add_colorbar=False, add_labels=False)
     ax[1].imshow(hillshade, cmap="Greys", alpha=0.5)
 
-    output_picks = pd.DataFrame({'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' : []})
-                        
+    output_picks = pd.DataFrame(
+        {
+            "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": [],
+        }
+    )
+
     for i in range(len(waveforms)):
         print(f"Plotting waveform {i+1}/{len(waveforms)}")
         current_P = p_phases[i]
         current_S = s_phases[i]
-        
-        x = [t0s[i] + pd.Timedelta(seconds=k/100) for k in np.linspace(0,6000,6000)]
+
+        x = [t0s[i] + pd.Timedelta(seconds=k / 100) for k in np.linspace(0, 6000, 6000)]
         x = mdates.date2num(x)
 
         # 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()
+        p_conf = 1 / (current_P.std() / p_max_confidence).item()
+        s_conf = 1 / (current_S.std() / s_max_confidence).item()
 
         delta_t = t0s[i].timestamp - obspy.UTCDateTime(timestamp).timestamp
 
-        ax[0].plot(x, waveforms[i][0, 0]*10+distances[i]*111.2, color='black', alpha=0.5, lw=1)
+        ax[0].plot(
+            x,
+            waveforms[i][0, 0] * 10 + distances[i] * 111.2,
+            color="black",
+            alpha=0.5,
+            lw=1,
+        )
+
+        if (current_P.std().item() * 60 < conf_thres_P) or (
+            current_S.std().item() * 60 < conf_thres_S
+        ):
+            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="|",
+            )
 
-        if (current_P.std().item()*60 < conf_thres_P) or (current_S.std().item()*60 < conf_thres_S):
-            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='|')
-        
-            velocity_p = (distances[i]*111.2)/(delta_t+current_P.mean()*60).item()
-            velocity_s = (distances[i]*111.2)/(delta_t+current_S.mean()*60).item()
+            velocity_p = (distances[i] * 111.2) / (
+                delta_t + current_P.mean() * 60
+            ).item()
+            velocity_s = (distances[i] * 111.2) / (
+                delta_t + current_S.mean() * 60
+            ).item()
 
             # Generate an array from st_lat to eq_lat and from st_lon to eq_lon
             x = np.linspace(st_lons[i], eq_lon, 50)
             y = np.linspace(st_lats[i], eq_lat, 50)
-            
+
             # Plot the array
-            ax[1].scatter(x, y, c=np.zeros_like(x)+velocity_p, alpha=0.1, vmin=0, vmax=8)
-            ax[2].scatter(x, y, c=np.zeros_like(x)+velocity_s, alpha=0.1, vmin=0, vmax=8)
+            ax[1].scatter(
+                x, y, c=np.zeros_like(x) + velocity_p, alpha=0.1, vmin=0, vmax=8
+            )
+            ax[2].scatter(
+                x, y, c=np.zeros_like(x) + velocity_s, alpha=0.1, vmin=0, vmax=8
+            )
 
         else:
             velocity_p = np.nan
             velocity_s = np.nan
-        
-        ax[0].set_ylabel('Z')
-        print(f"Station {st_lats[i]}, {st_lons[i]} has P velocity {velocity_p} and S velocity {velocity_s}")
-
-        output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]], 
-                                                        'st_lat' : [st_lats[i]], 'st_lon' : [st_lons[i]],
-                                                        'starttime' : [str(t0s[i])], 
-                                                        'p_phase, s' : [(delta_t+current_P.mean()*60).item()], 'p_uncertainty, s' : [current_P.std().item()*60], 
-                                                        's_phase, s' : [(delta_t+current_S.mean()*60).item()], 's_uncertainty, s' : [current_S.std().item()*60],
-                                                        'velocity_p, km/s' : [velocity_p], 'velocity_s, km/s' : [velocity_s]}))
-        
-        
+
+        ax[0].set_ylabel("Z")
+        print(
+            f"Station {st_lats[i]}, {st_lons[i]} has P velocity {velocity_p} and S velocity {velocity_s}"
+        )
+
+        output_picks = output_picks.append(
+            pd.DataFrame(
+                {
+                    "station_name": [names[i]],
+                    "st_lat": [st_lats[i]],
+                    "st_lon": [st_lons[i]],
+                    "starttime": [str(t0s[i])],
+                    "p_phase, s": [(delta_t + current_P.mean() * 60).item()],
+                    "p_uncertainty, s": [current_P.std().item() * 60],
+                    "s_phase, s": [(delta_t + current_S.mean() * 60).item()],
+                    "s_uncertainty, s": [current_S.std().item() * 60],
+                    "velocity_p, km/s": [velocity_p],
+                    "velocity_s, km/s": [velocity_s],
+                }
+            )
+        )
+
     # Add legend
-    ax[0].scatter(None, None, color='r', marker='|', label='P')
-    ax[0].scatter(None, None, color='b', marker='|', label='S')
-    ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))
+    ax[0].scatter(None, None, color="r", marker="|", label="P")
+    ax[0].scatter(None, None, color="b", marker="|", label="S")
+    ax[0].xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
     ax[0].xaxis.set_major_locator(mdates.SecondLocator(interval=20))
     ax[0].legend()
 
-    print('Plotting stations')
-    for i in range(1,3):
-        ax[i].scatter(st_lons, st_lats, color='b', label='Stations')
-        ax[i].scatter(eq_lon, eq_lat, color='r', marker='*', label='Earthquake')
-        ax[i].set_aspect('equal')
-        ax[i].set_xticklabels(ax[i].get_xticks(), rotation = 50)
+    print("Plotting stations")
+    for i in range(1, 3):
+        ax[i].scatter(st_lons, st_lats, color="b", label="Stations")
+        ax[i].scatter(eq_lon, eq_lat, color="r", marker="*", label="Earthquake")
+        ax[i].set_aspect("equal")
+        ax[i].set_xticklabels(ax[i].get_xticks(), rotation=50)
+
+    fig.subplots_adjust(
+        bottom=0.1, top=0.9, left=0.1, right=0.8, wspace=0.02, hspace=0.02
+    )
 
-    fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
-                    wspace=0.02, hspace=0.02)
-    
     cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
-    cbar = fig.colorbar(ax[2].scatter(None, None, c=velocity_p, alpha=0.5, vmin=0, vmax=8), cax=cb_ax)
+    cbar = fig.colorbar(
+        ax[2].scatter(None, None, c=velocity_p, alpha=0.5, vmin=0, vmax=8), cax=cb_ax
+    )
 
-    cbar.set_label('Velocity (km/s)')
-    ax[1].set_title('P Velocity')
-    ax[2].set_title('S Velocity')
+    cbar.set_label("Velocity (km/s)")
+    ax[1].set_title("P Velocity")
+    ax[2].set_title("S Velocity")
 
     for a in ax:
-        a.tick_params(axis='both', which='major', labelsize=8)
-        
-    plt.subplots_adjust(hspace=0., wspace=0.5)
-    fig.canvas.draw();
+        a.tick_params(axis="both", which="major", labelsize=8)
+
+    plt.subplots_adjust(hspace=0.0, wspace=0.5)
+    fig.canvas.draw()
     image = np.array(fig.canvas.renderer.buffer_rgba())
     plt.close(fig)
-    output_picks.to_csv(f'data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv', index=False)
-    output_csv = f'data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv'
+    output_picks.to_csv(
+        f"data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv", index=False
+    )
+    output_csv = f"data/velocity/{eq_lat}_{eq_lon}_{timestamp}_{len(waveforms)}.csv"
 
     return image, output_picks, output_csv
 
+
 model = torch.jit.load("model.pt")
 
 with gr.Blocks() as demo:
-    gr.HTML("""
+    gr.HTML(
+        """
 <div style="padding: 20px; border-radius: 10px;">
     <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>
 
@@ -394,132 +540,195 @@ with gr.Blocks() as demo:
         <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>
     </ul>
 </div>
-""")
+"""
+    )
 
     with gr.Tab("Try on a single station"):
-        with gr.Row(): 
+        with gr.Row():
             # Define the input and output types for Gradio
             inputs = gr.Dropdown(
-                ["data/sample/sample_0.npy", 
-                "data/sample/sample_1.npy", 
-                "data/sample/sample_2.npy"], 
-                label="Sample waveform", 
+                [
+                    "data/sample/sample_0.npy",
+                    "data/sample/sample_1.npy",
+                    "data/sample/sample_2.npy",
+                ],
+                label="Sample waveform",
                 info="Select one of the samples",
-                value = "data/sample/sample_0.npy"
+                value="data/sample/sample_0.npy",
             )
+            with gr.Column(scale=1):
+                P_thres_inputs = gr.Slider(
+                    minimum=0.01,
+                    maximum=1,
+                    value=0.1,
+                    label="P uncertainty threshold, s",
+                    step=0.01,
+                    info="Acceptable uncertainty for P picks expressed in std() seconds",
+                    interactive=True,
+                )
+
+                S_thres_inputs = gr.Slider(
+                    minimum=0.01,
+                    maximum=1,
+                    value=0.2,
+                    label="S uncertainty threshold, s",
+                    step=0.01,
+                    info="Acceptable uncertainty for S picks expressed in std() seconds",
+                    interactive=True,
+                )
 
             upload = gr.File(label="Or upload your own waveform")
 
         button = gr.Button("Predict phases")
-        outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)
-    
-        button.click(mark_phases, inputs=[inputs, upload], outputs=outputs)
-        
+        outputs = gr.Image(
+            label="Waveform with Phases Marked", type="numpy", interactive=False
+        )
+
+        button.click(
+            mark_phases,
+            inputs=[inputs, upload, P_thres_inputs, S_thres_inputs],
+            outputs=outputs,
+        )
+
     with gr.Tab("Select earthquake from catalogue"):
 
-        gr.HTML("""
+        gr.HTML(
+            """
         <div style="padding: 20px; border-radius: 10px; font-size: 16px;">
         <p style="font-weight: bold; font-size: 24px; margin-bottom: 20px;">Using PhaseHunter to Analyze Seismic Waveforms</p>
         <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>
         <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>
         <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>
         </div>
-        """)
-        with gr.Row(): 
+        """
+        )
+        with gr.Row():
             with gr.Column(scale=2):
                 client_inputs = gr.Dropdown(
-                    choices = list(URL_MAPPINGS.keys()), 
-                    label="FDSN Client", 
+                    choices=list(URL_MAPPINGS.keys()),
+                    label="FDSN Client",
                     info="Select one of the available FDSN clients",
-                    value = "IRIS",
-                    interactive=True
+                    value="IRIS",
+                    interactive=True,
                 )
 
                 velocity_inputs = gr.Dropdown(
-                    choices = ['1066a', '1066b', 'ak135', 
-                            'ak135f', 'herrin', 'iasp91', 
-                            'jb', 'prem', 'pwdk'], 
-                    label="1D velocity model", 
+                    choices=[
+                        "1066a",
+                        "1066b",
+                        "ak135",
+                        "ak135f",
+                        "herrin",
+                        "iasp91",
+                        "jb",
+                        "prem",
+                        "pwdk",
+                    ],
+                    label="1D velocity model",
                     info="Velocity model for station selection",
-                    value = "1066a",
-                    interactive=True
+                    value="1066a",
+                    interactive=True,
                 )
 
             with gr.Column(scale=2):
-                timestamp_inputs = gr.Textbox(value='2019-07-04 17:33:49',
-                                    placeholder='YYYY-MM-DD HH:MM:SS',
-                                    label="Timestamp",
-                                    info="Timestamp of the earthquake",
-                                    max_lines=1,
-                                    interactive=True)
-                
-                source_depth_inputs = gr.Number(value=10,
+                timestamp_inputs = gr.Textbox(
+                    value="2019-07-04 17:33:49",
+                    placeholder="YYYY-MM-DD HH:MM:SS",
+                    label="Timestamp",
+                    info="Timestamp of the earthquake",
+                    max_lines=1,
+                    interactive=True,
+                )
+
+                source_depth_inputs = gr.Number(
+                    value=10,
                     label="Source depth (km)",
                     info="Depth of the earthquake",
-                    interactive=True)
-                
+                    interactive=True,
+                )
+
             with gr.Column(scale=2):
-                eq_lat_inputs = gr.Number(value=35.766, 
-                                label="Latitude", 
-                                info="Latitude of the earthquake",
-                                interactive=True)
-                
-                eq_lon_inputs = gr.Number(value=-117.605,
-                                    label="Longitude",
-                                    info="Longitude of the earthquake",
-                                    interactive=True)
-                
+                eq_lat_inputs = gr.Number(
+                    value=35.766,
+                    label="Latitude",
+                    info="Latitude of the earthquake",
+                    interactive=True,
+                )
+
+                eq_lon_inputs = gr.Number(
+                    value=-117.605,
+                    label="Longitude",
+                    info="Longitude of the earthquake",
+                    interactive=True,
+                )
+
             with gr.Column(scale=2):
-                radius_inputs = gr.Slider(minimum=1, 
-                                        maximum=200, 
-                                        value=50, 
-                                        label="Radius (km)", 
-                                        step=10,
-                                        info="""Select the radius around the earthquake to download data from.\n 
+                radius_inputs = gr.Slider(
+                    minimum=1,
+                    maximum=200,
+                    value=50,
+                    label="Radius (km)",
+                    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)
-                
-                max_waveforms_inputs = gr.Slider(minimum=1,
-                                maximum=100,
-                                value=10,
-                                label="Max waveforms per section",
-                                step=1,
-                                info="Maximum number of waveforms to show per section\n (to avoid long prediction times)",
-                                interactive=True,
-                                )
+                    interactive=True,
+                )
+
+                max_waveforms_inputs = gr.Slider(
+                    minimum=1,
+                    maximum=100,
+                    value=10,
+                    label="Max waveforms per section",
+                    step=1,
+                    info="Maximum number of waveforms to show per section\n (to avoid long prediction times)",
+                    interactive=True,
+                )
             with gr.Column(scale=2):
-                P_thres_inputs = gr.Slider(minimum=0.01,
-                                maximum=1,
-                                value=0.1,
-                                label="P uncertainty threshold, s",
-                                step=0.01,
-                                info="Acceptable uncertainty for P picks expressed in std() seconds",
-                                interactive=True,
-                                )
-                S_thres_inputs = gr.Slider(minimum=0.01,
-                                maximum=1,
-                                value=0.2,
-                                label="S uncertainty threshold, s",
-                                step=0.01,
-                                info="Acceptable uncertainty for S picks expressed in std() seconds",
-                                interactive=True,
-                                )
-            
+                P_thres_inputs = gr.Slider(
+                    minimum=0.01,
+                    maximum=1,
+                    value=0.1,
+                    label="P uncertainty threshold, s",
+                    step=0.01,
+                    info="Acceptable uncertainty for P picks expressed in std() seconds",
+                    interactive=True,
+                )
+                S_thres_inputs = gr.Slider(
+                    minimum=0.01,
+                    maximum=1,
+                    value=0.2,
+                    label="S uncertainty threshold, s",
+                    step=0.01,
+                    info="Acceptable uncertainty for S picks expressed in std() seconds",
+                    interactive=True,
+                )
+
         button = gr.Button("Predict phases")
-        output_image = gr.Image(label='Waveforms with Phases Marked', type='numpy', interactive=False)
+        output_image = gr.Image(
+            label="Waveforms with Phases Marked", type="numpy", interactive=False
+        )
 
         with gr.Row():
-            output_picks = gr.Dataframe(label='Pick data', 
-                                        type='pandas', 
-                                        interactive=False)
+            output_picks = gr.Dataframe(
+                label="Pick data", type="pandas", interactive=False
+            )
             output_csv = gr.File(label="Output File", file_types=[".csv"])
 
-        button.click(predict_on_section, 
-                 inputs=[client_inputs, timestamp_inputs, 
-                         eq_lat_inputs, eq_lon_inputs, 
-                         radius_inputs, source_depth_inputs, 
-                         velocity_inputs, max_waveforms_inputs,
-                         P_thres_inputs, S_thres_inputs],
-                 outputs=[output_image, output_picks, output_csv])
-
-demo.launch()
+        button.click(
+            predict_on_section,
+            inputs=[
+                client_inputs,
+                timestamp_inputs,
+                eq_lat_inputs,
+                eq_lon_inputs,
+                radius_inputs,
+                source_depth_inputs,
+                velocity_inputs,
+                max_waveforms_inputs,
+                P_thres_inputs,
+                S_thres_inputs,
+            ],
+            outputs=[output_image, output_picks, output_csv],
+        )
+
+demo.launch()

data/velocity/35.766_-117.605_2019-07-04 17:33:49_3.csv

@@ -0,0 +1,4 @@
+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"
+CI.SRT,35.69235,-117.75051,2019-07-04T17:33:38.029990Z,4.52931022644043,0.014338254695758224,9.210612297058105,0.013868825335521251,3.417590792273917,1.6805971661817796
+CI.JRC2,35.98249,-117.80885,2019-07-04T17:33:39.947494Z,7.320904731750488,0.018777401419356465,13.390795707702637,0.037158007035031915,4.136213094741051,2.261316106809097
+CI.WMF,36.11758,-117.85486,2019-07-04T17:33:41.867962Z,9.504384994506836,0.01592034415807575,17.031570434570312,0.04673818708397448,4.745724852828504,2.6483286583912338

data/velocity/35.766_-117.605_2019-07-04 17:33:49_9.csv

@@ -0,0 +1,10 @@
+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"
+LB.DAC,36.277,-117.593697,2019-07-04T17:33:43.387184Z,34.49410629272461,0.22159690037369728,41.63465881347656,0.886770598590374,,
+CI.WMF,36.11758,-117.85486,2019-07-04T17:33:41.867962Z,9.505252838134766,0.01779856625944376,17.014007568359375,0.06349349627271295,4.7452915611377335,2.6510624200710167
+CI.ISA,35.66278,-118.47403,2019-07-04T17:33:46.297658Z,14.401352882385254,0.04201323492452502,26.648784637451172,0.040058575104922056,5.506237698660711,2.9756431008588833
+CI.JRC2,35.98249,-117.80885,2019-07-04T17:33:39.947494Z,7.321362495422363,0.018879775889217854,13.384379386901855,0.05111166858114302,4.135954480569838,2.2624001562934875
+CI.EDW2,34.8811,-117.993881,2019-07-04T17:33:49.567241Z,17.187454223632812,0.025604498223401606,31.232393264770508,0.16874447464942932,6.08203947635048,3.3469985537104074
+CI.RRX,34.875332,-116.996841,2019-07-04T17:33:50.712219Z,17.309837341308594,0.04084662417881191,30.36077117919922,0.091552734375,6.549756328808083,3.734266695918123
+CI.SRT,35.69235,-117.75051,2019-07-04T17:33:38.029990Z,4.529020309448242,0.010629930475261062,9.221219062805176,0.019179217633791268,3.4178095631283902,1.6786640486259041
+CI.CWC,36.439049,-118.080498,2019-07-04T17:33:47.189005Z,16.294597625732422,0.03656624467112124,28.510168075561523,0.047869981499388814,5.288708979140877,3.02268946806275
+NN.QSM,35.965,-116.869102,2019-07-04T17:33:45.081547Z,10.827579498291016,0.05653449450619519,19.081382751464844,0.02382224891334772,6.456704830678366,3.663806012475838

test.npy

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:32166b7b4d6bd898ebce5da148b1080ffd3f50682a34ec877c14a905aa441d91
+size 48128