Figure8_for_cpp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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.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\n",
    "from scipy.stats import linregress\n",
    "import ctypes\n",
    "from mpi4py import MPI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/yuantian/anaconda3/envs/obspy/lib/python3.7/site-packages/obspy/signal/filter.py:67: UserWarning: Selected high corner frequency (25.0) of bandpass is at or above Nyquist (25.0). Applying a high-pass instead.\n",
      "  warnings.warn(msg)\n"
     ]
    }
   ],
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "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'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>longitude</th>\n",
       "      <th>latitude</th>\n",
       "      <th>depth</th>\n",
       "      <th>magnitude</th>\n",
       "      <th>disD</th>\n",
       "      <th>event_type</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-146.8571</td>\n",
       "      <td>62.8881</td>\n",
       "      <td>44100.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.681260</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-09T19:12:06.535000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-146.8485</td>\n",
       "      <td>62.8784</td>\n",
       "      <td>45600.0</td>\n",
       "      <td>2.5</td>\n",
       "      <td>1.684795</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-09T19:02:13.756000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-147.7625</td>\n",
       "      <td>63.5276</td>\n",
       "      <td>10800.0</td>\n",
       "      <td>0.8</td>\n",
       "      <td>1.727153</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-09T15:29:05.744000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-148.3746</td>\n",
       "      <td>63.7534</td>\n",
       "      <td>84600.0</td>\n",
       "      <td>1.2</td>\n",
       "      <td>1.943564</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-09T11:12:10.065000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-149.0041</td>\n",
       "      <td>62.7214</td>\n",
       "      <td>62500.0</td>\n",
       "      <td>2.4</td>\n",
       "      <td>2.574220</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-09T10:00:10.896000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>-71.3062</td>\n",
       "      <td>-27.9686</td>\n",
       "      <td>21000.0</td>\n",
       "      <td>6.8</td>\n",
       "      <td>107.824510</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-01T04:09:28.470000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>204</th>\n",
       "      <td>-139.7911</td>\n",
       "      <td>64.8705</td>\n",
       "      <td>7600.0</td>\n",
       "      <td>1.2</td>\n",
       "      <td>2.021049</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-09T02:02:00.010000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>206</th>\n",
       "      <td>-147.7667</td>\n",
       "      <td>64.7974</td>\n",
       "      <td>11700.0</td>\n",
       "      <td>0.9</td>\n",
       "      <td>1.827002</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-07T21:47:51.243000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>217</th>\n",
       "      <td>-141.8599</td>\n",
       "      <td>65.8902</td>\n",
       "      <td>10800.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.115992</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-04T23:44:01.585000Z</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>234</th>\n",
       "      <td>-152.0282</td>\n",
       "      <td>61.0431</td>\n",
       "      <td>97100.0</td>\n",
       "      <td>1.9</td>\n",
       "      <td>4.723774</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>2020-09-08T11:59:14.009000Z</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>204 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     longitude  latitude    depth  magnitude        disD  event_type  \\\n",
       "0    -146.8571   62.8881  44100.0        2.0    1.681260  earthquake   \n",
       "1    -146.8485   62.8784  45600.0        2.5    1.684795  earthquake   \n",
       "2    -147.7625   63.5276  10800.0        0.8    1.727153  earthquake   \n",
       "3    -148.3746   63.7534  84600.0        1.2    1.943564  earthquake   \n",
       "4    -149.0041   62.7214  62500.0        2.4    2.574220  earthquake   \n",
       "..         ...       ...      ...        ...         ...         ...   \n",
       "199   -71.3062  -27.9686  21000.0        6.8  107.824510  earthquake   \n",
       "204  -139.7911   64.8705   7600.0        1.2    2.021049  earthquake   \n",
       "206  -147.7667   64.7974  11700.0        0.9    1.827002  earthquake   \n",
       "217  -141.8599   65.8902  10800.0        2.0    2.115992  earthquake   \n",
       "234  -152.0282   61.0431  97100.0        1.9    4.723774  earthquake   \n",
       "\n",
       "                            time  \n",
       "0    2020-09-09T19:12:06.535000Z  \n",
       "1    2020-09-09T19:02:13.756000Z  \n",
       "2    2020-09-09T15:29:05.744000Z  \n",
       "3    2020-09-09T11:12:10.065000Z  \n",
       "4    2020-09-09T10:00:10.896000Z  \n",
       "..                           ...  \n",
       "199  2020-09-01T04:09:28.470000Z  \n",
       "204  2020-09-09T02:02:00.010000Z  \n",
       "206  2020-09-07T21:47:51.243000Z  \n",
       "217  2020-09-04T23:44:01.585000Z  \n",
       "234  2020-09-08T11:59:14.009000Z  \n",
       "\n",
       "[204 rows x 7 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bh_catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)\n",
    "def cal_tauc_cpp(data, df, sat_t, end_t):\n",
    "    npts = len(data)\n",
    "    t = np.arange(npts, dtype=np.float32) / df\n",
    "    ind_range = (t < end_t) * (t > sat_t)\n",
    "    r_num = np.sum(np.diff(data[ind_range])**2 * df)\n",
    "    r_den = np.sum(data[ind_range]**2)\n",
    "    r = r_num / r_den\n",
    "    return 2 * np.pi / np.sqrt(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bh_catalog_nd=bh_catalog.reset_index()\n",
    "#trace=st_10day2[0]\n",
    "event_array=[]\n",
    "event_arlt_array=[]\n",
    "pre_min=0.5\n",
    "after_min=5\n",
    "for i in range(len(bh_catalog_nd)):\n",
    "\n",
    "    event=bh_catalog_nd.loc[i]\n",
    "    arrivals = model.get_travel_times(source_depth_in_km=event.depth/1000,distance_in_degree=event.disD)  \n",
    "    arl_time=event.time+arrivals[0].time\n",
    "    #arl_time=event_arlt_array[i]\n",
    "    trace_slc=trace.slice(arl_time-pre_min*60,arl_time+after_min*60)\n",
    "    event_array.append(trace_slc)\n",
    "    event_arlt_array.append(arl_time)\n",
    "pre_time=80\n",
    "after_time=10\n",
    "th_fold=10\n",
    "th_value=100\n",
    "window_len=50*60*1\n",
    "rest_time=30*60\n",
    "event_sel_array=[]\n",
    "event_arlt_sel_array=[]\n",
    "ind_slt=[]\n",
    "for i in range(len(bh_catalog_nd)):\n",
    "#for i in range(1):\n",
    "    event=bh_catalog_nd.loc[i]\n",
    "    tr=event_array[i]\n",
    "    df = tr.stats.sampling_rate\n",
    "    npts = tr.stats.npts\n",
    "    t = np.arange(npts, dtype=np.float32) / df\n",
    "    arl_time=event_arlt_array[i]\n",
    "    #n,ind_th=threshold_trigger_w_avg(th_fold,rest_time,tr,window_len)\n",
    "    n,ind_th=threshold_trigger(th_value,rest_time,tr)\n",
    "    for k in ind_th:\n",
    "        if arl_time-tr.stats.starttime>t[k]-pre_time and arl_time-tr.stats.starttime<t[k]+after_time:\n",
    "            event_sel_array.append(tr)\n",
    "            event_arlt_sel_array.append(arl_time)\n",
    "            ind_slt.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for detected events\n",
    "sat_t=20\n",
    "end_t=45\n",
    "tauc_over_events=[]\n",
    "mag_array=np.array(bh_catalog_nd.loc[ind_slt]['magnitude'])\n",
    "for i,trc in enumerate(event_sel_array):\n",
    "    #event=bh_cat_sort_mag_sel.loc[i]\n",
    "    tauc_over_events.append(cal_tauc(trc,sat_t,end_t))\n",
    "A = np.vstack([mag_array, np.ones(len(mag_array))]).T\n",
    "logtauc=np.log(tauc_over_events)\n",
    "m, c = np.linalg.lstsq(A, logtauc, rcond=None)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "#To run the compiled cpp cal_tauc, run: 'g++ -std=c++11 -shared -fPIC  -o cal_tauc cal_tauc.cpp' in the terminal\n",
    "# Define the start and end times for the tauc calculation\n",
    "sat_t = 20\n",
    "end_t = 45\n",
    "# Load the shared library\n",
    "lib = ctypes.CDLL('./cal_tauc_cpp.so')\n",
    "\n",
    "# Define the function arguments and return type\n",
    "lib.cal_tauc_cpp.argtypes = [np.ctypeslib.ndpointer(dtype=np.float32), ctypes.c_int, ctypes.c_float, ctypes.c_float, ctypes.c_float]\n",
    "lib.cal_tauc_cpp.restype = ctypes.c_float\n",
    "# Create an empty list to store the results\n",
    "tauc_over_events2 = []\n",
    "mag_array=np.array(bh_catalog_nd.loc[ind_slt]['magnitude'])\n",
    "# Loop over the event_sel_array and calculate the tauc for each trace\n",
    "for i, trc in enumerate(event_sel_array):\n",
    "    # Extract the data and sample rate from the trace\n",
    "    data = trc.data.astype(np.float32)\n",
    "    df = trc.stats.sampling_rate\n",
    "\n",
    "    # Calculate the tauc using the cal_tauc function\n",
    "    tauc = lib.cal_tauc_cpp(data, len(data), df, sat_t, end_t)\n",
    "\n",
    "    # Append the tauc value to the list\n",
    "    tauc_over_events2.append(tauc)\n",
    "A = np.vstack([mag_array, np.ones(len(mag_array))]).T\n",
    "logtauc=np.log(tauc_over_events2)\n",
    "m3, c3 = np.linalg.lstsq(A, logtauc, rcond=None)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for all events\n",
    "sat_t=20\n",
    "end_t=45\n",
    "tauc_over_events=[]\n",
    "mag_array=np.array(bh_catalog_nd['magnitude'])\n",
    "for i,trc in enumerate(event_array):\n",
    "    #event=bh_cat_sort_mag_sel.loc[i]\n",
    "    tauc_over_events.append(cal_tauc(trc,sat_t,end_t))\n",
    "A = np.vstack([mag_array, np.ones(len(mag_array))]).T\n",
    "logtauc=np.log(tauc_over_events)\n",
    "m, c = np.linalg.lstsq(A, logtauc, rcond=None)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f4616e09adc84d5fbbc3191c7baf12d0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'correlation coefficient = 0.791')"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y=m3*mag_array+c3\n",
    "msft=np.sqrt(np.sum((y-mag_array)**2))\n",
    "convar=np.sum((mag_array-mag_array.mean())*(logtauc-np.array(logtauc).mean()))/mag_array.size\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.scatter(mag_array,logtauc,label='Original tauc')\n",
    "plt.plot(mag_array,y,'b-',label='fitted tauc')\n",
    "plt.xlabel('magnitude (Mw)')\n",
    "plt.ylabel('log(tauc)')\n",
    "#plt.title(\"misfit= %f, m= %f, cov= %f\" %(msft,m,convar))\n",
    "corr_mat=np.corrcoef(mag_array,logtauc)\n",
    "plt.title('correlation coefficient = %4.3f' %(corr_mat[0][1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "55d9eb78b1424c7f8297300e1cd93689",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'correlation coefficient = 0.791')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y=m*mag_array+c\n",
    "msft=np.sqrt(np.sum((y-mag_array)**2))\n",
    "convar=np.sum((mag_array-mag_array.mean())*(logtauc-np.array(logtauc).mean()))/mag_array.size\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.scatter(mag_array,logtauc,label='Original tauc')\n",
    "plt.plot(mag_array,y,'b-',label='fitted tauc')\n",
    "plt.xlabel('magnitude (Mw)')\n",
    "plt.ylabel('log(tauc)')\n",
    "#plt.title(\"misfit= %f, m= %f, cov= %f\" %(msft,m,convar))\n",
    "corr_mat=np.corrcoef(mag_array,logtauc)\n",
    "plt.title('correlation coefficient = %4.3f' %(corr_mat[0][1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moving_average(x, w):\n",
    "    return np.convolve(x, np.ones(w), 'valid') / w\n",
    "def cal_spec_from_trc(trc,wid):\n",
    "    dt=trc.stats.delta\n",
    "    N=trc.stats.npts\n",
    "    T=N*dt\n",
    "    dw=2*np.pi/T\n",
    "    w=np.array(range(int(N/2)))*dw\n",
    "    spec_amp=np.abs(np.fft.fft(trc.data)[:int(N/2)])\n",
    "    spec_smth=moving_average(spec_amp,wid)\n",
    "    #max_ind=spec_amp.argmax()\n",
    "    max_ind=spec_smth.argmax()\n",
    "    if wid==1:\n",
    "        w_plot=w\n",
    "    else:\n",
    "        w_plot=w[:-wid+1]\n",
    "    return spec_smth, w_plot[max_ind],w_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "tauc_test_peakf=[]\n",
    "win_len=1200\n",
    "peak_f_a=[]\n",
    "for i,trc in enumerate(event_sel_array):\n",
    "    #event=bh_cat_sort_mag_sel.loc[i]\n",
    "    tauc_test_peakf.append(cal_tauc(trc,sat_t,end_t))\n",
    "    spec_amp,peak_w,w_plot=cal_spec_from_trc(trc,win_len)\n",
    "    peak_f_a.append(peak_w/(2*np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "tauc_test_peakf=[]\n",
    "peak_f_a=[]\n",
    "win_len=100\n",
    "for i,trc in enumerate(event_sel_array):\n",
    "    #event=bh_cat_sort_mag_sel.loc[i]\n",
    "    #tauc_test_peakf.append(cal_tauc(trc,sat_t,end_t))\n",
    "    spec_amp,peak_w,w_plot=cal_spec_from_trc(trc,win_len)\n",
    "    peak_f_a.append(peak_w/(2*np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1d56ec1e1ddb4d7fb7844012ff0f3310",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'correlation coefficient = 0.801')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "peak_inv=1/np.array(peak_f_a)\n",
    "x=mag_array\n",
    "A = np.vstack([x, np.ones(len(x))]).T\n",
    "#logtaucf=np.log(tauc_test_peakf)\n",
    "#logtaucf=mag_array\n",
    "m, c = np.linalg.lstsq(A, peak_inv, rcond=None)[0]\n",
    "y=m*x+c\n",
    "plt.scatter(x,peak_inv)\n",
    "plt.plot(x,y,'b-')\n",
    "plt.ylabel('1/(peak frequency) (s)')\n",
    "plt.xlabel('magnitude (Mw)')\n",
    "#plt.title('mag vs 1/f (freq peak from 25s segment)')\n",
    "corr_mat=np.corrcoef(x,peak_inv)\n",
    "plt.title('correlation coefficient = %4.3f' %(corr_mat[0][1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1 Trace(s) in Stream:\n",
       "AK.SCRK..BHZ | 2020-09-01T00:00:00.008400Z - 2020-09-09T23:59:59.988400Z | 50.0 Hz, 38880000 samples"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "st_10day_unfilter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "st_10day_unfilter=st_10day.copy()\n",
    "st_10day_unfilter_1=st_10day[0:3].merge(method=0,fill_value=0).detrend()\n",
    "st_10day_unfilter_2=st_10day[3:5].merge(method=0,fill_value=0).detrend()\n",
    "st_10day_unfilter_3=st_10day[5:8].merge(method=0,fill_value=0).detrend()\n",
    "st_10day_unfilter=st_10day_unfilter_1+st_10day_unfilter_2+st_10day_unfilter_3\n",
    "#st_10day_unfilter[2].data=np.append(st_10day_unfilter[2].data,0)\n",
    "# get unfilterd data slice\n",
    "trace=st_10day_unfilter[0]\n",
    "df = trace.stats.sampling_rate\n",
    "npts = trace.stats.npts\n",
    "t = np.arange(npts, dtype=np.float32) / df\n",
    "T=npts/df\n",
    "dw=2*np.pi/T\n",
    "event_array_unflt=[]\n",
    "event_arlt_array_unflt=[]\n",
    "pre_min=0.3\n",
    "after_min=2\n",
    "cat_bh_uft=bh_catalog_nd.loc[ind_slt].reset_index()\n",
    "for i in range(len(cat_bh_uft)):\n",
    "    event=cat_bh_uft.loc[i]\n",
    "    arrivals = model.get_travel_times(source_depth_in_km=event.depth/1000,distance_in_degree=event.disD)  \n",
    "    arl_time=event.time+arrivals[0].time\n",
    "    trace_slc=trace.slice(arl_time-pre_min*60,arl_time+after_min*60)\n",
    "    event_array_unflt.append(trace_slc)\n",
    "    event_arlt_array_unflt.append(arl_time)\n",
    "event_arlt_array_unflt_h0p1=[]\n",
    "for i in range(len(cat_bh_uft)):\n",
    "    flt_trc=event_array_unflt[i].detrend().filter('bandpass', freqmin=0.1,freqmax=25)\n",
    "    event_arlt_array_unflt_h0p1.append(flt_trc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "win_len=80\n",
    "f_divide=1\n",
    "f_highcut=5\n",
    "ratio_freq=[]\n",
    "for i in range(len(cat_bh_uft)):\n",
    "    trc=event_arlt_array_unflt_h0p1[i]\n",
    "    dt=trc.stats.delta\n",
    "    N=trc.stats.npts\n",
    "    T=N*dt\n",
    "    dw=2*np.pi/T\n",
    "    w=np.array(range(int(N/2)))*dw\n",
    "    spec_plot,peak_w,w_plot=cal_spec_from_trc(trc,win_len)\n",
    "    w_sel_a=(w_plot>f_divide*2*np.pi)*(w_plot<f_highcut*2*np.pi)\n",
    "    w_sel_an=w_plot<=f_divide*2*np.pi\n",
    "    ratio_freq.append(spec_plot[w_sel_an].sum()/spec_plot[w_sel_a].sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "82771e55c8cb4e1082b8e99260e4440b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'correlation coefficient = 0.731')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "x=mag_array\n",
    "log_ratio=np.log(ratio_freq)\n",
    "A = np.vstack([x, np.ones(len(x))]).T\n",
    "#logtaucf=np.log(tauc_test_peakf)\n",
    "#logtaucf=mag_array\n",
    "m, c = np.linalg.lstsq(A, log_ratio, rcond=None)[0]\n",
    "y=m*x+c\n",
    "plt.plot(x,y,'b-')\n",
    "plt.scatter(x,log_ratio)\n",
    "plt.xlabel('magnitude (Mw)')\n",
    "plt.ylabel('log(low/high freq content)')\n",
    "corr_mat=np.corrcoef(x,log_ratio)\n",
    "plt.title('correlation coefficient = %4.3f' %(corr_mat[0][1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testtrcn=5\n",
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().plot()\n",
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().filter('bandpass',freqmin=0.1, freqmax=24).plot()\n",
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().filter('bandpass',freqmin=2, freqmax=8).plot()\n",
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().filter('bandpass',freqmin=0.8, freqmax=2).plot()\n",
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().filter('bandpass',freqmin=0.8, freqmax=2).plot()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().filter('bandpass',freqmin=2, freqmax=4).plot()\n",
    "test_trc=event_arlt_array_unflt_h0p1[testtrcn].copy()\n",
    "test_trc.detrend().filter('bandpass',freqmin=2, freqmax=20).plot()"
   ]
  },
  {
   "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
}