seismic_petrophysics_interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating synthetic data, the INTERACTIVE notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This IPython notebook is intended as a companion to my [Seismic Petrophysics](http://nbviewer.ipython.org/github/aadm/geophysical_notes/blob/master/seismic_petrophysics.ipynb)  notebook to demonstrate how to use its [interactive cabalities](https://github.com/ipython/ipywidgets). I'm sure this feature was introduced for educational purposes (or maybe simply to show off) but I can see a future where geophysicists will build some sort of interactive reports for asset teams or clients to allow them to do exploratory data analysis.\n",
    "\n",
    "For all the details please refer to the above notebook and this github repo: <https://github.com/aadm/geophysical_notes>.\n",
    "\n",
    "_Note: To see the interactive features in action you have to download and run this notebook locally._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "from ipywidgets import interact, interactive\n",
    "from IPython.display import Audio, display\n",
    "\n",
    "%matplotlib inline\n",
    "# comment out the following if you're not on a Mac with HiDPI display\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.colors as colors\n",
    "#      0=undef   1=bri  2=oil   3=gas 4=shale\n",
    "ccc = ['#B3B3B3','blue','green','red','#996633',]\n",
    "cmap_facies = colors.ListedColormap(ccc[0:len(ccc)], 'indexed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "L=pd.read_csv('qsiwell2_frm_facies.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code builds a summary plot showing the following tracks, from left to right:\n",
    "\n",
    "1. petrophysical logs\n",
    "2. Ip log\n",
    "3. Vp/Vs log\n",
    "4. litho-facies class\n",
    "\n",
    "To the extreme right I have an interactive Ip-Vp/Vs crossplot; two sliders are used to select top and bottom limits (in meters) of the analysis window where data is dynamically taken to populate this crossplot.\n",
    "\n",
    "As a visual aid, the extend of the analysis window is displayed as a transparent rectangle in the first track."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "054b2a0aa73a43428b54884e19745cb1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=2150, description='z1', max=2360, min=2000), IntSlider(value=2250, descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Select top (z1) and bottom (z2) of analysis window to populate Ip-Vp/Vs crossplot.\n",
      "Analysis window will be highlighted in yellow in leftmost track.\n"
     ]
    }
   ],
   "source": [
    "@interact(z1=(2000,2360),z2=(2040,2400))\n",
    "def crossplot(z1=2150,z2=2250):\n",
    "    from matplotlib.patches import Rectangle\n",
    "    ll=L[(L.index>=2000) & (L.index<=2400)]\n",
    "    cluster=np.repeat(np.expand_dims(ll['LFC'].values,1),100,1)\n",
    "    \n",
    "    fig = plt.figure(figsize=(10,5))\n",
    "    ax0 = plt.subplot2grid((1,10), (0,0), colspan=2)\n",
    "    ax1 = plt.subplot2grid((1,10), (0,2), colspan=2,sharey=ax0)\n",
    "    ax2 = plt.subplot2grid((1,10), (0,4), colspan=2,sharey=ax0)\n",
    "    ax3 = plt.subplot2grid((1,10), (0,6))\n",
    "    ax4 = plt.subplot2grid((1,10), (0,7), colspan=3)\n",
    "    \n",
    "    ax0.plot(ll.VSH, ll.index, '-g', label='Vsh')\n",
    "    ax0.plot(ll.SWE, ll.index, '-b', label='Sw')\n",
    "    ax0.plot(ll.PHIE, ll.index, '-k', label='phi')\n",
    "    ax1.plot(ll.IP, ll.index, '-k')\n",
    "    ax2.plot(ll.VPVS, ll.index, '-k')\n",
    "    im=ax3.imshow(cluster, interpolation='none', aspect='auto',cmap=cmap_facies,vmin=0,vmax=4)  \n",
    "    ll=L[(L.index>=z1) & (L.index<=z2)]\n",
    "    ax4.scatter(ll.IP,ll.VPVS,20,ll.LFC,marker='o',edgecolors='none',alpha=0.7,cmap=cmap_facies,vmin=0,vmax=4)\n",
    "\n",
    "    ax0.invert_yaxis(), ax3.set_yticklabels([]), ax3.set_xticklabels([])\n",
    "    ax0.legend(fontsize='small', loc='lower right')\n",
    "    ax0.set_xlabel('Vcl/phi/Sw'),ax0.set_xlim(-.1,1.1)\n",
    "    ax1.set_xlabel('Ip [m/s*g/cc]'),ax1.set_xlim(3000,9000)\n",
    "    ax2.set_xlabel('Vp/Vs'),ax2.set_xlim(1.5,3)\n",
    "    ax3.set_xlabel('LFC')\n",
    "    ax4.set_xlabel('Ip'),ax4.set_ylabel('Vp/Vs'),ax4.set_xlim(3000,9000),ax4.set_ylim(1.5,3); \n",
    "    fig.tight_layout()     \n",
    "    cbar=plt.colorbar(im, ax=ax4)\n",
    "    cbar.set_label((12*' ').join(['undef', 'brine', 'oil', 'gas', 'shale']))\n",
    "    cbar.set_ticks(range(0,1)); cbar.set_ticklabels('')\n",
    "    ax0.add_patch(Rectangle((0, z1),1,z2-z1,edgecolor='None',facecolor=[0.9,0.9,0.0],alpha=0.5))\n",
    "\n",
    "print('Select top (z1) and bottom (z2) of analysis window to populate Ip-Vp/Vs crossplot.')\n",
    "print('Analysis window will be highlighted in yellow in leftmost track.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# comparing augmented well data and synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First I read in previously computed _augmented_ (i.e., a dataset containing original plus fluid-replaced logs simulating different pore fluid scenarios) and _synthetic_ (i.e., created through Monte Carlo simulation) data.\n",
    "\n",
    "Then I show side by side the two datasets, color-coded by litho-fluid class; a slider is used to select one of the 4 available classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ww=pd.read_csv('qsiwell2_augmented.csv')\n",
    "mc=pd.read_csv('qsiwell2_synthetic.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5b66e09f07cb4861b60aa110985d9964",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1, description='n', max=4, min=1), Output()), _dom_classes=('widget-inte…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Select Litho-Fluid Class: n=1 (brine), 2 (oil), 3 (gas), 4 (shale)\n"
     ]
    }
   ],
   "source": [
    "@interact(n=(1,4))\n",
    "def crossplot(n=1):\n",
    "    f, ax = plt.subplots(nrows=1, ncols=2, sharey=True, sharex=True, figsize=(8, 6))\n",
    "    scatt1=ax[0].scatter(ww.IP[ww.LFC==n],ww.VPVS[ww.LFC==n],20,ww.LFC[ww.LFC==n],marker='o',edgecolors='none', alpha=0.5, cmap=cmap_facies,vmin=0,vmax=4)\n",
    "    scatt2=ax[1].scatter(mc.IP[mc.LFC==n],mc.VPVS[mc.LFC==n],20,mc.LFC[mc.LFC==n],marker='o',edgecolors='none', alpha=0.5, cmap=cmap_facies,vmin=0,vmax=4)\n",
    "    ax[0].set_xlim(3000,9000), ax[0].set_ylim(1.5,3.0)\n",
    "    ax[0].set_title('augmented well data');\n",
    "    ax[1].set_title('synthetic data');\n",
    "    for i in range(len(ax)): ax[i].grid()  \n",
    "    uu=f.add_axes([0.95, 0.2, 0.015, 0.6])\n",
    "    cbar=f.colorbar(scatt2, cax=uu);\n",
    "    cbar.set_label((10*' ').join(['undef', 'brine', 'oil', 'gas', 'shale']))\n",
    "    cbar.set_ticks(range(0,1)); cbar.set_ticklabels('')\n",
    "    cbar.set_alpha(1)\n",
    "    cbar.draw_all()\n",
    "    return n\n",
    "print('Select Litho-Fluid Class: n=1 (brine), 2 (oil), 3 (gas), 4 (shale)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this is a second qc, showing a comparison between the Ip and Vp/Vs PDFs calculated for real (augmented) vs synthetic datasets, with the possibility to slide through the different classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6e89e8b1e5214e648474be9d9824737f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=2, description='n', max=4, min=1), Output()), _dom_classes=('widget-inte…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Select Litho-Fluid Class: n=1 (brine), 2 (oil), 3 (gas), 4 (shale)\n"
     ]
    }
   ],
   "source": [
    "@interact(n=(1,4))\n",
    "def histcompare(n):\n",
    "    f, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 4))\n",
    "\n",
    "    ww[ww.LFC==n]['IP'].plot(kind='density',style='k-',ax=ax[0],label='real',legend=True,title='Ip')\n",
    "    mc[mc.LFC==n]['IP'].plot(kind='density',style='k--',ax=ax[0],label='synthetic',legend=True)\n",
    "    ax[0].set_xlim(3000,9000)\n",
    "    ww[ww.LFC==n]['VPVS'].plot(kind='density',style='k-',ax=ax[1],label='real',legend=True,title='Vp/Vs')\n",
    "    mc[mc.LFC==n]['VPVS'].plot(kind='density',style='k--',ax=ax[1],label='synthetic',legend=True)\n",
    "    ax[1].set_xlim(1.5,3.0)\n",
    "    return n\n",
    "\n",
    "print('Select Litho-Fluid Class: n=1 (brine), 2 (oil), 3 (gas), 4 (shale)')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.5"
  },
  "widgets": {
   "state": {
    "585145afeb1a459593b3481e153b4069": {
     "views": [
      {
       "cell_index": 13
      }
     ]
    },
    "6d7b8cd98a3e401790c4bcd7051958ab": {
     "views": [
      {
       "cell_index": 7
      }
     ]
    },
    "f308a8c4272f4efbac605d48da564df7": {
     "views": [
      {
       "cell_index": 11
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}