Volcano_practical.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3def2d0c",
   "metadata": {},
   "source": [
    "<img src=\"figures/logos.png\" style=\"float: centre;\" width=\"800\"/>\n",
    "\n",
    "# Volcano practical\n",
    "\n",
    "## Written for the COMET InSAR Workshop 2023.\n",
    "\n",
    "Welcome to this COMET practical on measuring volcanic deformation with InSAR.\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 open the bender fresh if anything breaks.\n",
    "\n",
    "If no cells are running, check if the notebook says \"Not trusted\" in the top right corner. If so, click this box and accept the pop-up message.\n",
    "\n",
    "### Acknowledgements\n",
    "\n",
    "Detailed guidance on the content was provided by Susanna Ebmeier. Advice on the computational requirements of the practical was provided by Richard Rigby. Thanks to everyone for their contributions.\n",
    "\n",
    "Andrew Watson - 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2f5af21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required modules. Run this before continuing.\n",
    "\n",
    "# these are for performing the analysis\n",
    "import numpy as np\n",
    "import pyproj\n",
    "\n",
    "# these packages are required for 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",
    "# this is our own library of functions, and can be found in the same directory as this practical\n",
    "import volcano_lib as lib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5f14c1c",
   "metadata": {},
   "source": [
    "## 0. Introduction\n",
    "\n",
    "Deformation at volcano can be caused by many processes including the movement of magma, volume changes within a reservoir as well as flow subsidence or faulting. Most InSAR studies of volcanic deformation focus on pressure changes related to magmatic processes such as intrusion or eruption, and provide evidence of magma depth, source geometry and change in pressure. Not all erupting volcanoes deform in a manner that we can measure, and much deformation at volcanoes is not directly related to eruptions. However, unrest involving deformation is more likely to be detectable in pre-eruptive periods than other parameters measurable from satellite imagery such as temperature anomalies and gas emission. Deformation measurements are therefore important not only for research into magmatic processes, but also for monitoring volcanic activity.  \n",
    "\n",
    "Further reading:  \n",
    "- Pinel, V., Poland, M. P., & Hooper, A. (2014). Volcanology: Lessons learned from synthetic aperture radar imagery. Journal of Volcanology and Geothermal Research, 289, 81-113. \n",
    "- Segall, P. (2019). Magma chambers: what we can, and cannot, learn from volcano geodesy. Philosophical Transactions of the Royal Society A, 377(2139), 20180158 \n",
    "- Ebmeier, S.K., Andrews, B.J., Araya, M.C., Arnold, D.W.D., Biggs, J., Cooper, C., Cottrell, E., Furtney, M., Hickey, J., Jay, J.J.J.A.V. and Lloyd, R., 2018. Synthesis of global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring & the lateral extent of magmatic domains. Journal of Applied Volcanology, 7(1), pp.1-26. \n",
    "\n",
    "In this practical we will look at an example of a very simple but widely used analytical deformation source model (Mogi source) and explore the impact of topography and line-of-sight on the predicted displacements. We will then look at a real example from the COMET volcano portal and use forward models to try and match the observed displacements with our simple model.  \n",
    "\n",
    "A Mogi model is a point-source model used to explain surface deformation above an expanding magma chamber.\n",
    "\n",
    "<img src=\"figures/mogi.png\" style=\"float: centre;\" width=\"600\"/>\n",
    "\n",
    "*Figure 1: Schematic diagram of a Mogi source. From Lisowski (2007)*\n",
    "\n",
    "We make the following assumptions:\n",
    "- the source is in an isotropic elastic half space, which is defined by the Poisson's ratio $v$, and the Shear modulus $\\mu$.\n",
    "- the radius of the point source is lower than it depth (can't break the surface).\n",
    "- the magma itself is incompressible.\n",
    "\n",
    "We can express a mogi source as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix} ux \\\\ uy \\\\ uz \\end{pmatrix} = \\alpha^3 \\delta P \\frac{(1-v)}{\\mu} \\begin{pmatrix} x / R^3 \\\\ y / R^3 \\\\ z / R^3 \\end{pmatrix} \\qquad \\qquad (1)\n",
    "\\end{equation}\n",
    "\n",
    "where $ux$, $uy$, and $uz$ are the displacements in the x, y, and z axes, respecitvely. \\alpha is the radius of the source at depth, $d$, measured to the centre of the source. $\\delta P$ is the pressure change in the source. $R$ is the distance between the centre of the source and a position on the ground surface at which the displacements are calculated.\n",
    "\n",
    "By taking the following approximation:\n",
    "\\begin{equation}\n",
    "\\delta V \\approx \\frac{\\delta P}{\\mu} \\pi \\alpha^3\n",
    "\\end{equation}\n",
    "\n",
    "where $\\delta V$ is the change in volume of the source.\n",
    "\n",
    "We can rewrite Equation 1 as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix} ux \\\\ uy \\\\ uz \\end{pmatrix} = \\delta V \\frac{(1-v)}{\\pi} \\begin{pmatrix} x / R^3 \\\\ y / R^3 \\\\ z / R^3 \\end{pmatrix} \\qquad \\qquad (2)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "148837c4",
   "metadata": {},
   "source": [
    "## 1. Synthetic Mogi Source\n",
    "\n",
    "Let's start by generating a simple models for an inflating magma chamber - a Mogi source. You can read more about the assumptions behind the Mogi source and about its sensitivity to these assumptions being broken here:  \n",
    "- Mogi, K. 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Earthq Res Inst, 36, pp.99-134. \n",
    "- Masterlark, T., 2007. Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions. Journal of Geophysical Research: Solid Earth, 112(B6).\n",
    "\n",
    "An advantage of InSAR is that it allows measurements over large areas and at postings of every few metres to tens of metres.  Prior to this, most modelling of volcanic deformation sources was done on the basis of a few measurements from levelling, distance meters or GNSS sites along roughly radial lines – similar to a transect through an interferogram. Even if a full 3D displacement field cannot be retrieved, the relationship between vertical and East-West motion can give us some idea of source type and depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a8ccb3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup model params\n",
    "#####################################\n",
    "source_depth = 3 # depth to point (km)\n",
    "dV = 1e7 # change in volume (m^3)\n",
    "#####################################\n",
    "\n",
    "# setup coord limits in km\n",
    "xmin, xmax, xint = -20, 20, 0.1\n",
    "ymin, ymax, yint = -20, 20, 0.1\n",
    "\n",
    "# make grids of x and y coordinates\n",
    "x = np.arange(xmin, xmax+xint, xint)\n",
    "y = np.arange(ymin, ymax+yint, yint)\n",
    "xx, yy = np.meshgrid(x, y)\n",
    "\n",
    "# poisson's ratio for medium\n",
    "nu = 0.25 \n",
    "\n",
    "# generate displacements using Mogi source (note the unit conversions to metres)\n",
    "ux, uy, uz = lib.mogi_source(xx*1000, yy*1000, d=source_depth*1000, dV=dV, nu=nu)\n",
    "\n",
    "# set limits for colour bar, using different limits for z\n",
    "cmin, cmax = -100, 100\n",
    "cmin_z, cmax_z = 0, 250\n",
    "\n",
    "# plot the results\n",
    "fig, axs = plt.subplots(1,3,figsize=(20,10))\n",
    "fig.tight_layout()\n",
    "\n",
    "# displacement in x\n",
    "im = axs[0].imshow(ux*1000, extent=[xmin, xmax, ymin, ymax], vmin=cmin, vmax=cmax, cmap=cm.vik)\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0].set_title('Displacement in x', fontsize=20)\n",
    "axs[0].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[0].set_ylabel('y-coord (km)', fontsize=18)\n",
    "axs[0].tick_params(labelsize=16)\n",
    "\n",
    "# displacement in y\n",
    "im = axs[1].imshow(uy*1000, extent=[xmin, xmax, ymin, ymax], vmin=cmin, vmax=cmax, cmap=cm.vik)\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[1].set_title('Displacement in y', fontsize=20)\n",
    "axs[1].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[1].tick_params(labelsize=16)\n",
    "axs[1].set_yticklabels([])\n",
    "axs[1].set_yticks([])\n",
    "\n",
    "# displacement in z\n",
    "im = axs[2].imshow(uz*1000, extent=[xmin, xmax, ymin, ymax], vmin=cmin_z, vmax=cmax_z, cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[2])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[2].set_title('Displacement in z', fontsize=20)\n",
    "axs[2].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[2].set_ylabel('y-coord (km)', fontsize=18)\n",
    "axs[2].tick_params(labelsize=16)\n",
    "axs[2].yaxis.tick_right()\n",
    "axs[2].yaxis.set_label_position(\"right\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f692242",
   "metadata": {},
   "source": [
    "Let's also plot these as a single 3D surface displacement map, using colour to show vertical motion, and vectors to show horizontal motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3614780",
   "metadata": {},
   "outputs": [],
   "source": [
    "# downsample ux and uy so that we only draw arrows every 2 km\n",
    "arr_spacing = 2 # km\n",
    "xstep, ystep = round(arr_spacing/xint), round(arr_spacing/yint)\n",
    "ux_ds = ux[0::ystep,0::xstep]\n",
    "uy_ds = uy[0::ystep,0::xstep]\n",
    "xx_ds = xx[0::ystep,0::xstep]\n",
    "yy_ds = yy[0::ystep,0::xstep]\n",
    "\n",
    "# plot figure\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "im = plt.imshow(uz*1000, extent=[xmin, xmax, ymin, ymax], vmin=cmin_z, vmax=cmax_z, cmap=cm.batlow)\n",
    "plt.quiver(xx_ds, yy_ds, ux_ds, uy_ds, color='white')\n",
    "cbar = plt.colorbar(im, shrink=0.8)\n",
    "cbar.set_label(label='Vertical displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "plt.title('3D displacement', fontsize=20)\n",
    "plt.xlabel('x-coord (km)', fontsize=18)\n",
    "plt.ylabel('y-coord (km)', fontsize=18)\n",
    "plt.tick_params(labelsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e45fb76e",
   "metadata": {},
   "source": [
    "Let's also plot a profile through the $x$ and $u$ displacements at $y=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79b87fcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get indexes of y=0\n",
    "y0_ind = np.argmin(np.abs(y))\n",
    "\n",
    "# profile across x=0 using simple indexing\n",
    "prof_uz = uz[y0_ind,:]\n",
    "prof_ux = ux[y0_ind,:]\n",
    "\n",
    "# plot profiles\n",
    "fig = plt.figure(figsize=(14,5))\n",
    "\n",
    "plt.plot(x, prof_uz*1000, color='red', label='Displacement in z')\n",
    "plt.plot(x, prof_ux*1000, color='blue', label='Displacement in x')\n",
    "plt.legend(fontsize=16)\n",
    "plt.title('Profiles along y = 0 km', fontsize=20)\n",
    "plt.xlabel('Distance (km)', fontsize=18)\n",
    "plt.ylabel('Displacement (mm)', fontsize=18)\n",
    "plt.tick_params(labelsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a479be5",
   "metadata": {},
   "source": [
    "## 2. Impact of Satellite Line-of-Sight\n",
    "\n",
    "SAR satellite instruments are side looking, and can make measurements in just one dimension, the satellite’s line of sight (LOS). This means that the actual measurements made with InSAR look different to those predicted by our Mogi model.\n",
    "\n",
    "<img src=\"figures/line-of-sight.png\" style=\"float: centre;\" width=\"600\"/>\n",
    "\n",
    "*Figure 2: The apparent centre of displacements over a volcano will shift depending on the satellite look direction and incidence angle. From Magee et al. (2018)*\n",
    "\n",
    "Let's run the model again but include the effect of LOS, starting with a single value for both the satellite heading and the incidence angle. \n",
    "\n",
    "An incidence angle of 38° is reasonable for both ascending and descending. In real data, this incidence angle varies across the image, but we'll stick to a fixed value for simplicity. \n",
    "\n",
    "For heading, try -170° for ascending, and -10° for descending. See what impact this has on the shape of displacements predicted by the Mogi model by varying these values.\n",
    "\n",
    "While we don't perform this here, we could combine the different look directions of ascending and descending satellite passes to estimate the East-West and vertical displacements if we assume the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1012c3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# params\n",
    "#####################################\n",
    "heading = -10 # degrees\n",
    "incidence = 38 # degrees\n",
    "#####################################\n",
    "\n",
    "# create grids of these angles\n",
    "heading_grid = np.ones(xx.shape) * heading\n",
    "incidence_grid = np.ones(xx.shape) * incidence\n",
    "\n",
    "# generate conversions for each components to LOS\n",
    "x2los = np.cos(np.deg2rad(heading_grid)) * np.sin(np.deg2rad(incidence_grid))\n",
    "y2los = np.sin(np.deg2rad(heading_grid)) * np.sin(np.deg2rad(incidence_grid))\n",
    "z2los = np.cos(np.deg2rad(incidence_grid))\n",
    "\n",
    "# apply projection\n",
    "ux_los = ux * x2los\n",
    "uy_los = uy * y2los\n",
    "uz_los = uz * z2los\n",
    "\n",
    "# calculate total displacement in los \n",
    "ulos = ux_los + uy_los + uz_los\n",
    "\n",
    "# plot our new displacements\n",
    "fig, axs = plt.subplots(1,4,figsize=(15,10))\n",
    "fig.tight_layout(w_pad=0.1)\n",
    "\n",
    "# displacement in x\n",
    "im = axs[0].imshow(ux_los*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.vik)\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=11)\n",
    "cbar.ax.tick_params(labelsize=11)\n",
    "axs[0].set_title('Component of x in line-of-sight')\n",
    "axs[0].set_xlabel('x-coord (km)', fontsize=11)\n",
    "axs[0].set_ylabel('y-coord (km)', fontsize=11)\n",
    "\n",
    "# displacement in y\n",
    "im = axs[1].imshow(uy_los*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.vik)\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=11)\n",
    "cbar.ax.tick_params(labelsize=11)\n",
    "axs[1].set_title('Component of y in line-of-sight')\n",
    "axs[1].set_xlabel('x-coord (km)', fontsize=11)\n",
    "axs[1].set_yticklabels([])\n",
    "axs[1].set_yticks([])\n",
    "\n",
    "# displacement in z\n",
    "im = axs[2].imshow(uz_los*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[2])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=11)\n",
    "cbar.ax.tick_params(labelsize=11)\n",
    "axs[2].set_title('Component of z in line-of-sight')\n",
    "axs[2].set_xlabel('x-coord (km)', fontsize=11)\n",
    "axs[2].set_yticklabels([])\n",
    "axs[2].set_yticks([])\n",
    "\n",
    "# displacement in LOS\n",
    "im = axs[3].imshow(ulos*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.batlow)\n",
    "#axs[3].contour(y, x, ulos*1000)\n",
    "divider = make_axes_locatable(axs[3])\n",
    "cax = divider.append_axes(\"bottom\", size=\"7%\", pad=\"15%\")\n",
    "cbar = plt.colorbar(im, cax=cax, orientation=\"horizontal\")\n",
    "cbar.set_label(label='Displacement (mm)', fontsize=11)\n",
    "cbar.ax.tick_params(labelsize=11)\n",
    "axs[3].set_title('Total displacement in LOS')\n",
    "axs[3].set_xlabel('x-coord (km)', fontsize=11)\n",
    "axs[3].set_ylabel('y-coord (km)', fontsize=11)\n",
    "axs[3].yaxis.tick_right()\n",
    "axs[3].yaxis.set_label_position(\"right\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab0a533b",
   "metadata": {},
   "source": [
    "Let's compare our original 3D displacement figure with a wrapped version of our LOS displacements, which should make the change in the shape of the signal more apparent. We'll wrap the data into fringes of 10 mm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "881ebd4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,10))\n",
    "fig.tight_layout(w_pad=10)\n",
    "\n",
    "# 3D displacements\n",
    "im = axs[0].imshow(uz*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.batlow)\n",
    "axs[0].quiver(xx_ds, yy_ds, ux_ds, uy_ds, color='white')\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Vertical displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0].set_title('3D displacement', fontsize=20)\n",
    "axs[0].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[0].set_ylabel('y-coord (km)', fontsize=18)\n",
    "axs[0].tick_params(labelsize=16)\n",
    "\n",
    "# wrapped LOS\n",
    "im = axs[1].imshow(np.remainder(ulos*1000, 10), extent=[xmin, xmax, ymin, ymax], cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[1].set_title('Wrapped line-of-sight displacement (mm)', fontsize=20)\n",
    "axs[1].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[1].tick_params(labelsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2be5c1c9",
   "metadata": {},
   "source": [
    "Let's again profile the displacements and compare the new LOS displacements to the original x and z displacements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5576b43c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# profile los across y=0\n",
    "prof_los = ulos[y0_ind,:]\n",
    "\n",
    "# plot profiles\n",
    "fig = plt.figure(figsize=(14,5))\n",
    "\n",
    "plt.plot(x, prof_uz*1000, color='red', label='Displacement in z')\n",
    "plt.plot(x, prof_ux*1000, color='blue', label='Displacement in x')\n",
    "plt.plot(x, prof_los*1000, color='orange', label='Displacement in LOS')\n",
    "plt.legend(fontsize=16)\n",
    "plt.xlabel('Distance (km)', fontsize=18)\n",
    "plt.ylabel('Displacement (mm)', fontsize=18)\n",
    "plt.tick_params(labelsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c774ba2",
   "metadata": {},
   "source": [
    "## 3. Impact of topography\n",
    "\n",
    "Our Mogi source currently assumes that the Earth’s surface is entirely flat. In the real world, topography above the source will change the shape and magnitude of displacements at the Earth’s surface. The full impact of topography on the displacements can be assessed using numerical modelling (e.g., finite element or boundary element). \n",
    "\n",
    "However, we can make a simple approximation of the effects of topography on our Mogi source by correcting just for the changes in distance between source and the surface by varying the source depth, which is a good approximation of the impact of topography on vertical displacements (Figure 3, Williams and Wadge, 1998).   \n",
    "\n",
    "<img src=\"figures/topo_correction.png\" style=\"float: centre;\" width=\"400\"/>\n",
    "\n",
    "*Figure 3: Illustration of the vary-depth approximation for the impact of topography on surface displacements as described by Williams and Wadge (1998).  From Lisowski (2007), Figure 8.32.*\n",
    "\n",
    "- Williams, C.A. and Wadge, G., 1998. The effects of topography on magma chamber deformation models: Application to Mt. Etna and radar interferometry. Geophysical Research Letters, 25(10), pp.1549-1552. \n",
    "\n",
    "We'll load in elevations for Domuyo volcano (more on that below). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3d2cca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# source parameters\n",
    "#####################################\n",
    "source_depth = 3 # depth to point (km)\n",
    "dV = 1e7 # change in volume (m^3)\n",
    "nu = 0.25 # poisson's ratio for medium\n",
    "#####################################\n",
    "\n",
    "# load in a topography file\n",
    "synth_elev = np.load('data/synth_dem.npy')\n",
    "\n",
    "# shift the elevations so that the mean is at 0\n",
    "synth_elev = synth_elev - np.nanmean(synth_elev)\n",
    "\n",
    "# generate mogi source with topography\n",
    "ux_topo, uy_topo, uz_topo = lib.mogi_source_topo(xx*1000, yy*1000, synth_elev, d=source_depth*1000, dV=dV, nu=nu)\n",
    "\n",
    "# apply projection to LOS\n",
    "ux_topo_los = ux_topo * x2los\n",
    "uy_topo_los = uy_topo * y2los\n",
    "uz_topo_los = uz_topo * z2los\n",
    "\n",
    "# calculate total displacement in LOS \n",
    "ulos_topo = ux_topo_los + uy_topo_los + uz_topo_los\n",
    "\n",
    "# plot topography, original LOS dispalcements, and topographically corrected LOS displacements\n",
    "fig, axs = plt.subplots(1,3,figsize=(15,10))\n",
    "fig.tight_layout(w_pad=10)\n",
    "\n",
    "# topography\n",
    "im = axs[0].imshow(synth_elev/1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.bamako)\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Elevation (km)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0].set_title('Topography', fontsize=20)\n",
    "axs[0].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[0].set_ylabel('y-coord (km)', fontsize=18)\n",
    "axs[0].tick_params(labelsize=16)\n",
    "\n",
    "# original LOS\n",
    "im = axs[1].imshow(ulos*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[1].set_title('Original', fontsize=20)\n",
    "axs[1].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[1].tick_params(labelsize=16)\n",
    "\n",
    "# topography-corrected\n",
    "im = axs[2].imshow(ulos_topo*1000, extent=[xmin, xmax, ymin, ymax], cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[2])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[2].set_title('Topography-corrected', fontsize=20)\n",
    "axs[2].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[2].tick_params(labelsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f77d6be",
   "metadata": {},
   "source": [
    "Let's take profiles through y = 0 again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d871ecb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# profile across y=0\n",
    "prof_u_orig = ulos[y0_ind,:]\n",
    "prof_u_topo = ulos_topo[y0_ind,:]\n",
    "prof_elev = synth_elev[y0_ind,:]\n",
    "\n",
    "# plot profiles\n",
    "fig, ax1 = plt.subplots(figsize=(14,5))\n",
    "\n",
    "ax1.plot(x, prof_u_orig*1000, color='red', label='Original')\n",
    "ax1.plot(x, prof_u_topo*1000, color='blue', label='Topography-corrected')\n",
    "ax1.set_xlabel('Distance (km)', fontsize=18)\n",
    "ax1.set_ylabel('LOS displacement (mm)', fontsize=18)\n",
    "ax1.tick_params(labelsize=16)\n",
    "ax1.legend(loc='upper left',fontsize=16)\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(x, prof_elev/1000, color='green', label='Elevation')\n",
    "ax2.set_ylabel('Elevation (km)', fontsize=18)  # we already handled the x-label with ax1\n",
    "ax2.tick_params(labelsize=16)\n",
    "ax2.legend(fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a885afa1",
   "metadata": {},
   "source": [
    "## 4. Real world example - Domuyo\n",
    "\n",
    "We'll now apply our Mogi model to a real world example. We'll be looking at Domuyo volcano in Argentina, which showed relatively consistent uplift between 2015 and 2020. Since the end of 2020, the volcano has begun to subside again. We'll model just the 2015-2020 period using our Mogi source. \n",
    "\n",
    "Domuyo is a volcanic complex in Argentina with no historical eruptions, but with a geothermal reservoir that extends ~15 km to the southwest of the volcano. It uplifted between 2015 and 2020, coinciding with a period of surface warming (Lundgren et al., 2020).\n",
    "\n",
    "The time series for Domuyo can be found on the COMET Volcano Portal at:\n",
    "https://comet.nerc.ac.uk/comet-volcano-portal/volcano-index/South%20America/Argentina/Domuyo/S1_analysis\n",
    "\n",
    "<img src=\"figures/Domuyo_from_portal.png\" style=\"float: centre;\" width=\"1000\"/>\n",
    "\n",
    "*Figure 4: Total displacement (left) and the displacement time series (right) from frame 018A_12668_131313 for Domuyo volcano, taken from the COMET Volcano Portal.*\n",
    "\n",
    "This time series is produced from automatically processed Sentinel-1 imagery – you can read more about how it was constructed and how to interpret it here: https://comet.nerc.ac.uk/comet-volcano-portal/about-tools\n",
    "\n",
    "Background information on the deformation at Domuyo can be found at:\n",
    "https://www.nature.com/articles/s41598-020-67982-8\n",
    "\n",
    "The volcano portal allows us to download a csv showing the cumulative displacement, coherence, and elevation at any time step for every point in the grid. I've provided two of these csv files, which we will load in and regrid.\n",
    "\n",
    "The columns of these csv files are:\n",
    "**latitude, longitude, displacement (mm), coherence, elevation (m)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73e7f360",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# use numpy to open the two csv files\n",
    "input_asc = np.loadtxt('data/domuyo_018A_12668_131313_disp_raw.csv', delimiter=',', skiprows=1)\n",
    "input_desc = np.loadtxt('data/domuyo_083D_12636_131313_disp_raw.csv', delimiter=',', skiprows=1)\n",
    "\n",
    "# We can see that every grid position is included, with nan's used for padding,\n",
    "# meaning we can simply reshape the columns into grids.\n",
    "\n",
    "# get lat and lon coords by pulling the unique values from those columns\n",
    "# these data are on the same grid, so we can use the ascending coords for both\n",
    "lat = np.flip(np.unique(input_asc[:,0]))\n",
    "lon = np.unique(input_asc[:,1])\n",
    "\n",
    "# get axes limits\n",
    "lat_min, lat_max = np.min(lat), np.max(lat)\n",
    "lon_min, lon_max = np.min(lon), np.max(lon)\n",
    "\n",
    "# reshape other columns into grids based on size on lat and lon\n",
    "disp_asc = np.flipud(np.reshape(input_asc[:,2], (len(lat), len(lon))))\n",
    "disp_desc = np.flipud(np.reshape(input_desc[:,2], (len(lat), len(lon))))\n",
    "\n",
    "# same for coherence\n",
    "coh_asc = np.flipud(np.reshape(input_asc[:,3], (len(lat), len(lon))))\n",
    "coh_desc = np.flipud(np.reshape(input_desc[:,3], (len(lat), len(lon))))\n",
    "\n",
    "# and same for elevation\n",
    "elev_asc = np.flipud(np.reshape(input_asc[:,4], (len(lat), len(lon))))\n",
    "elev_desc = np.flipud(np.reshape(input_desc[:,4], (len(lat), len(lon))))\n",
    "\n",
    "# plot displacements, coherence, and elevation\n",
    "fig, axs = plt.subplots(2,2,figsize=(15,15))\n",
    "fig.tight_layout(w_pad=10, h_pad=2)\n",
    "\n",
    "# coherence (we'll just show the ascending, almost the same for both)\n",
    "im = axs[0,0].imshow(coh_asc, extent=[lon_min, lon_max, lat_min, lat_max], cmap=cm.grayC.reversed())\n",
    "divider = make_axes_locatable(axs[0,0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Coherence', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0,0].set_title('Coherence', fontsize=20)\n",
    "axs[0,0].set_xlabel('Longitude', fontsize=18)\n",
    "axs[0,0].set_ylabel('Latitude', fontsize=18)\n",
    "axs[0,0].tick_params(labelsize=16)\n",
    "\n",
    "# topography (we'll just show the ascending, identical for both)\n",
    "im = axs[0,1].imshow(elev_asc, extent=[lon_min, lon_max, lat_min, lat_max], cmap=cm.bamako)\n",
    "divider = make_axes_locatable(axs[0,1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Elevation (m)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0,1].set_title('Topography', fontsize=20)\n",
    "axs[0,1].set_xlabel('Longitude', fontsize=18)\n",
    "axs[0,1].tick_params(labelsize=16)\n",
    "\n",
    "# ascending displacement\n",
    "im = axs[1,0].imshow(disp_asc, extent=[lon_min, lon_max, lat_min, lat_max])\n",
    "divider = make_axes_locatable(axs[1,0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[1,0].set_title('018A_12668_131313', fontsize=20)\n",
    "axs[1,0].set_xlabel('Longitude', fontsize=18)\n",
    "axs[1,0].set_ylabel('Latitude', fontsize=18)\n",
    "axs[1,0].tick_params(labelsize=16)\n",
    "\n",
    "# descending displacement\n",
    "im = axs[1,1].imshow(disp_desc, extent=[lon_min, lon_max, lat_min, lat_max])\n",
    "divider = make_axes_locatable(axs[1,1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[1,1].set_title('083D_12636_131313', fontsize=20)\n",
    "axs[1,1].set_xlabel('Longitude', fontsize=18)\n",
    "axs[1,1].tick_params(labelsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf0a886e",
   "metadata": {},
   "source": [
    "Let's now take a profile across the displacements from the ascending frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb8ef38d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# start and end coordinates of the profile, and extra parameters\n",
    "#####################################\n",
    "prof_start = (-70.6, -36.75)\n",
    "prof_end = (-70.3, -36.45)\n",
    "prof_params = {\n",
    "    \"nbins\": 100, # number of bins to split the profile into\n",
    "    \"width\": 0.025 # 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,lat,disp_asc,prof_start,prof_end,prof_params)\n",
    "\n",
    "# plot the ascending displacements and the profiled data\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,5))\n",
    "plt.tight_layout(w_pad=2)\n",
    "\n",
    "# ascending map\n",
    "im = axs[0].imshow(disp_asc, extent=[lon_min, lon_max, lat_min, lat_max], cmap=cm.batlow)\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=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=16)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0].set_title('018A_12668_131313', fontsize=16)\n",
    "axs[0].set_xlabel('Longitude', fontsize=16)\n",
    "axs[0].set_ylabel('Latitude', fontsize=16)\n",
    "axs[0].tick_params(labelsize=14)\n",
    "axs[0].set_xlim(np.amin(lon), np.amax(lon))\n",
    "axs[0].set_ylim(np.amin(lat), np.amax(lat))\n",
    "\n",
    "# velocities projected onto profile\n",
    "axs[1].scatter(points_dist, points_val, s=10)\n",
    "axs[1].set_title('Profiled displacements', fontsize=16)\n",
    "axs[1].set_xlabel(\"Distance (degrees)\", fontsize=16)\n",
    "axs[1].set_ylabel(\"LOS displacement (mm)\", fontsize=16)\n",
    "axs[1].tick_params(labelsize=14)\n",
    "axs[1].yaxis.tick_right()\n",
    "axs[1].yaxis.set_label_position(\"right\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eb5fefc",
   "metadata": {},
   "source": [
    "The left plot shows our line-of-sight displacements for 018A, where motion towards the satellite is positive. The start and end of our profile are marked with the pink and black dots, respectively. The three red lines show the middle and outside edge of the profile. The right plot is every pixel within the outer edges of our profile projected onto the profile line.\n",
    "\n",
    "It would be useful if our data were in metres instead of degrees for modelling, so let's project from lat-long to UTM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87148c53",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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),\n",
    "        south_lat_degree = np.amin(lat),\n",
    "        east_lon_degree  = np.amax(lon),\n",
    "        north_lat_degree = np.amax(lat),\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",
    "x_utm, y_utm = transformer.transform(lat, lon)\n",
    "x_utm, y_utm = x_utm/1000, y_utm/1000 # convert to km\n",
    "xx_utm, yy_utm = np.meshgrid(x_utm, y_utm)\n",
    "x_utm_min, x_utm_max, y_utm_min, y_utm_max = np.amin(x_utm), np.amax(x_utm), np.amin(y_utm), np.amax(y_utm)\n",
    "\n",
    "# plot with new coords\n",
    "fig, axs = plt.subplots(1,2,figsize=(15,7))\n",
    "plt.tight_layout(w_pad=5)\n",
    "\n",
    "# ascending displacements\n",
    "im = axs[0].imshow(disp_asc, extent=[x_utm_min, x_utm_max, y_utm_min, y_utm_max], cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[0])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[0].set_title('018A_12668_131313', fontsize=20)\n",
    "axs[0].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[0].set_ylabel('y-coord (km)', fontsize=18)\n",
    "axs[0].tick_params(labelsize=16)\n",
    "axs[0].set_xlim(x_utm_min, x_utm_max)\n",
    "axs[0].set_ylim(y_utm_min, y_utm_max)\n",
    "\n",
    "# ascending displacements\n",
    "im = axs[1].imshow(disp_desc, extent=[x_utm_min, x_utm_max, y_utm_min, y_utm_max], cmap=cm.batlow)\n",
    "divider = make_axes_locatable(axs[1])\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"5%\")\n",
    "cbar = plt.colorbar(im, cax=cax)\n",
    "cbar.set_label(label='Line-of-sight displacement (mm)', fontsize=18)\n",
    "cbar.ax.tick_params(labelsize=16)\n",
    "axs[1].set_title('083D_12636_131313', fontsize=20)\n",
    "axs[1].set_xlabel('x-coord (km)', fontsize=18)\n",
    "axs[1].tick_params(labelsize=16)\n",
    "axs[1].set_xlim(x_utm_min, x_utm_max)\n",
    "axs[1].set_ylim(y_utm_min, y_utm_max)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecf75195",
   "metadata": {},
   "source": [
    "We'll now run our forward model, projecting the modelled displacements into the line-of-sight for both frames. We'll plot the modelled displacements, the real displacements, and the residual between the two.\n",
    "\n",
    "We can now use the Mogi forward model in this exercise to make a preliminary assessment of the depth and volume change of the geodetic source active at Domuyo. Edit the source location, volume (dV) and depth to reduce the magnitude of the residuals for both the ascending and descending data. Keep updating these values to get a small residual between the observed and modelled displacements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "947ad85c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup model params\n",
    "#####################################\n",
    "source_depth = 5 # depth to point (km)\n",
    "dV = 2e8 # change in volume (m^3)\n",
    "xcen, ycen = 375, 5945\n",
    "#####################################\n",
    "\n",
    "# poisson's ratio for medium\n",
    "nu = 0.25 \n",
    "\n",
    "# shift topography so mean is at 0 m\n",
    "elev_asc_shift = elev_asc - np.nanmean(elev_asc)\n",
    "\n",
    "# run model\n",
    "ux, uy, uz = lib.mogi_source_topo(xx_utm*1000, yy_utm*1000, elev_asc_shift, xcen=xcen*1000, ycen=ycen*1000, d=source_depth*1000, dV=dV, nu=nu)\n",
    "ux, uy, uz = ux*1000, uy*1000, uz*1000\n",
    "\n",
    "# heading and incidence angles\n",
    "heading_asc, heading_desc = -170, -10 # degrees\n",
    "incidence = 51 # degrees\n",
    "\n",
    "# create grids of these angles\n",
    "heading_grid_asc = np.ones(xx_utm.shape) * heading_asc\n",
    "heading_grid_desc = np.ones(xx_utm.shape) * heading_desc\n",
    "incidence_grid = np.ones(xx_utm.shape) * incidence\n",
    "\n",
    "# generate conversions for each components to LOS (asc)\n",
    "x2los_asc = np.cos(np.deg2rad(heading_grid_asc)) * np.sin(np.deg2rad(incidence_grid))\n",
    "y2los_asc = np.sin(np.deg2rad(heading_grid_asc)) * np.sin(np.deg2rad(incidence_grid))\n",
    "z2los_asc = np.cos(np.deg2rad(incidence_grid))\n",
    "\n",
    "# same for desc\n",
    "x2los_desc = np.cos(np.deg2rad(heading_grid_desc)) * np.sin(np.deg2rad(incidence_grid))\n",
    "y2los_desc = np.sin(np.deg2rad(heading_grid_desc)) * np.sin(np.deg2rad(incidence_grid))\n",
    "z2los_desc = np.cos(np.deg2rad(incidence_grid))\n",
    "\n",
    "# calculate total displacement in los \n",
    "ulos_asc = (ux * x2los_asc) + (uy * y2los_asc) + (uz * z2los_asc)\n",
    "ulos_desc = (ux * x2los_desc) + (uy * y2los_desc) + (uz * z2los_desc)\n",
    "\n",
    "# plot comparison of model, observed, and residual displacements for both frames\n",
    "extents = [x_utm_min, x_utm_max, y_utm_min, y_utm_max]\n",
    "clim = [-50, 600]\n",
    "lib.plot_comparison(ulos_asc, disp_asc, ulos_desc, disp_desc, xcen, ycen, extents, clim)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ad4d6d4",
   "metadata": {},
   "source": [
    "You are likely to have found that a Mogi model, even after an approximate topographic correction, only provides a very rough fit to these data.  There could be many reasons for this, including that the assumptions required by the Mogi model about crustal properties are too simple – e.g., the Earth is neither homogeneous nor purely elastic. It may also be that a point source is not the best geometry for modelling this uplift pattern (e.g., Figure 4).\n",
    "\n",
    "<img src=\"figures/source_models.png\" style=\"float: centre;\" width=\"1000\"/>\n",
    "\n",
    "*Figure 5: Displacement patterns for analytical elastic half space sources for different geometries: (i) A Mogi source, (ii) An opening dyke and (iii) inflating sill (Okada sources) and a penny-shaped crack.*\n",
    "\n",
    "In fact, the uplift at Domuyo between 2015 and 2019 is best fit by a slightly tilted, oblate tabular source with a centre depth of 6.5 km.   Lundgren et al., 2020 apply an inversion method that used compound dislocation models (Nikhoo et al, 2016) to find the best fit geometry for the source (see supplementary figures from Lundgren et al., 2020).   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22da6d65",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "Astort, A., Walter, T.R., Ruiz, F., Sagripanti, L., Nacif, A., Acosta, G. and Folguera, A., 2019. Unrest at Domuyo volcano, Argentina, detected by geophysical and geodetic data and morphometric analysis. Remote Sensing, 11(18), p.2175. \n",
    "\n",
    "Ebmeier, S.K., Andrews, B.J., Araya, M.C., Arnold, D.W.D., Biggs, J., Cooper, C., Cottrell, E., Furtney, M., Hickey, J., Jay, J.J.J.A.V. and Lloyd, R., 2018. Synthesis of global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring & the lateral extent of magmatic domains. Journal of Applied Volcanology, 7(1), pp.1-26.\n",
    "\n",
    "Lisowski, M., 2007. Analytical volcano deformation source models. In Volcano deformation (pp. 279-304). Springer, Berlin, Heidelberg.\n",
    "\n",
    "Lundgren, P., Girona, T., Bato, M.G., Realmuto, V.J., Samsonov, S., Cardona, C., Franco, L., Gurrola, E. and Aivazis, M., 2020. The dynamics of large silicic systems from satellite remote sensing observations: The intriguing case of Domuyo volcano, Argentina. Scientific reports, 10(1), pp.1-15. \n",
    "\n",
    "Magee, C., Stevenson, C.T., Ebmeier, S.K., Keir, D., Hammond, J.O., Gottsmann, J.H., Whaler, K.A., Schofield, N., Jackson, C.A., Petronis, M.S. and O’driscoll, B., 2018. Magma plumbing systems: a geophysical perspective. Journal of Petrology, 59(6), pp.1217-1251. \n",
    "\n",
    "Masterlark, T., 2007. Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions. Journal of Geophysical Research: Solid Earth, 112(B6).\n",
    "\n",
    "Mogi, K. 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Earthq Res Inst, 36, pp.99-134.\n",
    "\n",
    "Nikkhoo, M., Walter, T. R., Lundgren, P. R. & Prats-Iraola, P. Compound dislocation models (CDMs) for volcano deformation analyses. Geophys. J. Int. 208, 877–894 (2016).\n",
    "\n",
    "Pinel, V., Poland, M. P., & Hooper, A. (2014). Volcanology: Lessons learned from synthetic aperture radar imagery. Journal of Volcanology and Geothermal Research, 289, 81-113.\n",
    "\n",
    "Segall, P. (2019). Magma chambers: what we can, and cannot, learn from volcano geodesy. Philosophical Transactions of the Royal Society A, 377(2139), 20180158\n",
    "\n",
    "Williams, C.A. and Wadge, G., 1998. The effects of topography on magma chamber deformation models: Application to Mt. Etna and radar interferometry. Geophysical Research Letters, 25(10), pp.1549-1552."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ca57797",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}