channelFileReader.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5a703362",
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import math\n",
    "import numpy as np\n",
    "import datetime \n",
    "\n",
    "import obspy\n",
    "from obspy import UTCDateTime\n",
    "from obspy.taup import TauPyModel\n",
    "from obspy.geodetics.base import locations2degrees\n",
    "from obspy.clients.fdsn import Client\n",
    "\n",
    "client = Client(\"IRIS\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7b5d4315",
   "metadata": {},
   "outputs": [],
   "source": [
    "def open_lat_long_file(file, lat_long_list):\n",
    "    with open(file) as csv_file:\n",
    "        csv_reader = csv.reader(csv_file, delimiter=',')\n",
    "        line_count = 0\n",
    "        for row in csv_reader:\n",
    "            if line_count == 0:\n",
    "                print(f'{\", \".join(row)}')\n",
    "                line_count += 1\n",
    "            elif line_count == 1:\n",
    "                print(f'Column names are {\", \".join(row)}')\n",
    "                line_count += 1\n",
    "            else:\n",
    "                lat_long_list.append([(float(row[0]), float(row[1])), None])\n",
    "                line_count += 1\n",
    "        print(f'Processed {line_count} lines.')\n",
    "    return lat_long_list\n",
    "\n",
    "\n",
    "def distance_calc(origin, destination):\n",
    "    lat1, lon1 = origin\n",
    "    lat2, lon2 = destination\n",
    "    radius = 6371  # km\n",
    "\n",
    "    dlat = math.radians(lat2 - lat1)\n",
    "    dlon = math.radians(lon2 - lon1)\n",
    "    a = (math.sin(dlat / 2) * math.sin(dlat / 2) +\n",
    "         math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *\n",
    "         math.sin(dlon / 2) * math.sin(dlon / 2))\n",
    "    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))\n",
    "    d = radius * c\n",
    "\n",
    "    return d*1000 #To m\n",
    "\n",
    "def total_distance(lat_long_list):\n",
    "    total_dis = 0\n",
    "    for i, lat_long in enumerate(lat_long_list[:-1]):\n",
    "        total_dis += distance_calc(lat_long[0], lat_long_list[i + 1][0])\n",
    "    return total_dis\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3f1fd547",
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign_channel_number(lat_long_list, chan_spac):\n",
    "    optical_distance = channel_spacing * total_channel\n",
    "    dis = total_distance(lat_long_list)\n",
    "    ratio = dis/optical_distance\n",
    "    lat_long_list[0][1] = 0 #The first lat/long value being the origin\n",
    "    \n",
    "    for i, lat_long in enumerate(lat_long_list[:-1]):\n",
    "        dis = distance_calc(lat_long[0], lat_long_list[i + 1][0])\n",
    "        curr = lat_long[1]\n",
    "        channel_between = round((dis/chan_spac) * ratio)\n",
    "        lat_long_list[i + 1][1] = channel_between + curr\n",
    "    return lat_long_list\n",
    "\n",
    "def create_sim_channels(lat_long_list, chan_spac):\n",
    "    totalDistance = total_distance(lat_long_list)\n",
    "    new_channels = []\n",
    "    new_chan = []\n",
    "    \n",
    "    for i, lat_long in enumerate(lat_long_list[:-1]):\n",
    "        current_chan = lat_long\n",
    "        next_chan = lat_long_list[i + 1]\n",
    "        current_chan_num = current_chan[-1]\n",
    "        next_chan_num = next_chan[-1]\n",
    "        chan_num_list = np.arange(current_chan_num, next_chan_num)\n",
    "        num = len(chan_num_list)\n",
    "        if num == 0:\n",
    "            new_channels.append(next_chan)\n",
    "        else:\n",
    "            x_val = [current_chan_num, next_chan_num]\n",
    "            lat_val = [current_chan[0][-1], next_chan[0][-1]]\n",
    "            long_val = [current_chan[0][0], next_chan[0][0]]\n",
    "            new_lats = np.interp(chan_num_list, x_val, lat_val)\n",
    "            new_longs = np.interp(chan_num_list, x_val, long_val)\n",
    "            for i in range(len(chan_num_list)):\n",
    "                new_channels.append([(round(new_longs[i], 6), round(new_lats[i], 6)), current_chan_num + i])\n",
    "                \n",
    "    return new_channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "11aa1a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_arrival_times(lat_long_list, t, min_mag):\n",
    "    m = TauPyModel(model = 'ak135')\n",
    "    time = UTCDateTime(t.strftime(\"%Y-%m-%dT%H:%M:%S\"))\n",
    "    event_time = obspy.UTCDateTime(time)\n",
    "    cat = client.get_events(starttime = event_time - 10, endtime = event_time + 10, minmagnitude = 2)\n",
    "    origin = cat[0].preferred_origin()\n",
    "    depth = origin.depth\n",
    "    if depth < 0:\n",
    "        depth = 1 \n",
    "    for i, lat_long in enumerate(lat_long_list):\n",
    "        lat = lat_long_list[i][0][0]\n",
    "        long = lat_long_list[i][0][1]\n",
    "        distance = locations2degrees(origin.latitude, origin.longitude, lat, long)\n",
    "        arrival = m.get_ray_paths(distance_in_degree = distance, \n",
    "                         source_depth_in_km = depth / 1000, phase_list=[\"S\"])\n",
    "        arrival = arrival[0].time\n",
    "        lat_long.append(round(arrival, 6))\n",
    "        lat_long_list[i] = lat_long \n",
    "    return lat_long_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f6b7bf9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fiber Route\n",
      "Column names are Lat, Long, Elev, Description\n",
      "Processed 213 lines.\n"
     ]
    }
   ],
   "source": [
    "channel_spacing = 6.28 #channel spacing in meters\n",
    "total_channel = 4555\n",
    "\n",
    "eq_time = datetime.datetime(2022, 5, 25, 10, 38, 7) #min_mag 2\n",
    "lat_long_list = open_lat_long_file('fiberroute-seadasn.csv', [])\n",
    "\n",
    "#optical distance number of channels in file vs total distance of cable. \n",
    "\n",
    "lat_long_list = assign_channel_number(lat_long_list, channel_spacing);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e5b7cc14",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_time = datetime.datetime(2022, 5, 25, 10, 38, 7)\n",
    "lat_long_list = add_arrival_times(lat_long_list, eq_time, 2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1941c6cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "new_channels_list = create_sim_channels(lat_long_list, channel_spacing)\n",
    "new_channels_list = add_arrival_times(new_channels_list, eq_time, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f93b8566",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (Obspy)",
   "language": "python",
   "name": "python-obspy"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}