tensor_code_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# this code is written on google colab, with all required files stored in the user's google drive folder\n",
    "# please replace the file path with your own\n",
    "\n",
    "### Importing modules ###\n",
    "from google.colab import drive\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "!pip install ipympl\n",
    "import matplotlib.image as mpimg # base image\n",
    "!pip install basemap\n",
    "from mpl_toolkits.basemap import Basemap # matplotlib map overlay\n",
    "!pip install shapely\n",
    "from shapely.geometry import Point, Polygon, LineString # polygon data extraction\n",
    "from numpy.linalg import eig\n",
    "!pip install obspy\n",
    "from obspy.imaging.beachball import beachball\n",
    "\n",
    "!pip install geopy\n",
    "from geopy.distance import geodesic\n",
    "from pyproj import Geod # area for strain rates\n",
    "from shapely import wkt\n",
    "from math import radians, sin, cos, sqrt, atan2, degrees\n",
    "!pip install pyrocko\n",
    "from pyrocko import moment_tensor as pmt\n",
    "\n",
    "from google.colab import output\n",
    "output.enable_custom_widget_manager()\n",
    "\n",
    "\n",
    "\n",
    "drive.mount('/content/drive')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### read file with mw and m0 ###\n",
    "# read full CMT file from google drive\n",
    "file_path = \"/content/drive/My Drive/Colab Notebooks/cascadiafull3_mw_m0.xy\" # remember to mount drive\n",
    "# 2823 from file, but 2040 when summing up types of earthquakes (ie strikeslip, normal, subduction\n",
    "headers = [\"lon\", \"lat\", \"depth\", \"mrr\", \"mtt\", \"mpp\", \"mrt\", \"mrp\", \"mtp\", \"iexp\", \"X\", \"Y\", \"name\", \"mw_mttk\", \"m0_mttk\"]\n",
    "cascadia_full = pd.read_csv(file_path, delim_whitespace=True, names=headers, usecols=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])\n",
    "cascadia_full.drop(cascadia_full.columns[[10, 11, 12]], axis=1, inplace=True) # no name\n",
    "cascadia_full[\"iexp\"] = pd.to_numeric(cascadia_full[\"iexp\"], downcast='float')\n",
    "cascadia_full[\"m0_mttk\"] = pd.to_numeric(cascadia_full[\"m0_mttk\"], downcast='float')\n",
    "# converts 64bit data to 32bit float (stores exponent and decimal as opposed to entire number), courtesy to a friend pointing out why its not working\n",
    "\n",
    "\n",
    "### read files with CMT seperated earthquake types seperately (normal, strikeslip & subduction) ###\n",
    "# too lazy to append m0 into my files\n",
    "file_path = \"/content/drive/My Drive/Colab Notebooks/cascadianormal_mw.xy\"\n",
    "headers = [\"lon\", \"lat\", \"depth\", \"mrr\", \"mtt\", \"mpp\", \"mrt\", \"mrp\", \"mtp\", \"iexp\", \"X\", \"Y\", \"name\", \"mw_mttk\"]\n",
    "cascadia_normal = pd.read_csv(file_path, delim_whitespace=True, names=headers, usecols=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\n",
    "cascadia_normal.drop(cascadia_normal.columns[[10, 11, 12]], axis=1, inplace=True) # no name\n",
    "cascadia_normal[\"iexp\"] = pd.to_numeric(cascadia_normal[\"iexp\"], downcast='float')\n",
    "\n",
    "print('cascadia_normal len', len(cascadia_normal))\n",
    "\n",
    "file_path = \"/content/drive/My Drive/Colab Notebooks/cascadiastrikeslip_mw.xy\"\n",
    "headers = [\"lon\", \"lat\", \"depth\", \"mrr\", \"mtt\", \"mpp\", \"mrt\", \"mrp\", \"mtp\", \"iexp\", \"X\", \"Y\", \"name\", \"mw_mttk\"]\n",
    "cascadia_strikeslip = pd.read_csv(file_path, delim_whitespace=True, names=headers, usecols=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\n",
    "cascadia_strikeslip.drop(cascadia_strikeslip.columns[[10, 11, 12]], axis=1, inplace=True) # no name\n",
    "cascadia_strikeslip[\"iexp\"] = pd.to_numeric(cascadia_strikeslip[\"iexp\"], downcast='float')\n",
    "\n",
    "print('cascadia_strikeslip len', len(cascadia_strikeslip))\n",
    "\n",
    "file_path = \"/content/drive/My Drive/Colab Notebooks/cascadiasubduction_mw.xy\"\n",
    "headers = [\"lon\", \"lat\", \"depth\", \"mrr\", \"mtt\", \"mpp\", \"mrt\", \"mrp\", \"mtp\", \"iexp\", \"X\", \"Y\", \"name\", \"mw_mttk\"]\n",
    "cascadia_subduction = pd.read_csv(file_path, delim_whitespace=True, names=headers, usecols=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\n",
    "cascadia_subduction.drop(cascadia_subduction.columns[[10, 11, 12]], axis=1, inplace=True) # no name\n",
    "cascadia_subduction[\"iexp\"] = pd.to_numeric(cascadia_subduction[\"iexp\"], downcast='float')\n",
    "\n",
    "print('cascadia_subduction len', len(cascadia_subduction))\n",
    "\n",
    "\n",
    "# here for each row in each df, we check whether this earthquake exists in cascadia_full (it should!)\n",
    "# then we append a value (normal=1; strikeslip=2; subduction=3) in the new \"type\" column in cascadia_full\n",
    "\n",
    "# new column in cascadia_full called \"type\"\n",
    "cascadia_full['type'] = 0\n",
    "count = 0\n",
    "for row in range(len(cascadia_normal)):\n",
    "  match = (cascadia_normal.iloc[row, 0:10].values == cascadia_full.iloc[:, 0:10].values).all(axis=1) # outputs boolean\n",
    "  if match.any():\n",
    "    count += 1\n",
    "    cascadia_full.loc[match, 'type'] = 1\n",
    "\n",
    "print(\"count normal\", count)\n",
    "\n",
    "count = 0\n",
    "for row in range(len(cascadia_strikeslip)):\n",
    "  match = (cascadia_strikeslip.iloc[row, 0:10].values == cascadia_full.iloc[:, 0:10].values).all(axis=1) # outputs boolean\n",
    "  if match.any():\n",
    "    count += 1\n",
    "    cascadia_full.loc[match, 'type'] = 2\n",
    "\n",
    "print(\"count strikeslip\", count)\n",
    "\n",
    "count = 0\n",
    "for row in range(len(cascadia_subduction)):\n",
    "  match = (cascadia_subduction.iloc[row, 0:10].values == cascadia_full.iloc[:, 0:10].values).all(axis=1) # outputs boolean\n",
    "  if match.any():\n",
    "    count += 1\n",
    "    cascadia_full.loc[match, 'type'] = 3\n",
    "\n",
    "print(f\"count subduction {count}\\n\")\n",
    "\n",
    "\n",
    "print(cascadia_full)\n",
    "\n",
    "# Mw = moment magnitude; M0 = scaler moment = mu D A\n",
    "# \"type\": normal=1; strikeslip=2; subduction/thrust=3\n",
    "# 2823 from file, but 2040 when summing up types of eqs. this code overwrites eqs with percisely 45deg plunge with subduction"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### functions and stuff i stole from mttk ###\n",
    "# stolen and edited from mttk add_m0\n",
    "\n",
    "## function to get eigen values and vectors\n",
    "def eigen(mrr, mtt, mpp, mrt, mrp, mtp):\n",
    "  Mij = np.matrix([[ mtt, -mtp,  mrt],\n",
    "                         [ -mtp, mpp, -mrp],\n",
    "                         [ mrt, -mrp,  mrr]])\n",
    "\n",
    "  #Calculate eigenvalues and associated eigenvectors\n",
    "  eig_vals, eig_vecs = np.linalg.eig(Mij)\n",
    "\n",
    "  idx = eig_vals.argsort()[::-1]\n",
    "  eig_vals = eig_vals[idx]\n",
    "  eig_vecs = eig_vecs[:,idx]\n",
    "\n",
    "  return eig_vals, eig_vecs\n",
    "\n",
    "\n",
    "\n",
    "## function to get m0 for a moment tensor\n",
    "def get_m0(mrr, mtt, mpp, mrt, mrp, mtp, expon):\n",
    "    Mij = np.matrix([[ mtt, -mtp,  mrt],\n",
    "                         [ -mtp, mpp, -mrp],\n",
    "                         [ mrt, -mrp,  mrr]])\n",
    "\n",
    "    #Calculate eigenvalues and associated eigenvectors\n",
    "    eig_vals, eig_vecs = np.linalg.eig(Mij)\n",
    "\n",
    "    #Form diagonalised moment tensor (of elements with units dyne-cm)\n",
    "    eig_vals = eig_vals*(10**expon)\n",
    "    diag_Mij = np.diag(eig_vals)\n",
    "\n",
    "    #Calculate scalar seismic moment:\n",
    "    moment_iso   = np.trace(diag_Mij)/3.\n",
    "    moment_dev   = max([j - moment_iso for j in eig_vals])\n",
    "    moment_total = moment_iso + moment_dev\n",
    "\n",
    "    #Calculate moment magnitude from dyne-cm (e.g. Kanamori & Hanks 1979):\n",
    "    Mw = np.log10(moment_total)*(2./3) - 10.7\n",
    "\n",
    "    return moment_total\n",
    "\n",
    "'''\n",
    "## homemade m0 code THAT ALIGNS WITH HAVARD CMT for box m0\n",
    "\n",
    "def get_m0(mrr, mtt, mpp, mrt, mrp, mtp):\n",
    "  moment_tensor = np.array([[mtt, -mtp, mrt], [-mtp, mpp, -mrp], [mrt, -mrp, mrr]])\n",
    "  evals, evecs = eig(moment_tensor)\n",
    "  # https://www.geeksforgeeks.org/numpy-linalg-eig-method-in-python/\n",
    "  m0 = np.sqrt(0.5*(evals[0]**2 + evals[1]**2 + evals[2]**2)) # Bowers & Hudson (1999) Eqn 5\n",
    "  return m0\n",
    "'''"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### Loading GMT basemap and overlaying matplotlib map for interactive point picking ###\n",
    "## seting up figure and enabling interactive plot\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "# fig = plt.figure(figsize=(15.25, 12.3))\n",
    "fig = plt.figure(figsize=(10.5, 8.2)) #a4 portrait figure size\n",
    "ax1 = plt.subplot(111)\n",
    "\n",
    "\n",
    "## import GMT-generated basemap png\n",
    "gmtmap_image = mpimg.imread(\"/content/drive/My Drive/Colab Notebooks/focal_mechanisms_CMT_mercator_cropped.png\")\n",
    "\n",
    "\n",
    "## creating matplotlib basemap for overlay\n",
    "# https://matplotlib.org/basemap/index.html\n",
    "m = Basemap(projection='merc', resolution='l',\\\n",
    "                          llcrnrlat=38.5,urcrnrlat=66.3,\\\n",
    "                          llcrnrlon=-180, urcrnrlon=-118.5)\n",
    "\n",
    "# interactive_map.drawcoastlines() # drawing coastlines on matplotlib disabled\n",
    "\n",
    "\n",
    "## plotting parallels and meridians, annotating axis\n",
    "parallels = np.arange(30, 80, 10)\n",
    "meridians = np.arange(-190, -110, 10)\n",
    "m.drawparallels(parallels, labels=[1,0,0,0], fontsize=12, linewidth=0.5) # label parallels on right and top\n",
    "m.drawmeridians(meridians, labels=[0,0,0,1], fontsize=12, linewidth=0.5) # meridians on bottom and left\n",
    "\n",
    "\n",
    "## plotting image in bottom layer - https://stackoverflow.com/questions/11487797/python-matplotlib-basemap-overlay-small-image-on-map-plot\n",
    "m.imshow(gmtmap_image, origin='upper')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "\n",
    "### Interactive point picking ###\n",
    "# mainly follows this - https://sgcorcoran.github.io/Python-for-MSE/lessons_3114/06_selecting_pixel_coordinates.html\n",
    "# useful reference - https://matplotlib.org/stable/gallery/event_handling/pick_event_demo.html\n",
    "\n",
    "\n",
    "## conversion codes - between x y and lon lat\n",
    "# references this - https://matplotlib.org/basemap/users/mapcoords.html\n",
    "# lon,lat can be scalars, lists or numpy arrays\n",
    "\n",
    "def pos_to_lonlat(x, y):\n",
    "  lon, lat = m(x, y, inverse=True)\n",
    "  return lon, lat\n",
    "\n",
    "def lonlat_to_pos(lon, lat):\n",
    "  x, y = m(lon, lat)\n",
    "  return x, y\n",
    "\n",
    "\n",
    "## motified from https://sgcorcoran.github.io/Python-for-MSE/lessons_3114/06_selecting_pixel_coordinates.html\n",
    "pos = [] # in format [[None, None], [x, y], [x, y], ...]\n",
    "lonlat = []\n",
    "\n",
    "\n",
    "def onclick(event):\n",
    "    pos.append([event.xdata, event.ydata])\n",
    "\n",
    "    lon, lat = pos_to_lonlat(pos[-1][0], pos[-1][1]) # pos[-1] represent last click (list with x, y)\n",
    "    lonlat.append([lon, lat]) # converts x y to lon lat and appends\n",
    "\n",
    "    ax1.set_title(f'Click {len(pos)}: {lon}, {lat}')\n",
    "\n",
    "\n",
    "cid=fig.canvas.mpl_connect('button_press_event', onclick)\n",
    "     "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### manual data input from excel spreedsheet ###\n",
    "# REMEMBER TO PRESS RESET FIRST\n",
    "lon_excel = '-155.83456221 -156.36359447 -156.43917051 -156.43917051 -156.66589862 -156.89262673 -157.27050691 -157.0437788  -155.91013825 -155.68341014 -155.6078341  -155.30552995 -154.32304147 -153.340553   -152.50921659 -152.20691244 -151.82903226 -151.45115207 -150.84654378 -150.77096774 -149.78847926 -148.27695853 -147.14331797 -146.23640553 -145.93410138 -147.97465438 -148.73041475 -148.88156682 -149.41059908 -149.93963134 -150.31751152 -150.84654378 -151.90460829 -152.05576037 -152.35806452 -152.88709677 -153.340553   -153.71843318 -154.32304147 -154.24746544 -154.47419355 -154.47419355 -155.30552995 -155.38110599'\n",
    "lat_excel = '55.91381866 55.74402509 55.44509968 55.23018503 55.10067606 54.92734241 54.40283436 54.31475695 54.18228652 54.40283436 54.70962033 54.8403944 55.18706211 55.35927371 55.65885013 55.78654298 55.87143974 55.95615131 56.37693896 56.46054402 56.75171707 57.40886108 58.09438961 58.84527923 59.0790819  59.73295455 59.0790819  58.92338966 58.68852849 58.33324697 58.37290045 58.13431081 57.73308242 57.44954678 57.28653242 57.32735393 57.24566558 57.32735393 57.0406504  56.87581932 56.79313016 56.46054402 56.25118657 56.0828719'\n",
    "lon_excel = [float(lon) for lon in lon_excel.split()]\n",
    "lat_excel = [float(lat) for lat in lat_excel.split()]\n",
    "\n",
    "# check whether length is the same\n",
    "length_check = len(lon_excel) if len(lon_excel) == len(lat_excel) else \"Lengths are not the same\"\n",
    "print(\"Number of polygon bounding points:\", length_check)\n",
    "\n",
    "lonlat = list(zip(lon_excel, lat_excel))\n",
    "polygon_df = pd.DataFrame(lonlat)\n",
    "\n",
    "print(polygon_df)\n",
    "\n",
    "\n",
    "## conversion\n",
    "\n",
    "m = Basemap(projection='merc', resolution='l',\\\n",
    "                          llcrnrlat=38.5,urcrnrlat=66.3,\\\n",
    "                          llcrnrlon=-180, urcrnrlon=-118.5)\n",
    "\n",
    "def pos_to_lonlat(x, y):\n",
    "  lon, lat = m(x, y, inverse=True)\n",
    "  return lon, lat\n",
    "\n",
    "def lonlat_to_pos(lon, lat):\n",
    "  x, y = m(lon, lat)\n",
    "  return x, y\n",
    "\n",
    "\n",
    "pos = [lonlat_to_pos(lon, lat) for lon, lat in zip(lon_excel, lat_excel)]\n",
    "pts = np.array(pos[:])\n",
    "pts = np.append(pts,[pts[0]], axis=0) # to close our polygon when plotting. axis=0 specify along rows\n",
    "\n",
    "ax1.plot(pts[:,0],pts[:,1], 'b-')\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "file_path = \"/content/drive/My Drive/Colab Notebooks/Project_Cascadia_Kostrov.xlsx\"\n",
    "headers = [\"name\", \"info\", \"to_do\", \"polygon_pts\", \"lon\", \"lat\", \"area\", \"volume\", \"eq_no\", \"max_depth\", \"mmr\", \"mtt\", \"mpp\", \"mrt\", \"mrp\", \"mtp\", \"summed_m0\", \"individual_m0\", \"Cs\",\\\n",
    "           \"eigenP\", \"eigenI\", \"eigenT\", \"strike\", \"dip\", \"lon_l\", \"lat_l\", \"total_l\", \"w_length\", \"v_ms\", \"v_mmyr\", \"strainP\", 'strainI', 'strainT']\n",
    "excel_full = pd.read_excel(file_path, names=headers, usecols=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,\\\n",
    "                                                                                        21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33])\n",
    "\n",
    "print(excel_full)\n",
    "'''\n",
    "     "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### polygon drawing and data saving ###\n",
    "pts = np.array(pos[:])\n",
    "pts = np.append(pts,[pts[0]], axis=0) # to close our polygon when plotting. axis=0 specify along rows\n",
    "\n",
    "ax1.plot(pts[:,0],pts[:,1], 'b-')\n",
    "\n",
    "polygon_df = pd.DataFrame(lonlat)\n",
    "\n",
    "print(\"Number of polygon bounding points:\", len(lonlat))\n",
    "print(polygon_df)\n",
    "\n",
    "print(\"\\n for excel\")\n",
    "print(\"longitude:\")\n",
    "print(polygon_df[0].to_numpy())\n",
    "print(\"latitude:\")\n",
    "print(polygon_df[1].to_numpy())\n",
    "\n",
    "\n",
    "print(\"\\n for gmt; lon lat:\")\n",
    "for row in range(len(polygon_df)):\n",
    "  print(str(polygon_df.iloc[row, 0]), str(polygon_df.iloc[row, 1]))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### reset selection ###\n",
    "pos = []\n",
    "lonlat = []\n",
    "ax1.set_title(f'Clicks reset')\n",
    "\n",
    "for line in ax1.lines:\n",
    "        line.remove()\n",
    "plt.show()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### input moment tensor from CMT based on box chosen ###\n",
    "\n",
    "## check if a point is inside the polygon\n",
    "polygon = Polygon(lonlat)\n",
    "def is_inside_polygon(lon, lat):\n",
    "    point = Point(lon, lat)\n",
    "    return polygon.contains(point)\n",
    "\n",
    "## empty df for earthquakes in polygon\n",
    "headers = [\"lon\", \"lat\", \"depth\", \"mrr\", \"mtt\", \"mpp\", \"mrt\", \"mrp\", \"mtp\", \"iexp\", \"X\", \"Y\", \"name\", \"mw_mttk\", \"m0_mttk\", \"type\"]\n",
    "earthquakes_in_polygon_df = pd.DataFrame(columns=headers)\n",
    "earthquakes_in_polygon_df.drop(earthquakes_in_polygon_df.columns[[10, 11, 12]], axis=1, inplace=True) # no name\n",
    "\n",
    "\n",
    "## extract data within the polygon\n",
    "for index, row in cascadia_full.iterrows():\n",
    "    # iterrows() returns both the index and pd series\n",
    "    lon, lat = row['lon'], row['lat']\n",
    "    if is_inside_polygon(lon, lat):\n",
    "        earthquakes_in_polygon_df = earthquakes_in_polygon_df.append(row, ignore_index=True)\n",
    "\n",
    "earthquakes_in_polygon_df_backup = earthquakes_in_polygon_df.copy() # backup\n",
    "\n",
    "\n",
    "### defining depth range and earthquake type if applicable ###\n",
    "\n",
    "min_depth = 0.\n",
    "max_depth = 50\n",
    "# earthquake_type = 3   # normal=1; strikeslip=2; subduction/thrust=3\n",
    "\n",
    "\n",
    "earthquakes_in_polygon_df = earthquakes_in_polygon_df_backup.copy() # reset earthquake_in_polygon_df\n",
    "\n",
    "earthquakes_in_polygon_df_filtered = earthquakes_in_polygon_df[\\\n",
    "(earthquakes_in_polygon_df['depth'] >= min_depth) & (earthquakes_in_polygon_df['depth'] <= max_depth)] # depth filter\n",
    "\n",
    "# earthquakes_in_polygon_df_filtered = earthquakes_in_polygon_df_filtered[(earthquakes_in_polygon_df['type'] == 3)] # type filter\n",
    "# earthquakes_in_polygon_df_filtered = earthquakes_in_polygon_df_filtered[(earthquakes_in_polygon_df['type'] == 2) | (earthquakes_in_polygon_df['type'] == 1)] # 2 types filter\n",
    "\n",
    "earthquakes_in_polygon_df = earthquakes_in_polygon_df_filtered # confusing but so that I dont have to rename all of my variables\n",
    "\n",
    "print(\"\\n No. of eqs enclosed within polygon:\", len(earthquakes_in_polygon_df_backup))\n",
    "\n",
    "print(f'\\nMax earthquake depth WITHOUT depth filter: {earthquakes_in_polygon_df_backup[\"depth\"].max()} km')\n",
    "polygon_max_depth = earthquakes_in_polygon_df[\"depth\"].max()\n",
    "print(f'Max earthquake depth w depth filter: {polygon_max_depth} km')\n",
    "\n",
    "print(\"\\n No. of eqs enclosed within polygon AFTER depth and type filter:\", len(earthquakes_in_polygon_df))\n",
    "print(\"\\n Filtered df\", earthquakes_in_polygon_df)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### Kostrov tensor summation ###\n",
    "# empty 1x7 array [mrr, mtt, mpp, mrt, mrp, mtp] to store summation\n",
    "tensor_sum = np.zeros_like(earthquakes_in_polygon_df.iloc[0, 3:9]) # np.zeros_like allows the use of ref array\n",
    "\n",
    "## polygon summation\n",
    "for row in range(len(earthquakes_in_polygon_df)):\n",
    "  tensor_row = earthquakes_in_polygon_df.iloc[row, 3:9] * (10**earthquakes_in_polygon_df.iloc[row, 9]) #multiply by iexp ie normalized\n",
    "  tensor_sum += tensor_row\n",
    "\n",
    "# reduce exponent of tensor_sum by 27, add iexp at the end\n",
    "tensor_sum /= 10**27\n",
    "\n",
    "\n",
    "\n",
    "## m0 from INDIVIDUAL moment tensor summation\n",
    "'''\n",
    "sum_of_individual_m0 = 0.\n",
    "for row in range(len(earthquakes_in_polygon_df)):\n",
    "  sum_of_individual_m0 += earthquakes_in_polygon_df.iloc[row, 11]\n",
    "'''\n",
    "\n",
    "\n",
    "## prints\n",
    "print(\"\\n Summed moment tensor Mij (* 10^27):\")\n",
    "print(tensor_sum)\n",
    "\n",
    "\n",
    "m0_from_sum = get_m0(*tensor_sum, 27)\n",
    "print(\"\\n m0 from SUMMED moment tensor Mij:\")\n",
    "print(m0_from_sum)\n",
    "# func is 6 different inputs not an array --> * unpacks the array into seperate enteries\n",
    "# 1.5923\n",
    "\n",
    "sum_of_individual_m0 = np.sum(earthquakes_in_polygon_df['m0_mttk'])\n",
    "print(\"\\n m0 from INDIVIDUAL moment tensor summation:\")\n",
    "print(sum_of_individual_m0)\n",
    "\n",
    "\n",
    "Cs = m0_from_sum / sum_of_individual_m0\n",
    "print(\"\\n seismic consistency:\")\n",
    "print(Cs)\n",
    "\n",
    "\n",
    "eigen_vals = eigen(*tensor_sum)[1]\n",
    "print(\"\\nEigen values in descending order:\")\n",
    "print(eigen_vals)\n",
    "\n",
    "eigen_vec = eigen(*tensor_sum)[0]\n",
    "print(\"Eigen vectors:\")\n",
    "print(eigen_vec)\n",
    "\n",
    "\n",
    "## plotting beach ball\n",
    "# https://docs.obspy.org/tutorial/code_snippets/beachball_plot.html\n",
    "tensor_sum_list = tensor_sum.tolist()\n",
    "beachball(tensor_sum_list, size=200, linewidth=2, facecolor='c')\n",
    "fig.canvas.mpl_disconnect(cid) # disconnect interactive picking"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### plotting CUMULATIVE log(frequency)-magnitude plot ###\n",
    "# tested with small sample - ok\n",
    "\n",
    "fig, ax2 = plt.subplots(figsize =(3.5, 2.5))\n",
    "\n",
    "\n",
    "magnitude_bin = np.arange(4, 10, 0.5) # min, max (stops at 9.3), interval\n",
    "# barely any datapoints --> adjust histogram bin width\n",
    "bin_centers = 0.5*(magnitude_bin[:-1] + magnitude_bin[1:]) # midpoints, starts at 0 ends right before the end etc\n",
    "\n",
    "frequency_bin = np.zeros_like(bin_centers)\n",
    "\n",
    "\n",
    "for row in range(len(earthquakes_in_polygon_df)):\n",
    "    mw = earthquakes_in_polygon_df.iloc[row]['mw_mttk']\n",
    "    for bin in range(len(magnitude_bin) - 1):\n",
    "      if mw >= magnitude_bin[bin]:\n",
    "        frequency_bin[bin] += 1\n",
    "\n",
    "\n",
    "## plotting\n",
    "ax2.plot(bin_centers, frequency_bin, 'c.-')\n",
    "ax2.set_title(f'{len(earthquakes_in_polygon_df)} samples, bin width {magnitude_bin[1] - magnitude_bin[0]}')\n",
    "ax2.set_yscale('log') # log10N = a − bM\n",
    "\n",
    "# linear regression to find line of best fit; log10N = a - bM\n",
    "coefficients = np.polyfit(bin_centers[2:5], np.log10(frequency_bin[2:5] + 1e-7), 1) # + 1 to avoid dividing by 10\n",
    "a, b = coefficients[0], coefficients[1]\n",
    "print('best fit coefficients',a,b)\n",
    "\n",
    "x = np.linspace(4, 9.2, 100)\n",
    "best_fit_line = 10 ** (b + a * x)\n",
    "ax2.plot(x, best_fit_line, 'k-', label=f'log10N = {a} - {b}M')\n",
    "ax2.plot([3,9],[1,1])\n",
    "\n",
    "ax2.set_ylim(1, None) # adjusting limit\n",
    "ax2.set_xlim(4, 9.2)\n",
    "plt.show()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### interactive fault L picking - copy paste from earlier without comments ###\n",
    "## DOES NOT clip automatically within box - use eye\n",
    "%matplotlib inline\n",
    "%matplotlib widget\n",
    "\n",
    "# fig = plt.figure(figsize=(15.75, 12.3))\n",
    "fig = plt.figure(figsize=(10.5, 8.2))\n",
    "ax4 = plt.subplot(111)\n",
    "\n",
    "gmtmap_image = mpimg.imread(\"/content/drive/My Drive/Colab Notebooks/focal_mechanisms_CMT_mercator_cropped.png\")\n",
    "\n",
    "m = Basemap(projection='merc', resolution='l',\\\n",
    "                          llcrnrlat=38.5,urcrnrlat=66.3,\\\n",
    "                          llcrnrlon=-180, urcrnrlon=-118.5)\n",
    "\n",
    "\n",
    "parallels = np.arange(30, 80, 10)\n",
    "meridians = np.arange(-190, -110, 10)\n",
    "m.drawparallels(parallels, labels=[1,0,0,0], fontsize=12, linewidth=0.5) # label parallels on right and top\n",
    "m.drawmeridians(meridians, labels=[0,0,0,1], fontsize=12, linewidth=0.5) # meridians on bottom and left\n",
    "\n",
    "ax4.plot(pts[:,0],pts[:,1], 'b-') # showing original box\n",
    "\n",
    "m.imshow(gmtmap_image, origin='upper')\n",
    "ax4.set_title(f'Pick fault length')\n",
    "plt.show()\n",
    "\n",
    "## Interactive point picking\n",
    "## motified from https://sgcorcoran.github.io/Python-for-MSE/lessons_3114/06_selecting_pixel_coordinates.html\n",
    "pos_fault = [] # in format [[None, None], [x, y], [x, y], ...]\n",
    "lonlat_fault = []\n",
    "\n",
    "def onclick(event):\n",
    "    pos_fault.append([event.xdata, event.ydata])\n",
    "\n",
    "    lon, lat = pos_to_lonlat(pos_fault[-1][0], pos_fault[-1][1]) # pos[-1] represent last click (list with x, y)\n",
    "    lonlat_fault.append([lon, lat]) # converts x y to lon lat and appends\n",
    "\n",
    "    ax4.set_title(f'Click {len(pos_fault)}: {lon}, {lat}')\n",
    "\n",
    "\n",
    "cid=fig.canvas.mpl_connect('button_press_event', onclick)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "## fault L saving and printing\n",
    "fault_pts = np.array(pos_fault[:])\n",
    "\n",
    "ax4.plot(fault_pts[:,0],fault_pts[:,1], 'r-') # print fault_L\n",
    "\n",
    "fault_L_df = pd.DataFrame(lonlat_fault)\n",
    "print(fault_L_df)\n",
    "\n",
    "print(\"\\n for excel\")\n",
    "print(\"longitude fault L:\")\n",
    "print(fault_L_df[0].to_numpy())\n",
    "print(\"latitude fault L:\")\n",
    "print(fault_L_df[1].to_numpy())"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### reset selection fault_L ###\n",
    "pos_fault = []\n",
    "lonlat_fault = []\n",
    "ax4.set_title(f'Clicks reset')\n",
    "\n",
    "for line in ax4.lines:\n",
    "        line.remove()\n",
    "\n",
    "ax4.plot(pts[:,0],pts[:,1], 'b-') # showing original box\n",
    "\n",
    "plt.show()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### calculating velocity rate ###\n",
    "# expect values orders of mm to cm a year\n",
    "\n",
    "# calculating dustance between two points with haversine - https://en.wikipedia.org/wiki/Haversine_formula\n",
    "def haversine_distance(lat1, lon1, lat2, lon2):\n",
    "    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2]) # convert from deg to rad\n",
    "\n",
    "    # radius of the Earth in meters\n",
    "    R = 6371000.0\n",
    "\n",
    "    # Haversine formula\n",
    "    dlat = lat2 - lat1\n",
    "    dlon = lon2 - lon1\n",
    "\n",
    "\n",
    "    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2\n",
    "    distance = 2 * R * np.arcsin(sqrt(a))\n",
    "\n",
    "    return distance\n",
    "\n",
    "\n",
    "## fault length\n",
    "# distance between all consecutive points for fault L\n",
    "\n",
    "total_distance_fault_L = 0\n",
    "\n",
    "for i in range(len(lonlat_fault) - 1):\n",
    "    lat1, lon1 = lonlat_fault[i]\n",
    "    lat2, lon2 = lonlat_fault[i + 1]\n",
    "    distance = haversine_distance(lat1, lon1, lat2, lon2)\n",
    "    total_distance_fault_L += distance\n",
    "    # print(f\"Distance between point {i + 1} and point {i + 2}: {distance:.2f} m\")\n",
    "\n",
    "dip = 27.03 # dip in deg FROM MTTK\n",
    "polygon_max_depth = earthquakes_in_polygon_df_backup[\"depth\"].max()\n",
    "# polygon_max_depth = 20\n",
    "\n",
    "print(f'\\nMax earthquake depth WITHOUT depth filter: {earthquakes_in_polygon_df_backup[\"depth\"].max()} km')\n",
    "# polygon_max_depth\n",
    "print(f'Max depth used: {polygon_max_depth} km')\n",
    "print(\"Dip used:\", dip)\n",
    "print(f\"\\nTotal fault length: {total_distance_fault_L:.2f} m\")\n",
    "\n",
    "## fault W\n",
    "length_W = (polygon_max_depth * 1000) / np.sin(dip * (np.pi / 180)) # km to m; sin in radians\n",
    "print(f\"\\nlength_W: {length_W:.2f} m\")\n",
    "\n",
    "\n",
    "\n",
    "# other variables\n",
    "shear_modulus = 3.3 * (10 ** 10) # Stacey (1992)\n",
    "catalogue_duration_s = 17460 * 24 * 60 * 60\n",
    "# from 1st Jan 1976 to 20th Oct 2023; 17460 days\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "v_factor = shear_modulus * total_distance_fault_L * length_W * catalogue_duration_s\n",
    "v_rate = (1/ v_factor) * sum_of_individual_m0 * (1e-7) # either one m0 works; CONVERT FROM dyn cm to Nm\n",
    "\n",
    "print(f'\\nvelocity rate: {v_rate} m s^-1')\n",
    "\n",
    "## in mm per year (to compare with catalogues)\n",
    "v_rate_mmyr = v_rate * 31536000 * 1000 # 1 yr = 31536000 sec; convert from m to mm\n",
    "\n",
    "\n",
    "print(f'velocity rate: {v_rate_mmyr} mm year^-1')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### calculating average strain rate de / dt ###\n",
    "# expect values approx 10 ** -8 or even 10 ** -11\n",
    "print(f'\\nMax earthquake depth WITHOUT depth filter: {earthquakes_in_polygon_df_backup[\"depth\"].max()} km')\n",
    "\n",
    "polygon_max_depth = earthquakes_in_polygon_df[\"depth\"].max()\n",
    "polygon_max_depth = 30\n",
    "\n",
    "print(f'Max depth used: {polygon_max_depth} km')\n",
    "\n",
    "\n",
    "## variables\n",
    "shear_modulus = 3.3 * (10 ** 10) # Stacey (1992)\n",
    "catalogue_duration_s = 17460 * 24 * 60 * 60\n",
    "# from 1st Jan 1976 to 20th Oct 2023; 17460 days\n",
    "\n",
    "# calculating polygon area (geodesic)\n",
    "# https://stackoverflow.com/questions/23697374/calculate-polygon-area-in-planar-units-e-g-square-meters-in-shapely\n",
    "# first order approx checked from manually approximated google earth pro\n",
    "geod = Geod(ellps='WGS84') # around meter accuracy\n",
    "\n",
    "polygon_wkt = wkt.loads(f\"{Polygon(lonlat)}\")\n",
    "polygon_area = abs(geod.geometry_area_perimeter(polygon_wkt)[0])\n",
    "print('\\nGeodesic area: {:.3f} m^2'.format(polygon_area))\n",
    "\n",
    "# calculating volume\n",
    "\n",
    "polygon_volume = polygon_max_depth * 1000 * polygon_area # V = L * b * d; km to m\n",
    "print('Polygon volume: {:.3f} m^3'.format(polygon_volume))\n",
    "\n",
    "\n",
    "## calculating strain rate\n",
    "strain_factor = 2 * shear_modulus * polygon_volume * catalogue_duration_s\n",
    "strain_rate_eigen_vec = [0, 0, 0]\n",
    "\n",
    "print('\\nStrain eigen value')\n",
    "for i in range(len(eigen_vec)):\n",
    "  print(eigen_vec[i])\n",
    "  strain_rate_eigen_vec[i] = (1 / strain_factor) *  eigen_vec[i] * (10 ** 27) * (1e-7)\n",
    "\n",
    "# strain_rate_eigen_vec = (1 / strain_factor) *  eigen_vec * (10 ** 27) * (1e-7) # Mij doesnt have exponent information; CONVERT FROM dyn cm to Nm\n",
    "\n",
    "print(f'\\nstrain rate: {strain_rate_eigen_vec} s^-1')\n",
    "\n",
    "strain_rate_eigen_vec_yr = [0, 0, 0]\n",
    "for i in range(len(strain_rate_eigen_vec)):\n",
    "  strain_rate_eigen_vec_yr[i] = strain_rate_eigen_vec[i] * 31536000\n",
    "\n",
    "print(f'strain rate: {strain_rate_eigen_vec_yr} s^-1 year^-1')\n",
    "# takes\n",
    "     "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "### strain P T eigen value surface projection ###\n",
    "# normal box 2a\n",
    "strain_eigen_P_yr = strain_rate_eigen_vec_yr[0]\n",
    "strain_eigen_T_yr = strain_rate_eigen_vec_yr[2]\n",
    "strain_eigen_P_dip = 7.27\n",
    "strain_eigen_T_dip = 19.12\n",
    "\n",
    "\n",
    "eig_P_projection = strain_eigen_P_yr * np.cos(abs(strain_eigen_P_dip) * (np.pi / 180))\n",
    "eig_T_projection = strain_eigen_T_yr * np.cos(abs(strain_eigen_T_dip) * (np.pi / 180))\n",
    "\n",
    "print(\"Strain rate eigen vector P:\", strain_eigen_P_yr)\n",
    "print(\"Strain rate eigen vector T:\", strain_eigen_T_yr)\n",
    "print(\"Eigen P dip:\", strain_eigen_P_dip)\n",
    "print(\"Eigen T dip:\", strain_eigen_T_dip)\n",
    "\n",
    "\n",
    "print(\"\\nEigen P surface projection:\", eig_P_projection)\n",
    "print(\"Eigen T surface projection:\", eig_T_projection)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# generating a grid of lon lat points for unavco geodetic velocity vectors\n",
    "\n",
    "lonlat_geodetic_df = pd.DataFrame(columns=[\"lon\", \"lat\", \"depth\"])\n",
    "lon1 = 118\n",
    "lon2 = 180\n",
    "lat1 = 38\n",
    "lat2 = 66\n",
    "\n",
    "lon_range = lon2 - lon1\n",
    "print(\"lon range\", lon_range)\n",
    "lat_range = lat2 - lat1\n",
    "print(\"lat range\", lat_range)\n",
    "\n",
    "for lon in range(lon1, lon2+1, 2): # 2deg interval\n",
    "  for lat in range(lat1, lat2+1, 2):\n",
    "    lonlat_geodetic_df = lonlat_geodetic_df.append({'lon': -lon, 'lat': lat, 'depth': 0.}, ignore_index=True) # append returns a new df\n",
    "\n",
    "\n",
    "print(\"lon lat:\")\n",
    "for row in range(len(lonlat_geodetic_df)):\n",
    "  print(str(lonlat_geodetic_df.iloc[row, 0]), str(lonlat_geodetic_df.iloc[row, 1]), str(lonlat_geodetic_df.iloc[row, 2]), end=',')\n"
   ],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "orig_nbformat": 4,
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}