Interseis_practical.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/logos.png\" style=\"float: centre;\" width=\"800\"/>\n",
    "\n",
    "# Interseismic practical\n",
    "\n",
    "## Written for the COMET InSAR Workshop 2022.\n",
    "\n",
    "Welcome to this COMET practical on interseismic strain accumulation.\n",
    "\n",
    "This practical has been tested in the latest versions of Chrome and Firefox, but should work in most modern browsers.\n",
    "\n",
    "### Interacting with this practical\n",
    "\n",
    "To run a cell containing code, you can either press \"Run\" in the bar above, or press \"Shift+Enter\".\n",
    "\n",
    "A \"#\" indicates a comment within the code. \n",
    "\n",
    "##################################### <br />\n",
    "Blocks of code surrounded by multiple \"#\" (like this one) indicate variables that you can experiment changing the values to. While you can experiment with the full body of code further (and I encourage you to do so after the course), we are limited in our ability to help you debug any errors you may encounter. Changing just the values indicate should result in no major errors. <br />\n",
    "#####################################\n",
    "\n",
    "If you get a \"dead kernel\" error, either restart the kernel under \"Kernel\" on the top bar, or reopen the binder fresh. None of your changes here will affect the original binder, so feel free to reload it if anything breaks.\n",
    "\n",
    "### Acknowledgements\n",
    "Detailed guidance on the content was provided by [John Elliott](https://environment.leeds.ac.uk/see/staff/1248/dr-john-elliott). Advice on the computational requirements of the practical was provided by [Richard Rigby](https://environment.leeds.ac.uk/see/staff/2698/richard-rigby). This notebook was written using resources from the Tectonophysics ([Tim Wright](https://environment.leeds.ac.uk/see/staff/1613/professor-tim-wright)) and Inverse Theory ([Phil Livermore](https://environment.leeds.ac.uk/see/staff/1381/dr-phil-livermore)) undergraduate modules at the University of Leeds. Thanks to everyone for their contributions.\n",
    "\n",
    "[Andrew Watson](https://environment.leeds.ac.uk/see/pgr/2311/andrew-watson) - 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required modules. Run this before continuing.\n",
    "\n",
    "# these are for performing the analysis\n",
    "import numpy as np\n",
    "from scipy import interpolate\n",
    "from scipy import stats\n",
    "import statistics\n",
    "import pyproj\n",
    "import subprocess as subp\n",
    "from pathlib import Path\n",
    "\n",
    "# this is our own library of functions, which can be found in the same directory as this notebook\n",
    "import interseis_lib as lib\n",
    "\n",
    "# this package is for general plotting\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "# this is additional scientific colour maps, see \"https://www.fabiocrameri.ch/colourmaps/\"\n",
    "from cmcrameri import cm\n",
    "\n",
    "# these packages are only needed for the final multivariate plot\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "\n",
    "# this refreshes the import of interseis_lib if changes are made while the notebook is running\n",
    "import importlib\n",
    "importlib.reload(lib)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Introduction\n",
    "\n",
    "The earthquake cycle is a model for how earthquakes occur on faults. It consists of three phases:\n",
    "- **Interseismic** - the fault is locked towards the surface but is accumulating strain at depth driven by external forces (i.e. plate motion).\n",
    "- **Coseismic** - the strain is released seismically by slip on the previously-locked section of the fault.\n",
    "- **Postseismic** - mechanisms after the earthquake such as afterslip or viscoelastic relaxation.\n",
    "\n",
    "We'll be focusing on the interseismic period for this practical. Figure 1 shows a schematic diagram of how the ground deforms through the interseismic period (middle panel) and into the coseismic (right panel). Profile A-A' shows a characteristic interseismic signal in the middle panel (the shape of an arctan function), where the ground displacement decreases from the far-field to zero at the fault trace.\n",
    "\n",
    "<img src=\"figures/earthquakeCycle.png\" style=\"float: centre;\" width=\"600\"/>\n",
    "\n",
    "*Figure 1: Schematic map diagram of how expected ground deformation through the interseismic and coseismic periods. From Elliott (2009).*\n",
    "\n",
    "Geodesists are interested in this interseismic period as it provides information on the rate of strain accumulation and the size of the locked part of the fault, important parameters for understanding the seismic potential of a fault and the distribution of strain in the crust.\n",
    "\n",
    "We can model interseismic deformation by considering steady slip on a screw dislocation embedded within an elastic half-space (Savage and Burford, 1973, Wright et al., 2013). In this simple model, the fault-parallel surface velocities, $v$, are a function of the perpendicular distance from the fault, $x$, the slip rate, $s$, and the locking depth, $d$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\normalsize v(x) = \\frac{s}{\\pi}tan^{-1}\\left(\\frac{x}{d}\\right) + c \\qquad \\qquad (1)\n",
    "\\end{equation}\n",
    "\n",
    "Figure 2 shows a schematic diagram of this screw dislocation model.\n",
    "\n",
    "<img src=\"figures/screwDisc.png\" style=\"float: centre;\" width=\"600\"/>\n",
    "\n",
    "*Figure 2: Schematic diagram of a screw dislocation model. Reproduced from Vernant (2015).*\n",
    "\n",
    "Large volumes of InSAR data have been processed over the North Anatolian Fault (NAF) in Turkey by LiCS. We'll be using two of these LiCSBAS time series (087A_04904_121313 and 167D_04884_131212, shown in Figure 3) to estimate fault-parallel velocities, which we will then use to estimate the values of the parameters in Equation 1.\n",
    "\n",
    "<img src=\"figures/frameMap.png\" style=\"float: centre;\" width=\"600\"/>\n",
    "\n",
    "*Figure 3: Location of our two LiCS frames overlayed on faults (red lines) from the Global Earthquake Model (https://github.com/GEMScienceTools/gem-global-active-faults).*\n",
    "\n",
    "This practical is broken down into the following sections:\n",
    "\n",
    "**1)** Generating a simple synthetic screw dislocation model and investigating what we might expect InSAR to measure.\n",
    "\n",
    "**2)** Loading in our LiCSBAS velocities and other data files.\n",
    "\n",
    "**3)** Decomposing our two line-of-sight velocity maps into both East and vertical velocities, under the assumption of zero North velocity, and into fault-parallel and vertical velocities, under the assumption of zero fault-perpendicular velocity. This includes resampling our velocities onto a unified grid and setting a common reference point.\n",
    "\n",
    "**4)** Projecting the fault-parallel velocities onto a fault-perpendicular profile.\n",
    "\n",
    "**5)** Performing a forward model to explore our profiled velocities.\n",
    "\n",
    "**6)** Performing a Bayesian inversion to estimate the parameters in Equation 1, along with their associated uncertainties.\n",
    "\n",
    "**7)** Extension - incorporation of surface fault creep in our model.\n",
    "\n",
    "**8)** Extension - estimating shear strain rate.\n",
    "\n",
    "For an example of this general method, have a look at our pre-print on the Main Recent Fault, Iran:\n",
    "https://www.essoar.org/doi/abs/10.1002/essoar.10507719.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Synthetic screw dislocation model\n",
    "\n",
    "We'll start by generating some synthetic velocities using our screw dislocation model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set interseismic fault slip rate, s, and locking depth, d\n",
    "#####################################\n",
    "s = 20 # mm/yr          \n",
    "d = 20 # km             \n",
    "#####################################\n",
    "\n",
    "# set up x coords (fault-perpendicular distances)\n",
    "xmin, xmax, xint = -200, 200, 2 # km\n",
    "x = np.arange(xmin, xmax+xint, xint)\n",
    "\n",
    "# run the screw dislocation model. The 0 is just a static offset which we can ignore for now.\n",
    "# We convert all units into SI (metres) so that they are consistent.\n",
    "v = lib.screw_disc(x*1000, s/1000, d*1000, 0)\n",
    "\n",
    "# we'll also generate a model at s = 10 mm/yr and d = 10 km as a reference solution\n",
    "v_ref = lib.screw_disc(x*1000, 10/1000, 10*1000, 0)\n",
    "\n",
    "# expand into 2D, assuming our fault strikes East-West\n",
    "v_grid = np.tile(v, (len(x),1)).T\n",
    "\n",
    "\n",
    "# Plot the result\n",
    "# Let's break this first plotting section down\n",
    "\n",
    "# create a figure and axes for our plot\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,5))\n",
    "\n",
    "# plot the model velocities as a line to the first axes/subplot\n",
    "axs[0].plot(x, v*1000, color=\"blue\", label='Your model: \\ns = ' + str(s) + ' mm/yr, d = ' + str(d) + ' km') \n",
    "\n",
    "# plot the reference model\n",
    "axs[0].plot(x, v_ref*1000, color=\"red\", label='Reference model: \\ns = 10 mm/yr, d = 10 km')\n",
    "\n",
    "# draw dashed lines at x=0 (the fault trace) and at y=0\n",
    "axs[0].plot([0, 0],[axs[0].get_ylim()[0], axs[0].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[0].plot([axs[0].get_xlim()[0], axs[0].get_xlim()[1]], [0, 0], color='grey', linestyle='dashed')\n",
    "\n",
    "# add a legend\n",
    "axs[0].legend(fontsize=12)\n",
    "\n",
    "# set labels for the x and y axes\n",
    "axs[0].set_xlabel('Fault-perpendicular distance (km)')\n",
    "axs[0].set_ylabel('Fault-parallel velocity (mm/yr)')\n",
    "\n",
    "\n",
    "# colour plot of our 2D velocities, now on our second axes/subplot\n",
    "im = axs[1].imshow(v_grid*1000, extent=[xmin, xmax, xmin, xmax], cmap=cm.vik)\n",
    "\n",
    "# add a labeled colour bar\n",
    "plt.colorbar(im, ax=axs[1], label=\"Fault-parallel velocity (mm/yr)\")\n",
    "\n",
    "# add the fault as a solid line\n",
    "axs[1].plot([xmin, xmax], [0, 0], color='black', label='Fault trace') \n",
    "\n",
    "# add the profile as a dashed line\n",
    "axs[1].plot([0, 0],[xmin, xmax], color='black', linestyle='dashed', label=\"Profile line\")\n",
    "\n",
    "# add legend\n",
    "axs[1].legend(fontsize=12)\n",
    "\n",
    "# set labels for the x and y axes\n",
    "axs[1].set_xlabel('x-coord (km)')\n",
    "axs[1].set_ylabel('y-coord (km)')\n",
    "\n",
    "\n",
    "# this alters subplot spacing\n",
    "fig.tight_layout()\n",
    "\n",
    "# display the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left plot shows our simple 1D screw dislocation model. The right shows this model expanded to 2D, assuming that our fault strikes from East to West. The fault-parallel velocities are also East-West velocities in this setup.\n",
    "\n",
    "While these fault-parallel velocities are what we want for modelling a fault, they're not what we directly measure with InSAR. SAR satellites such as Sentinel-1 look down at the ground with an incidence angle, $\\theta$, of around 31-46&deg;. This angle changes through our scene. InSAR is also only capable of measuring displacements towards and away from the satellite. The result of these two factors is that InSAR measures only a component of the observed deformation projected into the satellite's line-of-sight (LOS), and a single look direction cannot differentiate between different components of displacement.\n",
    "\n",
    "<img src=\"figures/incAngle.png\" style=\"float: centre;\" width=\"300\"/>\n",
    "\n",
    "*Figure 4: Schematic diagram of satellite incidence angle, $\\theta$, measured between the vertical (dashed black line) and the satellite line-of-sight (blue arrow).*\n",
    "\n",
    "Let's apply a line-of-sight projection to our signals and see what effect this has on the velocities.\n",
    "\n",
    "The interseismic_lib includes a function to make a grid of incidence angles between a minimum and maximum value, rotated by the satellite heading. Remember that the Sentinel-1 satellites are side looking, so our incidence angles will change at 90&deg; to the satellite heading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# heading (-14 for ascending, -166 for descending), minimum incidence angle, and maximum incidence angle\n",
    "#####################################\n",
    "heading = -14 # degrees\n",
    "inc_min, inc_max = 31, 46 # degrees\n",
    "#####################################\n",
    "\n",
    "# create grids of these angles\n",
    "heading_grid = np.ones(v_grid.shape) * heading\n",
    "inc_grid = lib.gen_inc(inc_min, inc_max, heading, x, x)\n",
    "\n",
    "# projection from east-west to los\n",
    "e2los = np.cos(np.deg2rad(heading_grid)) * np.sin(np.deg2rad(inc_grid))\n",
    "v_grid_los = v_grid * e2los\n",
    "\n",
    "# plot \n",
    "fig, axs = plt.subplots(2,2,figsize=(15,12))\n",
    "\n",
    "# get colour limits for grids\n",
    "cmin, cmax = np.amin(v_grid*1000), np.amax(v_grid*1000)\n",
    "\n",
    "# incidence angle grid\n",
    "im = axs[0,0].imshow(inc_grid, extent=[xmin, xmax, xmin, xmax], cmap=cm.batlow)\n",
    "plt.colorbar(im, ax=axs[0,0], label=\"Incidence angle (degrees)\")\n",
    "\n",
    "# East-west velocity grid\n",
    "im = axs[0,1].imshow(v_grid*1000, extent=[xmin, xmax, xmin, xmax], cmap=cm.vik, vmin=cmin, vmax=cmax)\n",
    "plt.colorbar(im, ax=axs[0,1], label=\"East-west velocity (mm/yr)\")\n",
    "axs[0,1].plot([xmin, xmax], [0, 0], color='black', label='Fault trace') \n",
    "axs[0,1].plot([0, 0],[xmin, xmax], color='black', linestyle='dashed', label=\"Profile line\")\n",
    "axs[0,1].legend(fontsize=12)\n",
    "\n",
    "# 1D models\n",
    "axs[1,0].plot(x, v*1000, color='red', label='original (East-West) velocity')\n",
    "axs[1,0].plot(x, v_grid_los[:,np.argmin(np.abs(x))]*1000, color='blue', label='projected LOS velocity')\n",
    "axs[1,0].plot([0, 0],[axs[1,0].get_ylim()[0], axs[1,0].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[1,0].plot([axs[1,0].get_xlim()[0], axs[1,0].get_xlim()[1]], [0, 0], color='grey', linestyle='dashed')\n",
    "axs[1,0].set_xlabel('Fault-perpendicular distance (km)')\n",
    "axs[1,0].set_ylabel('Velocity (mm/yr)')\n",
    "axs[1,0].legend()\n",
    "\n",
    "# Line-of-sight velocity grid\n",
    "im = axs[1,1].imshow(v_grid_los*1000, extent=[xmin, xmax, xmin, xmax], cmap=cm.vik, vmin=cmin, vmax=cmax)\n",
    "plt.colorbar(im, ax=axs[1,1], label=\"Line-of-sight velocity (mm/yr)\")\n",
    "axs[1,1].plot([xmin, xmax], [0, 0], color='black')\n",
    "axs[1,1].plot([0, 0],[xmin, xmax], color='black', linestyle='dashed')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that not only does the overall signal get smaller, but it also changes across the scene because of the varying incidence angle of the satellite.\n",
    "\n",
    "Now that we've explored how a screw dislocation model might look in InSAR, we'll move on to applying the model to real world data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load the input time-series\n",
    "\n",
    "We'll begin by loading in our LiCSBAS outputs (you generated the equivalent outputs in the LiCSBAS practicals). These have been renamed to include the frame name for clarity. The _asc and _desc tags in the variable names indicate the ascending or descending frame.\n",
    "\n",
    "**vel_file** - file containing the LOS velocities (found in TS_GEOCml10/results/vel).\n",
    "\n",
    "**par_file** - parameter file which contains the geographic location of the data, along with the size in x and y (found in GEOCml##).\n",
    "\n",
    "**E/N/U_file** - unit vector components, which basically tell us how to go from the satellite LOS to East, North, and Up. We need these to perform the velocity decomposition later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup file names\n",
    "asc = '087A_04904_121313'\n",
    "desc = '167D_04884_131212'\n",
    "\n",
    "vel_file_asc = Path('data/' + asc + '.vel.mskd')\n",
    "par_file_asc = Path('data/' + asc + '.par')\n",
    "E_file_asc = Path('data/' + asc + '.E')\n",
    "N_file_asc = Path('data/' + asc + '.N')\n",
    "U_file_asc = Path('data/' + asc + '.U')\n",
    "\n",
    "vel_file_desc = Path('data/' + desc + '.vel.mskd')\n",
    "par_file_desc = Path('data/' + desc + '.par')\n",
    "E_file_desc = Path('data/' + desc + '.E')\n",
    "N_file_desc = Path('data/' + desc + '.N')\n",
    "U_file_desc = Path('data/' + desc + '.U')\n",
    "\n",
    "# read array dimensions from par file\n",
    "width_asc = int(lib.get_par(par_file_asc,'width'))\n",
    "length_asc = int(lib.get_par(par_file_asc,'nlines'))\n",
    "\n",
    "width_desc = int(lib.get_par(par_file_desc,'width'))\n",
    "length_desc = int(lib.get_par(par_file_desc,'nlines'))\n",
    "\n",
    "# get corner positions\n",
    "corner_lat_asc = float(lib.get_par(par_file_asc,'corner_lat'))\n",
    "corner_lon_asc = float(lib.get_par(par_file_asc,'corner_lon'))\n",
    "\n",
    "corner_lat_desc = float(lib.get_par(par_file_desc,'corner_lat'))\n",
    "corner_lon_desc = float(lib.get_par(par_file_desc,'corner_lon'))\n",
    "\n",
    "# get post spacing (distance between velocity measurements)\n",
    "post_lat_asc = float(lib.get_par(par_file_asc,'post_lat'))\n",
    "post_lon_asc = float(lib.get_par(par_file_asc,'post_lon'))\n",
    "\n",
    "post_lat_desc = float(lib.get_par(par_file_desc,'post_lat'))\n",
    "post_lon_desc = float(lib.get_par(par_file_desc,'post_lon'))\n",
    "\n",
    "# calculate grid spacings\n",
    "lat_asc = corner_lat_asc + post_lat_asc*np.arange(1,length_asc+1) - post_lat_asc/2\n",
    "lon_asc = corner_lon_asc + post_lon_asc*np.arange(1,width_asc+1) - post_lon_asc/2\n",
    "\n",
    "lat_desc = corner_lat_desc + post_lat_desc*np.arange(1,length_desc+1) - post_lat_desc/2\n",
    "lon_desc = corner_lon_desc + post_lon_desc*np.arange(1,width_desc+1) - post_lon_desc/2\n",
    "\n",
    "# load in velocities\n",
    "vel_asc = np.fromfile(vel_file_asc, dtype='float32').reshape((length_asc, width_asc))\n",
    "vel_desc = np.fromfile(vel_file_desc, dtype='float32').reshape((length_desc, width_desc))\n",
    "\n",
    "# load in unit vectors\n",
    "E_asc = np.fromfile(E_file_asc, dtype='float32').reshape((length_asc, width_asc))\n",
    "N_asc = np.fromfile(N_file_asc, dtype='float32').reshape((length_asc, width_asc))\n",
    "U_asc = np.fromfile(U_file_asc, dtype='float32').reshape((length_asc, width_asc))\n",
    "\n",
    "E_desc = np.fromfile(E_file_desc, dtype='float32').reshape((length_desc, width_desc))\n",
    "N_desc = np.fromfile(N_file_desc, dtype='float32').reshape((length_desc, width_desc))\n",
    "U_desc = np.fromfile(U_file_desc, dtype='float32').reshape((length_desc, width_desc))\n",
    "\n",
    "# load the frame polygons for plotting\n",
    "poly_asc = np.loadtxt(Path('data/' + asc + '.poly'))\n",
    "poly_desc = np.loadtxt(Path('data/' + desc + '.poly'))\n",
    "\n",
    "# load the naf fault trace\n",
    "fault_trace = np.loadtxt(Path('data/naf_trace.xy'))\n",
    "\n",
    "print('Files loaded.')\n",
    "print('Size of ascending velocity:')\n",
    "print(vel_asc.shape)\n",
    "print('Size of descending velocity:')\n",
    "print(vel_desc.shape)\n",
    "\n",
    "print('Limits of ascending (lon, lat):')\n",
    "print(np.amin(lon_asc), np.amax(lon_asc), np.amin(lat_asc), np.amax(lat_asc))\n",
    "print('Limits of descending (lon, lat):')\n",
    "print(np.amin(lon_desc), np.amax(lon_desc), np.amin(lat_desc), np.amax(lat_desc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll plot the velocities to make sure they've loaded in correctly. The trace of the NAF will be plotted as a red line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot input vels\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,5))\n",
    "\n",
    "# plot the ascending velocities, saturating the colour palette from -20 to 20 mm/yr.\n",
    "im = axs[0].imshow(vel_asc, extent=[np.amin(lon_asc), np.amax(lon_asc), np.amin(lat_asc), np.amax(lat_asc)], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "\n",
    "# add the NAF fault trace\n",
    "axs[0].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "\n",
    "# add the frame polygon\n",
    "axs[0].plot(poly_asc[:,0], poly_asc[:,1], color=\"black\", linestyle='dashed')\n",
    "\n",
    "# display the colorbar, this \"divider\" section is to scale the colorbar size to the plot\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='LOS velocity (mm/yr)')\n",
    "\n",
    "# set the title\n",
    "axs[0].set_title('087A_04904_121313')\n",
    "\n",
    "# set the x and y limits of the plot\n",
    "axs[0].set_xlim(np.amin(lon_asc), np.amax(lon_asc))\n",
    "axs[0].set_ylim(np.amin(lat_asc), np.amax(lat_asc))\n",
    "\n",
    "# label the x and y axes\n",
    "axs[0].set_xlabel('Longitude (degrees)')\n",
    "axs[0].set_ylabel('Latitude (degrees)')\n",
    "\n",
    "# repeat for the descending velocities\n",
    "im = axs[1].imshow(vel_desc, extent=[np.amin(lon_desc), np.amax(lon_desc), np.amin(lat_desc), np.amax(lat_desc)], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[1].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "axs[1].plot(poly_desc[:,0], poly_desc[:,1], color=\"black\", linestyle='dashed')\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='LOS velocity (mm/yr)')\n",
    "axs[1].set_title('167D_04884_131212')\n",
    "axs[1].set_xlim(np.amin(lon_desc), np.amax(lon_desc))\n",
    "axs[1].set_ylim(np.amin(lat_desc), np.amax(lat_desc))\n",
    "axs[1].set_xlabel('Longitude (degrees)')\n",
    "axs[1].set_ylabel('Latitude (degrees)')\n",
    "\n",
    "# show the plot\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, LiCSBAS sets motion towards the satellite as positive, and motion away from the satellite as negative. These velocities are relative to a LiCSBAS-selected reference point within each frame.\n",
    "\n",
    "We can already see a change a the LOS velocites across the fault in both maps. To investigate further, lets take a profile across the fault in the ascending velocities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start and end coordinates of the profile, and extra parameters\n",
    "#####################################\n",
    "prof_start = (33.5, 40.2)\n",
    "prof_end = (33.1, 41.8)\n",
    "prof_params = {\n",
    "    \"nbins\": 100, # number of bins to split the profile into\n",
    "    \"width\": 0.05 # total width of the profile in degrees (1 degree ~ 110 km)\n",
    "}\n",
    "#####################################\n",
    "\n",
    "# run the profiler, the outputs are as follows:\n",
    "# - bin_val = mean value of each bin\n",
    "# - prof_bin_mids = distance along the profile to the middle of each bin\n",
    "# - points_val = every velocity within the profile\n",
    "# - points_dist = distance along the profile to every point within the profile\n",
    "# - points_poly = polygon that defines the profile\n",
    "bin_val, prof_bin_mids, points_val, points_dist, points_poly \\\n",
    "    = lib.profile_data(lon_asc,lat_asc,vel_asc,prof_start,prof_end,prof_params)\n",
    "\n",
    "# calculate profile-fault intersection angle\n",
    "intersect_dist, intersect_angle = lib.profile_fault_intersection(prof_start,prof_end,fault_trace)\n",
    "\n",
    "# plot the ascending velocities and the profiled data\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,5))\n",
    "\n",
    "# ascending velocity map\n",
    "im = axs[0].imshow(vel_asc, extent=[np.amin(lon_asc), np.amax(lon_asc), np.amin(lat_asc), np.amax(lat_asc)], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[0].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "axs[0].plot([prof_start[0], prof_end[0]], [prof_start[1], prof_end[1]], color=\"red\")\n",
    "axs[0].plot(points_poly[:,0],points_poly[:,1], color=\"red\")\n",
    "axs[0].scatter(prof_start[0],prof_start[1], s=100, color='pink', edgecolor='black', zorder=3)\n",
    "axs[0].scatter(prof_end[0],prof_end[1], s=100, color='black', edgecolor='white', zorder=3)\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='LOS velocity (mm/yr)')\n",
    "axs[0].set_title('087A_04904_121313')\n",
    "axs[0].set_xlim(np.amin(lon_asc), np.amax(lon_asc))\n",
    "axs[0].set_ylim(np.amin(lat_asc), np.amax(lat_asc))\n",
    "\n",
    "# velocities projected onto profile\n",
    "axs[1].scatter(points_dist, points_val, s=2)\n",
    "axs[1].plot([intersect_dist, intersect_dist],[axs[1].get_ylim()[0], axs[1].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[1].plot([axs[1].get_xlim()[0], axs[1].get_xlim()[1]], [0, 0], color='grey', linestyle='dashed')\n",
    "axs[1].set_xlabel(\"Distance along profile (degrees)\")\n",
    "axs[1].set_ylabel(\"Line-of-sight velocity (mm/yr)\")\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left plot shows our velocity map. The pink dot marks the start of the profile, and the black dot marks the end. The three red line between these show the main profile line and the extents of the profile width.\n",
    "\n",
    "Our profiled velocities are shown on the right. These are projected onto the profile line, and so are given at distance along the profile from prof_start, which is located at 0 distance. The units of distance are the same as those for the profiled data set. The vertical dashed line marks where the profile line intersects the fault trace of the NAF.\n",
    "\n",
    "While the profiled velocities are a little noisy, you should be able to see an overall change in velocity across the fault."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Decompose the line-of-sight velocities into fault-parallel and vertical\n",
    "\n",
    "Equation 1 requires fault-parallel velocities, while ours are currently in the satellite LOS.\n",
    "We can estimate fault-parallel velocities by decomposing our LOS velocities into fault-parallel, fault-perpendicular, and vertical velocities.\n",
    "This is done by solving a linear system of equations between our InSAR velocities ($V_{asc}$ and $V_{desc}$) and our component velocities.\n",
    "We convert between the two using the radar incidence angle, $\\theta$, and the azimuth of the along-track satellite heading, $\\alpha$, and the strike of our fault, $\\gamma$. Both $\\theta$ and $\\alpha$ vary pixel-by-pixel.\n",
    "\n",
    "First, let's look at how we would decompose our LOS velocities into East ($V_E$), North ($V_N$), and vertical ($V_U$) velocities:\n",
    "\\begin{equation}\n",
    "    \\begin{bmatrix} V_{asc} \\\\ V_{desc} \\end{bmatrix} = \\begin{bmatrix} sin(\\theta_{asc})cos(\\alpha_{asc}) & sin(\\theta_{asc})sin(\\alpha_{asc}) & -cos(\\theta_{asc}) \\\\ sin(\\theta_{desc})cos(\\alpha_{desc}) & sin(\\theta_{desc})sin(\\alpha_{desc}) & -cos(\\theta_{desc}) \\end{bmatrix} \\begin{bmatrix} V_E \\\\ V_N \\\\ V_U \\end{bmatrix}\n",
    "    \\qquad \\qquad (2)\n",
    "\\end{equation}\n",
    "\n",
    "We calcualte fault-parallel ($V_{para}$) and fault-perpendicular ($V_{perp}$) velocities as a sum of relavent components of $V_E$ and $V_N$, including an extra term to account for the strike of the fault measured clockwise from North: \n",
    "\n",
    "\\begin{equation}\n",
    "    \\begin{bmatrix} V_{asc} \\\\ V_{desc} \\end{bmatrix} = \\begin{bmatrix} \\begin{pmatrix} sin(\\gamma)sin(\\theta_{asc})cos(\\alpha_{asc})\\\\+ cos(\\gamma)sin(\\theta_{asc})sin(\\alpha_{asc}) \\end{pmatrix} & \\begin{pmatrix} cos(\\gamma)sin(\\theta_{asc})cos(\\alpha_{asc})\\\\+ sin(\\gamma)sin(\\theta_{asc})sin(\\alpha_{asc}) \\end{pmatrix} & -cos(\\theta_{asc}) \\\\ \\begin{pmatrix} sin(\\gamma)sin(\\theta_{desc})cos(\\alpha_{desc})\\\\+ cos(\\gamma)sin(\\theta_{desc})sin(\\alpha_{desc}) \\end{pmatrix} & \\begin{pmatrix} cos(\\gamma)sin(\\theta_{desc})cos(\\alpha_{desc})\\\\+ sin(\\gamma)sin(\\theta_{desc})sin(\\alpha_{desc}) \\end{pmatrix} & -cos(\\theta_{desc}) \\end{bmatrix} \\begin{bmatrix} V_{para} \\\\ V_{perp} \\\\ V_U \\end{bmatrix}\n",
    "    \\qquad \\qquad (3)\n",
    "\\end{equation}\n",
    "\n",
    "In both cases, we have three unknowns and only two data sets, making this an underdetermined system.\n",
    "Given that the amount of data available to us is fixed, we need to reduce the number of component velocities we want to estimate by one so that we can solve for the remaining two.\n",
    "We'll assume that $V_{perp}$ is zero so that we can solve for $V_{para}$ and $V_U$. We choose $V_{perp}$ because we expect the dominant signal to be strike-slip motion on the East-West striking NAF, producing a small North-South component which the Sentinel-1 satellites are particularly insensitive to because of their near-polar orbit.\n",
    "\n",
    "We can now rewrite Equation 3:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\begin{bmatrix} V_{asc} \\\\ V_{desc} \\end{bmatrix} = \\begin{bmatrix} sin(\\gamma)sin(\\theta_{asc})cos(\\alpha_{asc})+cos(\\gamma)sin(\\theta_{asc})sin(\\alpha_{asc}) & -cos(\\theta_{asc}) \\\\ sin(\\gamma)sin(\\theta_{desc})cos(\\alpha_{desc})+cos(\\gamma)sin(\\theta_{desc})sin(\\alpha_{desc}) & -cos(\\theta_{desc}) \\end{bmatrix} \\begin{bmatrix} V_{para} \\\\ V_U \\end{bmatrix}\n",
    "    \\qquad \\qquad (4)\n",
    "\\end{equation}\n",
    "\n",
    "where $V_{asc}$ and $V_{desc}$ are the LOS velocities, $\\theta$ is the radar incidence angle, and $\\alpha$ is the azimuth of the along-track satellite heading. The conversion from LOS to fault-parallel can be simplified to $sin(\\theta)sin(\\alpha + \\gamma)$, however, we only have the unit component vectors available, not the individual values of $\\theta$ and $\\alpha$, so we'll stick with Equation 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently, these velocities are on two slightly different grids. We need to unify the grids so that we can perform the velocity decomposition for shared points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# limits and intervals for new grid\n",
    "lon_min = np.amin([lon_asc[0], lon_desc[0]])\n",
    "lon_max = np.amax([lon_asc[-1], lon_desc[-1]])\n",
    "lon_int = np.amin([post_lon_asc, post_lon_desc])\n",
    "lon_regrid = np.arange(lon_min,lon_max+lon_int,lon_int)\n",
    "\n",
    "lat_max = np.amax([lat_asc[0], lat_desc[0]])\n",
    "lat_min = np.amin([lat_asc[-1], lat_desc[-1]])\n",
    "lat_int = np.absolute(np.amin([post_lat_asc, post_lat_desc]))\n",
    "lat_regrid = np.arange(lat_min,lat_max+lat_int,lat_int)\n",
    "\n",
    "xx_regrid, yy_regrid = np.meshgrid(lon_regrid, lat_regrid[::-1])\n",
    "coords_regrid = np.transpose(np.vstack((xx_regrid.flatten(),yy_regrid.flatten())))\n",
    "\n",
    "# interpolate velocities onto this new grid\n",
    "xx_asc, yy_asc = np.meshgrid(lon_asc, lat_asc)\n",
    "vel_asc_regrid = interpolate.griddata((xx_asc.ravel(), yy_asc.ravel()), vel_asc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel()))\n",
    "vel_asc_regrid = vel_asc_regrid.reshape((len(lat_regrid),len(lon_regrid)))\n",
    "\n",
    "xx_desc, yy_desc = np.meshgrid(lon_desc, lat_desc)\n",
    "vel_desc_regrid = interpolate.griddata((xx_desc.ravel(), yy_desc.ravel()), vel_desc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel()))\n",
    "vel_desc_regrid = vel_desc_regrid.reshape((len(lat_regrid),len(lon_regrid)))\n",
    "\n",
    "# interpolate the component vectors into the new grid\n",
    "E_asc_regrid = interpolate.griddata((xx_asc.ravel(), yy_asc.ravel()), E_asc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel())).reshape((len(lat_regrid),len(lon_regrid)))\n",
    "N_asc_regrid = interpolate.griddata((xx_asc.ravel(), yy_asc.ravel()), N_asc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel())).reshape((len(lat_regrid),len(lon_regrid)))\n",
    "U_asc_regrid = interpolate.griddata((xx_asc.ravel(), yy_asc.ravel()), U_asc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel())).reshape((len(lat_regrid),len(lon_regrid)))\n",
    "\n",
    "E_desc_regrid = interpolate.griddata((xx_desc.ravel(), yy_desc.ravel()), E_desc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel())).reshape((len(lat_regrid),len(lon_regrid)))\n",
    "N_desc_regrid = interpolate.griddata((xx_desc.ravel(), yy_desc.ravel()), N_desc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel())).reshape((len(lat_regrid),len(lon_regrid)))\n",
    "U_desc_regrid = interpolate.griddata((xx_desc.ravel(), yy_desc.ravel()), U_desc.ravel(), (xx_regrid.ravel(), yy_regrid.ravel())).reshape((len(lat_regrid),len(lon_regrid)))\n",
    "\n",
    "print('Regridding complete.')\n",
    "print('Size of ascending velocity:')\n",
    "print(vel_asc_regrid.shape)\n",
    "print('Size of descending velocity:')\n",
    "print(vel_desc_regrid.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently the two velocities fields are relative to their own reference pixels. We need to make this reference the same so that the velocities can be decomposed correctly. We'll define a reference window and use the mean of all the velocities within that area as the reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define new reference pixel\n",
    "ref_xmin, ref_xmax = 32.65, 32.85\n",
    "ref_ymin, ref_ymax = 41.1, 41.2\n",
    "\n",
    "# ref poly for plotting\n",
    "ref_poly = np.array([[ref_xmin, ref_ymin],\n",
    "           [ref_xmin, ref_ymax],\n",
    "           [ref_xmax, ref_ymax],\n",
    "           [ref_xmax, ref_ymin],\n",
    "           [ref_xmin, ref_ymin]])\n",
    "\n",
    "# get index\n",
    "ind_xmin = np.argmin(np.absolute(lon_regrid-ref_xmin))\n",
    "ind_xmax = np.argmin(np.absolute(lon_regrid-ref_xmax))\n",
    "ind_ymin = np.argmin(np.absolute(lat_regrid-ref_ymin))\n",
    "ind_ymax = np.argmin(np.absolute(lat_regrid-ref_ymax))\n",
    "\n",
    "# get ref value\n",
    "ref_val_asc = np.nanmean(vel_asc_regrid[ind_ymin:ind_ymax+1,ind_xmin:ind_xmax+1])\n",
    "ref_val_desc = np.nanmean(vel_desc_regrid[ind_ymin:ind_ymax+1,ind_xmin:ind_xmax+1])\n",
    "\n",
    "# set as new reference\n",
    "vel_asc_regrid = vel_asc_regrid - ref_val_asc\n",
    "vel_desc_regrid = vel_desc_regrid - ref_val_desc\n",
    "\n",
    "print('Referencing complete.')\n",
    "print('Mean value of reference for ascending = ' + str(round(ref_val_asc,3)) + ' mm/yr')\n",
    "print('Mean value of reference for descending = ' + str(round(ref_val_desc,3)) + ' mm/yr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the regridded velocities to make sure they're fine\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,10))\n",
    "\n",
    "# ascending velocities\n",
    "im = axs[0].imshow(vel_asc_regrid, extent=[lon_regrid[0], lon_regrid[-1], lat_regrid[0], lat_regrid[-1]], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[0].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "axs[0].plot(ref_poly[:,0], ref_poly[:,1], color='black')\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='LOS velocity (mm/yr)')\n",
    "axs[0].set_title('087A_04904_121313')\n",
    "axs[0].set_xlim(np.amin(lon_regrid), np.amax(lon_regrid))\n",
    "axs[0].set_ylim(np.amin(lat_regrid), np.amax(lat_regrid))\n",
    "\n",
    "# descending velocities\n",
    "im = axs[1].imshow(vel_desc_regrid, extent=[lon_regrid[0], lon_regrid[-1], lat_regrid[0], lat_regrid[-1]], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[1].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "axs[1].plot(ref_poly[:,0], ref_poly[:,1], color='black')\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='LOS velocity (mm/yr)')\n",
    "axs[1].set_title('167D_04884_131212')\n",
    "axs[1].set_xlim(np.amin(lon_regrid), np.amax(lon_regrid))\n",
    "axs[1].set_ylim(np.amin(lat_regrid), np.amax(lat_regrid))\n",
    "\n",
    "fig.tight_layout(w_pad=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll run the decomposition for each pixel, assuming no correlation between pixels and a strike of 80&deg;."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set strike\n",
    "strike = 80\n",
    "strike_E, strike_N = np.sin(np.deg2rad(strike)), np.cos(np.deg2rad(strike))\n",
    "\n",
    "# pre-allocate\n",
    "vel_para = np.zeros((len(lat_regrid), len(lon_regrid)))\n",
    "vel_U = np.zeros((len(lat_regrid), len(lon_regrid)))\n",
    "\n",
    "# loop through every pixel\n",
    "for ii in np.arange(0,len(lat_regrid)):\n",
    "    for jj in np.arange(0,len(lon_regrid)):\n",
    "        \n",
    "        # create the design matrix\n",
    "        G = np.array([[strike_E*E_asc_regrid[ii,jj] + strike_N*N_asc_regrid[ii,jj], U_asc_regrid[ii,jj]], \n",
    "                      [strike_E*E_desc_regrid[ii,jj] + strike_N*N_desc_regrid[ii,jj], U_desc_regrid[ii,jj]]])\n",
    "        \n",
    "        # get the two velocities for this pixel\n",
    "        d = np.array([[vel_asc_regrid[ii,jj], vel_desc_regrid[ii,jj]]]).T\n",
    "        \n",
    "        # solve the linear system for the Up and East velocities\n",
    "        m = np.linalg.solve(G, d)\n",
    "        \n",
    "        # save to arrays\n",
    "        vel_para[ii,jj] = m[0]\n",
    "        vel_U[ii,jj] = m[1]\n",
    "        \n",
    "print('Velocity decomposition complete.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the decomposed velocities, including the frame polygons so that we can see how we only get the decomposed velocities in the overlap area, where we have both line-of-sights to work with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the results\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,10))\n",
    "\n",
    "# East velocities\n",
    "im = axs[0].imshow(vel_para, extent=[lon_regrid[0], lon_regrid[-1], lat_regrid[0], lat_regrid[-1]], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[0].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "axs[0].plot(poly_asc[:,0], poly_asc[:,1], color=\"black\", linestyle='dashed')\n",
    "axs[0].plot(poly_desc[:,0], poly_desc[:,1], color=\"black\", linestyle='dashed')\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='Fault-parallel velocity (mm/yr)')\n",
    "axs[0].set_title('Fault-parallel')\n",
    "axs[0].set_xlim(np.amin(lon_regrid), np.amax(lon_regrid))\n",
    "axs[0].set_ylim(np.amin(lat_regrid), np.amax(lat_regrid))\n",
    "\n",
    "# Up velocities\n",
    "im = axs[1].imshow(vel_U, extent=[lon_regrid[0], lon_regrid[-1], lat_regrid[0], lat_regrid[-1]], \\\n",
    "                   cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[1].plot(fault_trace[:,0], fault_trace[:,1], color=\"red\")\n",
    "axs[1].plot(poly_asc[:,0], poly_asc[:,1], color=\"black\", linestyle='dashed')\n",
    "axs[1].plot(poly_desc[:,0], poly_desc[:,1], color=\"black\", linestyle='dashed')\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='Vertical velocity (mm/yr)')\n",
    "axs[1].set_title('Up')\n",
    "axs[1].set_xlim(np.amin(lon_regrid), np.amax(lon_regrid))\n",
    "axs[1].set_ylim(np.amin(lat_regrid), np.amax(lat_regrid))\n",
    "\n",
    "plt.tight_layout(w_pad=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Profile the fault-parallel velocities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now profile our decomposed fault-parallel velocities to produce a dataset which we can then apply our screw dislocation model to. First, we need to project our lat-long coordinates into Universal Transverse Mercator (UTM) so that the profile distances are in metres and not degrees. After this projection, our velocities will no longer be on a regular grid, so we'll plot them using scatter instead of imshow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to UTM using gdal\n",
    "\n",
    "# get the utm zone projection\n",
    "utm_crs_list = pyproj.database.query_utm_crs_info(\n",
    "    datum_name=\"WGS 84\",\n",
    "    area_of_interest = pyproj.aoi.AreaOfInterest(\n",
    "        west_lon_degree  = np.amin(lon_regrid),\n",
    "        south_lat_degree = np.amin(lat_regrid),\n",
    "        east_lon_degree  = np.amax(lon_regrid),\n",
    "        north_lat_degree = np.amax(lat_regrid),\n",
    "    ),\n",
    ")\n",
    "utm_crs = pyproj.CRS.from_epsg(utm_crs_list[0].code)\n",
    "\n",
    "# create transformer for LL to UTM\n",
    "transformer = pyproj.Transformer.from_crs('epsg:4326', utm_crs)\n",
    "\n",
    "# apply transform to our grids of lat long coordinates\n",
    "xx_utm, yy_utm = transformer.transform(yy_regrid, xx_regrid)\n",
    "\n",
    "# apply transform to fault trace\n",
    "fault_trace_utm = fault_trace.copy()\n",
    "fault_trace_utm[:,0], fault_trace_utm[:,1] = transformer.transform(fault_trace[:,1], fault_trace[:,0])\n",
    "\n",
    "# convert utm metres to km (easier to work with for plotting)\n",
    "xx_utm = xx_utm / 1000\n",
    "yy_utm = yy_utm / 1000\n",
    "fault_trace_utm = fault_trace_utm / 1000\n",
    "\n",
    "# get new limits for plotting\n",
    "val_ind_row, val_ind_col = np.where(~np.isnan(vel_U).all(axis=0))[0], np.where(~np.isnan(vel_U).all(axis=1))[0]\n",
    "xlim = [xx_utm[-1,val_ind_row[0]], xx_utm[-1,val_ind_row[-1]]]\n",
    "ylim = [yy_utm[val_ind_col[-1],0], yy_utm[val_ind_col[0],0]]\n",
    "\n",
    "# replot with new coordinates\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,7))\n",
    "\n",
    "# East velocities\n",
    "im = axs[0].scatter(xx_utm.flatten(), yy_utm.flatten(), s=2, c=vel_para.flatten(), cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[0].plot(fault_trace_utm[:,0], fault_trace_utm[:,1], color=\"red\")\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='Fault-parallel velocity (mm/yr)')\n",
    "axs[0].set_aspect('equal', 'box')\n",
    "axs[0].set_title('Fault-parallel')\n",
    "axs[0].set_xlabel('x-coord (km)')\n",
    "axs[0].set_ylabel('y-coord (km)')\n",
    "axs[0].set_xlim(xlim)\n",
    "axs[0].set_ylim(ylim)\n",
    "\n",
    "# Up velocities\n",
    "im = axs[1].scatter(xx_utm.flatten(), yy_utm.flatten(), s=2, c=vel_U.flatten(), cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[1].plot(fault_trace_utm[:,0], fault_trace_utm[:,1], color=\"red\")\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='Vertical velocity (mm/yr)')\n",
    "axs[1].set_aspect('equal', 'box')\n",
    "axs[1].set_title('Up')\n",
    "axs[1].set_xlabel('x-coord (km)')\n",
    "axs[1].set_ylabel('y-coord (km)')\n",
    "axs[1].set_xlim(xlim)\n",
    "axs[1].set_ylim(ylim)\n",
    "\n",
    "fig.tight_layout(w_pad=5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the profile parameters\n",
    "\n",
    "# rough middle profile\n",
    "#####################################\n",
    "prof_start = (573, 4450)\n",
    "prof_end = (540, 4595)\n",
    "prof_params = {\n",
    "    \"nbins\": 100, # number of bins to split the profile into\n",
    "    \"width\": 5  # total width of the profile in km\n",
    "}\n",
    "#####################################\n",
    "\n",
    "# run the profiler\n",
    "bin_val, prof_bin_mids, points_val, points_dist, points_poly \\\n",
    "    = lib.profile_data(xx_utm,yy_utm,vel_para,prof_start,prof_end,prof_params)\n",
    "\n",
    "# calculate profile-fault intersection angle\n",
    "intersect_dist, intersect_angle = lib.profile_fault_intersection(prof_start,prof_end,fault_trace_utm)\n",
    "\n",
    "# plot\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,5))\n",
    "\n",
    "# East velocities\n",
    "im = axs[0].scatter(xx_utm.flatten(), yy_utm.flatten(), s=2, c=vel_para.flatten(), cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[0].plot(fault_trace_utm[:,0], fault_trace_utm[:,1], color=\"red\")\n",
    "axs[0].plot([prof_start[0], prof_end[0]], [prof_start[1], prof_end[1]], color=\"red\")\n",
    "axs[0].plot(points_poly[:,0],points_poly[:,1], color=\"red\")\n",
    "axs[0].scatter(prof_start[0],prof_start[1], s=100, color='pink', edgecolor='black', zorder=3)\n",
    "axs[0].scatter(prof_end[0],prof_end[1], s=100, color='black', edgecolor='white', zorder=3)\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='Fault-parallel velocity (mm/yr)')\n",
    "axs[0].set_aspect('equal', 'box')\n",
    "axs[0].set_title('Fault-parallel')\n",
    "axs[0].set_xlabel('x-coord (km)')\n",
    "axs[0].set_ylabel('y-coord (km)')\n",
    "axs[0].set_xlim(xlim)\n",
    "axs[0].set_ylim(ylim)\n",
    "\n",
    "# profile\n",
    "axs[1].scatter(points_dist,points_val)\n",
    "axs[1].plot(prof_bin_mids,bin_val,color=\"red\")\n",
    "axs[1].plot([intersect_dist, intersect_dist],[axs[1].get_ylim()[0], axs[1].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[1].plot([axs[1].get_xlim()[0], axs[1].get_xlim()[1]], [0, 0], color='grey', linestyle='dashed')\n",
    "axs[1].text(0.02, 0.90, 'intersection angle = ' + str(round(np.abs(intersect_angle))) + ' degrees', fontsize=14, transform = axs[1].transAxes)\n",
    "axs[1].set_xlabel(\"Distance along profile (km)\")\n",
    "axs[1].set_ylabel(\"Fault-parallel velocity (mm/yr)\")\n",
    "axs[1].yaxis.set_label_position(\"right\")\n",
    "axs[1].yaxis.tick_right()\n",
    "\n",
    "fig.tight_layout(w_pad=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I've included the intersection angle between the fault and the profile line on the right-hand plot.\n",
    "We want this angle to be 90&deg;, otherwise the signal will become smeared and effect our results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Forward model\n",
    "\n",
    "Now that we have our data, let's produce a simple forward model to get an idea of what model parameters produce a reasonable fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slip, locking depth, and offset\n",
    "#####################################\n",
    "s = 15 # mm/yr\n",
    "d = 15 # km\n",
    "c = -5 # mm/yr\n",
    "#####################################\n",
    "\n",
    "# shift profile so that the fault is at 0 km\n",
    "x_prof = points_dist - intersect_dist\n",
    "\n",
    "# Vector of positions every 1 km for -100 km to 100 km\n",
    "x = np.arange(-100,101,1)\n",
    "\n",
    "# run the forward model\n",
    "v = lib.screw_disc(x*1000, s/1000, d*1000, c/1000)\n",
    "\n",
    "# calculate rms misfit by running the forward model for x positions with velocities\n",
    "v_forward = lib.screw_disc(x_prof*1000, s/1000, d*1000, c/1000)\n",
    "rms_forward = lib.rms_misfit(points_val,v_forward*1000)\n",
    "\n",
    "# Plot comparison\n",
    "fig, axs = plt.subplots(1,1,figsize=(15,8))\n",
    "\n",
    "plt.scatter(x_prof,points_val)\n",
    "plt.plot(x, v*1000, c='r')\n",
    "plt.plot([0, 0],[axs.get_ylim()[0], axs.get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "plt.plot([axs.get_xlim()[0], axs.get_xlim()[1]], [0, 0], color='grey', linestyle='dashed')\n",
    "\n",
    "plt.text(0.15, 0.9, 'slip rate = ' + str(s) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.15, 0.84, 'locking depth = ' + str(d) + ' km', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.15, 0.78, 'offset = ' + str(c) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.15, 0.72, 'RMS misfit = ' + str(round(rms_forward,3)) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "\n",
    "plt.xlabel('Perpendicular distance from surface fault trace (km)')\n",
    "plt.ylabel('Fault-parallel velocity (mm/yr)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how well you can all model the fault, lets model the same profile and then share the results.\n",
    "\n",
    "Set your profile to the following:\n",
    "```\n",
    "prof_start = (573, 4450)\n",
    "prof_end = (540, 4595)\n",
    "```\n",
    "\n",
    "Once you have your best model, let me know your parameter values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Bayesian inversion\n",
    "\n",
    "Bayesian inversions are one way we can estimate model parameters and their uncertainties in geophysics. We take some prior knowledge about the model (e.g. limits of the model parameter values) and combine this with a likelihood function based on the observations to estimate the posterior probability (i.e. the likelihood of our model parameters being the true values given the data). By exploring a range of model parameter values, we can build up a a posterior probability distribution which tells us both our \"best fit\" values, known as the \"maximum a posteriori probability (MAP)\" solution, and allows us to estimate the uncertainties on this solution.\n",
    "\n",
    "We'll assume uniform priors for each model parameter, meaning all values between the upper and lower limits are equally likely. Values outside of these limits are \"impossible\" and will be rejected in our inversion.\n",
    "\n",
    "We use the following limits for our uniform priors:\n",
    "- $0 \\leq s \\leq 40$ mm/yr\n",
    "- $1 \\leq d \\leq 50$ km\n",
    "- $-10 \\leq c \\leq 10$ mm/yr\n",
    "\n",
    "Assuming our errors are Gaussian, we can find the maximum likelihood by minimising the weighted difference between the observed velocities, $\\mathbf d$, and the model velocities, $g(\\mathbf{m})$:\n",
    "\\begin{equation}\n",
    "    \\sum (\\mathbf d-g(\\mathbf{m}))*W*(\\mathbf d-g(\\mathbf{m}))\n",
    "\\end{equation}\n",
    "where $W$ contains the weight for every velocity. Weighting the data is useful when some velocities have larger errors than others. For simplicity, we'll assume a uniform uncertainty of 1 mm/yr (1 standard deviation) across all of our velocities.\n",
    "\n",
    "In order to explore the posterior distribution (i.e. test a large number of model parameter values) we'll use a Markov Chain Monte Carlo (MCMC) method.\n",
    "Our MCMC performs the following steps:\n",
    "- Take a starting model and add a random step to each parameters.\n",
    "- Test whether this new model is within the prior limits, rejecting the model if not.\n",
    "- Calculate the likelihoods for both models.\n",
    "- If the new model has a higher likelihood, then accept this model. Also accept a number of models that are worse to help the Monte Carlo explore the full parameter space.\n",
    "\n",
    "We need to allow the MCMC to sometimes accept less likely models to avoid local minima. These are models that are good compared to the surrounding models (in the parameter space), but are not the true best model (called the global minima)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "importlib.reload(lib)\n",
    "np.random.seed(1)\n",
    "\n",
    "# inversion setup\n",
    "#####################################\n",
    "n_iterations = 10000 # don't go above this value of 10,000 in the practical\n",
    "#####################################\n",
    "\n",
    "# take profiled velocities and shift the fault to be at zero\n",
    "x_prof = points_dist - intersect_dist\n",
    "v_prof = points_val\n",
    "\n",
    "# model params and their limits\n",
    "m_start = np.array([5, 5, 0], dtype=np.double) # slip in mm/yr, locking depth in km\n",
    "m_min = np.array([0, 1, -10], dtype=np.double)\n",
    "m_max = np.array([40, 50, 10], dtype=np.double)\n",
    "\n",
    "# variance-covariance matrix, currently 1 mm/yr uncertainties\n",
    "vcm = np.eye(len(x_prof))\n",
    "\n",
    "# convert all units to metres\n",
    "x = x_prof * 1000\n",
    "v = v_prof / 1000\n",
    "m_start = m_start * np.array([0.001, 1000, 0.001])\n",
    "m_min = m_min * np.array([0.001, 1000, 0.001])\n",
    "m_max = m_max * np.array([0.001, 1000, 0.001])\n",
    "vcm = vcm / 1e6\n",
    "\n",
    "# inversion setup\n",
    "burn_in = round(n_iterations/100)\n",
    "acceptance_coef = 0.001\n",
    "\n",
    "# weightings\n",
    "W = np.linalg.inv(vcm)\n",
    "\n",
    "# pre-allocate arrays\n",
    "models_saved=np.zeros((n_iterations,len(m_start))) * np.nan\n",
    "ll_saved=np.zeros((n_iterations,1)) * np.nan\n",
    "n_accept = 0 # number of accepted models\n",
    "n_reject = 0 # number of rejected models\n",
    "\n",
    "# pass starting model to current model\n",
    "m_current = m_start.copy()\n",
    "\n",
    "# run inversion\n",
    "for ii in range(n_iterations):\n",
    "    \n",
    "    # propose model using different step sizes for each parameter\n",
    "    m_trial = m_current.copy()\n",
    "    m_trial[0] = m_trial[0] + np.random.normal(loc=0, scale=2.5, size=1)/1000 # slip rate\n",
    "    m_trial[1] = m_trial[1] + np.random.normal(loc=0, scale=2.5, size=1)*1000 # locking depth\n",
    "    m_trial[2] = m_trial[2] + np.random.normal(loc=0, scale=1, size=1)/1000 # offset\n",
    "    \n",
    "    # check limits and skip the rest of the loop if any parameter is invalid\n",
    "    if not(lib.prior(m_trial,m_min,m_max)):\n",
    "        n_reject += 1\n",
    "        models_saved[ii,:] = m_current\n",
    "        continue\n",
    "    \n",
    "    # calculate likelihood for the current and trial models\n",
    "    ll_current = lib.likelihood(x, v, m_current, W)\n",
    "    ll_trial = lib.likelihood(x, v, m_trial, W)\n",
    "    \n",
    "    # test whether to keep trial model\n",
    "    #if np.exp(acceptance_coef*(ll_trial-ll_current)) > np.random.uniform(low=0, high=1,size=1):\n",
    "    if np.exp(acceptance_coef*(ll_current-ll_trial)) > np.random.uniform(low=0, high=1,size=1):\n",
    "        m_current = m_trial\n",
    "        ll_current = ll_trial\n",
    "        n_accept += 1\n",
    "    else:\n",
    "        n_reject += 1\n",
    "    models_saved[ii,:] = m_current\n",
    "    ll_saved[ii] = ll_current\n",
    "    \n",
    "# convert back from metres to mm/yr and km\n",
    "models_saved[:,0] = models_saved[:,0] * 1000 # slip rate\n",
    "models_saved[:,1] = models_saved[:,1] / 1000 # locking depth\n",
    "models_saved[:,2] = models_saved[:,2] * 1000 # offset\n",
    "\n",
    "# find best fit model using min of likelihood function\n",
    "best_model = models_saved[np.nanargmin(ll_saved),:]\n",
    "\n",
    "print('Inversion complete')\n",
    "print('Number of accepted models = ' + str(n_accept))\n",
    "print('Number of rejected models = ' + str(n_reject))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll plot the accepted trial models for the slip rate and locking depth to see how the algorithm explores the parameter space.\n",
    "We'll show the first 1% of iterations in red and the rest in blue.\n",
    "These early iterations are often termed a \"burn-in\" period, where the inversion is moving from the start position towards the (global) minimum, which it then begins to explore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot walker \n",
    "plt.figure()\n",
    "plt.plot(models_saved[:burn_in,0], models_saved[:burn_in,1], color='red', alpha=0.3)\n",
    "plt.plot(models_saved[burn_in:,0], models_saved[burn_in:,1], color='blue', alpha=0.3)\n",
    "plt.xlabel('Slip rate (mm/yr)')\n",
    "plt.ylabel('Locking depth (km)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at how the inversion explored the parameter space by plotting the likelihood function value for each iteration. The lowest likelihood will be located within a minima, while the higher likelihoods throughout show the inversion accepting worse models to help it escape local minima."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot likelihood value throughout monte carlo\n",
    "fig, axs = plt.subplots(1,1,figsize=(15,8))\n",
    "\n",
    "axs.scatter(np.arange(1,len(ll_saved)+1),ll_saved, s=3)\n",
    "axs.set_xlabel('Number of iterations')\n",
    "axs.set_ylabel('Weighted residual')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, lets plot our best model on top of the profiled velocities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate rms misfit\n",
    "v_final = lib.screw_disc(x_prof, best_model[0], best_model[1], best_model[2])\n",
    "best_rms = lib.rms_misfit(v_prof,v_final)\n",
    "\n",
    "# for plotting\n",
    "x = np.arange(-100,101,1)\n",
    "v_final_plot = lib.screw_disc(x, best_model[0], best_model[1], best_model[2])\n",
    "\n",
    "# Plot comparison\n",
    "fig, axs = plt.subplots(1,1,figsize=(15,8))\n",
    "\n",
    "plt.scatter(x_prof,v_prof)\n",
    "plt.plot(x, v_final_plot, c='r')\n",
    "plt.plot([0, 0],[axs.get_ylim()[0], axs.get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "\n",
    "plt.text(0.15, 0.8, 'slip = ' + str(round(best_model[0],2)) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.15, 0.74, 'locking depth = ' + str(round(best_model[1],2)) + ' km', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.15, 0.68, 'offset = ' + str(round(best_model[2],2)) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.15, 0.62, 'RMS misfit = ' + str(round(best_rms,3)) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "\n",
    "plt.xlabel('Perp. distance from fault (km)')\n",
    "plt.ylabel('Fault-parallel velocity (mm/yr)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll plot histograms for each model parameter. The overall shape of these histograms should be gaussian, although they may be skewed to one side. They may also be quite rough depending on the number of iterations that have been run. By fitting a gaussian to each histogram, we can estimate the standard deviation of each, which is a measure of the uncertainty on our estimates.\n",
    "\n",
    "We'll plot the best model values as vertical red lines. These should roughly line up with the peaks of our histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate mean and standard deviation for each parameter\n",
    "gauss_s = stats.norm.fit(models_saved[burn_in:,0])\n",
    "gauss_d = stats.norm.fit(models_saved[burn_in:,1])\n",
    "gauss_c = stats.norm.fit(models_saved[burn_in:,2])\n",
    "\n",
    "# create figure\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15,8))\n",
    "\n",
    "# plot histograms\n",
    "_, bins_s, _ = axs[0].hist(models_saved[burn_in:,0], bins=20, density=True)\n",
    "_, bins_d, _ = axs[1].hist(models_saved[burn_in:,1], bins=20, density=True)\n",
    "_, bins_c, _ = axs[2].hist(models_saved[burn_in:,2], bins=20, density=True)\n",
    "\n",
    "# create gaussian fit lines\n",
    "gauss_fit_s = stats.norm.pdf(bins_s, gauss_s[0], gauss_s[1])\n",
    "gauss_fit_d = stats.norm.pdf(bins_d, gauss_d[0], gauss_d[1])\n",
    "gauss_fit_c = stats.norm.pdf(bins_c, gauss_c[0], gauss_c[1])\n",
    "\n",
    "# plot gaussian fits\n",
    "axs[0].plot(bins_s, gauss_fit_s)\n",
    "axs[1].plot(bins_d, gauss_fit_d)\n",
    "axs[2].plot(bins_c, gauss_fit_c)\n",
    "\n",
    "# label mean and sd\n",
    "axs[0].text(0.02, 0.97, 'mean = ' + str(round(gauss_s[0],2)), fontsize=14, transform = axs[0].transAxes)\n",
    "axs[0].text(0.02, 0.94, 'sd = ' + str(round(gauss_s[1],2)), fontsize=14, transform = axs[0].transAxes)\n",
    "axs[1].text(0.02, 0.97, 'mean = ' + str(round(gauss_d[0],2)), fontsize=14, transform = axs[1].transAxes)\n",
    "axs[1].text(0.02, 0.94, 'sd = ' + str(round(gauss_d[1],2)), fontsize=14, transform = axs[1].transAxes)\n",
    "axs[2].text(0.02, 0.97, 'mean = ' + str(round(gauss_c[0],2)), fontsize=14, transform = axs[2].transAxes)\n",
    "axs[2].text(0.02, 0.94, 'sd = ' + str(round(gauss_c[1],2)), fontsize=14, transform = axs[2].transAxes)\n",
    "\n",
    "# plot best model values\n",
    "axs[0].plot([best_model[0], best_model[0]],[0, axs[0].get_ylim()[1]], color=\"red\")\n",
    "axs[1].plot([best_model[1], best_model[1]],[0, axs[1].get_ylim()[1]], color=\"red\")\n",
    "axs[2].plot([best_model[2], best_model[2]],[0, axs[2].get_ylim()[1]], color=\"red\")\n",
    "\n",
    "axs[0].set_title('Slip rate', fontsize=16)\n",
    "axs[1].set_title('Locking depth', fontsize=16)\n",
    "axs[2].set_title('Offset', fontsize=16)\n",
    "\n",
    "axs[0].set_xlabel('Slip rate (mm/yr)', fontsize=14)\n",
    "axs[1].set_xlabel('Locking depth (km)', fontsize=14)\n",
    "axs[2].set_xlabel('Offset (mm/yr)', fontsize=14)\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we saw in the \"walker\" plot, there is a trade off between slip rate and locking depth, with an increase in one producing an increase in the other.\n",
    "\n",
    "Lets replot the histograms above while also show the distribution of each model parameter against the other. This should make correlations between model parameters clearer to see. We'll use the Seaborn package for this, as it has a pre-made function for generate these multivariate plots. Its well worth an investigate if you want to use python to plot more scientific data in future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This plot can take a couple minutes to run (dependent on the number of iterations)\n",
    "\n",
    "# multivariate plot\n",
    "models_df = pd.DataFrame(data=models_saved, columns=[\"Slip rate (mm/yr)\", \"Locking depth (km)\", \"Offset (mm/yr)\"])\n",
    "g = sns.PairGrid(models_df, corner=True)\n",
    "g.map_diag(sns.histplot)\n",
    "g.map_lower(sns.kdeplot) # use this to plot as contours instead (runs slower)\n",
    "g.fig.set_size_inches(13,13)\n",
    "\n",
    "# add red dots for our best model\n",
    "g.axes[1,0].scatter(best_model[0], best_model[1], c='red')\n",
    "g.axes[2,0].scatter(best_model[0], best_model[2], c='red')\n",
    "g.axes[2,1].scatter(best_model[1], best_model[2], c='red')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare just slip rate and locking depth\n",
    "models_df = pd.DataFrame(data=models_saved[:,0:2], columns=[\"Slip rate (mm/yr)\", \"Locking depth (km)\"])\n",
    "ax = sns.jointplot(data=models_df, x=\"Slip rate (mm/yr)\", y=\"Locking depth (km)\", kind=\"hex\", height=8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Extension - Fault creep\n",
    "\n",
    "So far we've assumed that the NAF is locked from the surface down to a given depth, below which the fault is slipping aseismically. This gives a fault-parallel velocity of zero at the fault trace. Depending on where you have drawn your profile lines so far, you may have noticed a step in the velocities across the fault that is larger than the noise level. This may be fault creep, where the fault is also slipping aseismically at a shallow depth up to the surface. This fault creep occurs alongside the deeper interseismic slip.\n",
    "\n",
    "A schematic diagram of fault creep is shown in Figure 5.\n",
    "\n",
    "<img src=\"figures/faultCreep.png\" style=\"float: centre;\" width=\"500\"/>\n",
    "\n",
    "*Figure 5: (a) Schematic diagram of our fault creep model. (b) Forward model showing how the fault creep is expected to change our observed signal. From Hussian et al. (2016).*\n",
    "\n",
    "We can model a combination of interseismic strain accumulation and fault creep to calculate a fault parallel velocity, $v$, at fault-perpendicular distance, $x$, using Equation 4:\n",
    "\n",
    "\\begin{equation}\n",
    "v(x) = -\\frac{s}{\\pi}arctan\\left(\\frac{x+x_c}{d_1}\\right) + C\\left[\\frac{1}{\\pi}arctan\\left(\\frac{x+x_c}{d_2}\\right) - H(x+x_c)\\right] + c \\qquad \\qquad (4)\n",
    "\\end{equation}\n",
    "\n",
    "We retain the screw dislocation model from before, where slip, $s$ occurs below a depth, $d_1$. For the shallow creep, occuring at a rate of $C$ between the surface and depth $d_2$, we use a back slip approach. This models the creep as the sum of slip on the entire fault plane using a Heaviside function ($H$) plus a screw dislocation in the opposite sense to the plate motion at depth $d_2$. We also retain the static offset, $c$, and include a new offset to the fault location, $x_c$. Sometimes the observed signal may not exactly align with the mapped fault trace, and so it can be useful to include a small offset to the location of the fault.\n",
    "\n",
    "Let's take a new profile across a creeping part of the NAF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creep\n",
    "prof_start = (508, 4450)\n",
    "prof_end = (480, 4620)\n",
    "\n",
    "prof_params = {\n",
    "    \"nbins\": 100, # number of bins to split the profile into\n",
    "    \"width\": 5  # total width of the profile in km\n",
    "}\n",
    "\n",
    "# run the profiler\n",
    "bin_val, prof_bin_mids, points_val, points_dist, points_poly \\\n",
    "    = lib.profile_data(xx_utm,yy_utm,vel_para,prof_start,prof_end,prof_params)\n",
    "\n",
    "# calculate profile-fault intersection angle\n",
    "intersect_dist, intersect_angle = lib.profile_fault_intersection(prof_start,prof_end,fault_trace_utm)\n",
    "\n",
    "# plot\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,5))\n",
    "\n",
    "# East velocities\n",
    "im = axs[0].scatter(xx_utm.flatten(), yy_utm.flatten(), s=2, c=vel_para.flatten(), cmap=cm.vik, vmin=-20, vmax=20)\n",
    "axs[0].plot(fault_trace_utm[:,0], fault_trace_utm[:,1], color=\"red\")\n",
    "axs[0].plot([prof_start[0], prof_end[0]], [prof_start[1], prof_end[1]], color=\"red\")\n",
    "axs[0].plot(points_poly[:,0],points_poly[:,1], color=\"red\")\n",
    "axs[0].scatter(prof_start[0],prof_start[1], color='red')\n",
    "axs[0].scatter(prof_end[0],prof_end[1], color='black')\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"2%\")\n",
    "plt.colorbar(im, cax=cax, label='Fault-parallel velocity (mm/yr)')\n",
    "axs[0].set_aspect('equal', 'box')\n",
    "axs[0].set_title('Fault-parallel')\n",
    "axs[0].set_xlim(xlim)\n",
    "axs[0].set_ylim(ylim)\n",
    "\n",
    "# profile\n",
    "axs[1].scatter(points_dist,points_val)\n",
    "axs[1].plot(prof_bin_mids,bin_val,color=\"red\")\n",
    "axs[1].plot([intersect_dist, intersect_dist],[axs[1].get_ylim()[0], axs[1].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[1].text(0.02, 0.90, 'intersection angle = ' + str(round(np.abs(intersect_angle))) + ' degrees', fontsize=14, transform = axs[1].transAxes)\n",
    "axs[1].set_xlabel(\"Distance\")\n",
    "axs[1].set_ylabel(\"Fault-parallel velocity (mm/yr)\")\n",
    "axs[1].yaxis.set_label_position(\"right\")\n",
    "axs[1].yaxis.tick_right()\n",
    "\n",
    "plt.tight_layout(w_pad=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll run a forward model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slip and locking depth\n",
    "#####################################\n",
    "s1 = 20 # mm/yr\n",
    "s2 = 7 # mm/yr\n",
    "d1 = 15 # km\n",
    "d2 = 5 # km\n",
    "c = -7 # mm/yr\n",
    "xc = 0 # km\n",
    "#####################################\n",
    "\n",
    "# shift profile so thast fault is at 0 km\n",
    "x_prof = points_dist - intersect_dist\n",
    "\n",
    "# Vector of positions every 1 km for -100 km to 100 km\n",
    "x = np.arange(-80,81,0.1)\n",
    "\n",
    "# run the forward model\n",
    "v = lib.fault_creep(x*1000, s1/1000, s2/1000, d1*1000, d2*1000, c/1000, xc=xc*1000)\n",
    "\n",
    "# calculate rms misfit by running the forward model for x positions with velocities\n",
    "v_forward = lib.fault_creep(x_prof*1000, s1/1000, s2/1000, d1*1000, d2*1000, c/1000, xc=xc*1000)\n",
    "rms_forward = lib.rms_misfit(points_val,v_forward*1000)\n",
    "\n",
    "# Plot comparison\n",
    "fig, axs = plt.subplots(1,1,figsize=(15,8))\n",
    "\n",
    "plt.scatter(x_prof,points_val)\n",
    "plt.plot(x, v*1000, c='r')\n",
    "plt.plot([0, 0],[axs.get_ylim()[0], axs.get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "\n",
    "plt.text(0.02, 0.9, 'slip = ' + str(s1) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.25, 0.9, 'creep = ' + str(s2) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.02, 0.84, 'locking depth = ' + str(d1) + ' km', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.25, 0.84, 'creep depth = ' + str(d2) + ' km', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.02, 0.78, 'offset = ' + str(c) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "plt.text(0.02, 0.65, 'RMS misfit = ' + str(round(rms_forward,3)) + ' mm/yr', fontsize=14, transform = axs.transAxes)\n",
    "\n",
    "plt.xlabel('Perp. distance from fault (km)')\n",
    "plt.ylabel('Fault-parallel velocity (mm/yr)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Extension - Strain Rate\n",
    "\n",
    "Estimates of strain rate can highligh areas of concentrated deformation, an important product for understanding both seismic hazard the wider-scale kinematics. The state of strain in the lithosphere can be described using a strain rate tensor, $\\dot{\\varepsilon}$:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\dot{\\varepsilon} = \\begin{bmatrix}\n",
    "    \\frac{\\partial V_E}{\\partial x} & \\frac{1}{2}\\left(\\frac{\\partial V_E}{\\partial y}+\\frac{\\partial V_N}{\\partial x}\\right) & \\frac{1}{2}\\left(\\frac{\\partial V_E}{\\partial z}+\\frac{\\partial V_U}{\\partial x}\\right) \\\\    \n",
    "    \\frac{1}{2}\\left(\\frac{\\partial V_N}{\\partial x}+\\frac{\\partial V_E}{\\partial y}\\right) & \\frac{\\partial V_N}{\\partial y} & \\frac{1}{2}\\left(\\frac{\\partial V_N}{\\partial z}+\\frac{\\partial V_U}{\\partial y}\\right) \\\\    \n",
    "    \\frac{1}{2}\\left(\\frac{\\partial V_U}{\\partial x}+\\frac{\\partial V_E}{\\partial z}\\right) & \\frac{1}{2}\\left(\\frac{\\partial V_U}{\\partial y}+\\frac{\\partial V_N}{\\partial z}\\right) & \\frac{\\partial V_U}{\\partial z}    \n",
    "    \\end{bmatrix} \\qquad \\qquad (5)\n",
    "\\end{equation}\n",
    "\n",
    "where $v_E$, $v_N$, and $v_U$ are the velocities in the Cartesian coordinate direction, $x$, $y$, and $z$.\n",
    "\n",
    "In this instance, we have fault-parellel velocities at distances perpendicular to the fault, and so we can only estimate the shear strain across the fault. For our screw dislocation model, we can forward model the shear strain, $e_{shear}$, using (Savage & Burford, 1973):\n",
    "\n",
    "\\begin{equation}\n",
    "e_{shear} = \\left(\\frac{sd}{2\\pi}\\right)\\left(x^2+d^2\\right)^{-1} \\qquad \\qquad (6)\n",
    "\\end{equation}\n",
    "\n",
    "Let's start by calculating the shear rate for the screw dislocation model that we fit to our original profile, using both a forward model (Equation 6) and by calculating the gradient in the modelled velocities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slip rate and locking depth\n",
    "#####################################\n",
    "s = 19 # mm/yr\n",
    "d = 11 # km\n",
    "c = -5.7 # mm/yr\n",
    "#####################################\n",
    "\n",
    "# Original profile line\n",
    "prof_start = (573, 4450)\n",
    "prof_end = (540, 4595)\n",
    "prof_params = {\n",
    "    \"nbins\": 100, # number of bins to split the profile into\n",
    "    \"width\": 5  # total width of the profile in km\n",
    "}\n",
    "\n",
    "# run the profiler\n",
    "bin_val, prof_bin_mids, points_val, points_dist, points_poly \\\n",
    "    = lib.profile_data(xx_utm,yy_utm,vel_para,prof_start,prof_end,prof_params)\n",
    "\n",
    "# calculate profile-fault intersection angle\n",
    "intersect_dist, intersect_angle = lib.profile_fault_intersection(prof_start,prof_end,fault_trace_utm)\n",
    "\n",
    "# shift profile so that the fault is at 0 km\n",
    "x_prof = points_dist - intersect_dist\n",
    "\n",
    "# Vector of positions every 1 km for -100 km to 100 km\n",
    "x = np.arange(-100,101,1)\n",
    "\n",
    "# run the forward model\n",
    "v = lib.screw_disc(x*1000, s/1000, d*1000, c/1000)\n",
    "\n",
    "# calculate shear strain rate using Equation 6\n",
    "e_shear = lib.shear_strain(x*1000, s/1000, d*1000)\n",
    "\n",
    "# also calculate shear strain rate based purely on velocity gradient\n",
    "e_shear_grad = np.gradient(v,x*1000)\n",
    "\n",
    "# Plot the original profile and model, and strain rate\n",
    "fig, axs = plt.subplots(2,1,figsize=(15,15))\n",
    "\n",
    "axs[0].scatter(x_prof, points_val)\n",
    "axs[0].plot(x, v*1000, c='r')\n",
    "axs[0].plot([0, 0],[axs[0].get_ylim()[0], axs[0].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[0].text(0.02, 0.90, 'intersection angle = ' + str(round(np.abs(intersect_angle))) + ' degrees', fontsize=14, transform = axs[0].transAxes)\n",
    "axs[0].set_xlabel(\"Distance from fault (km)\")\n",
    "axs[0].set_ylabel(\"Fault-parallel velocity (mm/yr)\")\n",
    "axs[0].yaxis.set_label_position(\"right\")\n",
    "axs[0].yaxis.tick_right()\n",
    "\n",
    "axs[1].plot(x, e_shear, c='r', label='Forward model')\n",
    "axs[1].plot(x, e_shear_grad, c='b', label='Velocity gradient')\n",
    "axs[1].plot([0, 0],[axs[1].get_ylim()[0], axs[1].get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "axs[1].legend(fontsize='x-large')\n",
    "axs[1].set_xlabel(\"Distance from fault (km)\")\n",
    "axs[1].set_ylabel(\"Shear strain rate (/yr)\")\n",
    "axs[1].yaxis.set_label_position(\"right\")\n",
    "axs[1].yaxis.tick_right()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get nice smooth estimates of strain rate from our screw dislocation models. In real studies, we often want to estimate strain rates directly from our velocities, as this will capture signals that are not represented in our model. However, we still require some degree of smoothing to suppress short-wavelength noise in the InSAR. This noise, while small compared to the tectonic signal, can cause large gradients between adjacent points that then drown out the underying tectonic strain rate.\n",
    "\n",
    "Let's plot both the strain rate directly calculated from the binned profiled fault-parallel velocities, and one calculated from the same velocities but smoothed with a sliding window average. Try changing the size of the smoothing window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Window size for smoothing, must be an odd number\n",
    "#####################################\n",
    "wind_size = 9\n",
    "#####################################\n",
    "\n",
    "# shift fault location to 0 km\n",
    "x_prof = prof_bin_mids - intersect_dist\n",
    "\n",
    "# calcualte gradient of original binned velocities\n",
    "e_shear = np.gradient(bin_val/1000, x_prof*1000)\n",
    "\n",
    "# apply a sliding mean window to smooth out some of the noise\n",
    "x_smoothed, e_shear_smoothed = lib.sliding_window_mean(x_prof, bin_val, 11)\n",
    "\n",
    "# calculated the gradient from the smoothed velocities\n",
    "e_shear_smoothed = np.gradient(e_shear_smoothed/1000, x_smoothed*1000)\n",
    "\n",
    "# Plot comparison\n",
    "fig, axs = plt.subplots(1,1,figsize=(15,8))\n",
    "\n",
    "plt.plot(x_prof, e_shear, c='r', label='Original binned velocities')\n",
    "plt.plot(x_smoothed, e_shear_smoothed, c='b', label='Smoothed velocities')\n",
    "plt.plot([0, 0],[axs.get_ylim()[0], axs.get_ylim()[1]], color='grey', linestyle='dashed')\n",
    "plt.xlabel('Distance from fault (km)')\n",
    "plt.ylabel('Shear strain rate (/yr)')\n",
    "plt.legend(fontsize=\"x-large\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we calculate the strain rate directly from the binned velocities (which have themselves already been smoothed to some degree), then the strain signal we would expect to see based on Equation 5 is completely masked.\n",
    "\n",
    "By applying a smoothing window (a window size of 9 works quite well), we can bring out the higher strain across the fault trace. If the window size is increased too far, then we can actually smooth out the original interseismic signal, and this peak disappears again. How we filter out high-frequency noise in the velocities is integral to our estimates of strain rate.\n",
    "\n",
    "For further reading, I recommend Weiss et al. (2020), who apply the Velmap method to all of Turkey, and Ou et al. (2022), who use a median filter method to localise strain onto a number of faults in the Tibetan Plateau."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "Hopefully, you now have a basic understanding of how we measure and model interseismic deformation. This practical has been shared under an open-access licence (see repo), so feel free to pass it on and utilise any of the underlying code. \n",
    "\n",
    "If you are interested in performing velocity decompositions on a larger scale, have a look at our open-access Matlab package:\n",
    "https://github.com/andwatson/decompose_insar_velocities\n",
    "\n",
    "If you are interested in a more developed version of the bayesian inversion, have a look at both Hussain et al. (2016) and the GBIS inversion software, for interseismic and coseismic deformation, respectively:\n",
    "https://comet.nerc.ac.uk/gbis/\n",
    "\n",
    "If you found any issues in this practical, please log them under \"issues\" on the github repository.\n",
    "\n",
    "Andrew Watson - 2022"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "Elliott, J. (2009). Strain accumulation & release on the Tibetan Plateau measured using InSAR (Doctoral dissertation, Oxford University).\n",
    "\n",
    "Hussain, E., Wright, T. J., Walters, R. J., Bekaert, D., Hooper, A., & Houseman, G. A. (2016). Geodetic observations of postseismic creep in the decade after the 1999 Izmit earthquake, Turkey: Implications for a shallow slip deficit. Journal of Geophysical Research: Solid Earth, 121(4), 2980-3001.\n",
    "\n",
    "Ou, Q., Daout, S., Weiss, J. R., Shen, L., Lazecký, M., Wright, T. J., & Parsons, B. E. (2022). Large‐Scale Interseismic Strain Mapping of the NE Tibetan Plateau From Sentinel‐1 Interferometry. Journal of Geophysical Research: Solid Earth, 127(6), e2022JB024176.\n",
    "\n",
    "Savage, J. C., & Burford, R. O. (1973). Geodetic determination of relative plate motion in central California. Journal of Geophysical Research, 78(5), 832-845.\n",
    "\n",
    "Vernant, P. (2015). What can we learn from 20 years of interseismic GPS measurements across strike-slip faults?. Tectonophysics, 644, 22-39.\n",
    "\n",
    "Weiss, J. R., et al. (2020). High‐resolution surface velocities and strain for Anatolia from Sentinel‐1 InSAR and GNSS data. Geophysical Research Letters, 47(17), e2020GL087376."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}