Figure7.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from obspy import read_events\n",
    "from obspy import read_inventory\n",
    "from obspy.clients.fdsn import Client\n",
    "client_wm = Client(\"IRIS\")\n",
    "from obspy.clients.iris import Client  \n",
    "from obspy.taup import TauPyModel\n",
    "from obspy import Stream, Trace\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from obspy import UTCDateTime\n",
    "import matplotlib.pyplot as plt\n",
    "from obspy.geodetics import locations2degrees,degrees2kilometers\n",
    "from scipy import signal\n",
    "plt.rc('font', size=10) \n",
    "%matplotlib widget\n",
    "import datetime\n",
    "from obspy.signal.trigger import recursive_sta_lta,plot_trigger,coincidence_trigger\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv0 = client_wm.get_stations(\n",
    "    channel=\"BH*\", starttime=UTCDateTime(2020,9,1,0,0,0),endtime=UTCDateTime(2020,10,1,0,0,0),\n",
    "    network=\"AK\", sta=\"SCRK,J25K,K27K\",\n",
    "    level=\"response\")\n",
    "the_station=inv0[0][2]\n",
    "evid_dum='https://earthquake.usgs.gov/fdsnws/event/1/query?&starttime=2020-09-01&endtime=2020-9-10&latitude=63.98&longitude=-143.99&maxradius=4&minmagnitude=0'\n",
    "cat1 = read_events(evid_dum)\n",
    "evid_dum='https://earthquake.usgs.gov/fdsnws/event/1/query?&starttime=2020-09-01&endtime=2020-9-10&latitude=63.98&longitude=-143.99&minradius=4&maxradius=8&minmagnitude=1'\n",
    "cat2 = read_events(evid_dum)\n",
    "evid_dum='https://earthquake.usgs.gov/fdsnws/event/1/query?&starttime=2020-09-01&endtime=2020-9-10&latitude=63.98&longitude=-143.99&minradius=8&maxradius=20&minmagnitude=2'\n",
    "cat3 = read_events(evid_dum)\n",
    "evid_dum='https://earthquake.usgs.gov/fdsnws/event/1/query?&starttime=2020-09-01&endtime=2020-9-10&latitude=63.98&longitude=-143.99&minradius=20&maxradius=100&minmagnitude=3'\n",
    "cat4 = read_events(evid_dum)\n",
    "evid_dum='https://earthquake.usgs.gov/fdsnws/event/1/query?&starttime=2020-09-01&endtime=2020-9-10&latitude=63.98&longitude=-143.99&minradius=100&maxradius=180&minmagnitude=4'\n",
    "cat5 = read_events(evid_dum)\n",
    "cat_all=cat1+cat2+cat3+cat4+cat5\n",
    "channel = \"BHZ\"\n",
    "start_time=UTCDateTime(2020,9,1,0,0,0)\n",
    "end_time=UTCDateTime(2020,9,10,0,0,0)\n",
    "#st_2day = client_wm.get_waveforms(\"AK\", \"SCRK,K27K,J25K\", \"*\", \"BHZ\", start_time, end_time)\n",
    "st_10day = client_wm.get_waveforms(\"AK\", \"SCRK\", \"*\", \"BHZ\", start_time, end_time)\n",
    "EQ_data=pd.read_csv(\"detections/SCRK_outputs/X_prediction_results.csv\")\n",
    "eq_arv_t=EQ_data['event_start_time']\n",
    "eq_arv_s_t=EQ_data['s_arrival_time']\n",
    "st_10day2=st_10day.copy()\n",
    "st_10day2.merge(method=0,fill_value=0)\n",
    "st_10day2.detrend()\n",
    "st_10day2.filter('bandpass',freqmin=0.8, freqmax=25.0)\n",
    "trace=st_10day2[0]\n",
    "model = TauPyModel(model=\"prem\")\n",
    "df = trace.stats.sampling_rate\n",
    "npts = trace.stats.npts\n",
    "t = np.arange(npts, dtype=np.float32) / df\n",
    "cft = recursive_sta_lta(trace.data, int(5 * df), int(10 * df))\n",
    "eqt_array=[]\n",
    "for eqtm in eq_arv_t:\n",
    "    dumt=UTCDateTime(eqtm)\n",
    "    #ax.vlines(dumt-trace.stats.starttime, ii,jj, color='y',linestyles='dashdot', lw=3)\n",
    "    eqt_array.append(dumt-trace.stats.starttime)\n",
    "    picker_time_a=[]\n",
    "for event in cat_all:\n",
    "    disD=locations2degrees(the_station.latitude,the_station.longitude,event.origins[0].latitude,event.origins[0].longitude)\n",
    "    if event.origins[0].depth<0:\n",
    "        event.origins[0].depth=0\n",
    "    arrivals = model.get_travel_times(source_depth_in_km=event.origins[0].depth/1000,\n",
    "                                  distance_in_degree=disD)    \n",
    "    picker_time_a.append(event.origins[0].time+arrivals[0].time-trace.stats.starttime)\n",
    "\n",
    "max_search_range_b=10\n",
    "max_search_range_a=30\n",
    "cft_thr=1.35\n",
    "\n",
    "sel_cat_list=[]\n",
    "\n",
    "lij=-1\n",
    "for event in cat_all:\n",
    "    lij=lij+1\n",
    "    disD=locations2degrees(the_station.latitude,the_station.longitude,event.origins[0].latitude,event.origins[0].longitude)\n",
    "    picker_time=picker_time_a[lij]\n",
    "    t0=int(picker_time/trace.stats.delta)\n",
    "    tarray=range(int(t0-max_search_range_b*trace.stats.sampling_rate),int(t0+max_search_range_a*trace.stats.sampling_rate))\n",
    "    try:\n",
    "        if np.max(cft[tarray])>cft_thr:\n",
    "            sel_cat_list.append(event)\n",
    "            # if flag==0:\n",
    "            #     cat_all_sel=event\n",
    "            #     flag=1\n",
    "            # else:\n",
    "            #     cat_all_sel=cat_all_sel+event\n",
    "    except Exception:\n",
    "        continue\n",
    "    \n",
    "for i,event in enumerate(sel_cat_list):\n",
    "    if i==0:\n",
    "        dum=cat_all[0:1]\n",
    "        dum.clear()\n",
    "        all_sel_event=dum+event\n",
    "    else:\n",
    "        all_sel_event=all_sel_event+event\n",
    "lij=-1\n",
    "eqt_max_time=25\n",
    "for i,picker_time in enumerate(picker_time_a):\n",
    "    if np.abs(np.array(eqt_array)-picker_time).min()<eqt_max_time:\n",
    "        lij=lij+1\n",
    "        if lij==0:\n",
    "            dum=cat_all[0:1]\n",
    "            dum.clear()\n",
    "            eqt_sel_event=dum+cat_all[i]\n",
    "        else:\n",
    "            eqt_sel_event=eqt_sel_event+cat_all[i]\n",
    "eqtpall_sel_event=all_sel_event+eqt_sel_event\n",
    "event_in_list=[]\n",
    "for event in eqtpall_sel_event:\n",
    "    disD=locations2degrees(the_station.latitude,the_station.longitude,event.origins[0].latitude,event.origins[0].longitude)\n",
    "    event_list=[event.origins[0].longitude,event.origins[0].latitude,event.origins[0].depth,event.magnitudes[0].mag,disD,event.event_type,event.origins[0].time]\n",
    "    event_in_list.append(event_list)\n",
    "eqtpall_sel_event_pd=pd.DataFrame(event_in_list,columns=['longitude','latitude','depth','magnitude','disD','event_type','time'])\n",
    "bh_catalog=eqtpall_sel_event_pd.drop_duplicates(subset=['longitude','latitude','magnitude'])\n",
    "\n",
    "def threshold_trigger(th,record_time,trace):\n",
    "    record_number=record_time*trace.stats.sampling_rate\n",
    "    abs_st_data=np.abs(trace.data)\n",
    "    ind_th=[]\n",
    "    recording=False\n",
    "    for i,abs_st in enumerate(abs_st_data):\n",
    "        if abs_st>th and recording is False:\n",
    "            ind_th.append(i)\n",
    "            recording=True\n",
    "        if ind_th!=[] and (i-ind_th[-1])>record_number:\n",
    "            recording=False\n",
    "    return len(ind_th),ind_th\n",
    "\n",
    "def cal_bhcat_t(bh_catalog):\n",
    "    cat_arl_t=[]\n",
    "    cat_arl_i=[]\n",
    "    for i in bh_catalog.index:\n",
    "        arrivals = model.get_travel_times(source_depth_in_km=bh_catalog.loc[i,'depth']/1000,\n",
    "                                    distance_in_degree=bh_catalog.loc[i,'disD'])    \n",
    "        picker_time=bh_catalog.loc[i,'time']+arrivals[0].time-trace.stats.starttime\n",
    "        cat_arl_t.append(picker_time)\n",
    "        cat_arl_i.append(i)\n",
    "    return cat_arl_t,cat_arl_i\n",
    "\n",
    "\n",
    "def event_associate(pre_time,after_time,cat_arl_t,t,ind_th,cat_arl_i):\n",
    "    n_recorded=0\n",
    "    index_triggered=[]\n",
    "    index_cat=[]\n",
    "    for i in ind_th:\n",
    "        for j,cat_t in enumerate(cat_arl_t):\n",
    "            if cat_t>t[i]-pre_time and cat_t<t[i]+after_time:\n",
    "                n_recorded=n_recorded+1\n",
    "                index_triggered.append(cat_arl_i[j])\n",
    "                index_cat.append(j)\n",
    "                break\n",
    "    return n_recorded,index_triggered,index_cat\n",
    "def threshold_trigger_w_avg(th_fold,record_time,trace,window_len):\n",
    "    record_number=record_time*trace.stats.sampling_rate\n",
    "    abs_st_data=np.abs(trace.data)\n",
    "    ind_th=[]\n",
    "    recording=False\n",
    "    rtt=0\n",
    "    cm_sum=0\n",
    "    ltavg=np.sum(abs_st_data[0:window_len])/window_len\n",
    "    for i,abs_st in enumerate(abs_st_data):\n",
    "        rtt=rtt+1\n",
    "        if rtt>window_len:\n",
    "            rtt=0\n",
    "            ltavg=cm_sum/window_len\n",
    "            #print(ltavg)\n",
    "            cm_sum=abs_st\n",
    "        cm_sum=cm_sum+abs_st\n",
    "        if abs_st>th_fold*ltavg and recording is False:\n",
    "            ind_th.append(i)\n",
    "            recording=True\n",
    "        if ind_th!=[] and (i-ind_th[-1])>record_number:\n",
    "            recording=False\n",
    "    return len(ind_th),ind_th\n",
    "def cal_tauc(tr,sat_t,end_t):\n",
    "    df = tr.stats.sampling_rate\n",
    "    npts = tr.stats.npts\n",
    "    t = np.arange(npts, dtype=np.float32) / df\n",
    "    ind_range=(t<end_t) * (t>sat_t)\n",
    "    r=np.sum(np.diff(tr[ind_range])**2/tr.stats.delta)/np.sum(tr[ind_range]**2)\n",
    "    return 2*np.pi/np.sqrt(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat_arl_t,cat_arl_i=cal_bhcat_t(bh_catalog)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pre_time=80\n",
    "after_time=10\n",
    "th_fold=200\n",
    "th=200\n",
    "record_time=60*30\n",
    "fminay=np.arange(0.02,10,0.4)\n",
    "window_len=50*60*2\n",
    "freq_r_mat=np.zeros([25,75])\n",
    "freq_r_mat_ww=np.zeros([25,75])\n",
    "freq_r_r_mat_ww=np.zeros([25,75])\n",
    "freq_r_r_mat=np.zeros([25,75])\n",
    "for i,fmin in enumerate(fminay):\n",
    "    fmaxay=np.arange(fmin+0.4,30,0.4)\n",
    "    for j,fmax in enumerate(fmaxay):\n",
    "        print(fmin,fmax)\n",
    "        st_10day2=st_10day.copy()\n",
    "        st_10day2.merge(method=0,fill_value=0)\n",
    "        st_10day2.detrend()\n",
    "        st_10day2.filter('bandpass',freqmin=fmin, freqmax=fmax)\n",
    "        trace=st_10day2[0]\n",
    "\n",
    "        n_dum,ind_th=threshold_trigger(th,record_time,trace)\n",
    "        freq_r_mat[i][i+j]=n_dum\n",
    "        n_recorded,index_triggered,index_cat1=event_associate(pre_time,after_time,cat_arl_t,t,ind_th,cat_arl_i)\n",
    "        try:\n",
    "            freq_r_r_mat[i][i+j]=n_recorded/n_dum\n",
    "        except:\n",
    "            freq_r_r_mat[i][i+j]=0\n",
    "        n_dum,ind_th=threshold_trigger_w_avg(th_fold,record_time,trace,window_len)\n",
    "        freq_r_mat_ww[i][i+j]=n_dum\n",
    "        n_recorded,index_triggered,index_cat2=event_associate(pre_time,after_time,cat_arl_t,t,ind_th,cat_arl_i)\n",
    "        try:\n",
    "            freq_r_r_mat_ww[i][i+j]=n_recorded/n_dum\n",
    "        except:\n",
    "            freq_r_r_mat_ww[i][i+j]=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('font', size=15) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure(figsize=(9, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "#c = ax.contourf(range(60*5,60*60*2,60*5),range(100,1000,50),tri_matrix, 20,cmap='jet',norm=LogNorm(vmin=np.array(tri_matrix).min(), vmax=np.array(tri_matrix).max()))\n",
    "c= ax.pcolor(np.arange(0.02,30,0.4),np.arange(0.42,10,0.4),freq_r_mat[:][1:],cmap='jet')\n",
    "\n",
    "#ax.clabel(c, fmt='%2.1f', colors='w', fontsize=14)\n",
    "ax.set_title('Threshold method number of triggers')\n",
    "c.set_clim(0, 20)\n",
    "cbar=fig.colorbar(c, ax=ax)\n",
    "ax.set_xlabel('max frequency (Hz)')\n",
    "ax.set_ylabel('min frequency (Hz)')\n",
    "cbar.ax.set_ylabel('number of triggers')\n",
    "ax.set_xlim(0,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure(figsize=(9, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "#c = ax.contourf(range(60*5,60*60*2,60*5),range(100,1000,50),tri_matrix, 20,cmap='jet',norm=LogNorm(vmin=np.array(tri_matrix).min(), vmax=np.array(tri_matrix).max()))\n",
    "c= ax.pcolor(np.arange(0.02,30,0.4),np.arange(0.42,10,0.4),freq_r_mat_ww[:][1:],cmap='jet')\n",
    "\n",
    "#ax.clabel(c, fmt='%2.1f', colors='w', fontsize=14)\n",
    "ax.set_title('Segmented window method number of triggers')\n",
    "c.set_clim(0, 20)\n",
    "cbar=fig.colorbar(c, ax=ax)\n",
    "ax.set_xlabel('max frequency (Hz)')\n",
    "ax.set_ylabel('min frequency (Hz)')\n",
    "cbar.ax.set_ylabel('number of triggers')\n",
    "ax.set_xlim(0,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure(figsize=(9, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "#c = ax.contourf(range(60*5,60*60*2,60*5),range(100,1000,50),tri_matrix, 20,cmap='jet',norm=LogNorm(vmin=np.array(tri_matrix).min(), vmax=np.array(tri_matrix).max()))\n",
    "c= ax.pcolor(np.arange(0.02,30,0.4),np.arange(0.42,10,0.4),freq_r_r_mat[:][1:],cmap='jet')\n",
    "\n",
    "#ax.clabel(c, fmt='%2.1f', colors='w', fontsize=14)\n",
    "ax.set_title('Threshold method correct ratio')\n",
    "cbar=fig.colorbar(c, ax=ax)\n",
    "c.set_clim(0, 1)\n",
    "ax.set_xlabel('max frequency (Hz)')\n",
    "ax.set_ylabel('min frequency (Hz)')\n",
    "cbar.ax.set_ylabel('correct detection ratio')\n",
    "ax.set_xlim(0,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure(figsize=(9, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "#c = ax.contourf(range(60*5,60*60*2,60*5),range(100,1000,50),tri_matrix, 20,cmap='jet',norm=LogNorm(vmin=np.array(tri_matrix).min(), vmax=np.array(tri_matrix).max()))\n",
    "c= ax.pcolor(np.arange(0.02,30,0.4),np.arange(0.42,10,0.4),freq_r_r_mat_ww[:][1:],cmap='jet')\n",
    "\n",
    "#ax.clabel(c, fmt='%2.1f', colors='w', fontsize=14)\n",
    "ax.set_title('Segmented window method correct ratio')\n",
    "cbar=fig.colorbar(c, ax=ax)\n",
    "c.set_clim(0, 1)\n",
    "ax.set_xlabel('max frequency (Hz)')\n",
    "ax.set_ylabel('min frequency (Hz)')\n",
    "cbar.ax.set_ylabel('correct detection ratio')\n",
    "ax.set_xlim(0,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import savemat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savemat('freq_range_data.mat',{'freq_r_r_mat_ww':freq_r_r_mat_ww,'freq_r_r_mat':freq_r_r_mat,'freq_r_mat_ww':freq_r_mat_ww,'freq_r_mat':freq_r_mat})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bh_catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "c71a260edab3fcf64a811359c50653f22ebe325a9efcb1a4563b0bac5f87d06d"
  },
  "kernelspec": {
   "display_name": "Python 3.7.3 64-bit ('obspy': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}