notebooks/Getting_started.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Getting started\n",
    "\n",
    "## Importing used libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing GemPy\n",
    "import gempy as gp\n",
    "\n",
    "# Importing aux libraries\n",
    "from ipywidgets import interact\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "\n",
    "# Embedding matplotlib figures in the notebooks\n",
    "%matplotlib qt5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initializing the model:\n",
    "\n",
    "The first step to create a GemPy model is create a gempy.Model object that will contain all the other data structures and necessary functionality.\n",
    "\n",
    "In addition for this example we will define a regular grid since the beginning. This is the grid where we will interpolate the 3D geological model. GemPy comes with an array of different grids for different pourposes as we will see below. For visualization usually a regular grid is the one that makes more sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Active grids: ['regular']\n"
     ]
    }
   ],
   "source": [
    "geo_model = gp.create_model('Model1')\n",
    "geo_model = gp.init_data(geo_model, extent= [0, 791, 0, 200, -582, 0], resolution=[100, 10, 100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GemPy core code is written in Python. However for efficiency (and other reasons) most of heavy computations happend in optimize compile code, either C or CUDA for GPU. To do so, GemPy rely on the library theano. To guarantee maximum optimization theano requires to compile the code for every Python kernel. The compilation is done by calling the following line at any point (before computing the model):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setting kriging parameters to their default values.\n",
      "Compiling theano function...\n",
      "Level of Optimization:  fast_compile\n",
      "Device:  cpu\n",
      "Precision:  float64\n",
      "Number of faults:  0\n",
      "Compilation Done!\n",
      "Kriging values: \n",
      "                   values\n",
      "range             1002.2\n",
      "$C_o$            23914.4\n",
      "drift equations      [3]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<gempy.core.interpolator.InterpolatorModel at 0x7fec7956c990>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gp.set_interpolator(geo_model, theano_optimizer='fast_compile', verbose=[])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating figure:\n",
    "\n",
    "GemPy uses matplotlib and pyvista-vtk libraries for 2d and 3d visualization of the model respectively. One of the design decisions of GemPy is to allow real time construction of the model. What this means is that you can start adding input data and see in real time how the 3D surfaces evolve. Lets initialize the visualization windows.\n",
    "\n",
    "The first one is the 2d figure. Just place the window where you can see it (maybe move the jupyter notebook to half screen and use the other half for the renderers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt5\n",
    "import pyvista as pv\n",
    "pv.set_plot_theme(\"document\")\n",
    "\n",
    "p2d = gp.plot_2d(geo_model, section_names=[], direction=None, show=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Add model section\n",
    "\n",
    "In the 2d renderer we can add several cross section of the model. In this case, for simplicity sake we are just adding one perpendicular to y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this case perpendicular to the y axes\n",
    "ax = p2d.add_section(cell_number=1, direction='y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading cross-section image:\n",
    "\n",
    "Remember that gempy is simply using matplotlib and therofe the ax object created above is a standard matplotlib axes. This allow to manipulate it freely. Lets load an image with the information of couple of boreholes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reading image\n",
    "img = mpimg.imread('../examples/getting_started/wells.png')\n",
    "# Plotting it inplace\n",
    "ax.imshow(img, origin='upper', alpha=.8, extent = (0, 791, -582,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do the same in 3D through pyvista and vtk rendering. Click the qt5 button Back (+Y) to have the same view as in the 2D viwer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3d = gp.plot_3d(geo_model, plotter_type='background', notebook=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the model\n",
    "\n",
    "Now that we have everything initialize we can start the construction of the geological model. \n",
    "\n",
    "### Surfaces\n",
    "\n",
    "GemPy is a surface based interpolator. This means that all the input data we add has to be refered to a surface. The surfaces always mark the bottom of a unit. By default GemPy surfaces are empty:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.surfaces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we do not care about the names and we just want to interpolate a surface we can use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Default surfaces:\n",
    "geo_model.set_default_surfaces()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can start adding data. GemPy input data consist on surface points and orientations (perpendicular to the layers). The 2D plot gives you the X and Z coordinates when hovering the mouse over. We can add a surface point as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a point\n",
    "geo_model.add_surface_points(X=223, Y=0.01, Z=-94, surface='surface1')\n",
    "\n",
    "# Plot in 2D\n",
    "p2d.plot_data(ax, cell_number=11)\n",
    "\n",
    "# Plot in 3D\n",
    "p3d.plot_surface_points()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can add the other two points of the layer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add points\n",
    "geo_model.add_surface_points(X=458, Y=0, Z=-107, surface='surface1')\n",
    "geo_model.add_surface_points(X=612, Y=0, Z=-14, surface='surface1')\n",
    "\n",
    "# Plotting\n",
    "p2d.plot_data(ax, cell_number=11)\n",
    "p3d.plot_surface_points()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum amount of data to interpolate anything in gempy is:\n",
    "a) 2 surface points per surface\n",
    "b) One orientation per series.\n",
    "\n",
    "Lets add an orientation anywhere in space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding orientation\n",
    "geo_model.add_orientations(X=350, Y=0, Z=-300, surface='surface1', pole_vector= (0,0,1))\n",
    "p2d.plot_data(ax, cell_number=5)\n",
    "p3d.plot_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have enough data for finally interpolate!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.compute_model(geo_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.additional_data.kriging_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, we have interpolated the 3D surface. We can visualize with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In 2D\n",
    "p2d.plot_contacts(ax, cell_number=5)\n",
    "\n",
    "# In 3D\n",
    "p3d.plot_surfaces()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding more layers:\n",
    "\n",
    "So far we only need 2 units defined. The cross-section image that we load have 4 however. Lets add two layers more:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.surfaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.add_surfaces(['surface3', 'basement'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Layer 2\n",
    "\n",
    "Add the layer next layers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here:\n",
    "geo_model.add_surface_points(X=225, Y=0, Z=-269, surface='surface2')\n",
    "geo_model.add_surface_points(X=459, Y=0, Z=-279, surface='surface2')\n",
    "\n",
    "\n",
    "#--------------------\n",
    "# Plot data\n",
    "p2d.remove(ax)\n",
    "p2d.plot_data(ax, cell_number=5, legend='force')\n",
    "p3d.plot_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute model\n",
    "gp.compute_model(geo_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.additional_data.kriging_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot 2D\n",
    "p2d.plot_lith(ax, cell_number=5)\n",
    "p2d.plot_contacts(ax, cell_number=5)\n",
    "\n",
    "# Plot 3D\n",
    "p3d.plot_surfaces()\n",
    "p3d.plot_structured_grid(opacity=.2, annotations = {1: 'surface1', 2:'surface2', 3:'surface3'})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Layer 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here:\n",
    "geo_model.add_surface_points(X=225, Y=0, Z=-439, surface='surface3')\n",
    "geo_model.add_surface_points(X=464, Y=0, Z=-456, surface='surface3')\n",
    "geo_model.add_surface_points(X=619, Y=0, Z=-433, surface='surface3')\n",
    "\n",
    "# ------------------\n",
    "# Plotting data\n",
    "p2d.remove(ax)\n",
    "p2d.plot_data(ax, cell_number=5)\n",
    "p3d.plot_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computing and plotting 3D\n",
    "gp.compute_model(geo_model)\n",
    "p2d.plot_lith(ax, cell_number=5)\n",
    "p2d.plot_contacts(ax, cell_number=5)\n",
    "\n",
    "p3d.plot_surfaces()\n",
    "p3d.plot_structured_grid(opacity=.2, annotations = {1: 'surface1', 2:'surface2', 3:'surface3', 4:'basement'})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.additional_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Faults:\n",
    "\n",
    "So far the model is simply a depositional unit. GemPy allows for unconformities and faults to build complex models. This input is given by categorical data. In general:\n",
    "\n",
    "input data (surface points/ orientations) <belong to< surface <belong to< series\n",
    "\n",
    "And series can be a fault---i.e. offset the rest of surface--- or not. We are going to show how to add a fault as an example.\n",
    "\n",
    "First we need to add a series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.add_series('Fault1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.reorder_series(['Fault1', 'Default series'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then define that is a fault:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.set_is_fault('Fault1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we also need to add a new surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.add_surfaces('fault1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally assign the new surface to the new series/fault"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.map_series_to_surfaces(geo_model, {'Fault1':'fault1'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can just add input data as before (remember the minimum amount of input data to compute a model):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add input data of the fault\n",
    "geo_model.add_surface_points(X=550, Y=0, Z=-30, surface='fault1')\n",
    "geo_model.add_surface_points(X=650, Y=0, Z=-200, surface='fault1')\n",
    "geo_model.add_orientations(X=600, Y=0, Z= -100, surface='fault1', pole_vector=(.3,0,.3))\n",
    "\n",
    "# Plotting Inpute data\n",
    "p2d.remove(ax)\n",
    "p2d.plot_data(ax, cell_number=5)\n",
    "p3d.plot_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now is computing as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute\n",
    "gp.compute_model(geo_model)\n",
    "\n",
    "# Plot\n",
    "p2d.plot_lith(ax, cell_number=5)\n",
    "p2d.plot_contacts(ax, cell_number=5)\n",
    "p3d.plot_surfaces()\n",
    "p3d.plot_structured_grid(opacity=.2, annotations = {2: 'surface1', 3:'surface2', 4:'surface3', 5:'basement'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see now instead of having folding layers we have a sharp jump. Building on this you can pretty much any model you can imagine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional features:\n",
    "\n",
    "Over the years we have built a bunch of assets integrate with gempy. Here we will show some of them:\n",
    "\n",
    "### Topogragpy\n",
    "\n",
    "GemPy has a built-in capanilities to read and manipulate topographic data (through gdal). To show an example we can just create a random topography:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Adding random topography\n",
    "geo_model.set_topography(source='random',fd=1.9, d_z=np.array([-150, 0]), resolution=np.array([200,200]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The topography can we visualize in both renderers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2d.plot_topography(ax, cell_number=5)\n",
    "p3d.plot_topography(scalars='topography')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But also allows us to compute the geological map of an area:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.compute_model(geo_model)\n",
    "p3d.plot_surfaces()\n",
    "p3d.plot_topography()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3d.plot_structured_grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gravity inversion\n",
    "\n",
    "GemPy also allows for inversions (in production only gravity so far). We can see a small demo how this works.\n",
    "\n",
    "The first thing to do is to assign densities to each of the units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.add_surface_values([0, 2.6, 2.4, 3.2, 3.6], ['density'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also we can create a centered grid around a device for precision:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_model.set_centered_grid(centers= [[400, 0, 0]],  resolution = [10, 10, 100], radius=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to modify the compile code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.set_interpolator(geo_model, output=['gravity'], theano_optimizer='fast_run', )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But now additionally to the interpolation we also compute the forward gravity of the model (at the point XYZ = 400, 0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.compute_model(geo_model)\n",
    "geo_model.solutions.fw_gravity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize it it better in the following figure. The aim of an inversion is to find the set of parameters that fit a measured point the better. In this example the red x symbolize the measured gravity while the blue dots are the current gravity (and previous) fw gravity values. The widget moves up and down the surface 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating new figure\n",
    "p_grav = gp.plot_2d(geo_model, direction=None, show=False, figsize=(6, 9))\n",
    "\n",
    "# Adding section for model\n",
    "ax2 = p_grav.add_section(cell_number=1, direction='y', ax_pos=211)\n",
    "\n",
    "# Adding section for gravity \n",
    "ax3 = p_grav.add_section(ax_pos=414)\n",
    "\n",
    "# Plotting model in section\n",
    "p2d.plot_data(ax2, cell_number=5, legend='force')\n",
    "p2d.plot_lith(ax2, cell_number=5)\n",
    "p2d.plot_contacts(ax2, cell_number=5)\n",
    "ax2.plot(400, 0, '^', markersize=40, c='red')\n",
    "\n",
    "# Plotting initial values of the gravity axes\n",
    "target_grav = -810\n",
    "ax3.tick_params( bottom=False)\n",
    "ax3.spines['top'].set_visible(False)\n",
    "ax3.spines['bottom'].set_visible(False)\n",
    "ax3.spines['right'].set_visible(False)\n",
    "ax3.plot(0, target_grav, 'X', label='Target Gravity', markersize=4, c='red')\n",
    "ax3.plot(1, geo_model.solutions.fw_gravity, 'o', label='Current Gravity', markersize=4, c='blue')\n",
    "ax3.set_ylabel('grav')\n",
    "\n",
    "# We store the original values of z for the surface 3\n",
    "Z_ori = geo_model.surface_points.df.loc[[5,6,7], 'Z'].copy()\n",
    "# init a list to store grav\n",
    "grav_ = []\n",
    "\n",
    "# Function that modify the model, compute it and plot\n",
    "def gravity_invert(dz):\n",
    "    new_z = Z_ori + dz\n",
    "    geo_model.modify_surface_points(indices=[5, 6, 7], Z= new_z)\n",
    "    gp.compute_model(geo_model)\n",
    "    grav_.append(geo_model.solutions.fw_gravity[0])\n",
    "    \n",
    "    p2d.remove(ax2)\n",
    "    p2d.plot_data(ax2, cell_number=5)\n",
    "    p2d.plot_lith(ax2, cell_number=5)\n",
    "    p2d.plot_contacts(ax2, cell_number=5)\n",
    "    ax3.plot(np.arange(len(grav_))+1, grav_, 'o', label='Current Gravity', markersize=4, c='blue')\n",
    "    ax3.set_xlim(-1, len(grav_)+1)\n",
    "    \n",
    "interact(gravity_invert, dz=(-200, 200, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Creating new figure\n",
    "p_grav = gp.plot_2d(geo_model, direction=None, show=False, figsize=(6, 9))\n",
    "\n",
    "# Adding section for model\n",
    "ax2 = p_grav.add_section(cell_number=1, direction='y', ax_pos=211)\n",
    "\n",
    "# Adding section for gravity \n",
    "ax3 = p_grav.add_section(ax_pos=414)\n",
    "\n",
    "# Plotting model in section\n",
    "p2d.plot_data(ax2, cell_number=5, legend='force')\n",
    "p2d.plot_lith(ax2, cell_number=5)\n",
    "p2d.plot_contacts(ax2, cell_number=5)\n",
    "ax2.plot(400, 0, '^', markersize=40, c='red')\n",
    "\n",
    "# Plotting initial values of the gravity axes\n",
    "target_grav = -81\n",
    "ax3.tick_params( bottom=False)\n",
    "ax3.spines['top'].set_visible(False)\n",
    "ax3.spines['bottom'].set_visible(False)\n",
    "ax3.spines['right'].set_visible(False)\n",
    "ax3.plot(0, target_grav, 'X', label='Target Gravity', markersize=4, c='red')\n",
    "ax3.plot(1, geo_model.solutions.fw_gravity, 'o', label='Current Gravity', markersize=4, c='blue')\n",
    "ax3.set_ylabel('grav')\n",
    "\n",
    "# We store the original values of z for the surface 3\n",
    "Z_ori = geo_model.surface_points.df.loc[[5,6,7], 'Z'].copy()\n",
    "# init a list to store grav\n",
    "grav_ = []\n",
    "\n",
    "# Function that modify the model, compute it and plot\n",
    "def gravity_invert(dz):\n",
    "    new_z = Z_ori + dz\n",
    "    geo_model.modify_surface_points(indices=[5, 6, 7], Z= new_z)\n",
    "    gp.compute_model(geo_model)\n",
    "    grav_.append(geo_model.solutions.fw_gravity[0])\n",
    "    \n",
    "    p2d.remove(ax2)\n",
    "    p2d.plot_data(ax2, cell_number=5)\n",
    "    p2d.plot_lith(ax2, cell_number=5)\n",
    "    p2d.plot_contacts(ax2, cell_number=5)\n",
    "    ax3.plot((np.arange(len(grav_))+1), grav_, 'o', label='Current Gravity', markersize=4, c='blue')\n",
    "    \n",
    "interact(gravity_invert, dz=(-200, 200, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "split_at_heading": true
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}